numbers-aux.lsp 11.3 KB
 pfdietz committed Apr 09, 2003 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) `````` pfdietz committed Apr 30, 2005 8 ``````(eval-when (:compile-toplevel :load-toplevel :execute) `````` pfdietz committed Sep 12, 2003 9 10 `````` (compile-and-load "random-aux.lsp")) `````` pfdietz committed Jun 22, 2004 11 12 13 14 15 16 ``````;;; 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." `````` sds committed Aug 12, 2009 17 `````` `````` pfdietz committed Jun 22, 2004 18 19 20 21 22 23 24 `````` (assert (functionp fn)) (assert (floatp lo)) (assert (floatp hi)) (assert (<= lo hi)) (assert (funcall fn hi)) (loop while (<= lo hi) `````` sds committed Aug 12, 2009 25 26 27 28 29 30 31 32 `````` 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)))))) `````` pfdietz committed Jun 22, 2004 33 `````` `````` pfdietz committed Dec 05, 2004 34 35 36 37 38 ``````(defun integer-binary-search (fn lo hi) "FN is a function that, if true for X, is true for all Y < X. Find the largest integer in [lo,hi) for which the function return true." `````` sds committed Aug 12, 2009 39 `````` `````` pfdietz committed Dec 05, 2004 40 41 42 43 44 45 46 `````` (assert (functionp fn)) (assert (integerp lo)) (assert (integerp hi)) (assert (<= lo hi)) (assert (funcall fn lo)) (loop while (< lo hi) `````` sds committed Aug 12, 2009 47 48 49 50 51 52 53 `````` do (let ((mid (ceiling (+ lo hi) 2))) (if (funcall fn mid) (setq lo mid) (if (= mid hi) (return lo) (setq hi mid)))) finally (return lo))) `````` pfdietz committed Dec 05, 2004 54 55 56 57 58 `````` (defun find-largest-exactly-floatable-integer (upper-bound) (integer-binary-search #'(lambda (i) (let* ((f (float i)) `````` sds committed Aug 12, 2009 59 60 61 `````` (i- (1- i)) (f- (float i-))) (and (= f i) (= f- i-)))) `````` pfdietz committed Dec 05, 2004 62 63 `````` 0 upper-bound)) `````` pfdietz committed Sep 05, 2003 64 65 66 67 68 ``````(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) `````` sds committed Aug 12, 2009 69 70 `````` (eqlzt (realpart x) (realpart y)) (eqlzt (imagpart x) (imagpart y)))) `````` pfdietz committed Sep 05, 2003 71 72 73 74 `````` ((zerop x) (eqlt (abs x) (abs y))) (t (eqlt x y)))) `````` pfdietz committed Feb 05, 2004 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 ``````(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)) `````` pfdietz committed Aug 23, 2003 104 105 106 107 ``````(defun numbers-are-compatible (x y) (cond ((complexp x) (and (numbers-are-compatible (realpart x) y) `````` sds committed Aug 12, 2009 108 `````` (numbers-are-compatible (imagpart x) y))) `````` pfdietz committed Aug 23, 2003 109 110 `````` ((complexp y) (and (numbers-are-compatible x (realpart y)) `````` sds committed Aug 12, 2009 111 `````` (numbers-are-compatible x (imagpart y)))) `````` pfdietz committed Aug 23, 2003 112 113 114 `````` (t (when (floatp x) (rotatef x y)) (or (floatp x) `````` sds committed Aug 12, 2009 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 `````` (not (floatp y)) (etypecase y (short-float (<= +rational-most-negative-short-float+ x +rational-most-positive-short-float+)) (single-float (<= +rational-most-negative-single-float+ x +rational-most-positive-single-float+)) (double-float (<= +rational-most-negative-double-float+ x +rational-most-positive-double-float+)) (long-float (<= +rational-most-negative-long-float+ x +rational-most-positive-long-float+))))))) `````` pfdietz committed Sep 06, 2003 133 134 135 136 137 `````` ;;; 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. `````` pfdietz committed Aug 23, 2003 138 `````` `````` pfdietz committed Apr 09, 2003 139 140 ``````(defun =.4-fn () (loop for x in *numbers* `````` sds committed Aug 12, 2009 141 142 143 144 145 `````` append (loop for y in *numbers* unless (or ;; (not (numbers-are-compatible x y)) (if (= x y) (= y x) (not (= y x)))) collect (list x y)))) `````` pfdietz committed Aug 03, 2003 146 147 148 `````` (defun /=.4-fn () (loop for x in *numbers* `````` sds committed Aug 12, 2009 149 150 151 152 153 `````` append (loop for y in *numbers* unless (or ;; (not (numbers-are-compatible x y)) (if (/= x y) (/= y x) (not (/= y x)))) collect (list x y)))) `````` pfdietz committed Aug 03, 2003 154 155 156 `````` (defun /=.4a-fn () (loop for x in *numbers* `````` sds committed Aug 12, 2009 157 158 159 160 161 162 163 `````` append (loop for y in *numbers* when (and ;; (numbers-are-compatible x y) (if (= x y) (/= x y) (not (/= x y)))) collect (list x y)))) `````` pfdietz committed Aug 03, 2003 164 165 `````` (defun <.8-fn () `````` pfdietz committed Aug 04, 2003 166 `````` (loop for x in *reals* `````` sds committed Aug 12, 2009 167 168 169 170 171 172 173 174 `````` nconc (loop for y in *reals* when (handler-case (and ;; (numbers-are-compatible x y) (and (< x y) (> x y))) (arithmetic-error () nil)) collect (list x y)))) `````` pfdietz committed Aug 03, 2003 175 176 `````` (defun <.9-fn () `````` pfdietz committed Aug 04, 2003 177 `````` (loop for x in *reals* `````` sds committed Aug 12, 2009 178 179 180 181 182 183 184 185 186 `````` nconc (loop for y in *reals* when (handler-case (and ;; (numbers-are-compatible x y) (if (< x y) (not (> y x)) (> y x))) (arithmetic-error () nil)) collect (list x y)))) `````` pfdietz committed Aug 03, 2003 187 188 `````` (defun <.10-fn () `````` pfdietz committed Aug 04, 2003 189 `````` (loop for x in *reals* `````` sds committed Aug 12, 2009 190 191 192 193 194 195 196 197 198 `````` nconc (loop for y in *reals* when (handler-case (and ;; (numbers-are-compatible x y) (if (< x y) (>= x y) (not (>= x y)))) (arithmetic-error () nil)) collect (list x y)))) `````` pfdietz committed Aug 03, 2003 199 200 `````` (defun <=.8-fn () `````` pfdietz committed Aug 04, 2003 201 `````` (loop for x in *reals* `````` sds committed Aug 12, 2009 202 203 204 205 206 207 208 209 210 211 `````` nconc (loop for y in *reals* when (handler-case (and ;; (numbers-are-compatible x y) (if (<= x y) (not (>= y x)) (>= y x))) (arithmetic-error () nil)) collect (list x y)))) `````` pfdietz committed Aug 03, 2003 212 ``````(defun <=.9-fn () `````` pfdietz committed Aug 04, 2003 213 `````` (loop for x in *reals* `````` sds committed Aug 12, 2009 214 215 216 217 218 219 220 221 222 `````` nconc (loop for y in *reals* 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)) collect (list x y)))) `````` pfdietz committed Aug 03, 2003 223 224 `````` (defun >.8-fn () `````` pfdietz committed Aug 04, 2003 225 `````` (loop for x in *reals* `````` sds committed Aug 12, 2009 226 227 228 229 230 231 232 233 234 `````` nconc (loop for y in *reals* when (handler-case (and ;; (numbers-are-compatible x y) (if (> x y) (<= x y) (not (<= x y)))) (arithmetic-error () nil)) collect (list x y)))) `````` pfdietz committed Aug 03, 2003 235 236 `````` (defun >=.8-fn () `````` pfdietz committed Aug 04, 2003 237 `````` (loop for x in *reals* `````` sds committed Aug 12, 2009 238 239 240 241 242 243 244 245 246 `````` nconc (loop for y in *reals* 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)) collect (list x y)))) `````` pfdietz committed Aug 03, 2003 247 248 249 250 251 `````` ;;; Comparison of rationsls (defun compare-random-rationals (n m rep) (loop for a = (- (random n) m) `````` sds committed Aug 12, 2009 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 `````` 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))) `````` pfdietz committed Aug 03, 2003 273 `````` `````` pfdietz committed Aug 04, 2003 274 275 ``````(defun max.2-fn () (loop for x in *reals* `````` sds committed Aug 12, 2009 276 277 278 279 280 281 282 283 284 285 286 `````` nconc (loop for y in *reals* when (numbers-are-compatible x y) 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)) collect (list x y (max x y))))) `````` pfdietz committed Aug 05, 2003 287 288 289 `````` (defun min.2-fn () (loop for x in *reals* `````` sds committed Aug 12, 2009 290 291 292 293 294 295 296 297 298 299 300 `````` nconc (loop for y in *reals* when (numbers-are-compatible x y) 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)) collect (list x y (min x y))))) `````` pfdietz committed Aug 09, 2003 301 `````` `````` pfdietz committed Aug 31, 2003 302 303 304 305 306 307 308 309 ``````;;; 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 > -epsilon. (defun find-epsilon (x) (assert (floatp x)) (let* ((one (float 1 x)) `````` sds committed Aug 12, 2009 310 311 `````` (radix (float-radix one)) (eps (/ 1 radix))) `````` pfdietz committed Aug 31, 2003 312 313 314 315 316 317 `````` (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)))))) `````` pfdietz committed Sep 07, 2003 318 `````` `````` 319 ``````(defun test-log-op-with-decls (op xlo xhi ylo yhi niters `````` sds committed Aug 12, 2009 320 321 322 `````` &optional (decls '((optimize (speed 3) (safety 1) (debug 1))))) `````` 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 `````` "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 `````` sds committed Aug 12, 2009 338 339 340 341 342 343 `````` `(lambda (x y) (declare (type (integer ,xlo ,xhi) x) (type (integer ,ylo ,yhi) y) ,@ decls) (,op x y))) (fn (compile nil source))) `````` 344 `````` (loop for i below niters `````` sds committed Aug 12, 2009 345 346 347 348 349 `````` for x = (random-from-interval (1+ xhi) xlo) for y = (random-from-interval (1+ yhi) ylo) unless (eql (funcall (the symbol op) x y) (funcall fn x y)) collect (list x y)))) `````` 350 351 352 `````` (defun test-log-op (op n1 n2) (flet ((%r () (let ((r (random 33))) `````` sds committed Aug 12, 2009 353 `````` (- (random (ash 1 (1+ r))) (ash 1 r))))) `````` 354 `````` (loop for x1 = (%r) `````` sds committed Aug 12, 2009 355 356 357 358 359 360 361 362 363 `````` 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)))) `````` pfdietz committed Feb 10, 2004 364 365 366 367 368 369 ``````(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)))``````