numbers-aux.lsp 9.54 KB
Newer Older
1
2
3
4
5
6
7
;-*- Mode:     Lisp -*-
;;;; Author:   Paul Dietz
;;;; Created:  Mon Apr  7 07:24:43 2003
;;;; Contains: Auxiliary functions for number tests

(in-package :cl-test)

8
9
10
(eval-when (load eval)
  (compile-and-load "random-aux.lsp"))

11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
;;; Binary search on reals

(defun float-binary-search (fn lo hi)
  "FN is a function that, if true for X, is true for all Y > X.
   Find the smallest float in [lo,hi] for which the function
   return true."
  
  (assert (functionp fn))
  (assert (floatp lo))
  (assert (floatp hi))
  (assert (<= lo hi))
  (assert (funcall fn hi))

  (loop while (<= lo hi)
	do (let ((mid (/ (+ lo hi) 2)))
	     (if (funcall fn mid)
		 (if (= mid hi)
		     (return hi)
		   (setq hi mid))
	       (if (= mid lo)
		   (return hi)
		 (setq lo mid))))))

34
35
36
37
38
39
40
41
42
43
44
(defun eqlzt (x y)
  "Return T if (eql x y) or if both are zero of the same type."
  (cond
   ((complexp x)
    (and (complexp y)
	 (eqlzt (realpart x) (realpart y))
	 (eqlzt (imagpart x) (imagpart y))))
   ((zerop x)
    (eqlt (abs x) (abs y)))
   (t (eqlt x y))))

45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
(defconstant +rational-most-negative-short-float+
  (rational-safely most-negative-short-float))

(defconstant +rational-most-negative-single-float+
  (rational-safely most-negative-single-float))

(defconstant +rational-most-negative-double-float+
  (rational-safely most-negative-double-float))

(defconstant +rational-most-negative-long-float+
  (rational-safely most-negative-long-float))

(defconstant +rational-most-positive-short-float+
  (rational-safely most-positive-short-float))

(defconstant +rational-most-positive-single-float+
  (rational-safely most-positive-single-float))

(defconstant +rational-most-positive-double-float+
  (rational-safely most-positive-double-float))

(defconstant +rational-most-positive-long-float+
  (rational-safely most-positive-long-float))

(defun float-exponent (x)
  (if (floatp x)
      (nth-value 1 (decode-float x))
    0))

74
75
76
77
78
79
80
81
82
83
84
85
86
87
(defun numbers-are-compatible (x y)
  (cond
   ((complexp x)
    (and (numbers-are-compatible (realpart x) y)
	 (numbers-are-compatible (imagpart x) y)))
   ((complexp y)
    (and (numbers-are-compatible x (realpart y))
	 (numbers-are-compatible x (imagpart y))))
   (t
    (when (floatp x) (rotatef x y))
    (or (floatp x)
	(not (floatp y))
	(etypecase y
	  (short-float
88
	   (<= +rational-most-negative-short-float+
89
	       x
90
	       +rational-most-positive-short-float+))
91
	  (single-float
92
	   (<= +rational-most-negative-single-float+
93
	       x
94
	       +rational-most-positive-single-float+))
95
	  (double-float
96
	   (<= +rational-most-negative-double-float+
97
	       x
98
	       +rational-most-positive-double-float+))
99
	  (long-float
100
	   (<= +rational-most-negative-long-float+
101
	       x
102
	       +rational-most-positive-long-float+)))))))
103
104
105
106
107

;;; NOTE!  According to section 12.1.4.1, when a rational is compared
;;; to a float, the effect is as if the float is convert to a rational
;;; (by RATIONAL), not as if the rational is converted to a float.
;;; This means the calls to numbers-are-compatible are not necessary.
108

109
110
111
112
(defun =.4-fn ()
  (loop for x in *numbers*
	append
	(loop for y in *numbers*
113
	      unless (or ;; (not (numbers-are-compatible x y))
114
			 (if (= x y) (= y x) (not (= y x))))
115
	      collect (list x y))))
pfdietz's avatar
pfdietz committed
116
117
118
119
120

(defun /=.4-fn ()
  (loop for x in *numbers*
	append
	(loop for y in *numbers*
121
	      unless (or ;; (not (numbers-are-compatible x y))
122
			 (if (/= x y) (/= y x) (not (/= y x))))
pfdietz's avatar
pfdietz committed
123
124
125
126
127
128
	      collect (list x y))))

