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))) |
---|