From 69a917778bad1b3c82a8cdd511097adf11a1531a Mon Sep 17 00:00:00 2001 From: Lutz Euler Date: Wed, 28 Aug 2013 23:09:05 +0200 Subject: [PATCH] Improve the test float.impure.lisp / (RANGE-REDUCTION PRECISE-PI). The way the test calculated its expected values was flawed and worked correctly only accidentally due to the specific test values used and to allowing a relatively large margin of error. This commit corrects these calculations, removes some test values and adds others and tightens the error margin. I do not expect this to cause the test's outcome on any platform to change. The flaw was to reduce the arguments by taking the remainder of truncating modulo 2 pi. This allows precise calculations only of the sine and the tangent of values slightly above even multiples of pi, but not for example for the sine of an argument near an odd multiple of pi. Instead the reduction is now done by taking the remainder of rounding to the nearest multiple of pi/2 so that all arguments near the zeroes of both sine and cosine reduce to values near zero. This change was prompted when the test unexpectedly failed with some values from gcc bug 43490 which I tried when investigating lp #1137924. --- tests/float.impure.lisp | 69 ++++++++++++++++++++++++++++++++++------------- 1 file changed, 51 insertions(+), 18 deletions(-) diff --git a/tests/float.impure.lisp b/tests/float.impure.lisp index 29ca23b..621520e9 100644 --- a/tests/float.impure.lisp +++ b/tests/float.impure.lisp @@ -303,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)) @@ -333,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)))))) -- 1.7.10.4