(defun /=.4a-fn ()
  (loop for x in *numbers*
	append
	(loop for y in *numbers*
129
	      when (and ;; (numbers-are-compatible x y)
130
131
132
			(if (= x y)
			    (/= x y)
			  (not (/= x y))))
pfdietz's avatar
pfdietz committed
133
134
135
	      collect (list x y))))

(defun <.8-fn ()
136
  (loop for x in *reals*
pfdietz's avatar
pfdietz committed
137
	nconc
138
	(loop for y in *reals*
139
140
141
142
143
	      when
	      (handler-case
	       (and ;; (numbers-are-compatible x y)
		(and (< x y) (> x y)))
	       (arithmetic-error () nil))
pfdietz's avatar
pfdietz committed
144
145
146
	      collect (list x y))))

(defun <.9-fn ()
147
  (loop for x in *reals*
pfdietz's avatar
pfdietz committed
148
	nconc
149
	(loop for y in *reals*
150
151
152
153
154
155
	      when
	      (handler-case
	       (and ;; (numbers-are-compatible x y)
		(if (< x y) (not (> y x))
		  (> y x)))
	       (arithmetic-error () nil))
pfdietz's avatar
pfdietz committed
156
157
158
	      collect (list x y))))

(defun <.10-fn ()
159
  (loop for x in *reals*
pfdietz's avatar
pfdietz committed
160
	nconc
161
	(loop for y in *reals*
162
163
164
165
166
167
	      when
	      (handler-case
	       (and ;; (numbers-are-compatible x y)
		(if (< x y) (>= x y)
		  (not (>= x y))))
	       (arithmetic-error () nil))
pfdietz's avatar
pfdietz committed
168
169
170
	      collect (list x y))))

(defun <=.8-fn ()
171
  (loop for x in *reals*
pfdietz's avatar
pfdietz committed
172
	nconc
173
	(loop for y in *reals*
174
175
176
177
178
179
	      when
	      (handler-case
	       (and ;; (numbers-are-compatible x y)
		(if (<= x y) (not (>= y x))
		  (>= y x)))
	       (arithmetic-error () nil))
pfdietz's avatar
pfdietz committed
180
181
182
	      collect (list x y))))
 
(defun <=.9-fn ()
183
  (loop for x in *reals*
pfdietz's avatar
pfdietz committed
184
	nconc
185
	(loop for y in *reals*
186
187
188
189
190
191
	      when
	      (handler-case
	       (and ;; (numbers-are-compatible x y)
		(if (<= x y) (not (or (= x y) (< x y)))
		  (or (= x y) (< x y))))
	       (arithmetic-error () nil))
pfdietz's avatar
pfdietz committed
192
193
194
	      collect (list x y))))

(defun >.8-fn ()
195
  (loop for x in *reals*
pfdietz's avatar
pfdietz committed
196
	nconc
197
	(loop for y in *reals*
198
199
200
201
202
203
	      when
	      (handler-case
	       (and ;; (numbers-are-compatible x y)
		(if (> x y) (<= x y)
		  (not (<= x y))))
	       (arithmetic-error () nil))
pfdietz's avatar
pfdietz committed
204
205
206
	      collect (list x y))))

(defun >=.8-fn ()
207
  (loop for x in *reals*
pfdietz's avatar
pfdietz committed
208
	nconc
209
	(loop for y in *reals*
210
211
212
213
214
215
	      when
	      (handler-case
	       (and ;; (numbers-are-compatible x y)
		(if (>= x y) (not (or (= x y) (> x y)))
		  (or (= x y) (> x y))))
	       (arithmetic-error () nil))
pfdietz's avatar
pfdietz committed
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
	      collect (list x y))))

;;; Comparison of rationsls

(defun compare-random-rationals (n m rep)
  (loop for a = (- (random n) m)
	for b = (- (random n) m)
	for c = (- (random n) m)
	for d = (- (random n) m)
	repeat rep
	when
	(and (/= b 0)
	     (/= d 0)
	     (let ((q1 (/ a b))
		   (q2 (/ c d))
		   (ad (* a d))
		   (bc (* b c)))
	       (when (< (* b d) 0)
		 (setq ad (- ad))
		 (setq bc (- bc)))
	       (or (if (< q1 q2) (not (< ad bc)) (< ad bc))
		   (if (<= q1 q2) (not (<= ad bc)) (<= ad bc))
		   (if (> q1 q2) (not (> ad bc)) (> ad bc))
		   (if (>= q1 q2) (not (>= ad bc)) (>= ad bc))
		   (if (= q1 q2) (not (= ad bc)) (= ad bc))
		   (if (/= q1 q2) (not (/= ad bc)) (/= ad bc)))))
	collect (list a b c d)))

