Commit 408aa78a authored by Raymond Toy's avatar Raymond Toy
Browse files

Merge branch 'master' into rtoy-simp-dd-trig

parents 00bd409b 01a3f47b
...@@ -187,30 +187,6 @@ ...@@ -187,30 +187,6 @@
(%sqrt x)) (%sqrt x))
) )
;;; The standard libm routines for sin, cos, and tan on x86 (Linux,
;;; 32-bit. 64-bit is apparently ok) and ppc are not very accurate
;;; for large arguments when compared to sparc (and maxima). This is
;;; basically caused by the fact that those libraries do not do an
;;; accurate argument reduction. The following functions use some
;;; routines Sun's free fdlibm library to do accurate reduction. Then
;;; we call the standard C functions (or vops for x86) on the reduced
;;; argument. This produces much more accurate values.
;;;
;;; You can test this by computing (cos (scale-float 1d0 120)). The
;;; true answer is -0.9258790228548379d0.
#+(or ppc x86)
(progn
(declaim (inline %%ieee754-rem-pi/2))
;; Basic argument reduction routine. It returns two values: n and y
;; such that (n + 8*k)*pi/2+y = x where |y|<pi/4 and n indicates in
;; which octant the arg lies. Y is actually computed in two parts,
;; y[0] and y[1] such that the sum is y, for accuracy.
(alien:def-alien-routine ("__ieee754_rem_pio2" %%ieee754-rem-pi/2) c-call:int
(x double-float)
(y (* double-float)))
;; Same as above, but instead of needing to pass an array in, the ;; Same as above, but instead of needing to pass an array in, the
;; output array is broken up into two output values instead. This is ;; output array is broken up into two output values instead. This is
;; easier for the user, and we don't have to wrap calls with ;; easier for the user, and we don't have to wrap calls with
...@@ -221,148 +197,438 @@ ...@@ -221,148 +197,438 @@
(y0 double-float :out) (y0 double-float :out)
(y1 double-float :out)) (y1 double-float :out))
) ;; Implement sin/cos/tan in Lisp. These are based on the routines
;; from fdlibm.
;; If the C library is accurate, use %trig as the Lisp name. ;; Block compile so the trig routines don't cons their args when
#-(or ppc (and sse2 (not darwin))) ;; calling the kernel trig routines.
(progn (declaim (ext:start-block kernel-sin kernel-cos kernel-tan
(declaim (inline %sin %cos %tan)) %sin %cos %tan
(macrolet ((frob (alien-name lisp-name) %sincos))
`(alien:def-alien-routine (,alien-name ,lisp-name) double-float
(x double-float))))
(frob "sin" %sin)
(frob "cos" %cos)
(frob "tan" %tan))
)
;; Make %%trig be the C library routines that don't do accurate
;; reduction. This is for PPC and for any SSE2 build except on
;; Darwin. Darwin has accurate C library routines.
#+(or ppc (and sse2 (not darwin)))
(progn
(declaim (inline %%sin %%cos %%tan))
(macrolet ((frob (alien-name lisp-name)
`(alien:def-alien-routine (,alien-name ,lisp-name) double-float
(x double-float))))
(frob "sin" %%sin)
(frob "cos" %%cos)
(frob "tan" %%tan))
)
;; When the C library is not accurate, define %trig to do accurate
;; argument reduction and call the appropriate C function on the
;; reduced arg. For x87, we can use the x87 FPU trig instructions.
#+(or ppc (and x86 (not darwin)))
(macrolet
((frob (sin cos tan)
`(progn
;; In all of the routines below, we just compute the sum of
;; y0 and y1 and use that as the (reduced) argument for the
;; trig functions. This is slightly less accurate than what
;; fdlibm does, which calls special functions using y0 and
;; y1 separately, for greater accuracy. This isn't
;; implemented, and some spot checks indicate that what we
;; have here is accurate.
;;
;; For x86 with an fsin/fcos/fptan instruction, the pi/4 is
;; probably too restrictive.
(defun %sin (x)
(declare (double-float x))
(if (< (abs x) (/ pi 4))
(,sin x)
;; Argument reduction needed
(multiple-value-bind (n y0 y1)
(%ieee754-rem-pi/2 x)
(let ((reduced (+ y0 y1)))
(case (logand n 3)
(0 (,sin reduced))
(1 (,cos reduced))
(2 (- (,sin reduced)))
(3 (- (,cos reduced))))))))
(defun %cos (x)
(declare (double-float x))
(if (< (abs x) (/ pi 4))
(,cos x)
;; Argument reduction needed
(multiple-value-bind (n y0 y1)
(%ieee754-rem-pi/2 x)
(let ((reduced (+ y0 y1)))
(case (logand n 3)
(0 (,cos reduced))
(1 (- (,sin reduced)))
(2 (- (,cos reduced)))
(3 (,sin reduced)))))))
(defun %tan (x)
(declare (double-float x))
(if (< (abs x) (/ pi 4))
(,tan x)
;; Argument reduction needed
(multiple-value-bind (n y0 y1)
(%ieee754-rem-pi/2 x)
(let ((reduced (+ y0 y1)))
(if (evenp n)
(,tan reduced)
(- (/ (,tan reduced)))))))))))
;; Don't want %sin-quick and friends with sse2.
#+(and x86 (not sse2))
(frob %sin-quick %cos-quick %tan-quick)
#+(or ppc sse2)
(frob %%sin %%cos %%tan))
;; Linux and sparc have a sincos function in the C library. Use it.
;; But on linux we need to do pi reduction ourselves because the C
;; library doesn't do accurate reduction. Sparc does accurate pi
;; reduction, so we don't need to do it ourselves.
#+(or (and linux x86) sparc)
(progn
(declaim (inline %%sincos))
(export '%%sincos)
(alien:def-alien-routine ("sincos" %%sincos) c-call:void
(x double-float)
(sin double-float :out)
(cos double-float :out))
#+(and linux x86)
(defun %sincos (theta)
(declare (double-float theta))
;; Accurately reduce theta.
(multiple-value-bind (n y0 y1)
(%ieee754-rem-pi/2 theta)
(multiple-value-bind (ignore s c)
(%%sincos y0)
(declare (ignore ignore))
;; Figure out which quadrant to use, and finish out the
;; computation using y1. This is done by using a 1st-order
;; Taylor expansion about y0.
(flet ((sin2 (s c y)
;; sin(x+y) = sin(x) + cos(x)*y
(+ s (* c y)))
(cos2 (s c y)
;; cos(x+y) = cos(x) - sin(x)*y
(- c (* s y))))
(case (logand n 3)
(0
(values (sin2 s c y1)
(cos2 s c y1)))
(1
(values (cos2 s c y1)
(- (sin2 s c y1))))
(2
(values (- (sin2 s c y1))
(- (cos2 s c y1))))
(3
(values (- (cos2 s c y1))
(sin2 s c y1))))))))
#+sparc
(declaim (inline %sinccos))
#+sparc
(defun %sincos (theta)
(multiple-value-bind (ignore s c)
(%%sincos theta)
(declare (ignore ignore))
(values s c)))
)
;; kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
;; Input x is assumed to be bounded by ~pi/4 in magnitude.
;; Input y is the tail of x.
;; Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
;;
;; Algorithm
;; 1. Since sin(-x) = -sin(x), we need only to consider positive x.
;; 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
;; 3. sin(x) is approximated by a polynomial of degree 13 on
;; [0,pi/4]
;; 3 13
;; sin(x) ~ x + S1*x + ... + S6*x
;; where
;;
;; |sin(x) 2 4 6 8 10 12 | -58
;; |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2
;; | x |
;;
;; 4. sin(x+y) = sin(x) + sin'(x')*y
;; ~ sin(x) + (1-x*x/2)*y
;; For better accuracy, let
;; 3 2 2 2 2
;; r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
;; then 3 2
;; sin(x) = x + (S1*x + (x *(r-y/2)+y))
(declaim (ftype (function (double-float double-float fixnum)
double-float)
kernel-sin))
(defun kernel-sin (x y iy)
(declare (type (double-float -1d0 1d0) x y)
(fixnum iy)
(optimize (speed 3) (safety 0)))
(let ((ix (ldb (byte 31 0) (kernel:double-float-high-bits x))))
(when (< ix #x3e400000)
;; |x| < 2^-27
;; Signal inexact if x /= 0
(if (zerop (truncate x))
(return-from kernel-sin x)
(return-from kernel-sin x)))
(let* ((s1 -1.66666666666666324348d-01) ; #xBFC55555 #x55555549
(s2 8.33333333332248946124d-03) ; #x3F811111 #x1110F8A6
(s3 -1.98412698298579493134d-04) ; #xBF2A01A0 #x19C161D5
(s4 2.75573137070700676789d-06) ; #x3EC71DE3 #x57B1FE7D
(s5 -2.50507602534068634195d-08) ; #xBE5AE5E6 #x8A2B9CEB
(s6 1.58969099521155010221d-10) ; #x3DE5D93A #x5ACFD57C
(z (* x x))
(v (* z x))
(r (+ s2
(* z
(+ s3
(* z
(+ s4
(* z
(+ s5
(* z s6))))))))))
(if (zerop iy)
(+ x (* v (+ s1 (* z r))))
(- x (- (- (* z (- (* .5 y)
(* v r)))
y)
(* v s1)))))))
;; kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
;; Input x is assumed to be bounded by ~pi/4 in magnitude.
;; Input y is the tail of x.
;;
;; Algorithm
;; 1. Since cos(-x) = cos(x), we need only to consider positive x.
;; 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
;; 3. cos(x) is approximated by a polynomial of degree 14 on
;; [0,pi/4]
;; 4 14
;; cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
;; where the remez error is
;;
;; | 2 4 6 8 10 12 14 | -58
;; |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
;; | |
;;
;; 4 6 8 10 12 14
;; 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
;; cos(x) = 1 - x*x/2 + r
;; since cos(x+y) ~ cos(x) - sin(x)*y
;; ~ cos(x) - x*y,
;; a correction term is necessary in cos(x) and hence
;; cos(x+y) = 1 - (x*x/2 - (r - x*y))
;; For better accuracy when x > 0.3, let qx = |x|/4 with
;; the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
;; Then
;; cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
;; Note that 1-qx and (x*x/2-qx) is EXACT here, and the
;; magnitude of the latter is at least a quarter of x*x/2,
;; thus, reducing the rounding error in the subtraction.
(declaim (ftype (function (double-float double-float)
double-float)
kernel-cos))
(defun kernel-cos (x y)
(declare (type (double-float -1d0 1d0) x y)
(optimize (speed 3) (safety 0)))
;; cos(-x) = cos(x), so we just compute cos(|x|).
(let ((ix (ldb (byte 31 0) (kernel:double-float-high-bits x))))
;; cos(x) = 1 when |x| < 2^-27
(when (< ix #x3e400000)
;; Signal inexact if x /= 0
(if (zerop (truncate x))
(return-from kernel-cos 1d0)
(return-from kernel-cos 1d0)))
(let* ((c1 4.16666666666666019037d-02)
(c2 -1.38888888888741095749d-03)
(c3 2.48015872894767294178d-05)
(c4 -2.75573143513906633035d-07)
(c5 2.08757232129817482790d-09)
(c6 -1.13596475577881948265d-11)
(z (* x x))
(r (* z
(+ c1
(* z
(+ c2
(* z
(+ c3
(* z
(+ c4
(* z
(+ c5
(* z c6)))))))))))))
(cond ((< ix #x3fd33333)
;; \x| < 0.3
(- 1 (- (* .5 z)
(- (* z r)
(* x y)))))
(t
;; qx = 0.28125 if |x| > 0.78125, else x/4 dropping the
;; least significant 32 bits.
(let* ((qx (if (> ix #x3fe90000)
0.28125d0
;; x/4, exactly, and also dropping the
;; least significant 32 bits of the
;; fraction.
(make-double-float (- ix #x00200000)
0)))
(hz (- (* 0.5 z) qx))
(a (- 1 qx)))
(- a (- hz (- (* z r)
(* x y))))))))))
(declaim (type (simple-array double-float (*)) tan-coef))
(defconstant tan-coef
(make-array 13 :element-type 'double-float
:initial-contents
'(3.33333333333334091986d-01
1.33333333333201242699d-01
5.39682539762260521377d-02
2.18694882948595424599d-02
8.86323982359930005737d-03
3.59207910759131235356d-03
1.45620945432529025516d-03
5.88041240820264096874d-04
2.46463134818469906812d-04
7.81794442939557092300d-05
7.14072491382608190305d-05
-1.85586374855275456654d-05
2.59073051863633712884d-05)))
;; kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
;; Input x is assumed to be bounded by ~pi/4 in magnitude.
;; Input y is the tail of x.
;; Input k indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned.
;;
;; Algorithm
;; 1. Since tan(-x) = -tan(x), we need only to consider positive x.
;; 2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
;; 3. tan(x) is approximated by a odd polynomial of degree 27 on
;; [0,0.67434]
;; 3 27
;; tan(x) ~ x + T1*x + ... + T13*x
;; where
;;
;; |tan(x) 2 4 26 | -59.2
;; |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2
;; | x |
;;
;; Note: tan(x+y) = tan(x) + tan'(x)*y
;; ~ tan(x) + (1+x*x)*y
;; Therefore, for better accuracy in computing tan(x+y), let
;; 3 2 2 2 2
;; r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
;; then
;; 3 2
;; tan(x+y) = x + (T1*x + (x *(r+y)+y))
;;
;; 4. For x in [0.67434,pi/4], let y = pi/4 - x, then
;; tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
;; = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
(declaim (ftype (function (double-float double-float fixnum)
double-float)
kernel-tan))
(defun kernel-tan (x y iy)
(declare (type (double-float -1d0 1d0) x y)
(type (member -1 1) iy)
(optimize (speed 3) (safety 0)))
(let* ((hx (kernel:double-float-high-bits x))
(ix (logand hx #x7fffffff))
(w 0d0)
(z 0d0)
(v 0d0)
(s 0d0)
(r 0d0))
(declare (double-float w z v s r))
(when (< ix #x3e300000)
;; |x| < 2^-28
(when (zerop (truncate x))
(cond ((zerop (logior (logior ix (kernel:double-float-low-bits x))
(+ iy 1)))
;; x = 0 (because hi and low bits are 0) and iy = -1
;; (cot)
(return-from kernel-tan (/ (abs x))))
((= iy 1)
(return-from kernel-tan x))
(t
;; x /= 0 and iy = -1 (cot)
;; Compute -1/(x+y) carefully
(let ((a 0d0)
(tt 0d0))
(setf w (+ x y))
(setf z (make-double-float (double-float-high-bits w) 0))
(setf v (- y (- z x)))
(setf a (/ -1 w))
(setf tt (make-double-float (double-float-high-bits a) 0))
(setf s (+ 1 (* tt z)))
(return-from kernel-tan (+ tt
(* a (+ s (* tt v))))))))))
(when (>= ix #x3FE59428)
;; |x| > .6744
(when (minusp hx)
(setf x (- x))
(setf y (- y)))
;; The two constants below are such that pi/4 + pi/4_lo is pi/4
;; to twice the accuracy of a double float.
;;
;; z = pi/4-x
(setf z (- (make-double-float #x3FE921FB #x54442D18) x))
;; w = pi/4_lo - y.
(setf w (- (make-double-float #x3C81A626 #x33145C07) y))
(setf x (+ z w))
(setf y 0d0))
(setf z (* x x))
(setf w (* z z))
;; Break x^5*(T[1]+x^2*T[2]+...) into
;; x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
;; x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
(setf r (+ (aref tan-coef 1)
(* w
(+ (aref tan-coef 3)
(* w
(+ (aref tan-coef 5)
(* w
(+ (aref tan-coef 7)
(* w
(+ (aref tan-coef 9)
(* w (aref tan-coef 11))))))))))))
(setf v (* z
(+ (aref tan-coef 2)
(* w
(+ (aref tan-coef 4)
(* w
(+ (aref tan-coef 6)
(* w
(+ (aref tan-coef 8)
(* w
(+ (aref tan-coef 10)
(* w (aref tan-coef 12)))))))))))))
(setf s (* z x))
(setf r (+ y (* z (+ (* s (+ r v))
y))))
(incf r (* s (aref tan-coef 0)))
(setf w (+ x r))
(when (>= ix #x3FE59428)
(let ((v (float iy 1d0)))
(return-from kernel-tan
(* (- 1 (logand 2 (ash hx -30)))
(- v
(* 2
(- x (- (/ (* w w)
(+ w v))
r))))))))
(when (= iy 1)
(return-from kernel-tan w))
;; Compute 1/w=1/(x+r) carefully
(let ((a 0d0)
(tt 0d0))
(setf z (kernel:make-double-float (kernel:double-float-high-bits w) 0))
(setf v (- r (- z x))) ; z + v = r + x
(setf a (/ -1 w))
(setf tt (kernel:make-double-float (kernel:double-float-high-bits a) 0))
(setf s (+ 1 (* tt z)))
(+ tt
(* a
(+ s (* tt v)))))))
;; Return sine function of x.
;;
;; kernel function:
;; __kernel_sin ... sine function on [-pi/4,pi/4]
;; __kernel_cos ... cose function on [-pi/4,pi/4]
;; __ieee754_rem_pio2 ... argument reduction routine
;;
;; Method.
;; Let S,C and T denote the sin, cos and tan respectively on
;; [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
;; in [-pi/4 , +pi/4], and let n = k mod 4.
;; We have
;;
;; n sin(x) cos(x) tan(x)
;; ----------------------------------------------------------
;; 0 S C T
;; 1 C -S -1/T
;; 2 -S -C T
;; 3 -C S -1/T
;; ----------------------------------------------------------
;;
;; Special cases:
;; Let trig be any of sin, cos, or tan.
;; trig(+-INF) is NaN, with signals;
;; trig(NaN) is that NaN;
;;
;; Accuracy:
;; TRIG(x) returns trig(x) nearly rounded
(defun %sin (x)
(declare (double-float x)
(optimize (speed 3)))
(let ((ix (ldb (byte 31 0) (kernel:double-float-high-bits x))))
(cond
((<= ix #x3fe921fb)
;; |x| < pi/4, approx
(kernel-sin x 0d0 0))
((>= ix #x7ff00000)
;; sin(Inf or NaN) is NaN
(- x x))
(t
;; Argument reduction needed
(multiple-value-bind (n y0 y1)
(%ieee754-rem-pi/2 x)
(case (logand n 3)
(0
(kernel-sin y0 y1 1))
(1
(kernel-cos y0 y1))
(2
(- (kernel-sin y0 y1 1)))
(3
(- (kernel-cos y0 y1)))))))))
(defun %cos (x)
(declare (double-float x)
(optimize (speed 3)))
(let ((ix (ldb (byte 31 0) (kernel:double-float-high-bits x))))
(cond
((< ix #x3fe921fb)
;;|x| < pi/4, approx
(kernel-cos x 0d0))
((>= ix #x7ff00000)
;; cos(Inf or NaN) is NaN
(- x x))
(t
;; Argument reduction needed
(multiple-value-bind (n y0 y1)
(%ieee754-rem-pi/2 x)
(ecase (logand n 3)
(0
(kernel-cos y0 y1))
(1
(- (kernel-sin y0 y1 1)))
(2
(- (kernel-cos y0 y1)))
(3
(kernel-sin y0 y1 1))))))))
(defun %tan (x)
(declare (double-float x)
(optimize (speed 3)))
(let ((ix (logand #x7fffffff (kernel:double-float-high-bits x))))
(cond ((<= ix #x3fe921fb)
;; |x| < pi/4
(kernel-tan x 0d0 1))
((>= ix #x7ff00000)
;; tan(Inf or Nan) is NaN
(- x x))
(t
(multiple-value-bind (n y0 y1)
(%ieee754-rem-pi/2 x)
(let ((flag (- 1 (ash (logand n 1) 1))))
;; flag = 1 if n even, -1 if n odd
(kernel-tan y0 y1 flag)))))))
;; Compute sin and cos of x, simultaneously.
(defun %sincos (x)
(declare (double-float x)
(optimize (speed 3)))
(cond ((<= (abs x) (/ pi 4))
(values (kernel-sin x 0d0 0)
(kernel-cos x 0d0)))
(t
;; Argument reduction needed
(multiple-value-bind (n y0 y1)
(%ieee754-rem-pi/2 x)
(case (logand n 3)
(0
(values (kernel-sin y0 y1 1)
(kernel-cos y0 y1)))
(1
(values (kernel-cos y0 y1)
(- (kernel-sin y0 y1 1))))
(2
(values (- (kernel-sin y0 y1 1))
(- (kernel-cos y0 y1))))
(3
(values (- (kernel-cos y0 y1))
(kernel-sin y0 y1 1))))))))
(declaim (ext:end-block))
;;;; Power functions. ;;;; Power functions.
...@@ -1018,9 +1284,6 @@ ...@@ -1018,9 +1284,6 @@
"Return cos(Theta) + i sin(Theta), AKA exp(i Theta)." "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
(if (complexp theta) (if (complexp theta)
(error (intl:gettext "Argument to CIS is complex: ~S") theta) (error (intl:gettext "Argument to CIS is complex: ~S") theta)
#-(or (