X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=tests%2Ffloat.impure.lisp;h=29ca23b5c0e766ef680100ab86b3da21d94c27c9;hb=1e786e1a23a1f2276ec2dbe197dcc31a53b43738;hp=aaf7eb8c57abbda00634373e064d35aa88c572a1;hpb=175c318c892b0627b36fa3c4db66f59680242204;p=sbcl.git diff --git a/tests/float.impure.lisp b/tests/float.impure.lisp index aaf7eb8..29ca23b 100644 --- a/tests/float.impure.lisp +++ b/tests/float.impure.lisp @@ -114,8 +114,227 @@ (assert (= (test 1.0d0) 2.0d0)) -(deftype myarraytype (&optional (length '*)) +(deftype myarraytype (&optional (length '*)) `(simple-array double-float (,length))) (defun new-pu-label-from-pu-labels (array) - (setf (aref (the myarraytype array) 0) + (setf (aref (the myarraytype array) 0) sb-ext:double-float-positive-infinity)) + +;;; bug 407 +;;; +;;; FIXME: it may be that TYPE-ERROR is wrong, and we should +;;; instead signal an overflow or coerce into an infinity. +(defun bug-407a () + (loop for n from (expt 2 1024) upto (+ 10 (expt 2 1024)) + do (handler-case + (coerce n 'single-float) + (simple-type-error () + (return-from bug-407a :type-error))))) +(assert (eq :type-error (bug-407a))) +(defun bug-407b () + (loop for n from (expt 2 1024) upto (+ 10 (expt 2 1024)) + do (handler-case + (format t "~E~%" (coerce n 'single-float)) + (simple-type-error () + (return-from bug-407b :type-error))))) +(assert (eq :type-error (bug-407b))) + +;; 1.0.29.44 introduces a ton of changes for complex floats +;; on x86-64. Huge test of doom to help catch weird corner +;; cases. +;; Abuse the framework to also test some float arithmetic +;; changes wrt constant arguments in 1.0.29.54. +(defmacro def-compute (name real-type + &optional (complex-type `(complex ,real-type))) + `(defun ,name (x y r) + (declare (type ,complex-type x y) + (type ,real-type r)) + (flet ((reflections (x) + (values x + (conjugate x) + (complex (- (realpart x)) (imagpart x)) + (- x))) + (compute (x y r) + (declare (type ,complex-type x y) + (type ,real-type r)) + (list (1+ x) (* 2 x) (/ x 2) (= 1 x) + (+ x y) (+ r x) (+ x r) + (- x y) (- r x) (- x r) + (* x y) (* x r) (* r x) + (unless (zerop y) + (/ x y)) + (unless (zerop r) + (/ x r)) + (unless (zerop x) + (/ r x)) + (conjugate x) (conjugate r) + (abs r) (- r) (= 1 r) + (- x) (1+ r) (* 2 r) (/ r 2) + (complex r) (complex r r) (complex 0 r) + (= x y) (= r x) (= y r) (= x (complex 0 r)) + (= r (realpart x)) (= (realpart x) r) + (> r (realpart x)) (< r (realpart x)) + (> (realpart x) r) (< (realpart x) r) + (eql x y) (eql x (complex r)) (eql y (complex r)) + (eql x (complex r r)) (eql y (complex 0 r)) + (eql r (realpart x)) (eql (realpart x) r)))) + (declare (inline reflections)) + (multiple-value-bind (x1 x2 x3 x4) (reflections x) + (multiple-value-bind (y1 y2 y3 y4) (reflections y) + #.(let ((form '(list))) + (dolist (x '(x1 x2 x3 x4) (reverse form)) + (dolist (y '(y1 y2 y3 y4)) + (push `(list ,x ,y r + (append (compute ,x ,y r) + (compute ,x ,y (- r)))) + form))))))))) + +(def-compute compute-number real number) +(def-compute compute-single single-float) +(def-compute compute-double double-float) + +(labels ((equal-enough (x y) + (cond ((eql x y)) + ((or (complexp x) + (complexp y)) + (or (eql (coerce x '(complex double-float)) + (coerce y '(complex double-float))) + (and (equal-enough (realpart x) (realpart y)) + (equal-enough (imagpart x) (imagpart y))))) + ((numberp x) + (or (eql (coerce x 'double-float) (coerce y 'double-float)) + (< (abs (- x y)) 1d-5)))))) + (let* ((reals '(0 1 2)) + (complexes '#.(let ((reals '(0 1 2)) + (cpx '())) + (dolist (x reals (nreverse cpx)) + (dolist (y reals) + (push (complex x y) cpx)))))) + (declare (notinline every)) + (dolist (r reals) + (dolist (x complexes) + (dolist (y complexes) + (let ((value (compute-number x y r)) + (single (compute-single (coerce x '(complex single-float)) + (coerce y '(complex single-float)) + (coerce r 'single-float))) + (double (compute-double (coerce x '(complex double-float)) + (coerce y '(complex double-float)) + (coerce r 'double-float)))) + (assert (every (lambda (pos ref single double) + (declare (ignorable pos)) + (every (lambda (ref single double) + (or (and (equal-enough ref single) + (equal-enough ref double)) + (and (not (numberp single)) ;; -ve 0s + (equal-enough single double)))) + (fourth ref) (fourth single) (fourth double))) + '((0 0) (0 1) (0 2) (0 3) + (1 0) (1 1) (1 2) (1 3) + (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-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." + (declare (type double-float x)) + (coerce (mod (rational x) rational-2pi) + 'double-float)) + (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 + ;; measure relative instead of absolute error. + (unless (or (= actual expected 0) + (and (= (signum actual) (signum expected)) + (< (abs (/ (- actual expected) + (+ actual expected))) + (/ 1d-12 2)))) + (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) + 6.543554061677196d28 + 1.5334254929660437d43 + 1.837213298702053d93 + 4.913896894631919d229)) + (test op val))))))