777b370a65487c2ab7b231c510a0925a1a3d005f
[sbcl.git] / arith.pure.lisp
1 ;;;; arithmetic tests with no side effects
2
3 ;;;; This software is part of the SBCL system. See the README file for
4 ;;;; more information.
5 ;;;;
6 ;;;; While most of SBCL is derived from the CMU CL system, the test
7 ;;;; files (like this one) were written from scratch after the fork
8 ;;;; from CMU CL.
9 ;;;;
10 ;;;; This software is in the public domain and is provided with
11 ;;;; absolutely no warranty. See the COPYING and CREDITS files for
12 ;;;; more information.
13
14 (cl:in-package :cl-user)
15
16 ;;; Once upon a time, in the process of porting CMUCL's SPARC backend
17 ;;; to SBCL, multiplications were excitingly broken.  While it's
18 ;;; unlikely that anything with such fundamental arithmetic errors as
19 ;;; these are going to get this far, it's probably worth checking.
20 (macrolet ((test (op res1 res2)
21              `(progn
22                (assert (= (,op 4 2) ,res1))
23                (assert (= (,op 2 4) ,res2))
24                (assert (= (funcall (compile nil '(lambda (x y) (,op x y))) 4 2)
25                         ,res1))
26                (assert (= (funcall (compile nil '(lambda (x y) (,op x y))) 2 4)
27                         ,res2)))))
28   (test + 6 6)
29   (test - 2 -2)
30   (test * 8 8)
31   (test / 2 1/2)
32   (test expt 16 16))
33
34 ;;; In a bug reported by Wolfhard Buss on cmucl-imp 2002-06-18 (BUG
35 ;;; 184), sbcl didn't catch all divisions by zero, notably divisions
36 ;;; of bignums and ratios by 0.  Fixed in sbcl-0.7.6.13.
37 (assert (raises-error? (/ 2/3 0) division-by-zero))
38 (assert (raises-error? (/ (1+ most-positive-fixnum) 0) division-by-zero))
39
40 ;;; In a bug reported by Raymond Toy on cmucl-imp 2002-07-18, (COERCE
41 ;;; <RATIONAL> '(COMPLEX FLOAT)) was failing to return a complex
42 ;;; float; a patch was given by Wolfhard Buss cmucl-imp 2002-07-19.
43 (assert (= (coerce 1 '(complex float)) #c(1.0 0.0)))
44 (assert (= (coerce 1/2 '(complex float)) #c(0.5 0.0)))
45 (assert (= (coerce 1.0d0 '(complex float)) #c(1.0d0 0.0d0)))
46
47 ;;; (COERCE #c(<RATIONAL> <RATIONAL>) '(complex float)) resulted in
48 ;;; an error up to 0.8.17.31
49 (assert (= (coerce #c(1 2) '(complex float)) #c(1.0 2.0)))
50
51 ;;; COERCE also sometimes failed to verify that a particular coercion
52 ;;; was possible (in particular coercing rationals to bounded float
53 ;;; types.
54 (assert (raises-error? (coerce 1 '(float 2.0 3.0)) type-error))
55 (assert (raises-error? (coerce 1 '(single-float -1.0 0.0)) type-error))
56 (assert (eql (coerce 1 '(single-float -1.0 2.0)) 1.0))
57
58 ;;; ANSI says MIN and MAX should signal TYPE-ERROR if any argument
59 ;;; isn't REAL. SBCL 0.7.7 didn't in the 1-arg case. (reported as a
60 ;;; bug in CMU CL on #lisp IRC by lrasinen 2002-09-01)
61 (assert (null (ignore-errors (min '(1 2 3)))))
62 (assert (= (min -1) -1))
63 (assert (null (ignore-errors (min 1 #(1 2 3)))))
64 (assert (= (min 10 11) 10))
65 (assert (null (ignore-errors (min (find-package "CL") -5.0))))
66 (assert (= (min 5.0 -3) -3))
67 (assert (null (ignore-errors (max #c(4 3)))))
68 (assert (= (max 0) 0))
69 (assert (null (ignore-errors (max "MIX" 3))))
70 (assert (= (max -1 10.0) 10.0))
71 (assert (null (ignore-errors (max 3 #'max))))
72 (assert (= (max -3 0) 0))
73
74 ;;; (CEILING x 2^k) was optimized incorrectly
75 (loop for divisor in '(-4 4)
76       for ceiler = (compile nil `(lambda (x)
77                                    (declare (fixnum x))
78                                    (declare (optimize (speed 3)))
79                                    (ceiling x ,divisor)))
80       do (loop for i from -5 to 5
81                for exact-q = (/ i divisor)
82                do (multiple-value-bind (q r)
83                       (funcall ceiler i)
84                     (assert (= (+ (* q divisor) r) i))
85                     (assert (<= exact-q q))
86                     (assert (< q (1+ exact-q))))))
87
88 ;;; (TRUNCATE x 2^k) was optimized incorrectly
89 (loop for divisor in '(-4 4)
90       for truncater = (compile nil `(lambda (x)
91                                       (declare (fixnum x))
92                                       (declare (optimize (speed 3)))
93                                       (truncate x ,divisor)))
94       do (loop for i from -9 to 9
95                for exact-q = (/ i divisor)
96                do (multiple-value-bind (q r)
97                       (funcall truncater i)
98                     (assert (= (+ (* q divisor) r) i))
99                     (assert (<= (abs q) (abs exact-q)))
100                     (assert (< (abs exact-q) (1+ (abs q)))))))
101
102 ;;; CEILING had a corner case, spotted by Paul Dietz
103 (assert (= (ceiling most-negative-fixnum (1+ most-positive-fixnum)) -1))
104
105 ;;; give any optimizers of constant multiplication a light testing.
106 ;;; 100 may seem low, but (a) it caught CSR's initial errors, and (b)
107 ;;; before checking in, CSR tested with 10000.  So one hundred
108 ;;; checkins later, we'll have doubled the coverage.
109 (dotimes (i 100)
110   (let* ((x (random most-positive-fixnum))
111          (x2 (* x 2))
112          (x3 (* x 3)))
113     (let ((fn (handler-bind ((sb-ext:compiler-note
114                               (lambda (c)
115                                 (when (<= x3 most-positive-fixnum)
116                                   (error c)))))
117                 (compile nil
118                          `(lambda (y)
119                             (declare (optimize speed) (type (integer 0 3) y))
120                             (* y ,x))))))
121       (unless (and (= (funcall fn 0) 0)
122                    (= (funcall fn 1) x)
123                    (= (funcall fn 2) x2)
124                    (= (funcall fn 3) x3))
125         (error "bad results for ~D" x)))))
126
127 ;;; Bugs reported by Paul Dietz:
128
129 ;;; (GCD 0 x) must return (abs x)
130 (dolist (x (list -10 (* 3 most-negative-fixnum)))
131   (assert (= (gcd 0 x) (abs x))))
132 ;;; LCM returns a non-negative number
133 (assert (= (lcm 4 -10) 20))
134 (assert (= (lcm 0 0) 0))
135
136 ;;; PPC bignum arithmetic bug:
137 (multiple-value-bind (quo rem)
138     (truncate 291351647815394962053040658028983955 10000000000000000000000000)
139   (assert (= quo 29135164781))
140   (assert (= rem 5394962053040658028983955)))
141
142 ;;; x86 LEA bug:
143 (assert (= (funcall
144             (compile nil '(lambda (x) (declare (bit x)) (+ x #xf0000000)))
145             1)
146            #xf0000001))
147
148 ;;; LOGBITP on bignums:
149 (dolist (x '(((1+ most-positive-fixnum) 1 nil)
150              ((1+ most-positive-fixnum) -1 t)
151              ((1+ most-positive-fixnum) (1+ most-positive-fixnum) nil)
152              ((1+ most-positive-fixnum) (1- most-negative-fixnum) t)
153              (1 (ash most-negative-fixnum 1) nil)
154              (#.(- sb-vm:n-word-bits sb-vm:n-fixnum-tag-bits 1) most-negative-fixnum t)
155              (#.(1+ (- sb-vm:n-word-bits sb-vm:n-fixnum-tag-bits 1)) (ash most-negative-fixnum 1) t)
156              (#.(+ 2 (- sb-vm:n-word-bits sb-vm:n-fixnum-tag-bits 1)) (ash most-negative-fixnum 1) t)
157              (#.(+ sb-vm:n-word-bits 32) (ash most-negative-fixnum #.(+ 32 sb-vm:n-fixnum-tag-bits 2)) nil)
158              (#.(+ sb-vm:n-word-bits 33) (ash most-negative-fixnum #.(+ 32 sb-vm:n-fixnum-tag-bits 2)) t)))
159   (destructuring-bind (index int result) x
160     (assert (eq (eval `(logbitp ,index ,int)) result))))
161
162 ;;; off-by-1 type inference error for %DPB and %DEPOSIT-FIELD:
163 (let ((f (compile nil '(lambda (b)
164                         (integer-length (dpb b (byte 4 28) -1005))))))
165   (assert (= (funcall f 1230070) 32)))
166 (let ((f (compile nil '(lambda (b)
167                         (integer-length (deposit-field b (byte 4 28) -1005))))))
168   (assert (= (funcall f 1230070) 32)))
169
170 ;;; type inference leading to an internal compiler error:
171 (let ((f (compile nil '(lambda (x)
172                         (declare (type fixnum x))
173                         (ldb (byte 0 0) x)))))
174   (assert (= (funcall f 1) 0))
175   (assert (= (funcall f most-positive-fixnum) 0))
176   (assert (= (funcall f -1) 0)))
177
178 ;;; Alpha bignum arithmetic bug:
179 (assert (= (* 966082078641 419216044685) 404997107848943140073085))
180
181 ;;; Alpha smallnum arithmetic bug:
182 (assert (= (ash -129876 -1026) -1))
183
184 ;;; Alpha middlenum (yes, really! Affecting numbers between 2^32 and
185 ;;; 2^64 :) arithmetic bug
186 (let ((fn (compile nil '(LAMBDA (A B C D)
187           (DECLARE (TYPE (INTEGER -1621 -513) A)
188                    (TYPE (INTEGER -3 34163) B)
189                    (TYPE (INTEGER -9485132993 81272960) C)
190                    (TYPE (INTEGER -255340814 519943) D)
191                    (IGNORABLE A B C D)
192                    (OPTIMIZE (SPEED 3) (SAFETY 1) (DEBUG 1)))
193           (TRUNCATE C (MIN -100 4149605))))))
194   (assert (= (funcall fn -1332 5864 -6963328729 -43789079) 69633287)))
195
196 ;;; Here's another fantastic Alpha backend bug: the code to load
197 ;;; immediate 64-bit constants into a register was wrong.
198 (let ((fn (compile nil '(LAMBDA (A B C D)
199           (DECLARE (TYPE (INTEGER -3563 2733564) A)
200                    (TYPE (INTEGER -548947 7159) B)
201                    (TYPE (INTEGER -19 0) C)
202                    (TYPE (INTEGER -2546009 0) D)
203                    (IGNORABLE A B C D)
204                    (OPTIMIZE (SPEED 3) (SAFETY 1) (DEBUG 1)))
205           (CASE A
206             ((89 125 16) (ASH A (MIN 18 -706)))
207             (T (DPB -3 (BYTE 30 30) -1)))))))
208   (assert (= (funcall fn 1227072 -529823 -18 -792831) -2147483649)))
209
210 ;;; ASH of a negative bignum by a bignum count would erroneously
211 ;;; return 0 prior to sbcl-0.8.4.4
212 (assert (= (ash (1- most-negative-fixnum) (1- most-negative-fixnum)) -1))
213
214 ;;; Whoops.  Too much optimization in division operators for 0
215 ;;; divisor.
216 (macrolet ((frob (name)
217              `(let ((fn (compile nil '(lambda (x)
218                                        (declare (optimize speed) (fixnum x))
219                                        (,name x 0)))))
220                (assert (raises-error? (funcall fn 1) division-by-zero)))))
221   (frob mod)
222   (frob truncate)
223   (frob rem)
224   (frob /)
225   (frob floor)
226   (frob ceiling))
227
228 ;; Check that the logic in SB-KERNEL::BASIC-COMPARE for doing fixnum/float
229 ;; comparisons without rationalizing the floats still gives the right anwers
230 ;; in the edge cases (had a fencepost error).
231 (macrolet ((test (range type sign)
232              `(let (ints
233                     floats
234                     (start (- ,(find-symbol (format nil
235                                                     "MOST-~A-EXACTLY-~A-FIXNUM"
236                                                     sign type)
237                                             :sb-kernel)
238                               ,range)))
239                 (dotimes (i (1+ (* ,range 2)))
240                   (let* ((x (+ start i))
241                          (y (coerce x ',type)))
242                     (push x ints)
243                     (push y floats)))
244                 (dolist (i ints)
245                   (dolist (f floats)
246                     (dolist (op '(< <= = >= >))
247                       (unless (eq (funcall op i f)
248                                   (funcall op i (rationalize f)))
249                         (error "(not (eq (~a ~a ~f) (~a ~a ~a)))~%"
250                                op i f
251                                op i (rationalize f)))
252                       (unless (eq (funcall op f i)
253                                   (funcall op (rationalize f) i))
254                         (error "(not (eq (~a ~f ~a) (~a ~a ~a)))~%"
255                                op f i
256                                op (rationalize f) i))))))))
257   (test 32 double-float negative)
258   (test 32 double-float positive)
259   (test 32 single-float negative)
260   (test 32 single-float positive))
261
262 ;; x86-64 sign-extension bug found using pfdietz's random tester.
263 (assert (= 286142502
264            (funcall (lambda ()
265                       (declare (notinline logxor))
266                       (min (logxor 0 0 0 286142502))))))
267
268 ;; Small bugs in LOGCOUNT can still allow SBCL to be built and thus go
269 ;; unnoticed, so check more thoroughly here.
270 (with-test (:name :logcount)
271   (flet ((test (x n)
272            (unless (= (logcount x) n)
273              (error "logcount failure for ~a" x))))
274     ;; Test with some patterns with well known number of ones/zeroes ...
275     (dotimes (i 128)
276       (let ((x (ash 1 i)))
277         (test x 1)
278         (test (- x) i)
279         (test (1- x) i)))
280     ;; ... and with some random integers of varying length.
281     (flet ((test-logcount (x)
282              (declare (type integer x))
283              (do ((result 0 (1+ result))
284                   (x (if (minusp x)
285                          (lognot x)
286                          x)
287                      (logand x (1- x))))
288                  ((zerop x) result))))
289       (dotimes (i 200)
290         (let ((x (random (ash 1 i))))
291           (test x (test-logcount x))
292           (test (- x) (test-logcount (- x))))))))
293
294 ;; 1.0 had a broken ATANH on win32
295 (with-test (:name :atanh)
296   (assert (= (atanh 0.9d0) 1.4722194895832204d0)))
297
298 ;; Test some cases of integer operations with constant arguments
299 (with-test (:name :constant-integers)
300   (labels ((test-forms (op x y header &rest forms)
301              (let ((val (funcall op x y)))
302                (dolist (form forms)
303                  (let ((new-val (funcall (compile nil (append header form)) x y)))
304                    (unless (eql val new-val)
305                      (error "~S /= ~S: ~S ~S ~S~%" val new-val (append header form) x y))))))
306            (test-case (op x y type)
307              (test-forms op x y `(lambda (x y &aux z)
308                                    (declare (type ,type x y)
309                                             (ignorable x y z)
310                                             (notinline identity)
311                                             (optimize speed (safety 0))))
312                          `((,op x ,y))
313                          `((setf z (,op x ,y))
314                            (identity x)
315                            z)
316                          `((values (,op x ,y) x))
317                          `((,op ,x y))
318                          `((setf z (,op ,x y))
319                            (identity y)
320                            z)
321                          `((values (,op ,x y) y))
322
323                          `((identity x)
324                            (,op x ,y))
325                          `((identity x)
326                            (setf z (,op x ,y))
327                            (identity x)
328                            z)
329                          `((identity x)
330                            (values (,op x ,y) x))
331                          `((identity y)
332                            (,op ,x y))
333                          `((identity y)
334                            (setf z (,op ,x y))
335                            (identity y)
336                            z)
337                          `((identity y)
338                            (values (,op ,x y) y))))
339            (test-op (op)
340              (let ((ub `(unsigned-byte ,sb-vm:n-word-bits))
341                    (sb `(signed-byte ,sb-vm:n-word-bits)))
342                (loop for (x y type)
343                      in `((2 1 fixnum)
344                           (2 1 ,ub)
345                           (2 1 ,sb)
346                           (,(1+ (ash 1 28)) ,(1- (ash 1 28)) fixnum)
347                           (,(+ 3 (ash 1 30)) ,(+ 2 (ash 1 30)) ,ub)
348                           (,(- -2 (ash 1 29)) ,(- 3 (ash 1 29)) ,sb)
349                           ,@(when (> sb-vm:n-word-bits 32)
350                               `((,(1+ (ash 1 29)) ,(1- (ash 1 29)) fixnum)
351                                 (,(1+ (ash 1 31)) ,(1- (ash 1 31)) ,ub)
352                                 (,(- -2 (ash 1 31)) ,(- 3 (ash 1 30)) ,sb)
353                                 (,(ash 1 40) ,(ash 1 39) fixnum)
354                                 (,(ash 1 40) ,(ash 1 39) ,ub)
355                                 (,(ash 1 40) ,(ash 1 39) ,sb)))
356                           ;; fixnums that can be represented as 32-bit
357                           ;; sign-extended immediates on x86-64
358                           ,@(when (and (> sb-vm:n-word-bits 32)
359                                        (< sb-vm:n-fixnum-tag-bits 3))
360                               `((,(1+ (ash 1 (- 31 sb-vm:n-fixnum-tag-bits)))
361                                  ,(1- (ash 1 (- 32 sb-vm:n-fixnum-tag-bits)))
362                                  fixnum))))
363                      do
364                   (test-case op x y type)
365                   (test-case op x x type)))))
366     (mapc #'test-op '(+ - * truncate
367                       < <= = >= >
368                       eql
369                       eq))))
370
371 ;; GCD used to sometimes return negative values. The following did, on 32 bit
372 ;; builds.
373 (with-test (:name :gcd)
374   ;; from lp#413680
375   (assert (plusp (gcd 20286123923750474264166990598656
376                       680564733841876926926749214863536422912)))
377   ;; from lp#516750
378   (assert (plusp (gcd 2596102012663483082521318626691873
379                       2596148429267413814265248164610048))))
380
381 (with-test (:name :expt-zero-zero)
382   ;; Check that (expt 0.0 0.0) and (expt 0 0.0) signal error, but (expt 0.0 0)
383   ;; returns 1.0
384   (assert (raises-error? (expt 0.0 0.0) sb-int:arguments-out-of-domain-error))
385   (assert (raises-error? (expt 0 0.0) sb-int:arguments-out-of-domain-error))
386   (assert (eql (expt 0.0 0) 1.0)))
387
388 (with-test (:name :multiple-constant-folding)
389   (let ((*random-state* (make-random-state t)))
390     (flet ((make-args ()
391              (let (args vars)
392                (loop repeat (1+ (random 12))
393                      do (if (zerop (random 2))
394                             (let ((var (gensym)))
395                               (push var args)
396                               (push var vars))
397                             (push (- (random 21) 10) args)))
398                (values args vars))))
399       (dolist (op '(+ * logior logxor logand logeqv gcd lcm - /))
400         (loop repeat 10
401               do (multiple-value-bind (args vars) (make-args)
402                    (let ((fast (compile nil `(lambda ,vars
403                                                (,op ,@args))))
404                          (slow (compile nil `(lambda ,vars
405                                                (declare (notinline ,op))
406                                                (,op ,@args)))))
407                      (loop repeat 3
408                            do (let* ((call-args (loop repeat (length vars)
409                                                       collect (- (random 21) 10)))
410                                      (fast-result (handler-case
411                                                       (apply fast call-args)
412                                                     (division-by-zero () :div0)))
413                                      (slow-result (handler-case
414                                                       (apply slow call-args)
415                                                     (division-by-zero () :div0))))
416                                 (if (eql fast-result slow-result)
417                                     (print (list :ok `(,op ,@args) :=> fast-result))
418                                     (error "oops: ~S, ~S" args call-args)))))))))))
419
420 ;;; (TRUNCATE <unsigned-word> <constant unsigned-word>) is optimized
421 ;;; to use multiplication instead of division. This propagates to FLOOR,
422 ;;; MOD and REM. Test that the transform is indeed triggered and test
423 ;;; several cases for correct results.
424 (with-test (:name (:integer-division-using-multiplication :used)
425                   :skipped-on '(not (or :x86-64 :x86)))
426   (dolist (fun '(truncate floor ceiling mod rem))
427     (let* ((foo (compile nil `(lambda (x)
428                                 (declare (optimize (speed 3)
429                                                    (space 1)
430                                                    (compilation-speed 0))
431                                          (type (unsigned-byte
432                                                 ,sb-vm:n-word-bits) x))
433                                 (,fun x 9))))
434            (disassembly (with-output-to-string (s)
435                           (disassemble foo :stream s))))
436       ;; KLUDGE copied from test :float-division-using-exact-reciprocal
437       ;; in compiler.pure.lisp.
438       (assert (and (not (search "DIV" disassembly))
439                    (search "MUL" disassembly))))))
440
441 (with-test (:name (:integer-division-using-multiplication :correctness))
442   (let ((*random-state* (make-random-state t)))
443     (dolist (dividend-type `((unsigned-byte ,sb-vm:n-word-bits)
444                              (and fixnum unsigned-byte)
445                              (integer 10000 10100)))
446       (dolist (divisor `(;; Some special cases from the paper
447                          7 10 14 641 274177
448                          ;; Range extremes
449                          3
450                          ,most-positive-fixnum
451                          ,(1- (expt 2 sb-vm:n-word-bits))
452                          ;; Some random values
453                          ,@(loop for i from 8 to sb-vm:n-word-bits
454                                  for r = (random (expt 2 i))
455                                  ;; We don't want 0, 1 and powers of 2.
456                                  when (not (zerop (logand r (1- r))))
457                                  collect r)))
458         (dolist (fun '(truncate ceiling floor mod rem))
459           (let ((foo (compile nil `(lambda (x)
460                                      (declare (optimize (speed 3)
461                                                         (space 1)
462                                                         (compilation-speed 0))
463                                               (type ,dividend-type x))
464                                      (,fun x ,divisor)))))
465             (dolist (dividend `(0 1 ,most-positive-fixnum
466                                 ,(1- divisor) ,divisor
467                                 ,(1- (* divisor 2)) ,(* divisor 2)
468                                 ,@(loop repeat 4
469                                         collect (+ 10000 (random 101)))
470                                 ,@(loop for i from 4 to sb-vm:n-word-bits
471                                         for pow = (expt 2 (1- i))
472                                         for r = (+ pow (random pow))
473                                         collect r)))
474               (when (typep dividend dividend-type)
475                 (multiple-value-bind (q1 r1)
476                     (funcall foo dividend)
477                   (multiple-value-bind (q2 r2)
478                       (funcall fun dividend divisor)
479                     (unless (and (= q1 q2)
480                                  (eql r1 r2))
481                       (error "bad results for ~s with dividend type ~s"
482                              (list fun dividend divisor)
483                              dividend-type))))))))))))
484
485 ;; The fast path for logbitp underestimated sb!vm:n-positive-fixnum-bits
486 ;; for > 61 bit fixnums.
487 (with-test (:name :logbitp-wide-fixnum)
488   (assert (not (logbitp (1- (integer-length most-positive-fixnum))
489                         most-negative-fixnum))))
490
491 ;; EXPT dispatches in a complicated way on the types of its arguments.
492 ;; Check that all possible combinations are covered.
493 (with-test (:name (:expt :argument-type-combinations))
494   (let ((numbers '(2                 ; fixnum
495                    3/5               ; ratio
496                    1.2f0             ; single-float
497                    2.0d0             ; double-float
498                    #c(3/5 1/7)       ; complex rational
499                    #c(1.2f0 1.3f0)   ; complex single-float
500                    #c(2.0d0 3.0d0))) ; complex double-float
501         (bignum (expt 2 64))
502         results)
503     (dolist (base (cons bignum numbers))
504       (dolist (power numbers)
505         (format t "(expt ~s ~s) => " base power)
506         (let ((result (expt base power)))
507           (format t "~s~%" result)
508           (push result results))))
509     (assert (every #'numberp results))))
510
511 (with-test (:name :bug-741564)
512   ;; The bug was that in (expt <fixnum> <(complex double-float)>) the
513   ;; calculation was partially done only to single-float precision,
514   ;; making the complex double-float result too unprecise. Some other
515   ;; combinations of argument types were affected, too; test that all
516   ;; of them are good to double-float precision.
517   (labels ((nearly-equal-p (x y)
518              "Are the arguments equal to nearly double-float precision?"
519              (declare (type double-float x y))
520              (< (/ (abs (- x y)) (abs y))
521                 (* double-float-epsilon 4))) ; Differences in the two least
522                                              ; significant mantissa bits
523                                              ; are OK.
524            (test-complex (x y)
525              (and (nearly-equal-p (realpart x) (realpart y))
526                   (nearly-equal-p (imagpart x) (imagpart y))))
527            (print-result (msg base power got expected)
528              (format t "~a (expt ~s ~s)~%got      ~s~%expected ~s~%"
529                      msg base power got expected)))
530     (let ((n-broken 0))
531       (flet ((test (base power coerce-to-type)
532                (let* ((got (expt base power))
533                       (expected (expt (coerce base coerce-to-type) power))
534                       (result (test-complex got expected)))
535                  (print-result (if result "Good:" "Bad:")
536                                base power got expected)
537                  (unless result
538                    (incf n-broken)))))
539         (dolist (base (list 2                    ; fixnum
540                             (expt 2 64)          ; bignum
541                             3/5                  ; ratio
542                             2.0f0))              ; single-float
543           (let ((power #c(-2.5d0 -4.5d0)))       ; complex double-float
544             (test base power 'double-float)))
545         (dolist (base (list #c(2.0f0 3.0f0)      ; complex single-float
546                             #c(2 3)              ; complex fixnum
547                             (complex (expt 2 64) (expt 2 65))
548                                                  ; complex bignum
549                             #c(3/5 1/7)))        ; complex ratio
550           (dolist (power (list #c(-2.5d0 -4.5d0) ; complex double-float
551                                -2.5d0))          ; double-float
552             (test base power '(complex double-float)))))
553       (when (> n-broken 0)
554         (error "Number of broken combinations: ~a" n-broken)))))
555
556 (with-test (:name (:ldb :rlwinm :ppc))
557   (let ((one (compile nil '(lambda (a) (ldb (byte 9 27) a))))
558         (two (compile nil '(lambda (a)
559                             (declare (type (integer -3 57216651) a))
560                             (ldb (byte 9 27) a)))))
561     (assert (= 0 (- (funcall one 10) (funcall two 10))))))
562
563 ;; The ISQRT implementation is sufficiently complicated that it should
564 ;; be tested.
565 (with-test (:name :isqrt)
566   (labels ((test (x)
567              (let* ((r (isqrt x))
568                     (r2 (expt r 2))
569                     (s2 (expt (1+ r) 2)))
570                (unless (and (<= r2 x)
571                             (> s2 x))
572                  (error "isqrt failure for ~a" x))))
573            (tests (x)
574              (test x)
575              (let ((x2 (expt x 2)))
576                (test x2)
577                (test (1+ x2))
578                (test (1- x2)))))
579     (loop for i from 1 to 200
580           for pow = (expt 2 (1- i))
581           for j = (+ pow (random pow))
582           do
583           (tests i)
584           (tests j))
585     (dotimes (i 10)
586       (tests (random (expt 2 (+ 1000 (random 10000))))))))