# Diff of /src/code/irrat.lisp

revision 1.48 by rtoy, Wed Jul 19 03:52:38 2006 UTC revision 1.49 by rtoy, Wed Jul 19 13:02:45 2006 UTC
# Line 171  Line 171
171      (%sqrt x))      (%sqrt x))
172    )    )
173
174    ;;; The standard libm routines for sin, cos, and tan on x86 (Linux)
175    ;;; and ppc are not very accurate for large arguments when compared to
176    ;;; sparc (and maxima).  This is basically caused by the fact that
177    ;;; those libraries do not do an accurate argument reduction.  The
178    ;;; following functions use some routines Sun's free fdlibm library to
179    ;;; do accurate reduction.  Then we call the standard C functions (or
180    ;;; vops for x86) on the reduced argument.  This produces much more
181    ;;; accurate values.
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
187    ;; 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,
189    ;; 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)))
# Line 194  Line 208
208  (macrolet  (macrolet
209      ((frob (sin cos tan)      ((frob (sin cos tan)
210         `(let ((y (make-array 2 :element-type 'double-float)))         `(let ((y (make-array 2 :element-type 'double-float)))
211              ;; The array y holds the result for %ieee754-rem-pi/2
212              ;;
213              ;; In all of the routines below, we just compute the sum of
214              ;; y[0] and y[1] and use that as the (reduced) argument for
215              ;; the trig functions.  This is slightly less accurate than
216              ;; what fdlibm does, which calls special functions using
217              ;; y[0] and y[1] separately, for greater accuracy.  This
218              ;; isn't implemented, and some spot checks indicate that
219              ;; what we have here is accurate.
220            (defun %sin (x)            (defun %sin (x)
221              (declare (double-float x))              (declare (double-float x))
222              (if (< (abs x) (/ pi 4))              (if (< (abs x) (/ pi 4))

Legend:
 Removed from v.1.48 changed lines Added in v.1.49