From: Lutz Euler <leuler@us...>  20120607 12:26:04

The branch "master" has been updated in SBCL: via 4de15256dd6387e52c8b9f5588c08044b68f6817 (commit) from d7c9385e5a7b27d38884ebed9df6e9e5c3ac5cbf (commit)  Log  commit 4de15256dd6387e52c8b9f5588c08044b68f6817 Author: Lutz Euler <lutz.euler@...> Date: Thu Jun 7 14:23:11 2012 +0200 Add a test for range reduction of trigonometric functions on nonx86. Trigonometric functions of large float arguments yield differing results depending on how range reduction is done. There is already a test for the results expected on x86, float.pure.lisp/(:rangereduction :x87). Add a test for the results expected from GNU libm and enable it for all platforms except x86. The new test is impure and therefore added to "float.impure.lisp". Move the x86 test there, too, to reduce the potential for future confusion.  tests/float.impure.lisp  101 +++++++++++++++++++++++++++++++++++++++++++++++ tests/float.pure.lisp  32  2 files changed, 101 insertions(+), 32 deletions() diff git a/tests/float.impure.lisp b/tests/float.impure.lisp index 1cb8c22..daef1f2 100644  a/tests/float.impure.lisp +++ b/tests/float.impure.lisp @@ 234,3 +234,104 @@ (2 0) (2 1) (2 2) (2 3) (3 0) (3 1) (3 2) (3 3)) value single double)))))))) + +;; The x86 port used not to reduce the arguments of transcendentals +;; correctly. +;; This test is valid only for x86: The x86 port uses the builtin x87 +;; FPU instructions to implement the trigonometric functions; other +;; ports rely on the system's math library. These two differ in the +;; precision of pi used for the range reduction and so yield results +;; that can differ by arbitrarily large amounts for large inputs. +;; The test expects the x87 results. +(withtest (:name (:rangereduction :x87) + :skippedon '(not :x86)) + (flet ((almost= (x y) + (< (abs ( x y)) 1d5))) + (macrolet ((foo (op value) + `(assert (almost= (,op (mod ,value (* 2 pi))) + (,op ,value))))) + (let ((big (* pi (expt 2d0 70))) + (mid (coerce mostpositivefixnum 'doublefloat)) + (odd (* pi mostpositivefixnum))) + (foo sin big) + (foo sin mid) + (foo sin odd) + (foo sin (/ odd 2d0)) + + (foo cos big) + (foo cos mid) + (foo cos odd) + (foo cos (/ odd 2d0)) + + (foo tan big) + (foo tan mid) + (foo tan odd))))) + +;; To test the range reduction of trigonometric functions we need a much +;; more accurate approximation of pi than CL:PI is. Calculating this is +;; more fun than copypasting a constant and GaussLegendre converges +;; extremely fast. +(defun pigausslegendre (nbits) + "Return a rational approximation to pi using the GaussLegendre +algorithm. The calculations are done with integers, representing +multiples of (expt 2 ( NBITS)), and the result is an integral multiple +of this number. The result is accurate to a few less than NBITS many +fractional bits." + (let ((a (ash 1 nbits)) ; scaled 1 + (b (isqrt (expt 2 (1 (* nbits 2))))) ; scaled (sqrt 1/2) + (c (ash 1 ( nbits 2))) ; scaled 1/4 + (d 0)) + (loop + (when (<= ( a b) 1) + (return)) + (let ((a1 (ash (+ a b) 1))) + (psetf a a1 + b (isqrt (* a b)) + c ( c (ash (expt ( a a1) 2) ( d nbits))) + d (1+ d)))) + (/ (round (expt (+ a b) 2) (* 4 c)) + (ash 1 nbits)))) + +;; Test that the range reduction of trigonometric functions is done +;; with a sufficiently accurate value of pi that the reduced argument +;; is correct to nearly doublefloat precision even for arguments of +;; very large absolute value. +;; This test is skipped on x86; as to why see the comment at the test +;; (:rangereduction :x87) above. +(withtest (:name (:rangereduction :precisepi) + :skippedon :x86) + (let ((rational2pi (* 2 (pigausslegendre 2200)))) + (labels ((mod2pi (x) + "Return X modulo 2 pi, where pi is precise enough that the + result is exact to doublefloat precision for all possible + doublefloat arguments." + (declare (type doublefloat x)) + (coerce (mod (rational x) rational2pi) + 'doublefloat)) + (test (op x) + (let ((actual (funcall op x)) + (expected (funcall op (mod2pi x)))) + ;; Some of the test values are chosen to reduce modulo + ;; 2 pi to small numbers (between 1d10 and 1d7), + ;; making their sine and tangent this small, too. + ;; For other test values the absolute value of the + ;; tangent may be much larger than 1. Therefore we + ;; measure relative instead of absolute error. + (unless (or (= actual expected 0) + (and (= (signum actual) (signum expected)) + (< (abs (/ ( actual expected) + (+ actual expected))) + (/ 1d12 2)))) + (error "Inaccurate result for ~a: expected ~a, got ~a" + (list op x) expected actual))))) + (dolist (op '(sin cos tan)) + (dolist (val `(,(coerce mostpositivefixnum 'doublefloat) + ,@(loop for v = mostpositivedoublefloat + then (expt v 4/5) + while (> v (expt 2 50)) + collect v) + 6.543554061677196d28 + 1.5334254929660437d43 + 1.837213298702053d93 + 4.913896894631919d229)) + (test op val)))))) diff git a/tests/float.pure.lisp b/tests/float.pure.lisp index 62a54b0..5612a73 100644  a/tests/float.pure.lisp +++ b/tests/float.pure.lisp @@ 269,38 +269,6 @@ (the (eql #c(1.0 2.0)) x)))))))) ;; The x86 port used not to reduce the arguments of transcendentals ;; correctly. ;; This test is valid only for x86: The x86 port uses the builtin x87 ;; FPU instructions to implement the trigonometric functions; other ;; ports rely on the system's math library. These two differ in the ;; precision of pi used for the range reduction and so yield results ;; that can differ by arbitrarily large amounts for large inputs. ;; The test expects the x87 results. (withtest (:name (:rangereduction :x87)  :skippedon '(not :x86))  (flet ((almost= (x y)  (< (abs ( x y)) 1d5)))  (macrolet ((foo (op value)  `(assert (almost= (,op (mod ,value (* 2 pi)))  (,op ,value)))))  (let ((big (* pi (expt 2d0 70)))  (mid (coerce mostpositivefixnum 'doublefloat))  (odd (* pi mostpositivefixnum)))  (foo sin big)  (foo sin mid)  (foo sin odd)  (foo sin (/ odd 2d0))   (foo cos big)  (foo cos mid)  (foo cos odd)  (foo cos (/ odd 2d0))   (foo tan big)  (foo tan mid)  (foo tan odd)))))  ;; Leakage from the host could result in wrong values for truncation. (withtest (:name :truncate) (assert (plusp (sbkernel:%unarytruncate/singlefloat (expt 2f0 33))))  hooks/postreceive  SBCL 