1 ;;;; This file is for floating-point-related tests which have side
2 ;;;; effects (e.g. executing DEFUN).
4 ;;;; This software is part of the SBCL system. See the README file for
7 ;;;; While most of SBCL is derived from the CMU CL system, the test
8 ;;;; files (like this one) were written from scratch after the fork
11 ;;;; This software is in the public domain and is provided with
12 ;;;; absolutely no warranty. See the COPYING and CREDITS files for
13 ;;;; more information.
15 (cl:in-package :cl-user)
17 ;;; Hannu Rummukainen reported a CMU CL bug on cmucl-imp@cons.org 26
18 ;;; Jun 2000. This is the test case for it.
20 ;;; The bug was listed as "39: .. Probably the same bug exists in
21 ;;; SBCL" for a while until Martin Atzmueller showed that it's not
22 ;;; present after all, presumably because the bug was introduced into
23 ;;; CMU CL after the fork. But we'll test for it anyway, in case
24 ;;; e.g. someone inadvertently ports the bad code.
27 :element-type 'double-float
28 :initial-contents (list x y)))
30 (declaim (inline point39-x point39-y))
32 (declare (type (simple-array double-float (2)) p))
35 (declare (type (simple-array double-float (2)) p))
37 (defun order39 (points)
38 (sort points (lambda (p1 p2)
39 (let* ((y1 (point39-y p1))
46 (order39 (make-array 4
47 :initial-contents (list (point39 0.0d0 0.0d0)
50 (point39 3.0d0 3.0d0)))))
51 (assert (equalp (test39)
57 (defun complex-double-float-ppc (x y)
58 (declare (type (complex double-float) x y))
59 (declare (optimize speed))
61 (compile 'complex-double-float-ppc)
62 (assert (= (complex-double-float-ppc #c(0.0d0 1.0d0) #c(2.0d0 3.0d0))
65 (defun single-float-ppc (x)
66 (declare (type (signed-byte 32) x) (optimize speed))
68 (compile 'single-float-ppc)
69 (assert (= (single-float-ppc -30) -30f0))
71 ;;; constant-folding irrational functions
74 ;; do not remove the ECASE here: the bug this checks for indeed
75 ;; depended on this configuration
76 (ecase x (1 least-positive-double-float)))
77 (macrolet ((test (fun)
78 (let ((name (intern (format nil "TEST-CONSTANT-~A" fun))))
80 (defun ,name () (,fun (df 1)))
98 ;;; Broken move-arg-double-float for non-rsp frame pointers on x86-64
100 (declare (optimize speed))
101 (multiple-value-bind (x)
103 (declare (double-float x))
110 (format t "y=~s~%" y)))
115 (assert (= (test 1.0d0) 2.0d0))
117 (deftype myarraytype (&optional (length '*))
118 `(simple-array double-float (,length)))
119 (defun new-pu-label-from-pu-labels (array)
120 (setf (aref (the myarraytype array) 0)
121 sb-ext:double-float-positive-infinity))
125 ;;; FIXME: it may be that TYPE-ERROR is wrong, and we should
126 ;;; instead signal an overflow or coerce into an infinity.
128 (loop for n from (expt 2 1024) upto (+ 10 (expt 2 1024))
130 (coerce n 'single-float)
131 (simple-type-error ()
132 (return-from bug-407a :type-error)))))
133 (assert (eq :type-error (bug-407a)))
135 (loop for n from (expt 2 1024) upto (+ 10 (expt 2 1024))
137 (format t "~E~%" (coerce n 'single-float))
138 (simple-type-error ()
139 (return-from bug-407b :type-error)))))
140 (assert (eq :type-error (bug-407b)))
142 ;; 1.0.29.44 introduces a ton of changes for complex floats
143 ;; on x86-64. Huge test of doom to help catch weird corner
145 ;; Abuse the framework to also test some float arithmetic
146 ;; changes wrt constant arguments in 1.0.29.54.
147 (defmacro def-compute (name real-type
148 &optional (complex-type `(complex ,real-type)))
149 `(defun ,name (x y r)
150 (declare (type ,complex-type x y)
152 (flet ((reflections (x)
155 (complex (- (realpart x)) (imagpart x))
158 (declare (type ,complex-type x y)
160 (list (1+ x) (* 2 x) (/ x 2) (= 1 x)
161 (+ x y) (+ r x) (+ x r)
162 (- x y) (- r x) (- x r)
163 (* x y) (* x r) (* r x)
170 (conjugate x) (conjugate r)
171 (abs r) (- r) (= 1 r)
172 (- x) (1+ r) (* 2 r) (/ r 2)
173 (complex r) (complex r r) (complex 0 r)
174 (= x y) (= r x) (= y r) (= x (complex 0 r))
175 (= r (realpart x)) (= (realpart x) r)
176 (> r (realpart x)) (< r (realpart x))
177 (> (realpart x) r) (< (realpart x) r)
178 (eql x y) (eql x (complex r)) (eql y (complex r))
179 (eql x (complex r r)) (eql y (complex 0 r))
180 (eql r (realpart x)) (eql (realpart x) r))))
181 (declare (inline reflections))
182 (multiple-value-bind (x1 x2 x3 x4) (reflections x)
183 (multiple-value-bind (y1 y2 y3 y4) (reflections y)
184 #.(let ((form '(list)))
185 (dolist (x '(x1 x2 x3 x4) (reverse form))
186 (dolist (y '(y1 y2 y3 y4))
188 (append (compute ,x ,y r)
189 (compute ,x ,y (- r))))
192 (def-compute compute-number real number)
193 (def-compute compute-single single-float)
194 (def-compute compute-double double-float)
196 (labels ((equal-enough (x y)
200 (or (eql (coerce x '(complex double-float))
201 (coerce y '(complex double-float)))
202 (and (equal-enough (realpart x) (realpart y))
203 (equal-enough (imagpart x) (imagpart y)))))
205 (or (eql (coerce x 'double-float) (coerce y 'double-float))
206 (< (abs (- x y)) 1d-5))))))
207 (let* ((reals '(0 1 2))
208 (complexes '#.(let ((reals '(0 1 2))
210 (dolist (x reals (nreverse cpx))
212 (push (complex x y) cpx))))))
213 (declare (notinline every))
215 (dolist (x complexes)
216 (dolist (y complexes)
217 (let ((value (compute-number x y r))
218 (single (compute-single (coerce x '(complex single-float))
219 (coerce y '(complex single-float))
220 (coerce r 'single-float)))
221 (double (compute-double (coerce x '(complex double-float))
222 (coerce y '(complex double-float))
223 (coerce r 'double-float))))
224 (assert (every (lambda (pos ref single double)
225 (declare (ignorable pos))
226 (every (lambda (ref single double)
227 (or (and (equal-enough ref single)
228 (equal-enough ref double))
229 (and (not (numberp single)) ;; -ve 0s
230 (equal-enough single double))))
231 (fourth ref) (fourth single) (fourth double)))
232 '((0 0) (0 1) (0 2) (0 3)
233 (1 0) (1 1) (1 2) (1 3)
234 (2 0) (2 1) (2 2) (2 3)
235 (3 0) (3 1) (3 2) (3 3))
236 value single double))))))))
238 ;; The x86 port used not to reduce the arguments of transcendentals
240 ;; This test is valid only for x86: The x86 port uses the builtin x87
241 ;; FPU instructions to implement the trigonometric functions; other
242 ;; ports rely on the system's math library. These two differ in the
243 ;; precision of pi used for the range reduction and so yield results
244 ;; that can differ by arbitrarily large amounts for large inputs.
245 ;; The test expects the x87 results.
246 (with-test (:name (:range-reduction :x87)
247 :skipped-on '(not :x86))
248 (flet ((almost= (x y)
249 (< (abs (- x y)) 1d-5)))
250 (macrolet ((foo (op value)
251 `(let ((actual (,op ,value))
252 (expected (,op (mod ,value (* 2 pi)))))
253 (unless (almost= actual expected)
254 (error "Inaccurate result for ~a: expected ~a, got ~a"
255 (list ',op ,value) expected actual)))))
256 (let ((big (* pi (expt 2d0 70)))
257 (mid (coerce most-positive-fixnum 'double-float))
258 (odd (* pi most-positive-fixnum)))
262 (foo sin (/ odd 2d0))
267 (foo cos (/ odd 2d0))
273 ;; To test the range reduction of trigonometric functions we need a much
274 ;; more accurate approximation of pi than CL:PI is. Calculating this is
275 ;; more fun than copy-pasting a constant and Gauss-Legendre converges
277 (defun pi-gauss-legendre (n-bits)
278 "Return a rational approximation to pi using the Gauss-Legendre
279 algorithm. The calculations are done with integers, representing
280 multiples of (expt 2 (- N-BITS)), and the result is an integral multiple
281 of this number. The result is accurate to a few less than N-BITS many
283 (let ((a (ash 1 n-bits)) ; scaled 1
284 (b (isqrt (expt 2 (1- (* n-bits 2))))) ; scaled (sqrt 1/2)
285 (c (ash 1 (- n-bits 2))) ; scaled 1/4
290 (let ((a1 (ash (+ a b) -1)))
293 c (- c (ash (expt (- a a1) 2) (- d n-bits)))
295 (/ (round (expt (+ a b) 2) (* 4 c))
298 ;; Test that the range reduction of trigonometric functions is done
299 ;; with a sufficiently accurate value of pi that the reduced argument
300 ;; is correct to nearly double-float precision even for arguments of
301 ;; very large absolute value.
302 ;; This test is skipped on x86; as to why see the comment at the test
303 ;; (:range-reduction :x87) above.
304 (with-test (:name (:range-reduction :precise-pi)
306 (let ((rational-pi-half (/ (pi-gauss-legendre 2200) 2)))
307 (labels ((round-pi-half (x)
308 "Return two values as if (ROUND X (/ PI 2)) was called
309 but where PI is precise enough that for all possible
310 double-float arguments the quotient is exact and the
311 remainder is exact to double-float precision."
312 (declare (type double-float x))
313 (multiple-value-bind (q r)
314 (round (rational x) rational-pi-half)
315 (values q (coerce r 'double-float))))
317 "Calculate (OP X) precisely by shifting the argument by
318 an integral multiple of (/ PI 2) into the range from
319 (- (/ PI 4)) to (/ PI 4) and applying the phase-shift
320 formulas for the trigonometric functions. PI here is
321 precise enough that the result is exact to double-float
323 (labels ((precise-val (op q r)
325 (sin (let ((x (if (zerop (mod q 2))
331 (cos (precise-val 'sin (1+ q) r))
332 (tan (if (zerop (mod q 2))
335 (multiple-value-bind (q r)
337 (precise-val op q r))))
339 (let ((actual (funcall op x))
340 (expected (expected-val op x)))
341 ;; Some of the test values are chosen to lie very near
342 ;; to an integral multiple of pi/2 (within a distance of
343 ;; between 1d-11 and 1d-8), making the absolute value of
344 ;; their sine or cosine this small, too. The absolute
345 ;; value of the tangent is then either similarly small or
346 ;; as large as the reciprocal of this value. Therefore we
347 ;; measure relative instead of absolute error.
348 (unless (or (= actual expected 0)
349 (and (= (signum actual) (signum expected))
350 (< (abs (/ (- actual expected)
351 (+ actual expected)))
352 (* 8 double-float-epsilon))))
353 (error "Inaccurate result for ~a: expected ~a, got ~a"
354 (list op x) expected actual)))))
355 (dolist (op '(sin cos tan))
356 (dolist (val `(,(coerce most-positive-fixnum 'double-float)
357 ,@(loop for v = most-positive-double-float
359 while (> v (expt 2 50))
361 ;; The following values cover all eight combinations
362 ;; of values slightly below or above integral
363 ;; multiples of pi/2 with the integral factor
364 ;; congruent to 0, 1, 2 or 3 modulo 4.
366 4.913896894631919d229
368 3.8335637324151093d42
369 1.8178427396473695d155
371 4.2766818550391727d188
372 1.635888515419299d28))