244
245
246
247
(defun max.2-fn ()
  (loop for x in *reals*
	nconc
	(loop for y in *reals*
248
	      when (numbers-are-compatible x y)
249
250
251
252
253
254
255
	      unless
	      (handler-case
	       (let ((m (max x y)))
		 (and (>= m x) (>= m y)
		      (or (= m x) (= m y))))
	       (floating-point-underflow () t)
	       (floating-point-overflow () t))
256
	      collect (list x y (max x y)))))
pfdietz's avatar
pfdietz committed
257
258
259
260
261

(defun min.2-fn ()
  (loop for x in *reals*
	nconc
	(loop for y in *reals*
262
	      when (numbers-are-compatible x y)
263
264
265
266
267
268
269
	      unless
	      (handler-case
	       (let ((m (min x y)))
		 (and (<= m x) (<= m y)
		      (or (= m x) (= m y))))
	       (floating-point-underflow () t)
	       (floating-point-overflow () t))
270
	      collect (list x y (min x y)))))
271
272
273

(defun epsilon (number)
  (etypecase number
274
    (complex (* 2 (epsilon (realpart number)))) ;; crude
275
276
277
278
279
280
281
282
    (short-float short-float-epsilon)
    (single-float single-float-epsilon)
    (double-float double-float-epsilon)
    (long-float long-float-epsilon)
    (rational 0)))

(defun negative-epsilon (number)
  (etypecase number
283
    (complex (* 2 (negative-epsilon (realpart number)))) ;; crude
284
285
286
287
288
289
    (short-float short-float-negative-epsilon)
    (single-float single-float-negative-epsilon)
    (double-float double-float-negative-epsilon)
    (long-float long-float-negative-epsilon)
    (rational 0)))

290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
;;; Compute the number of digits that can be added to 1.0 in the appropriate
;;; float type, a rational representation of the smallest radix^(-k) s.t.
;;; 1.0 + radix^(-k) /= 1.0, and the float representation of that value.
;;; Note that this will in general be > <float-type>-epsilon.

(defun find-epsilon (x)
  (assert (floatp x))
  (let* ((one (float 1 x))
	 (radix (float-radix one))
	 (eps (/ 1 radix)))
    (loop
     for next-eps = (/ eps radix)
     for i from 1
     until (eql one (+ one next-eps))
     do (setq eps next-eps)
     finally (return (values i eps (float eps one))))))
306
307
308
309

;;; Approximate equality function
(defun approx= (x y &optional (eps (epsilon x)))
  (<= (abs (/ (- x y) (max (abs x) 1))) eps))
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336

(defun test-log-op-with-decls (op xlo xhi ylo yhi niters
				  &optional
				  (decls '((optimize (speed 3) (safety 1)
						     (debug 1)))))
  "Test that a compiled form of the LOG* function OP computes
   the expected result on two random integers drawn from the
   types `(integer ,xlo ,xhi) and `(integer ,ylo ,yhi).  Try
   niters choices.  Return a list of pairs on which the test fails."

  (assert (symbolp op))
  (assert (integerp xlo))
  (assert (integerp xhi))
  (assert (integerp ylo))
  (assert (integerp yhi))
  (assert (integerp niters))
  (assert (<= xlo xhi))
  (assert (<= ylo yhi))

  (let* ((source
	  `(lambda (x y)
	     (declare (type (integer ,xlo ,xhi) x)
		      (type (integer ,ylo ,yhi) y)
		      ,@ decls)
	     (,op x y)))
	 (fn (compile nil source)))
    (loop for i below niters
337
338
	  for x = (random-from-interval (1+ xhi) xlo)
	  for y = (random-from-interval (1+ yhi) ylo)
339
340
341
	  unless (eql (funcall (the symbol op) x y)
		      (funcall fn x y))
	  collect (list x y))))
342
343
344
345
346
347
348
349
350
351
352
353
354
355

(defun test-log-op (op n1 n2)
  (flet ((%r () (let ((r (random 33)))
		  (- (random (ash 1 (1+ r))) (ash 1 r)))))
    (loop for x1 = (%r)
	  for x2 = (%r)
	  for y1 = (%r)
	  for y2 = (%r)
	  repeat n1
	  nconc
	  (test-log-op-with-decls op
				  (min x1 x2) (max x1 x2)
				  (min y1 y2) (max y1 y2)
				  n2))))
356
357
358
359
360
361
(defun safe-tan (x &optional (default 0.0))
  (handler-case
   (let ((result (multiple-value-list (tan x))))
     (assert (null (cdr result)))
     (car result))
   (arithmetic-error () default)))