source: trunk/source/tests/ansi-tests/numbers-aux.lsp @ 8991

Last change on this file since 8991 was 8991, checked in by gz, 11 years ago

Check in the gcl ansi test suite (original, in preparation for making local changes)

File size: 9.6 KB
Line 
1;-*- Mode:     Lisp -*-
2;;;; Author:   Paul Dietz
3;;;; Created:  Mon Apr  7 07:24:43 2003
4;;;; Contains: Auxiliary functions for number tests
5
6(in-package :cl-test)
7
8(eval-when (:compile-toplevel :load-toplevel :execute)
9  (compile-and-load "random-aux.lsp"))
10
11;;; Binary search on reals
12
13(defun float-binary-search (fn lo hi)
14  "FN is a function that, if true for X, is true for all Y > X.
15   Find the smallest float in [lo,hi] for which the function
16   return true."
17 
18  (assert (functionp fn))
19  (assert (floatp lo))
20  (assert (floatp hi))
21  (assert (<= lo hi))
22  (assert (funcall fn hi))
23
24  (loop while (<= lo hi)
25        do (let ((mid (/ (+ lo hi) 2)))
26             (if (funcall fn mid)
27                 (if (= mid hi)
28                     (return hi)
29                   (setq hi mid))
30               (if (= mid lo)
31                   (return hi)
32                 (setq lo mid))))))
33
34(defun integer-binary-search (fn lo hi)
35
36  "FN is a function that, if true for X, is true for all Y < X.
37   Find the largest integer in [lo,hi) for which the function
38   return true."
39 
40  (assert (functionp fn))
41  (assert (integerp lo))
42  (assert (integerp hi))
43  (assert (<= lo hi))
44  (assert (funcall fn lo))
45
46  (loop while (< lo hi)
47        do (let ((mid (ceiling (+ lo hi) 2)))
48             (if (funcall fn mid)
49                 (setq lo mid)
50               (if (= mid hi)
51                   (return lo)
52                 (setq hi mid))))
53        finally (return lo)))
54
55(defun find-largest-exactly-floatable-integer (upper-bound)
56  (integer-binary-search
57   #'(lambda (i)
58       (let* ((f  (float i))
59              (i- (1- i))
60              (f- (float i-)))
61         (and (= f i) (= f- i-))))
62   0 upper-bound))
63
64(defun eqlzt (x y)
65  "Return T if (eql x y) or if both are zero of the same type."
66  (cond
67   ((complexp x)
68    (and (complexp y)
69         (eqlzt (realpart x) (realpart y))
70         (eqlzt (imagpart x) (imagpart y))))
71   ((zerop x)
72    (eqlt (abs x) (abs y)))
73   (t (eqlt x y))))
74
75(defconstant +rational-most-negative-short-float+
76  (rational-safely most-negative-short-float))
77
78(defconstant +rational-most-negative-single-float+
79  (rational-safely most-negative-single-float))
80
81(defconstant +rational-most-negative-double-float+
82  (rational-safely most-negative-double-float))
83
84(defconstant +rational-most-negative-long-float+
85  (rational-safely most-negative-long-float))
86
87(defconstant +rational-most-positive-short-float+
88  (rational-safely most-positive-short-float))
89
90(defconstant +rational-most-positive-single-float+
91  (rational-safely most-positive-single-float))
92
93(defconstant +rational-most-positive-double-float+
94  (rational-safely most-positive-double-float))
95
96(defconstant +rational-most-positive-long-float+
97  (rational-safely most-positive-long-float))
98
99(defun float-exponent (x)
100  (if (floatp x)
101      (nth-value 1 (decode-float x))
102    0))
103
104(defun numbers-are-compatible (x y)
105  (cond
106   ((complexp x)
107    (and (numbers-are-compatible (realpart x) y)
108         (numbers-are-compatible (imagpart x) y)))
109   ((complexp y)
110    (and (numbers-are-compatible x (realpart y))
111         (numbers-are-compatible x (imagpart y))))
112   (t
113    (when (floatp x) (rotatef x y))
114    (or (floatp x)
115        (not (floatp y))
116        (etypecase y
117          (short-float
118           (<= +rational-most-negative-short-float+
119               x
120               +rational-most-positive-short-float+))
121          (single-float
122           (<= +rational-most-negative-single-float+
123               x
124               +rational-most-positive-single-float+))
125          (double-float
126           (<= +rational-most-negative-double-float+
127               x
128               +rational-most-positive-double-float+))
129          (long-float
130           (<= +rational-most-negative-long-float+
131               x
132               +rational-most-positive-long-float+)))))))
133
134;;; NOTE!  According to section 12.1.4.1, when a rational is compared
135;;; to a float, the effect is as if the float is convert to a rational
136;;; (by RATIONAL), not as if the rational is converted to a float.
137;;; This means the calls to numbers-are-compatible are not necessary.
138
139(defun =.4-fn ()
140  (loop for x in *numbers*
141        append
142        (loop for y in *numbers*
143              unless (or ;; (not (numbers-are-compatible x y))
144                         (if (= x y) (= y x) (not (= y x))))
145              collect (list x y))))
146
147(defun /=.4-fn ()
148  (loop for x in *numbers*
149        append
150        (loop for y in *numbers*
151              unless (or ;; (not (numbers-are-compatible x y))
152                         (if (/= x y) (/= y x) (not (/= y x))))
153              collect (list x y))))
154
155(defun /=.4a-fn ()
156  (loop for x in *numbers*
157        append
158        (loop for y in *numbers*
159              when (and ;; (numbers-are-compatible x y)
160                        (if (= x y)
161                            (/= x y)
162                          (not (/= x y))))
163              collect (list x y))))
164
165(defun <.8-fn ()
166  (loop for x in *reals*
167        nconc
168        (loop for y in *reals*
169              when
170              (handler-case
171               (and ;; (numbers-are-compatible x y)
172                (and (< x y) (> x y)))
173               (arithmetic-error () nil))
174              collect (list x y))))
175
176(defun <.9-fn ()
177  (loop for x in *reals*
178        nconc
179        (loop for y in *reals*
180              when
181              (handler-case
182               (and ;; (numbers-are-compatible x y)
183                (if (< x y) (not (> y x))
184                  (> y x)))
185               (arithmetic-error () nil))
186              collect (list x y))))
187
188(defun <.10-fn ()
189  (loop for x in *reals*
190        nconc
191        (loop for y in *reals*
192              when
193              (handler-case
194               (and ;; (numbers-are-compatible x y)
195                (if (< x y) (>= x y)
196                  (not (>= x y))))
197               (arithmetic-error () nil))
198              collect (list x y))))
199
200(defun <=.8-fn ()
201  (loop for x in *reals*
202        nconc
203        (loop for y in *reals*
204              when
205              (handler-case
206               (and ;; (numbers-are-compatible x y)
207                (if (<= x y) (not (>= y x))
208                  (>= y x)))
209               (arithmetic-error () nil))
210              collect (list x y))))
211 
212(defun <=.9-fn ()
213  (loop for x in *reals*
214        nconc
215        (loop for y in *reals*
216              when
217              (handler-case
218               (and ;; (numbers-are-compatible x y)
219                (if (<= x y) (not (or (= x y) (< x y)))
220                  (or (= x y) (< x y))))
221               (arithmetic-error () nil))
222              collect (list x y))))
223
224(defun >.8-fn ()
225  (loop for x in *reals*
226        nconc
227        (loop for y in *reals*
228              when
229              (handler-case
230               (and ;; (numbers-are-compatible x y)
231                (if (> x y) (<= x y)
232                  (not (<= x y))))
233               (arithmetic-error () nil))
234              collect (list x y))))
235
236(defun >=.8-fn ()
237  (loop for x in *reals*
238        nconc
239        (loop for y in *reals*
240              when
241              (handler-case
242               (and ;; (numbers-are-compatible x y)
243                (if (>= x y) (not (or (= x y) (> x y)))
244                  (or (= x y) (> x y))))
245               (arithmetic-error () nil))
246              collect (list x y))))
247
248;;; Comparison of rationsls
249
250(defun compare-random-rationals (n m rep)
251  (loop for a = (- (random n) m)
252        for b = (- (random n) m)
253        for c = (- (random n) m)
254        for d = (- (random n) m)
255        repeat rep
256        when
257        (and (/= b 0)
258             (/= d 0)
259             (let ((q1 (/ a b))
260                   (q2 (/ c d))
261                   (ad (* a d))
262                   (bc (* b c)))
263               (when (< (* b d) 0)
264                 (setq ad (- ad))
265                 (setq bc (- bc)))
266               (or (if (< q1 q2) (not (< ad bc)) (< ad bc))
267                   (if (<= q1 q2) (not (<= ad bc)) (<= ad bc))
268                   (if (> q1 q2) (not (> ad bc)) (> ad bc))
269                   (if (>= q1 q2) (not (>= ad bc)) (>= ad bc))
270                   (if (= q1 q2) (not (= ad bc)) (= ad bc))
271                   (if (/= q1 q2) (not (/= ad bc)) (/= ad bc)))))
272        collect (list a b c d)))
273
274(defun max.2-fn ()
275  (loop for x in *reals*
276        nconc
277        (loop for y in *reals*
278              when (numbers-are-compatible x y)
279              unless
280              (handler-case
281               (let ((m (max x y)))
282                 (and (>= m x) (>= m y)
283                      (or (= m x) (= m y))))
284               (floating-point-underflow () t)
285               (floating-point-overflow () t))
286              collect (list x y (max x y)))))
287
288(defun min.2-fn ()
289  (loop for x in *reals*
290        nconc
291        (loop for y in *reals*
292              when (numbers-are-compatible x y)
293              unless
294              (handler-case
295               (let ((m (min x y)))
296                 (and (<= m x) (<= m y)
297                      (or (= m x) (= m y))))
298               (floating-point-underflow () t)
299               (floating-point-overflow () t))
300              collect (list x y (min x y)))))
301
302;;; Compute the number of digits that can be added to 1.0 in the appropriate
303;;; float type, a rational representation of the smallest radix^(-k) s.t.
304;;; 1.0 + radix^(-k) /= 1.0, and the float representation of that value.
305;;; Note that this will in general be > <float-type>-epsilon.
306
307(defun find-epsilon (x)
308  (assert (floatp x))
309  (let* ((one (float 1 x))
310         (radix (float-radix one))
311         (eps (/ 1 radix)))
312    (loop
313     for next-eps = (/ eps radix)
314     for i from 1
315     until (eql one (+ one next-eps))
316     do (setq eps next-eps)
317     finally (return (values i eps (float eps one))))))
318
319(defun test-log-op-with-decls (op xlo xhi ylo yhi niters
320                                  &optional
321                                  (decls '((optimize (speed 3) (safety 1)
322                                                     (debug 1)))))
323  "Test that a compiled form of the LOG* function OP computes
324   the expected result on two random integers drawn from the
325   types `(integer ,xlo ,xhi) and `(integer ,ylo ,yhi).  Try
326   niters choices.  Return a list of pairs on which the test fails."
327
328  (assert (symbolp op))
329  (assert (integerp xlo))
330  (assert (integerp xhi))
331  (assert (integerp ylo))
332  (assert (integerp yhi))
333  (assert (integerp niters))
334  (assert (<= xlo xhi))
335  (assert (<= ylo yhi))
336
337  (let* ((source
338          `(lambda (x y)
339             (declare (type (integer ,xlo ,xhi) x)
340                      (type (integer ,ylo ,yhi) y)
341                      ,@ decls)
342             (,op x y)))
343         (fn (compile nil source)))
344    (loop for i below niters
345          for x = (random-from-interval (1+ xhi) xlo)
346          for y = (random-from-interval (1+ yhi) ylo)
347          unless (eql (funcall (the symbol op) x y)
348                      (funcall fn x y))
349          collect (list x y))))
350
351(defun test-log-op (op n1 n2)
352  (flet ((%r () (let ((r (random 33)))
353                  (- (random (ash 1 (1+ r))) (ash 1 r)))))
354    (loop for x1 = (%r)
355          for x2 = (%r)
356          for y1 = (%r)
357          for y2 = (%r)
358          repeat n1
359          nconc
360          (test-log-op-with-decls op
361                                  (min x1 x2) (max x1 x2)
362                                  (min y1 y2) (max y1 y2)
363                                  n2))))
364(defun safe-tan (x &optional (default 0.0))
365  (handler-case
366   (let ((result (multiple-value-list (tan x))))
367     (assert (null (cdr result)))
368     (car result))
369   (arithmetic-error () default)))
Note: See TracBrowser for help on using the repository browser.