[git] CMU Common Lisp branch rtoy-simp-dd-trig created. snapshot-2013-12-a-25-g712df0b
Raymond Toy
rtoy at common-lisp.net
Sat Dec 21 01:01:09 UTC 2013
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 at 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 at 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 at 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 at 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 at 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 at 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 at 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
--
CMU Common Lisp
More information about the cmucl-cvs
mailing list