numbers-aux.lsp 11.3 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
(eval-when (:compile-toplevel :load-toplevel :execute)
9 10
  (compile-and-load "random-aux.lsp"))

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's avatar
sds committed
17

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's avatar
sds committed
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))))))
33

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's avatar
sds committed
39

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's avatar
sds committed
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)))
54 55 56 57 58

(defun find-largest-exactly-floatable-integer (upper-bound)
  (integer-binary-search
   #'(lambda (i)
       (let* ((f  (float i))
sds's avatar
sds committed
59 60 61
              (i- (1- i))
              (f- (float i-)))
         (and (= f i) (= f- i-))))
62 63
   0 upper-bound))

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's avatar
sds committed
69 70
         (eqlzt (realpart x) (realpart y))
         (eqlzt (imagpart x) (imagpart y))))
71 72 73 74
   ((zerop x)
    (eqlt (abs x) (abs y)))
   (t (eqlt x y))))

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

104 105 106 107
(defun numbers-are-compatible (x y)
  (cond
   ((complexp x)
    (and (numbers-are-compatible (realpart x) y)
sds's avatar
sds committed
108
         (numbers-are-compatible (imagpart x) y)))
109 110
   ((complexp y)
    (and (numbers-are-compatible x (realpart y))
sds's avatar
sds committed
111
         (numbers-are-compatible x (imagpart y))))
112 113 114
   (t
    (when (floatp x) (rotatef x y))
    (or (floatp x)
sds's avatar
sds committed
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+)))))))
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.
138

139 140
(defun =.4-fn ()
  (loop for x in *numbers*
sds's avatar
sds committed
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's avatar
pfdietz committed
146 147 148

(defun /=.4-fn ()
  (loop for x in *numbers*
sds's avatar
sds committed
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's avatar
pfdietz committed
154 155 156

(defun /=.4a-fn ()
  (loop for x in *numbers*
sds's avatar
sds committed
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's avatar
pfdietz committed
164 165

(defun <.8-fn ()
166
  (loop for x in *reals*
sds's avatar
sds committed
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's avatar
pfdietz committed
175 176

(defun <.9-fn ()
177
  (loop for x in *reals*
sds's avatar
sds committed
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's avatar
pfdietz committed
187 188

(defun <.10-fn ()
189
  (loop for x in *reals*
sds's avatar
sds committed
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's avatar
pfdietz committed
199 200

(defun <=.8-fn ()
201
  (loop for x in *reals*
sds's avatar
sds committed
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's avatar
pfdietz committed
212
(defun <=.9-fn ()
213
  (loop for x in *reals*
sds's avatar
sds committed
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's avatar
pfdietz committed
223 224

(defun >.8-fn ()
225
  (loop for x in *reals*
sds's avatar
sds committed
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's avatar
pfdietz committed
235 236

(defun >=.8-fn ()
237
  (loop for x in *reals*
sds's avatar
sds committed
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's avatar
pfdietz committed
247 248 249 250 251

;;; Comparison of rationsls

(defun compare-random-rationals (n m rep)
  (loop for a = (- (random n) m)
sds's avatar
sds committed
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's avatar
pfdietz committed
273

274 275
(defun max.2-fn ()
  (loop for x in *reals*
sds's avatar
sds committed
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's avatar
pfdietz committed
287 288 289

(defun min.2-fn ()
  (loop for x in *reals*
sds's avatar
sds committed
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)))))
301

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 > <float-type>-epsilon.

(defun find-epsilon (x)
  (assert (floatp x))
  (let* ((one (float 1 x))
sds's avatar
sds committed
310 311
         (radix (float-radix one))
         (eps (/ 1 radix)))
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))))))
318

319
(defun test-log-op-with-decls (op xlo xhi ylo yhi niters
sds's avatar
sds committed
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's avatar
sds committed
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's avatar
sds committed
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's avatar
sds committed
353
                  (- (random (ash 1 (1+ r))) (ash 1 r)))))
354
    (loop for x1 = (%r)
sds's avatar
sds committed
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))))
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)))