X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=tests%2Ffloat.impure.lisp;h=621520e9d47cbfdc8807c1c593af031990c501db;hb=31f68584d0732dc0d17f379773e5f87f1e5a78ad;hp=daef1f2c92d6e017faee4bf132099427e8ce63b5;hpb=4de15256dd6387e52c8b9f5588c08044b68f6817;p=sbcl.git

diff --git a/tests/float.impure.lisp b/tests/float.impure.lisp
index daef1f2..621520e9 100644
--- a/tests/float.impure.lisp
+++ b/tests/float.impure.lisp
@@ -248,8 +248,11 @@
   (flet ((almost= (x y)
            (< (abs (- x y)) 1d-5)))
     (macrolet ((foo (op value)
-                 `(assert (almost= (,op (mod ,value (* 2 pi)))
-                                   (,op ,value)))))
+                 `(let ((actual (,op ,value))
+                        (expected (,op (mod ,value (* 2 pi)))))
+                    (unless (almost= actual expected)
+                      (error "Inaccurate result for ~a: expected ~a, got ~a"
+                             (list ',op ,value) expected actual)))))
       (let ((big (* pi (expt 2d0 70)))
             (mid (coerce most-positive-fixnum 'double-float))
             (odd (* pi most-positive-fixnum)))
@@ -300,28 +303,53 @@ fractional bits."
 ;; (:range-reduction :x87) above.
 (with-test (:name (:range-reduction :precise-pi)
             :skipped-on :x86)
-  (let ((rational-2pi (* 2 (pi-gauss-legendre 2200))))
-    (labels ((mod-2pi (x)
-               "Return X modulo 2 pi, where pi is precise enough that the
-                result is exact to double-float precision for all possible
-                double-float arguments."
+  (let ((rational-pi-half (/ (pi-gauss-legendre 2200) 2)))
+    (labels ((round-pi-half (x)
+               "Return two values as if (ROUND X (/ PI 2)) was called
+                but where PI is precise enough that for all possible
+                double-float arguments the quotient is exact and the
+                remainder is exact to double-float precision."
                (declare (type double-float x))
-               (coerce (mod (rational x) rational-2pi)
-                       'double-float))
+               (multiple-value-bind (q r)
+                   (round (rational x) rational-pi-half)
+                 (values q (coerce r 'double-float))))
+             (expected-val (op x)
+               "Calculate (OP X) precisely by shifting the argument by
+                an integral multiple of (/ PI 2) into the range from
+                (- (/ PI 4)) to (/ PI 4) and applying the phase-shift
+                formulas for the trigonometric functions. PI here is
+                precise enough that the result is exact to double-float
+                precision."
+               (labels ((precise-val (op q r)
+                          (ecase op
+                            (sin (let ((x (if (zerop (mod q 2))
+                                              (sin r)
+                                              (cos r))))
+                                   (if (<= (mod q 4) 1)
+                                       x
+                                       (- x))))
+                            (cos (precise-val 'sin (1+ q) r))
+                            (tan (if (zerop (mod q 2))
+                                     (tan r)
+                                     (/ (- (tan r))))))))
+                 (multiple-value-bind (q r)
+                     (round-pi-half x)
+                   (precise-val op q r))))
              (test (op x)
                (let ((actual (funcall op x))
-                     (expected (funcall op (mod-2pi x))))
-                 ;; Some of the test values are chosen to reduce modulo
-                 ;; 2 pi to small numbers (between 1d-10 and 1d-7),
-                 ;; 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
+                     (expected (expected-val op x)))
+                 ;; Some of the test values are chosen to lie very near
+                 ;; to an integral multiple of pi/2 (within a distance of
+                 ;; between 1d-11 and 1d-8), making the absolute value of
+                 ;; their sine or cosine this small, too. The absolute
+                 ;; value of the tangent is then either similarly small or
+                 ;; as large as the reciprocal of this value. Therefore we
                  ;; measure relative instead of absolute error.
                  (unless (or (= actual expected 0)
                              (and (= (signum actual) (signum expected))
                                   (< (abs (/ (- actual expected)
                                              (+ actual expected)))
-                                     (/ 1d-12 2))))
+                                     (* 8 double-float-epsilon))))
                    (error "Inaccurate result for ~a: expected ~a, got ~a"
                           (list op x) expected actual)))))
       (dolist (op '(sin cos tan))
@@ -330,8 +358,16 @@ fractional bits."
                                then (expt v 4/5)
                                while (> v (expt 2 50))
                                collect v)
-                       6.543554061677196d28
-                       1.5334254929660437d43
-                       1.837213298702053d93
-                       4.913896894631919d229))
+                       ;; The following values cover all eight combinations
+                       ;; of values slightly below or above integral
+                       ;; multiples of pi/2 with the integral factor
+                       ;; congruent to 0, 1, 2 or 3 modulo 4.
+                       5.526916451564098d71
+                       4.913896894631919d229
+                       7.60175752894437d69
+                       3.8335637324151093d42
+                       1.8178427396473695d155
+                       9.41634760758887d89
+                       4.2766818550391727d188
+                       1.635888515419299d28))
           (test op val))))))