# Diff of /src/code/irrat.lisp

revision 1.53 by rtoy, Fri Aug 3 14:28:22 2007 UTC revision 1.54 by rtoy, Mon Jan 28 18:21:03 2008 UTC
# Line 182  Line 182
182
183  #+(or ppc x86)  #+(or ppc x86)
184  (progn  (progn
185  (declaim (inline %ieee754-rem-pi/2))  (declaim (inline %%ieee754-rem-pi/2))
186  ;; Basic argument reduction routine.  It returns two values: n and y  ;; Basic argument reduction routine.  It returns two values: n and y
187  ;; such that (n + 8*k)*pi/2+y = x where |y|<pi/4 and n indicates in  ;; such that (n + 8*k)*pi/2+y = x where |y|<pi/4 and n indicates in
188  ;; which octant the arg lies.  Y is actually computed in two parts,  ;; which octant the arg lies.  Y is actually computed in two parts,
189  ;; y[0] and y[1] such that the sum is y, for accuracy.  ;; y[0] and y[1] such that the sum is y, for accuracy.
190
191  (alien:def-alien-routine ("__ieee754_rem_pio2" %ieee754-rem-pi/2) c-call:int  (alien:def-alien-routine ("__ieee754_rem_pio2" %%ieee754-rem-pi/2) c-call:int
192    (x double-float)    (x double-float)
193    (y (* double-float)))    (y (* double-float)))
194
195    ;; Same as above, but instead of needing to pass an array in, the
196    ;; output array is broken up into two output values instead.  This is
197    ;; easier for the user, and we don't have to wrap calls with
198    ;; without-gcing.
199    (declaim (inline %ieee754-rem-pi/2))
200    (alien:def-alien-routine ("ieee754_rem_pio2" %ieee754-rem-pi/2) c-call:int
201      (x double-float)
202      (y0 double-float :out)
203      (y1 double-float :out))
204
205  )  )
206
207  #+ppc  #+ppc
# Line 207  Line 218
218  #+(or ppc x86)  #+(or ppc x86)
219  (macrolet  (macrolet
220      ((frob (sin cos tan)      ((frob (sin cos tan)
221         `(let ((y (make-array 2 :element-type 'double-float)))         `(progn
222            ;; The array y holds the result for %ieee754-rem-pi/2            ;; The array y holds the result for %ieee754-rem-pi/2
223            ;;            ;;
224            ;; In all of the routines below, we just compute the sum of            ;; In all of the routines below, we just compute the sum of
# Line 222  Line 233
233              (if (< (abs x) (/ pi 4))              (if (< (abs x) (/ pi 4))
234                  (,sin x)                  (,sin x)
235                  ;; Argument reduction needed                  ;; Argument reduction needed
236                  (let* ((n (%ieee754-rem-pi/2 x (vector-sap y)))                  (multiple-value-bind (n y0 y1)
237                         (reduced (+ (aref y 0) (aref y 1))))                      (%ieee754-rem-pi/2 x)
238                    (case (logand n 3)                    (let ((reduced (+ y0 y1)))
239                      (0 (,sin reduced))                      (case (logand n 3)
240                      (1 (,cos reduced))                        (0 (,sin reduced))
241                      (2 (- (,sin reduced)))                        (1 (,cos reduced))
242                      (3 (- (,cos reduced)))))))                        (2 (- (,sin reduced)))
243                          (3 (- (,cos reduced))))))))
244            (defun %cos (x)            (defun %cos (x)
245              (declare (double-float x))              (declare (double-float x))
246              (if (< (abs x) (/ pi 4))              (if (< (abs x) (/ pi 4))
247                  (,cos x)                  (,cos x)
248                  ;; Argument reduction needed                  ;; Argument reduction needed
249                  (let* ((n (%ieee754-rem-pi/2 x (vector-sap y)))                  (multiple-value-bind (n y0 y1)
250                         (reduced (+ (aref y 0) (aref y 1))))                      (%ieee754-rem-pi/2 x)
251                    (case (logand n 3)                    (let ((reduced (+ y0 y1)))
252                      (0 (,cos reduced))                      (case (logand n 3)
253                      (1 (- (,sin reduced)))                        (0 (,cos reduced))
254                      (2 (- (,cos reduced)))                        (1 (- (,sin reduced)))
255                      (3 (,sin reduced))))))                        (2 (- (,cos reduced)))
256                          (3 (,sin reduced)))))))
257            (defun %tan (x)            (defun %tan (x)
258              (declare (double-float x))              (declare (double-float x))
259              (if (< (abs x) (/ pi 4))              (if (< (abs x) (/ pi 4))
260                  (,tan x)                  (,tan x)
261                  ;; Argument reduction needed                  ;; Argument reduction needed
262                  (let* ((n (%ieee754-rem-pi/2 x (vector-sap y)))                  (multiple-value-bind (n y0 y1)
263                         (reduced (+ (aref y 0) (aref y 1))))                      (%ieee754-rem-pi/2 x)
264                    (if (evenp n)                    (let ((reduced (+ y0 y1)))
265                        (,tan reduced)                      (if (evenp n)
266                        (- (/ (,tan reduced))))))))))                          (,tan reduced)
267                            (- (/ (,tan reduced)))))))))))
268    #+x86    #+x86
269    (frob %sin-quick %cos-quick %tan-quick)    (frob %sin-quick %cos-quick %tan-quick)
270    #+ppc    #+ppc

Legend:
 Removed from v.1.53 changed lines Added in v.1.54