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-simp-dd-trig has been created at 712df0bc4e655226bc5c9ed91aa9c875b4a5eb0d (commit)
- Log ----------------------------------------------------------------- commit 712df0bc4e655226bc5c9ed91aa9c875b4a5eb0d Author: Raymond Toy toy.raymond@gmail.com Date: Fri Dec 20 16:47:21 2013 -0800
Add tests for dd-%sincos.
diff --git a/src/tests/trig.lisp b/src/tests/trig.lisp index 911b623..b46c904 100644 --- a/src/tests/trig.lisp +++ b/src/tests/trig.lisp @@ -388,3 +388,46 @@ -4.080663888418042385451434945255951177650840227682488471558860153w-1 1.888w-33)))
+(define-test dd-sincos.signed-zeroes + "Test sincos at 0d0, -0d0" + (:tag :sincos :signed-zeroes :double-double) + (assert-equal '(0w0 1w0) + (multiple-value-list (kernel::dd-%sincos 0w0))) + (assert-equal '(-0w0 1w0) + (multiple-value-list (kernel::dd-%sincos -0w0)))) + +;; Test sincos at a bunch of random points and compare the result from +;; sin and cos. If they differ, save the result in a list to be +;; returned. +(defun dd-sincos-test (limit n) + (let (results) + (dotimes (k n) + (let* ((x (random limit)) + (s-exp (sin x)) + (c-exp (cos x))) + (multiple-value-bind (s c) + (kernel::dd-%sincos x) + (unless (and (eql s s-exp) + (eql c c-exp)) + (push (list x + (list s s-exp) + (list c c-exp)) + results))))) + results)) + +(define-test dd-sincos.consistent + "Test sincos is consistent with sin and cos" + (:tag :sincos :double-double) + ;; Small values + (assert-eql nil + (dd-sincos-test (/ kernel:dd-pi 4) 1000)) + ;; Medium + (assert-eql nil + (dd-sincos-test 16w0 1000)) + ;; Large + (assert-eql nil + (dd-sincos-test (scale-float 1w0 120) 1000)) + ;; Very large + (assert-eql nil + (dd-sincos-test (scale-float 1w0 1023) 1000))) +
commit bf84dbc8c5bd5478fd36b55f99e119cfff11ca6d Author: Raymond Toy toy.raymond@gmail.com Date: Fri Dec 20 16:47:07 2013 -0800
Add dd-%sincos and use it as needed instead of calling sin and cos separately.
diff --git a/src/code/irrat-dd.lisp b/src/code/irrat-dd.lisp index 6661f2b..381d678 100644 --- a/src/code/irrat-dd.lisp +++ b/src/code/irrat-dd.lisp @@ -1191,6 +1191,29 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 (dd-%%tan reduced) (- (/ (dd-%%tan reduced))))))))
+(defun dd-%sincos (x) + (declare (double-double-float x)) + (cond ((< (abs x) (/ pi 4)) + (values (dd-%%sin x) + (dd-%%cos x))) + (t + ;; Argument reduction needed + (multiple-value-bind (n reduced) + (reduce-arg x) + (case (logand n 3) + (0 + (values (dd-%%sin reduced) + (dd-%%cos reduced))) + (1 + (values (dd-%%cos reduced) + (- (dd-%%sin reduced)))) + (2 + (values (- (dd-%%sin reduced)) + (- (dd-%%cos reduced)))) + (3 + (values (- (dd-%%cos reduced)) + (dd-%%sin reduced)))))))) + ;;; dd-%log2 ;;; Base 2 logarithm.
diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp index 4ccf80a..078e56f 100644 --- a/src/code/irrat.lisp +++ b/src/code/irrat.lisp @@ -1298,7 +1298,9 @@ (coerce s '(dispatch-type theta))))) #+double-double ((double-double-float) - (complex (cos theta) (sin theta)))))) + (multiple-value-bind (s c) + (dd-%sincos theta) + (complex c s))))))
(defun asin (number) "Return the arc sine of NUMBER." diff --git a/src/compiler/float-tran.lisp b/src/compiler/float-tran.lisp index d123d18..34acdeb 100644 --- a/src/compiler/float-tran.lisp +++ b/src/compiler/float-tran.lisp @@ -748,8 +748,9 @@
#+double-double (deftransform cis ((z) (double-double-float) *) - ;; Cis. - '(complex (cos z) (sin z))) + `(multiple-value-bind (s c) + (kernel::dd-%sincos x) + (complex c s)))
;;; The argument range is limited on the x86 FP trig. functions. A
commit 2e3e48d466c67a09ea7aeb23106fdc50143be3b5 Author: Raymond Toy toy.raymond@gmail.com Date: Fri Dec 20 16:07:12 2013 -0800
Add tests for double-double trig functions.
diff --git a/src/tests/trig.lisp b/src/tests/trig.lisp index 58d10b6..911b623 100644 --- a/src/tests/trig.lisp +++ b/src/tests/trig.lisp @@ -215,3 +215,176 @@ (assert-eql nil (sincos-test (scale-float 1d0 1023) 1000)))
+;; Compute the relative error between actual and expected if expected +;; is not zero. Otherwise, return absolute error between actual and +;; expected. If the error is less than the threshold, return T. +;; Otherwise return the actual (relative or absolute) error. +(defun rel-or-abs-error (actual expected &optional (threshold double-float-epsilon)) + (let ((err (if (zerop expected) + (abs (- actual expected)) + (/ (abs (- actual expected)) + (abs expected))))) + (if (<= err threshold) + t + err))) + +(define-test dd-sin.signed-zeroes + "Test sin for 0w0 and -0w0" + (:tag :sin :double-double :signed-zeroes) + (assert-eql 0w0 (sin 0w0)) + (assert-equal -0w0 (sin -0w0))) + +(define-test dd-sin.no-reduction + "Test sin for small args without reduction" + (:tag :sin :double-double) + (assert-eq t (rel-or-abs-error + (sin .5w0) + 4.794255386042030002732879352155713880818033679406006751886166131w-1 + 1w-32)) + (assert-eq t (rel-or-abs-error + (sin -0.5w0) + -4.794255386042030002732879352155713880818033679406006751886166131w-1 + 1w-32))) + +(define-test dd-sin.pi/2 + "Test for arg near pi/2" + (:tag :sin :double-double) + (assert-eq t (rel-or-abs-error + (sin (/ kernel:dd-pi 2)) + 1w0 + 1w-50))) + +;; The reference value were computed using maxima. Here's how to +;; compute the reference value. Set fpprec:64 to tell maxima to use +;; 64 digits of precision. For 7/4*pi, do (integer-decode-float (* 7/4 +;; kernel:dd-pi)) to get the exact rational representation of the +;; desired double-double-float. Then bfloat(sin(<rational>)). +(define-test dd-sin.arg-reduction + "Test for sin with arg reduction" + (:tag :sin :double-double) + ;; Test for argument reduction with n mod 4 = 0 + (assert-eq t (rel-or-abs-error + (sin (* 7/4 kernel:dd-pi)) + -7.07106781186547524400844362104849691328261037289050238659653433w-1 + 0w0)) + ;; Test for argument reduction with n mod 4 = 1 + (assert-eq t (rel-or-abs-error + (sin (* 9/4 kernel:dd-pi)) + 7.07106781186547524400844362104858161816423215627023442400880643w-1 + 0w0)) + ;; Test for argument reduction with n mod 4 = 2 + (assert-eq t (rel-or-abs-error + (sin (* 11/4 kernel:dd-pi)) + 7.071067811865475244008443621048998682901731241858306822215522497w-1 + 8.716w-33)) + ;; Test for argument reduction with n mod 4 = 3 + (assert-eq t (rel-or-abs-error + (sin (* 13/4 kernel:dd-pi)) + -7.071067811865475244008443621048777109664479707052746581685893187w-1 + 8.716w-33)) + ;; Test for argument reduction, big value + (assert-eq t (rel-or-abs-error + (sin (scale-float 1w0 120)) + 3.778201093607520226555484700569229919605866976512306642257987199w-1 + 8.156w-33))) + +(define-test dd-cos.signed-zeroes + "Test cos for 0w0 and -0w0" + (:tag :cos :double-double :signed-zeroes) + (assert-eql 1w0 (cos 0w0)) + (assert-equal 1w0 (cos -0w0))) + +(define-test dd-cos.no-reduction + "Test cos for small args without reduction" + (:tag :cos :double-double) + (assert-eq t (rel-or-abs-error + (cos .5w0) + 8.775825618903727161162815826038296519916451971097440529976108683w-1 + 0w0)) + (assert-eq t (rel-or-abs-error + (cos -0.5w0) + 8.775825618903727161162815826038296519916451971097440529976108683w-1 + 0w0))) + +(define-test dd-cos.pi/2 + "Test for arg near pi/2" + (:tag :cos :double-double) + (assert-eq t (rel-or-abs-error + (cos (/ kernel:dd-pi 2)) + -1.497384904859169777320797133937725094986669701841027904483071358w-33 + 0w0))) + +(define-test dd-cos.arg-reduction + "Test for cos with arg reduction" + (:tag :cos :double-double) + ;; Test for argument reduction with n mod 4 = 0 + (assert-eq t (rel-or-abs-error + (cos (* 7/4 kernel:dd-pi)) + 7.07106781186547524400844362104849691328261037289050238659653433w-1 + 0w0)) + ;; Test for argument reduction with n mod 4 = 1 + (assert-eq t (rel-or-abs-error + (cos (* 9/4 kernel:dd-pi)) + 7.07106781186547524400844362104858161816423215627023442400880643w-1 + 3.487w-32)) + ;; Test for argument reduction with n mod 4 = 2 + (assert-eq t (rel-or-abs-error + (cos (* 11/4 kernel:dd-pi)) + -7.071067811865475244008443621048998682901731241858306822215522497w-1 + 1.482w-31)) + ;; Test for argument reduction with n mod 4 = 3 + (assert-eq t (rel-or-abs-error + (cos (* 13/4 kernel:dd-pi)) + -7.071067811865475244008443621048777109664479707052746581685893187w-1 + 7.845w-32)) + ;; Test for argument reduction, big value + (assert-eq t (rel-or-abs-error + (cos (scale-float 1w0 120)) + -9.258790228548378673038617641074149467308332099286564602360493726w-1 + 0w0))) + +(define-test dd-tan.signed-zeroes + "Test tan for 0w0 and -0w0" + (:tag :tan :double-double :signed-zeroes) + (assert-eql 0w0 (tan 0w0)) + (assert-equal -0w0 (tan -0w0))) + +(define-test dd-tan.no-reduction + "Test tan for small args without reduction" + (:tag :tan :double-double) + (assert-eq t (rel-or-abs-error + (tan .5w0) + 5.463024898437905132551794657802853832975517201797912461640913859w-1 + 0w0)) + (assert-eq t (rel-or-abs-error + (tan -0.5w0) + -5.463024898437905132551794657802853832975517201797912461640913859w-1 + 0w0))) + +(define-test dd-tan.pi/2 + "Test for arg near pi/2" + (:tag :tan :double-double) + (assert-eq t (rel-or-abs-error + (tan (/ kernel:dd-pi 2)) + -6.67830961000672557834948096545679895621313886078988606234681001w32 + 0w0))) + +(define-test dd-tan.arg-reduction + "Test for tan with arg reduction" + (:tag :tan :double-double) + ;; Test for argument reduction with n even + (assert-eq t (rel-or-abs-error + (tan (* 7/4 kernel:dd-pi)) + -1.000000000000000000000000000000001844257310064121018312678894979w0 + 6.467w-33)) + ;; Test for argument reduction with n odd + (assert-eq t (rel-or-abs-error + (tan (* 9/4 kernel:dd-pi)) + 1.000000000000000000000000000000025802415787810837455445433037983w0 + 5.773w-33)) + ;; Test for argument reduction, big value + (assert-eq t (rel-or-abs-error + (tan (scale-float 1w0 120)) + -4.080663888418042385451434945255951177650840227682488471558860153w-1 + 1.888w-33))) +
commit 949acab5719c22e51a45a936271b46b04edaa8ac Author: Raymond Toy toy.raymond@gmail.com Date: Fri Dec 20 16:06:40 2013 -0800
For dd-%%sin, return x if x is small enough. (Makes sin(-0w0) be -0w0).
diff --git a/src/code/irrat-dd.lisp b/src/code/irrat-dd.lisp index 170bc06..6661f2b 100644 --- a/src/code/irrat-dd.lisp +++ b/src/code/irrat-dd.lisp @@ -1000,8 +1000,11 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 (declare (type (double-double-float -1w0 1w0) x) (optimize (speed 2) (space 0) (inhibit-warnings 3))) - (let ((x2 (* x x))) - (+ x (* x (* x2 (poly-eval x2 sincof)))))) + (if (< (abs (double-double-hi x)) + (scale-float 1d0 -52)) + x + (let ((x2 (* x x))) + (+ x (* x (* x2 (poly-eval x2 sincof)))))))
;; cos(x) = 1 - .5 x^2 + x^2 (x^2 P(x^2)) ;; Theoretical peak relative error = 2.1e-37,
commit 6f25e2e894bdf13c00a29651031ad0bbedb50f0e Merge: 408aa78 82d0a77 Author: Raymond Toy toy.raymond@gmail.com Date: Fri Dec 20 12:39:48 2013 -0800
Merge branch 'master' into rtoy-simp-dd-trig
commit 408aa78aa3947f2c8f8a5b2a03429d5c05e93fbe Merge: 00bd409 01a3f47 Author: Raymond Toy toy.raymond@gmail.com Date: Fri Dec 20 07:44:09 2013 -0800
Merge branch 'master' into rtoy-simp-dd-trig
commit 00bd409b8d9de4f5f6223bf4d017e6f1c0826e48 Author: Raymond Toy toy.raymond@gmail.com Date: Thu Dec 12 18:48:41 2013 -0800
Simplify dd-%%sin, dd-%%cos, and dd-%%tan.
These routines did argument reduction, but since we use __kernel_rem_pio2 to do accurate argument reduction, the argument reduction in these routines is a waste of time. This greatly simplifies the routines to just the polynomial (or rational) approximations.
diff --git a/src/code/irrat-dd.lisp b/src/code/irrat-dd.lisp index 4c57165..170bc06 100644 --- a/src/code/irrat-dd.lisp +++ b/src/code/irrat-dd.lisp @@ -995,6 +995,14 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 -1.666666666666666666666666666666666647199w-1 )))
+;; Compute sin(x) for |x| < pi/4 (approx). +(defun dd-%%sin (x) + (declare (type (double-double-float -1w0 1w0) x) + (optimize (speed 2) (space 0) + (inhibit-warnings 3))) + (let ((x2 (* x x))) + (+ x (* x (* x2 (poly-eval x2 sincof)))))) + ;; cos(x) = 1 - .5 x^2 + x^2 (x^2 P(x^2)) ;; Theoretical peak relative error = 2.1e-37, ;; relative peak error spread = 1.4e-8 @@ -1016,101 +1024,17 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 4.166666666666666666666666666666459301466w-2 )))
-(defconstant dp1 - (scale-float (float #b1100100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111000000011 1w0) -106)) - -(defconstant dp2 - (scale-float (float #b0111000001110011010001001010010000001001001110000010001000101001100111110011000111010000000010000010111011 1w0) (* 2 -106))) - -(defconstant dp3 - (scale-float (float #b1110101001100011101100010011100110110010001001010001010010100000100001111001100011100011010000000100110111 1w0) (* 3 -106))) - -(defconstant dp4 - (scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106))) - -(defun dd-%%sin (x) - (declare (type double-double-float x) - (optimize (speed 2) (space 0) - (inhibit-warnings 3))) - (when (minusp x) - (return-from dd-%%sin (- (the double-double-float (dd-%%sin (- x)))))) - ;; y = integer part of x/(pi/4). - (let* ((y (float (floor (/ x dd-pi/4)) 1w0)) - (z (scale-float y -4))) - (declare (type double-double-float y z)) - (setf z (float (floor z) 1w0)) ; integer part of y/8 - (setf z (- y (scale-float z 4))) ; y - 16*(y/16) - - (let ((j (truncate z)) - (sign 1)) - (declare (type (integer -1 1) sign)) - (unless (zerop (logand j 1)) - (incf j) - (incf y)) - (setf j (logand j 7)) - - (when (> j 3) - (setf sign (- sign)) - (decf j 4)) - - ;; Extended precision modular arithmetic - (setf z (- (- (- x (* y dp1)) - (* y dp2)) - (* y dp3))) - (let ((zz (* z z))) - (if (or (= j 1) - (= j 2)) - (setf y (+ (- 1 (scale-float zz -1)) - (* zz zz (poly-eval zz coscof)))) - (setf y (+ z (* z (* zz (poly-eval zz sincof)))))) - (if (< sign 0) - (- y) - y))))) - +;; Compue cos(x) for |x| < pi/4 (approx) (defun dd-%%cos (x) - (declare (type double-double-float x) + (declare (type (double-double-float -1w0 1w0) x) (optimize (speed 2) (space 0) (inhibit-warnings 3))) - (when (minusp x) - (return-from dd-%%cos (dd-%%cos (- x)))) - ;; y = integer part of x/(pi/4). - (let* ((y (float (floor (/ x dd-pi/4)) 1w0)) - (z (scale-float y -4))) - (declare (type double-double-float y z)) - (setf z (float (floor z) 1w0)) ; integer part of y/8 - (setf z (- y (scale-float z 4))) ; y - 16*(y/16) - - (let ((i (truncate z)) - (j 0) - (sign 1)) - (declare (type (integer 0 7) j) - (type (integer -1 1) sign)) - (unless (zerop (logand i 1)) - (incf i) - (incf y)) - (setf j (logand i 7)) - - (when (> j 3) - (setf sign (- sign)) - (decf j 4)) - (when (> j 1) - (setf sign (- sign))) - - ;; Extended precision modular arithmetic. This is basically - ;; computing x - y*(pi/4) accurately so that |z| < pi/4. - (setf z (- (- (- x (* y dp1)) - (* y dp2)) - (* y dp3))) - (let ((zz (* z z))) - (if (or (= j 1) - (= j 2)) - (setf y (+ z (* z (* zz (poly-eval zz sincof))))) - (setf y (+ (- 1 (scale-float zz -1)) - (* zz (poly-eval zz coscof) zz)))) - (if (< sign 0) - (- y) - y))))) + (let ((x2 (* x x))) + (+ (- 1 (scale-float x2 -1)) + (* x2 (poly-eval x2 coscof) x2))))
+;; Compute tan(x) or cot(x) for |x| < pi/4 (approx). If cotflag is +;; non-nil, cot(x) is returned. Otherwise, return tan(x). (let ((P (make-array 6 :element-type 'double-double-float :initial-contents '( @@ -1132,50 +1056,18 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 -4.152206921457208101480801635640958361612w10 8.650244186622719093893836740197250197602w10 )))) - (defun dd-tancot (xx cotflag) - (declare (type double-double-float xx) - (optimize (speed 2) (space 0))) - (let ((x 0w0) - (sign 1)) - (declare (type double-double-float x) - (type (integer -1 1) sign)) - (cond ((minusp xx) - (setf x (- xx)) - (setf sign -1)) - (t - (setf x xx))) - (let* ((y (float (floor (/ x dd-pi/4)) 1w0)) - (z (scale-float y -4)) - (j 0)) - (declare (type double-double-float y z) - (type fixnum j)) - (setf z (float (floor z) 1w0)) - (setf z (- y (scale-float z 4))) - - (setf j (truncate z)) - - (unless (zerop (logand j 1)) - (incf j) - (incf y)) - - (setf z (- (- (- x (* y dp1)) - (* y dp2)) - (* y dp3))) - (let ((zz (* z z))) - (if (> zz 1w-40) - (setf y (+ z - (* z (* zz (/ (poly-eval zz p) - (poly-eval-1 zz q)))))) - (setf y z)) - (if (not (zerop (logand j 2))) - (if cotflag - (setf y (- y)) - (setf y (/ -1 y))) - (if cotflag - (setf y (/ y)))) - (if (< sign 0) - (- y) - y)))))) + (defun dd-tancot (x cotflag) + (declare (type (double-double-float -1w0 1w0) x) + (optimize (speed 2) (space 0) (inhibit-warnings 3))) + (let* ((xx (* x x)) + (y (if (> xx 1w-40) + (+ x + (* x (* xx (/ (poly-eval xx p) + (poly-eval-1 xx q))))) + x))) + (if cotflag + (/ y) + y))))
(defun dd-%%tan (x) (declare (type double-double-float x)) @@ -1254,9 +1146,7 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 dd-%sin)) (defun dd-%sin (x) (declare (double-double-float x)) - (cond ((minusp (float-sign x)) - (- (dd-%sin (- x)))) - ((< (abs x) (/ pi 4)) + (cond ((< (abs x) (/ pi 4)) (dd-%%sin x)) (t ;; Argument reduction needed @@ -1272,9 +1162,7 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 dd-%cos)) (defun dd-%cos (x) (declare (double-double-float x)) - (cond ((minusp x) - (dd-%cos (- x))) - ((< (abs x) (/ pi 4)) + (cond ((< (abs x) (/ pi 4)) (dd-%%cos x)) (t ;; Argument reduction needed @@ -1290,9 +1178,7 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010 dd-%tan)) (defun dd-%tan (x) (declare (double-double-float x)) - (cond ((minusp (float-sign x)) - (- (dd-%tan (- x)))) - ((< (abs x) (/ pi 4)) + (cond ((< (abs x) (/ pi 4)) (dd-%%tan x)) (t ;; Argument reduction needed
-----------------------------------------------------------------------
hooks/post-receive