This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "CMU Common Lisp". The branch, rtoy-lisp-trig has been created at 7190b61cf97c8320d6a218c430471c0fb0bf518e (commit) - Log ----------------------------------------------------------------- commit 7190b61cf97c8320d6a218c430471c0fb0bf518e Author: Raymond Toy <toy.raymond@gmail.com> Date: Sat Dec 14 21:21:50 2013 -0800 Add test for sincos(-0d0). diff --git a/src/tests/trig.lisp b/src/tests/trig.lisp index a0c8c15..9555b11 100644 --- a/src/tests/trig.lisp +++ b/src/tests/trig.lisp @@ -142,6 +142,10 @@ -4.08066388841804238545143494525595117765084022768d-1) +(rt:deftest sincos.0 + (multiple-value-list (kernel::%sincos -0d0)) + (-0d0 1d0)) + (rt:deftest sincos.1 (let (results) (dotimes (k 1000) commit b79c28727f40f1dd3cdd035eb86fd929594b0d64 Author: Raymond Toy <toy.raymond@gmail.com> Date: Sat Dec 14 21:20:00 2013 -0800 Implement sincos using the new Lisp trig routines. This can now be used for all platforms. code/irrat.lisp:: * Implement %SINCOS compiler/float-tran.lisp:: * Update deftransforms for CIS. %SINCOS can be used on any platform. tests/trig.lisp: * Add tests to verify %sincos returns exactly the same values as for sin and cos. diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 660c519..c23321d 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -202,7 +202,9 @@ ;; Block compile so the trig routines don't cons their args when ;; calling the kernel trig routines. -(declaim (ext:start-block kernel-sin kernel-cos kernel-tan %sin %cos %tan)) +(declaim (ext:start-block kernel-sin kernel-cos kernel-tan + %sin %cos %tan + %sincos)) ;; kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854 ;; Input x is assumed to be bounded by ~pi/4 in magnitude. @@ -592,63 +594,32 @@ ;; flag = 1 if n even, -1 if n odd (kernel-tan y0 y1 flag))))))) +(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)) -;; 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))) -) - -  ;;;; Power functions. @@ -1303,9 +1274,6 @@ "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)." (if (complexp theta) (error (intl:gettext "Argument to CIS is complex: ~S") theta) - #-(or (and linux x86) sparc) - (complex (cos theta) (sin theta)) - #+(or (and linux x86) sparc) (number-dispatch ((theta real)) ((rational) (let ((arg (coerce theta 'double-float))) diff --git a/src/compiler/float-tran.lisp b/src/compiler/float-tran.lisp index a8147d9..d123d18 100644 --- a/src/compiler/float-tran.lisp +++ b/src/compiler/float-tran.lisp @@ -731,8 +731,6 @@ (deftransform name ((x) '(double-float) rtype :eval-name t :when :both) `(,prim x)))) -#+(or (and linux x86) sparc) -(progn (defknown (kernel::%sincos) (double-float) (values double-float double-float) (movable foldable flushable)) @@ -752,7 +750,7 @@ (deftransform cis ((z) (double-double-float) *) ;; Cis. '(complex (cos z) (sin z))) -) + ;;; The argument range is limited on the x86 FP trig. functions. A ;;; post-test can detect a failure (and load a suitable result), but diff --git a/src/tests/trig.lisp b/src/tests/trig.lisp index 49f16fb..a0c8c15 100644 --- a/src/tests/trig.lisp +++ b/src/tests/trig.lisp @@ -140,4 +140,89 @@ (rt:deftest tan.misc.1 (tan (scale-float 1d0 120)) -4.08066388841804238545143494525595117765084022768d-1) - + + +(rt:deftest sincos.1 + (let (results) + (dotimes (k 1000) + (let* ((x (random (/ pi 4))) + (s-exp (sin x)) + (c-exp (cos x))) + (multiple-value-bind (s c) + (kernel::%sincos x) + (unless (and (= s s-exp) + (= c c-exp)) + (push (list x + (list s s-exp) + (list c c-exp)) + results))))) + results) + nil) + +(rt:deftest sincos.2 + (let (results) + (dotimes (k 1000) + (let* ((x (random 16d0)) + (s-exp (sin x)) + (c-exp (cos x))) + (multiple-value-bind (s c) + (kernel::%sincos x) + (unless (and (= s s-exp) + (= c c-exp)) + (push (list x + (list s s-exp) + (list c c-exp)) + results))))) + results) + nil) + +(rt:deftest sincos.3 + (let (results) + (dotimes (k 1000) + (let* ((x (random (scale-float 1d0 120))) + (s-exp (sin x)) + (c-exp (cos x))) + (multiple-value-bind (s c) + (kernel::%sincos x) + (unless (and (= s s-exp) + (= c c-exp)) + (push (list x + (list s s-exp) + (list c c-exp)) + results))))) + results) + nil) + +(rt:deftest sincos.3a + (let (results) + (dotimes (k 1000) + (let* ((x (- (random (scale-float 1d0 120)))) + (s-exp (sin x)) + (c-exp (cos x))) + (multiple-value-bind (s c) + (kernel::%sincos x) + (unless (and (= s s-exp) + (= c c-exp)) + (push (list x + (list s s-exp) + (list c c-exp)) + results))))) + results) + nil) + +(rt:deftest sincos.4 + (let (results) + (dotimes (k 1000) + (let* ((x (random (scale-float 1d0 1023))) + (s-exp (sin x)) + (c-exp (cos x))) + (multiple-value-bind (s c) + (kernel::%sincos x) + (unless (and (= s s-exp) + (= c c-exp)) + (push (list x + (list s s-exp) + (list c c-exp)) + results))))) + results) + nil) commit e6a9577f0093b72d5d5e0c90cb0930df6a16bb8b Author: Raymond Toy <toy.raymond@gmail.com> Date: Sat Dec 14 20:35:56 2013 -0800 Implement trig functions in Lisp code/irrat.lisp:: * Add Lisp implementation for sin, cos, and tan, based on code from fdlibm. Requires the C reduction routines. Only working so far on systems that already include the reduction routies. tests/trig.lisp:: * Tests for the new sin, cos, and tan functions. Tests pass on x86/darwin. diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 270f1dc..660c519 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -187,30 +187,6 @@ (%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 ;; 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 @@ -221,93 +197,402 @@ (y0 double-float :out) (y1 double-float :out)) -) - -;; If the C library is accurate, use %trig as the Lisp name. -#-(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)) -) +;; Implement sin/cos/tan in Lisp. These are based on the routines +;; from fdlibm. -;; 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)) -) +;; Block compile so the trig routines don't cons their args when +;; calling the kernel trig routines. +(declaim (ext:start-block kernel-sin kernel-cos kernel-tan %sin %cos %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)) +;; 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) + (if (zerop (truncate x)) + (return-from kernel-sin x) + (return-from kernel-sin x))) + (let* ((s1 -1.66666666666666324348d-01) + (s2 8.33333333332248946124d-03) + (s3 -1.98412698298579493134d-04) + (s4 2.75573137070700676789d-06) + (s5 -2.50507602534068634195d-08) + (s6 1.58969099521155010221d-10) + (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 + (let* ((qx (if (> ix #x3fe90000) + 0.28125d0 + ;; x/4, exactly, and also dropping the + ;; least significant 32 bits of the + ;; fraction. (Why?) + (kernel: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 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 (kernel:make-double-float (kernel:double-float-high-bits w) 0)) + (setf v (- y (- z x))) + (setf a (/ -1 w)) + (setf tt (kernel:make-double-float (kernel: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))) + ;; z = pi/4-x + (setf z (- (kernel:make-double-float #x3FE921FB #x54442D18) x)) + ;; w = pi/4_lo - y + (setf w (- (kernel: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)) + ;; + (let ((a 0d0) + (tt 0d0)) + (setf z (kernel:make-double-float (kernel:double-float-high-bits w) 0)) + (setf v (- r (- r 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) + (kernel::%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) + (kernel::%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) + (kernel-tan x 0d0 1)) + ((>= ix #x7ff00000) + (- x x)) + (t + (multiple-value-bind (n y0 y1) + (kernel::%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))))))) + +(declaim (ext:end-block)) ;; 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 diff --git a/src/tests/trig.lisp b/src/tests/trig.lisp new file mode 100644 index 0000000..49f16fb --- /dev/null +++ b/src/tests/trig.lisp @@ -0,0 +1,143 @@ +(rt:deftest sin.1 + (sin 0d0) + 0d0) + +(rt:deftest sin.2 + (sin -0d0) + -0d0) + +(rt:deftest sin.3 + ;; Tests the case for |x| < 2^-27, but not 0. + (sin (scale-float 1d0 -28)) + #.(scale-float 1d0 -28)) + +(rt:deftest sin.4 + ;; Just a random test, without argument reduction + (sin .5d0) + 0.479425538604203d0) + +(rt:deftest sin.5 + ;; Test for arg near pi/2 + (sin (/ pi 2)) + 1d0) + +(rt:deftest sin.red.0 + ;; Test for argument reduction with n mod 4 = 0 + (sin (* 7/4 pi)) + -7.07106781186547675943154203316156531867416581156d-1) + +(rt:deftest sin.red.1 + ;; Test for argument reduction with n mod 4 = 1 + (sin (* 9/4 pi)) + 7.07106781186547329560731709118834541043171055432d-1) + +(rt:deftest sin.red.2 + ;; Test for argument reduction with n mod 4 = 2 + (sin (* 11/4 pi)) + 7.07106781186548390575743300374993861263439430213d-1) + +(rt:deftest sin.red.3 + ;; Test for argument reduction with n mod 4 = 3 + (sin (* 13/4 pi)) + -7.07106781186547871002109559079472349116005337743d-1) + +(rt:deftest sin.misc.1 + ;; Test for argument reduction + (sin (scale-float 1d0 120)) + 0.377820109360752d0) + +(rt:deftest cos.1 + (cos 0d0) + 1d0) + +(rt:deftest cos.2 + (cos -0d0) + 1d0) + +(rt:deftest cos.3 + ;; Test for |x| < 2^-27 + (cos (scale-float 1d0 -28)) + 1d0) + +(rt:deftest cos.4 + ;; Test for branch |x| < .3 + (cos 0.25d0) + 0.9689124217106447d0) + +(rt:deftest cos.5 + ;; Test for branch |x| > .3 and \x| < .78125 + (cos 0.5d0) + 8.7758256189037271611628158260382965199164519711d-1) + +(rt:deftest cos.6 + ;; Test for branch |x| > .3 and |x| > .78125 + (cos 0.785d0) + 0.7073882691671998d0) + +(rt:deftest cos.7 + ;; Random test near pi/2 + (cos (/ pi 2)) + 6.123233995736766d-17) + +(rt:deftest cos.misc.1 + ;; Test for argument reduction + (cos (scale-float 1d0 120)) + -0.9258790228548379d0) + +(rt:deftest cos.red.0 + ;; Test for argument reduction with n mod 4 = 0 + (cos (* 7/4 pi)) + 7.07106781186547372858534520893509069186435867941d-1) + +(rt:deftest cos.red.1 + ;; Test for argument reduction with n mod 4 = 1 + (cos (* 9/4 pi)) + 7.0710678118654771924095701509080985020443197242d-1) + +(rt:deftest cos.red.2 + ;; Test for argument reduction with n mod 4 = 2 + (cos (* 11/4 pi)) + -7.07106781186546658225945423833643190916000739026d-1) + +(rt:deftest cos.red.3 + ;; Test for argument reduction with n mod 4 = 3 + (cos (* 13/4 pi)) + -7.07106781186547177799579165130055836531929091466d-1) + +(rt:deftest tan.1 + (tan 0d0) + 0d0) + +(rt:deftest tan.2 + (tan -0d0) + -0d0) + +(rt:deftest tan.3 + ;; |x| < 2^-28 + (tan (scale-float 1d0 -29)) + #.(scale-float 1d0 -29)) + +(rt:deftest tan.4 + ;; |x| < .6744 + (tan 0.5d0) + 5.4630248984379051325517946578028538329755172018d-1) + +(rt:deftest tan.5 + ;; |x = 11/16 = 0.6875 > .6744 + (tan (float 11/16 1d0)) + 8.21141801589894121911423965374711700875371645309d-1) + +(rt:deftest tan.red.0 + ;; Test for argument reduction with n even + (tan (* 7/4 pi)) + -1.00000000000000042862637970157370388940976433505d0) + +(rt:deftest tan.red.1 + ;; Test for argument reduction with n odd + (tan (* 9/4 pi)) + 9.99999999999999448908940383691222098948324989275d-1) + +(rt:deftest tan.misc.1 + (tan (scale-float 1d0 120)) + -4.08066388841804238545143494525595117765084022768d-1) + commit 32bdd53bf002fca1c9ad6b543d522a6558cae768 Author: Raymond Toy <toy.raymond@gmail.com> Date: Sat Dec 14 20:22:40 2013 -0800 Add RT. src/contrib/rt:: * Add RT code, including asdf. src/code/module.lisp:: * Add RT as a module diff --git a/src/code/module.lisp b/src/code/module.lisp index 70ccba7..42b0ac2 100644 --- a/src/code/module.lisp +++ b/src/code/module.lisp @@ -148,6 +148,12 @@ (defmodule "asdf" "modules:asdf/asdf") +(defmodule :rt + "modules:rt/rt") + +(defmodule "rt" + "modules:rt/rt") + ;; Allow user to specify "cmu-contribs" or :cmu-contribs. (defmodule "cmu-contribs" "modules:contrib") diff --git a/src/contrib/rt/rt.asd b/src/contrib/rt/rt.asd new file mode 100644 index 0000000..718e965 --- /dev/null +++ b/src/contrib/rt/rt.asd @@ -0,0 +1,33 @@ +;;;; -*- Mode: LISP; Syntax: ANSI-Common-Lisp; Base: 10 -*- +;;;; ************************************************************************* +;;;; FILE IDENTIFICATION +;;;; +;;;; Name: rt.asd +;;;; Purpose: ASDF definition file for Rt +;;;; Programmer: Kevin M. Rosenberg +;;;; Date Started: Sep 2002 +;;;; +;;;; $Id$ +;;;; +;;;; This file, part of cl-rt, is Copyright (c) 2002 by Kevin M. Rosenberg +;;;; +;;;; cl-rt users are granted the rights to distribute and use this software +;;;; as governed by the terms of the GNU Lesser General Public License +;;;; (http://www.gnu.org/licenses/lgpl.html) +;;;; ************************************************************************* + +(in-package :asdf) + +(defsystem :rt + :name "cl-rt" + :version "1990.12.19" + :maintainer "Kevin M. Rosenberg <kmr@debian.org>" + :licence "MIT" + :description "MIT Regression Tester" + :long-description "RT provides a framework for writing regression test suites" + :perform (load-op :after (op rt) + (pushnew :rt cl:*features*)) + :components + ((:file "rt"))) + + diff --git a/src/contrib/rt/rt.lisp b/src/contrib/rt/rt.lisp new file mode 100644 index 0000000..3df87c4 --- /dev/null +++ b/src/contrib/rt/rt.lisp @@ -0,0 +1,409 @@ +;-*-syntax:COMMON-LISP;Package:(RT :use "COMMON-LISP" :colon-mode :external)-*- + +#|----------------------------------------------------------------------------| + | Copyright 1990 by the Massachusetts Institute of Technology, Cambridge MA. | + | | + | Permission to use, copy, modify, and distribute this software and its | + | documentation for any purpose and without fee is hereby granted, provided | + | that this copyright and permission notice appear in all copies and | + | supporting documentation, and that the name of M.I.T. not be used in | + | advertising or publicity pertaining to distribution of the software | + | without specific, written prior permission. M.I.T. makes no | + | representations about the suitability of this software for any purpose. | + | It is provided "as is" without express or implied warranty. | + | | + | M.I.T. DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING | + | ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL | + | M.I.T. BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR | + | ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, | + | WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, | + | ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS | + | SOFTWARE. | + |----------------------------------------------------------------------------|# + +(defpackage #:regression-test + (:nicknames #:rtest #-lispworks #:rt) + (:use #:cl) + (:export #:*do-tests-when-defined* #:*test* #:continue-testing + #:deftest #:do-test #:do-tests #:get-test #:pending-tests + #:rem-all-tests #:rem-test) + (:documentation "The MIT regression tester with pfdietz's modifications")) + +;;This was the December 19, 1990 version of the regression tester, but +;;has since been modified. + +(in-package :regression-test) + +(declaim (ftype (function (t) t) get-entry expanded-eval do-entries)) +(declaim (type list *entries*)) +(declaim (ftype (function (t &rest t) t) report-error)) +(declaim (ftype (function (t &optional t) t) do-entry)) + +(defvar *test* nil "Current test name") +(defvar *do-tests-when-defined* nil) +(defvar *entries* '(nil) "Test database. Has a leading dummy cell that does not contain an entry.") +(defvar *entries-tail* *entries* "Tail of the *entries* list") +(defvar *entries-table* (make-hash-table :test #'equal) + "Map the names of entries to the cons cell in *entries* that precedes the one whose car is the entry.") +(defvar *in-test* nil "Used by TEST") +(defvar *debug* nil "For debugging") +(defvar *catch-errors* t "When true, causes errors in a test to be caught.") +(defvar *print-circle-on-failure* nil + "Failure reports are printed with *PRINT-CIRCLE* bound to this value.") + +(defvar *compile-tests* nil "When true, compile the tests before running them.") +(defvar *expanded-eval* nil "When true, convert the tests into a form that is less likely to have compiler optimizations.") +(defvar *optimization-settings* '((safety 3))) + +(defvar *expected-failures* nil + "A list of test names that are expected to fail.") + +(defvar *notes* (make-hash-table :test 'equal) + "A mapping from names of notes to note objects.") + +(defstruct (entry (:conc-name nil)) + pend name props form vals) + +;;; Note objects are used to attach information to tests. +;;; A typical use is to mark tests that depend on a particular +;;; part of a set of requirements, or a particular interpretation +;;; of the requirements. + +(defstruct note + name + contents + disabled ;; When true, tests with this note are considered inactive + ) + +;; (defmacro vals (entry) `(cdddr ,entry)) + +(defmacro defn (entry) + (let ((var (gensym))) + `(let ((,var ,entry)) + (list* (name ,var) (form ,var) (vals ,var))))) + +(defun entry-notes (entry) + (let* ((props (props entry)) + (notes (getf props :notes))) + (if (listp notes) + notes + (list notes)))) + +(defun has-disabled-note (entry) + (let ((notes (entry-notes entry))) + (loop for n in notes + for note = (if (note-p n) n + (gethash n *notes*)) + thereis (and note (note-disabled note))))) + +(defun pending-tests () + (loop for entry in (cdr *entries*) + when (and (pend entry) (not (has-disabled-note entry))) + collect (name entry))) + +(defun rem-all-tests () + (setq *entries* (list nil)) + (setq *entries-tail* *entries*) + (clrhash *entries-table*) + nil) + +(defun rem-test (&optional (name *test*)) + (let ((pred (gethash name *entries-table*))) + (when pred + (if (null (cddr pred)) + (setq *entries-tail* pred) + (setf (gethash (name (caddr pred)) *entries-table*) pred)) + (setf (cdr pred) (cddr pred)) + (remhash name *entries-table*) + name))) + +(defun get-test (&optional (name *test*)) + (defn (get-entry name))) + +(defun get-entry (name) + (let ((entry ;; (find name (the list (cdr *entries*)) + ;; :key #'name :test #'equal) + (cadr (gethash name *entries-table*)) + )) + (when (null entry) + (report-error t + "~%No test with name ~:@(~S~)." + name)) + entry)) + +(defmacro deftest (name &rest body) + (let* ((p body) + (properties + (loop while (keywordp (first p)) + unless (cadr p) + do (error "Poorly formed deftest: ~A~%" + (list* 'deftest name body)) + append (list (pop p) (pop p)))) + (form (pop p)) + (vals p)) + `(add-entry (make-entry :pend t + :name ',name + :props ',properties + :form ',form + :vals ',vals)))) + +(defun add-entry (entry) + (setq entry (copy-entry entry)) + (let* ((pred (gethash (name entry) *entries-table*))) + (cond + (pred + (setf (cadr pred) entry) + (report-error nil + "Redefining test ~:@(~S~)" + (name entry))) + (t + (setf (gethash (name entry) *entries-table*) *entries-tail*) + (setf (cdr *entries-tail*) (cons entry nil)) + (setf *entries-tail* (cdr *entries-tail*)) + ))) + (when *do-tests-when-defined* + (do-entry entry)) + (setq *test* (name entry))) + +(defun report-error (error? &rest args) + (cond (*debug* + (apply #'format t args) + (if error? (throw '*debug* nil))) + (error? (apply #'error args)) + (t (apply #'warn args))) + nil) + +(defun do-test (&optional (name *test*)) + #-sbcl (do-entry (get-entry name)) + #+sbcl (handler-bind ((sb-ext:code-deletion-note #'muffle-warning)) + (do-entry (get-entry name)))) + +(defun my-aref (a &rest args) + (apply #'aref a args)) + +(defun my-row-major-aref (a index) + (row-major-aref a index)) + +(defun equalp-with-case (x y) + "Like EQUALP, but doesn't do case conversion of characters. + Currently doesn't work on arrays of dimension > 2." + (cond + ((eq x y) t) + ((consp x) + (and (consp y) + (equalp-with-case (car x) (car y)) + (equalp-with-case (cdr x) (cdr y)))) + ((and (typep x 'array) + (= (array-rank x) 0)) + (equalp-with-case (my-aref x) (my-aref y))) + ((typep x 'vector) + (and (typep y 'vector) + (let ((x-len (length x)) + (y-len (length y))) + (and (eql x-len y-len) + (loop + for i from 0 below x-len + for e1 = (my-aref x i) + for e2 = (my-aref y i) + always (equalp-with-case e1 e2)))))) + ((and (typep x 'array) + (typep y 'array) + (not (equal (array-dimensions x) + (array-dimensions y)))) + nil) + + ((typep x 'array) + (and (typep y 'array) + (let ((size (array-total-size x))) + (loop for i from 0 below size + always (equalp-with-case (my-row-major-aref x i) + (my-row-major-aref y i)))))) + + (t (eql x y)))) + +(defun do-entry (entry &optional + (s *standard-output*)) + (catch '*in-test* + (setq *test* (name entry)) + (setf (pend entry) t) + (let* ((*in-test* t) + ;; (*break-on-warnings* t) + (aborted nil) + r) + ;; (declare (special *break-on-warnings*)) + + (block aborted + (setf r + (flet ((%do + () + (cond + (*compile-tests* + (multiple-value-list + (funcall (compile + nil + `(lambda () + (declare + (optimize ,@*optimization-settings*)) + ,(form entry)))))) + (*expanded-eval* + (multiple-value-list + (expanded-eval (form entry)))) + (t + (multiple-value-list + (eval (form entry))))))) + (if *catch-errors* + (handler-bind + (#-ecl (style-warning #'muffle-warning) + (error #'(lambda (c) + (setf aborted t) + (setf r (list c)) + (return-from aborted nil)))) + (%do)) + (%do))))) + + (setf (pend entry) + (or aborted + (not (equalp-with-case r (vals entry))))) + + (when (pend entry) + (let ((*print-circle* *print-circle-on-failure*)) + (format s "~&Test ~:@(~S~) failed~ + ~%Form: ~S~ + ~%Expected value~P: ~ + ~{~S~^~%~17t~}~%" + *test* (form entry) + (length (vals entry)) + (vals entry)) + (handler-case + (let ((st (format nil "Actual value~P: ~ + ~{~S~^~%~15t~}.~%" + (length r) r))) + (format s "~A" st)) + (error () (format s "Actual value: #<error during printing>~%") + )) + (finish-output s) + )))) + (when (not (pend entry)) *test*)) + +(defun expanded-eval (form) + "Split off top level of a form and eval separately. This reduces the chance that + compiler optimizations will fold away runtime computation." + (if (not (consp form)) + (eval form) + (let ((op (car form))) + (cond + ((eq op 'let) + (let* ((bindings (loop for b in (cadr form) + collect (if (consp b) b (list b nil)))) + (vars (mapcar #'car bindings)) + (binding-forms (mapcar #'cadr bindings))) + (apply + (the function + (eval `(lambda ,vars ,@(cddr form)))) + (mapcar #'eval binding-forms)))) + ((and (eq op 'let*) (cadr form)) + (let* ((bindings (loop for b in (cadr form) + collect (if (consp b) b (list b nil)))) + (vars (mapcar #'car bindings)) + (binding-forms (mapcar #'cadr bindings))) + (funcall + (the function + (eval `(lambda (,(car vars) &aux ,@(cdr bindings)) ,@(cddr form)))) + (eval (car binding-forms))))) + ((eq op 'progn) + (loop for e on (cdr form) + do (if (null (cdr e)) (return (eval (car e))) + (eval (car e))))) + ((and (symbolp op) (fboundp op) + (not (macro-function op)) + (not (special-operator-p op))) + (apply (symbol-function op) + (mapcar #'eval (cdr form)))) + (t (eval form)))))) + +(defun continue-testing () + (if *in-test* + (throw '*in-test* nil) + (do-entries *standard-output*))) + +(defun do-tests (&optional + (out *standard-output*)) + (dolist (entry (cdr *entries*)) + (setf (pend entry) t)) + (if (streamp out) + (do-entries out) + (with-open-file + (stream out :direction :output) + (do-entries stream)))) + +(defun do-entries* (s) + (format s "~&Doing ~A pending test~:P ~ + of ~A tests total.~%" + (count t (the list (cdr *entries*)) :key #'pend) + (length (cdr *entries*))) + (finish-output s) + (dolist (entry (cdr *entries*)) + (when (and (pend entry) + (not (has-disabled-note entry))) + (format s "~@[~<~%~:; ~:@(~S~)~>~]" + (do-entry entry s)) + (finish-output s) + )) + (let ((pending (pending-tests)) + (expected-table (make-hash-table :test #'equal))) + (dolist (ex *expected-failures*) + (setf (gethash ex expected-table) t)) + (let ((new-failures + (loop for pend in pending + unless (gethash pend expected-table) + collect pend))) + (if (null pending) + (format s "~&No tests failed.") + (progn + (format s "~&~A out of ~A ~ + total tests failed: ~ + ~:@(~{~<~% ~1:;~S~>~ + ~^, ~}~)." + (length pending) + (length (cdr *entries*)) + pending) + (if (null new-failures) + (format s "~&No unexpected failures.") + (when *expected-failures* + (format s "~&~A unexpected failures: ~ + ~:@(~{~<~% ~1:;~S~>~ + ~^, ~}~)." + (length new-failures) + new-failures))) + )) + (finish-output s) + (null pending)))) + +(defun do-entries (s) + #-sbcl (do-entries* s) + #+sbcl (handler-bind ((sb-ext:code-deletion-note #'muffle-warning)) + (do-entries* s))) + +;;; Note handling functions and macros + +(defmacro defnote (name contents &optional disabled) + `(eval-when (:load-toplevel :execute) + (let ((note (make-note :name ',name + :contents ',contents + :disabled ',disabled))) + (setf (gethash (note-name note) *notes*) note) + note))) + +(defun disable-note (n) + (let ((note (if (note-p n) n + (setf n (gethash n *notes*))))) + (unless note (error "~A is not a note or note name." n)) + (setf (note-disabled note) t) + note)) + +(defun enable-note (n) + (let ((note (if (note-p n) n + (setf n (gethash n *notes*))))) + (unless note (error "~A is not a note or note name." n)) + (setf (note-disabled note) nil) + note)) ----------------------------------------------------------------------- hooks/post-receive -- CMU Common Lisp