X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=tests%2Ffloat.impure.lisp;h=621520e9d47cbfdc8807c1c593af031990c501db;hb=26e7568488a46369198e336808b4aba57bbe7a63;hp=1cb8c222c12a11a8ad06497bb2dc20a16ac9af91;hpb=ff8359b556436696f9e8b94195b073da636c719f;p=sbcl.git diff --git a/tests/float.impure.lisp b/tests/float.impure.lisp index 1cb8c22..621520e9 100644 --- a/tests/float.impure.lisp +++ b/tests/float.impure.lisp @@ -234,3 +234,140 @@ (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. +(with-test (:name (:range-reduction :x87) + :skipped-on '(not :x86)) + (flet ((almost= (x y) + (< (abs (- x y)) 1d-5))) + (macrolet ((foo (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))) + (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 copy-pasting a constant and Gauss-Legendre converges +;; extremely fast. +(defun pi-gauss-legendre (n-bits) + "Return a rational approximation to pi using the Gauss-Legendre +algorithm. The calculations are done with integers, representing +multiples of (expt 2 (- N-BITS)), and the result is an integral multiple +of this number. The result is accurate to a few less than N-BITS many +fractional bits." + (let ((a (ash 1 n-bits)) ; scaled 1 + (b (isqrt (expt 2 (1- (* n-bits 2))))) ; scaled (sqrt 1/2) + (c (ash 1 (- n-bits 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 n-bits))) + d (1+ d)))) + (/ (round (expt (+ a b) 2) (* 4 c)) + (ash 1 n-bits)))) + +;; 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 double-float precision even for arguments of +;; very large absolute value. +;; This test is skipped on x86; as to why see the comment at the test +;; (:range-reduction :x87) above. +(with-test (:name (:range-reduction :precise-pi) + :skipped-on :x86) + (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)) + (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 (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))) + (* 8 double-float-epsilon)))) + (error "Inaccurate result for ~a: expected ~a, got ~a" + (list op x) expected actual))))) + (dolist (op '(sin cos tan)) + (dolist (val `(,(coerce most-positive-fixnum 'double-float) + ,@(loop for v = most-positive-double-float + then (expt v 4/5) + while (> v (expt 2 50)) + collect v) + ;; 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))))))