floating-point.lisp

   1 ;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp -*-
   2 ;;;
   3 ;;; Floating tests and assertions for LISP-UNIT
   4 ;;;
   5 ;;; Copyright (c) 2009 Thomas M. Hermann
   6 ;;;
   7 ;;; Permission is hereby granted, free of charge, to any person obtaining
   8 ;;; a copy of this software and associated documentation files (the "Software"),
   9 ;;; to deal in the Software without restriction, including without limitation
  10 ;;; the rights to use, copy, modify, merge, publish, distribute, sublicense,
  11 ;;; and/or sell copies of the Software, and to permit persons to whom the
  12 ;;; Software is furnished to do so, subject to the following conditions:
  13 ;;;
  14 ;;; The above copyright notice and this permission notice shall be included
  15 ;;; in all copies or substantial portions of the Software.
  16 ;;;
  17 ;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  18 ;;; OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19 ;;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  20 ;;; THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
  21 ;;; OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  22 ;;; ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  23 ;;; OTHER DEALINGS IN THE SOFTWARE.
  24 ;;;
  25 ;;; References
  26 ;;; [NumAlgoC] Gisela Engeln-Mullges and Frank Uhlig "Numerical
  27 ;;;            Algorithms with C", Springer, 1996
  28 ;;;            ISBN: 3-540-60530-4
  29
  30 (common-lisp:in-package :lisp-unit)
  31
  32 (defparameter *epsilon* nil
  33   "Set the error epsilon if the defaults are not acceptable.")
  34
  35 (defparameter *significant-figures* 4
  36   "Default to 4 significant figures.")
  37
  38 ;;; (RELATIVE-ERROR x y) => number
  39 ;;; [NumAlgoC] : Definition 1.3, pg. 2
  40 ;;;              modified with Definition 1.1, pg. 1
  41 ;;;
  42 ;;; The definition of relative error in this routine is modified from
  43 ;;; the Definition 1.3 in [NumAlgoC] for cases when either the exact
  44 ;;; or the approximate value equals zero. According to Definition 1.3,
  45 ;;; the relative error is identically equal to 1 in those cases. This
  46 ;;; function returns the absolue error in those cases. This is more
  47 ;;; useful for testing.
  48 (defun relative-error (exact approximate)
  49   "Return the error delta between the exact and approximate floating
  50 point value."
  51   (abs (if (or (zerop exact) (zerop approximate))
  52            (- exact approximate)
  53            (/ (- exact approximate) exact))))
  54
  55 ;;; (FLOAT-EQUAL float1 float2 &optional epsilon) => true or false
  56 (defun %float-equal (float1 float2 epsilon)
  57   "Internal version that does not verify arguments as floats."
  58   (cond
  59     ((and (zerop float1) (zerop float2)))
  60     (epsilon
  61      (< (relative-error float1 float2) epsilon))
  62     ((or (typep float1 'short-float) (typep float2 'short-float))
  63      (< (relative-error float1 float2) (* 2.0 short-float-epsilon)))
  64     ((or (typep float1 'single-float) (typep float2 'single-float))
  65      (< (relative-error float1 float2) (* 2.0 single-float-epsilon)))
  66     ((or (typep float1 'double-float) (typep float2 'double-float))
  67      (< (relative-error float1 float2) (* 2.0 double-float-epsilon)))
  68     (t (< (relative-error float1 float2) (* 2.0 long-float-epsilon)))))
  69
  70 (defun float-equal (float1 float2 &optional (epsilon *epsilon*))
  71   "Return true if the relative error between float1 and float2 is less
  72 than some epsilon."
  73   (when (and (floatp float1) (floatp float2))
  74     (%float-equal float1 float2 epsilon)))
  75
  76 (defmacro assert-float-equal (expected form &rest extras)
  77   (expand-assert :equal form form expected extras :test #'float-equal))
  78
  79 ;;; (COMPLEX-EQUAL complex1 complex2 &optional epsilon) => true or false
  80 (defun %complex-equal (complex1 complex2 epsilon)
  81   "Internal version that does not verify arguments as (complex float)."
  82   (and
  83    (%float-equal (realpart complex1) (realpart complex2) epsilon)
  84    (%float-equal (imagpart complex1) (imagpart complex2) epsilon)))
  85
  86 (defun complex-equal (complex1 complex2 &optional (epsilon *epsilon*))
  87   "Return true if the relative error between Re(complex1),
  88 Re(complex2) and between Im(complex1), Im(complex2) are each less than
  89 epsilon."
  90   (when (and (typep complex1 '(complex float))
  91              (typep complex2 '(complex float)))
  92     (%complex-equal complex1 complex2 epsilon)))
  93
  94 (defmacro assert-complex-equal (expected form &rest extras)
  95   (expand-assert :equal form form expected extras :test #'complex-equal))
  96
  97 ;;; (NUMBER-EQUAL number1 number2) => true or false
  98 (defun number-equal (number1 number2 &optional (epsilon *epsilon*))
  99   "Return true if the numbers are equal using the appropriate
 100 comparison."
 101   (cond
 102     ((and (floatp number1) (floatp number2))
 103      (%float-equal number1 number2 epsilon))
 104     ((and (typep number1 '(complex float)) (typep number2 '(complex float)))
 105      (%complex-equal number1 number2 epsilon))
 106     ((and (numberp number1) (numberp number2))
 107      (= number1 number2))
 108     (t (error "~A and ~A are not numbers." number1 number2))))
 109
 110 (defmacro assert-number-equal (expected form &rest extras)
 111   (expand-assert :equal form form expected extras :test #'number-equal))
 112
 113 ;;; (NORMALIZE-FLOAT significand &optional exponent) => significand,exponent
 114 ;;; [NumAlgoC] : Definition 1.7, pg. 4
 115 ;;;
 116 ;;; To avoid using 0.1, first 1.0 <= significand < 10. On the final
 117 ;;; return, scale 0.1 <= significand < 1.
 118 (defun normalize-float (significand &optional (exponent 0))
 119   "Return the normalized floating point number and exponent."
 120   (cond
 121     ((zerop significand)
 122      (values significand 0))
 123     ((>= (abs significand) 10)
 124      (normalize-float (/ significand 10.0) (1+ exponent)))
 125     ((< (abs significand) 1)
 126      (normalize-float (* significand 10.0) (1- exponent)))
 127     (t (values (/ significand 10.0) (1+ exponent)))))
 128
 129 ;;; (SIGFIG-EQUAL float1 float2 significant-figures) => true or false
 130 (defun sigfig-equal (float1 float2 &optional (significant-figures *significant-figures*))
 131   "Return true if the floating point numbers have equal significant
 132 figures."
 133   ;; Convert 0.5 to precision of FLOAT1 and 10 to precision of FLOAT2.
 134   ;; Then, rely on Rule of Float and Rational Contagion, CLHS 12.1.4.1,
 135   ;; to obtain a DELTA of the proper precision.
 136   (let ((delta (* (float 0.5 float1) (expt (float 10 float2) (- significant-figures)))))
 137     (if (or (zerop float1) (zerop float2))
 138         (< (abs (+ float1 float2)) delta)
 139         (multiple-value-bind (sig1 exp1) (normalize-float float1)
 140           (multiple-value-bind (sig2 exp2) (normalize-float float2)
 141             (and (= exp1 exp2)
 142                  (< (abs (- sig1 sig2)) delta)))))))
 143
 144 (defmacro assert-sigfig-equal (expected form &rest extras)
 145   (expand-assert :equal form form expected extras :test #'sigfig-equal))
 146
 147 ;;; (SEQ-EQUAL seq1 seq2) => true or false
 148 (defun seq-equal (seq1 seq2 &key (test #'number-equal))
 149   "Return true if the elements of the sequences are equal."
 150   (and
 151    (typep seq1 'sequence) (typep seq2 'sequence)
 152    (= (length seq1) (length seq2))
 153    (every test seq1 seq2)))
 154
 155 (defmacro assert-seq-equal (test expected form &rest extras)
 156   (expand-assert :equal form form expected extras
 157                  :test `(lambda (s1 s2) (seq-equal s1 s2 :test ,test))))
 158
 159 ;;; (ARRAY-EQUAL array1 array2) => true or false
 160 (defun array-equal (array1 array2 &key (test #'number-equal))
 161   "Return true if the elements of the arrays are equal."
 162   (when (equal (array-dimensions array1) (array-dimensions array2))
 163     (every test
 164      (make-array (array-total-size array1)
 165                  :element-type (array-element-type array1)
 166                  :displaced-to array1)
 167      (make-array (array-total-size array2)
 168                  :element-type (array-element-type array2)
 169                  :displaced-to array2))))
 170
 171 (defmacro assert-array-equal (test expected form &rest extras)
 172   (expand-assert :equal form form expected extras
 173                  :test `(lambda (a1 a2) (array-equal a1 a2 :test ,test))))
 174
 175 ;;; (NUMERICAL-EQUAL result1 result2) => true or false
 176 ;;;
 177 ;;; This is a universal wrapper created by Liam Healy. It is
 178 ;;; implemented to support testing in GSLL. The interface is expanded,
 179 ;;; but backwards compatible with previous versions.
 180 ;;;
 181 (defun numerical-equal (result1 result2 &key (test #'number-equal))
 182   (cond
 183     ((and (numberp result1) (numberp result2))
 184      (funcall test result1 result2))
 185     ((and (typep result1 'sequence) (typep result2 'sequence))
 186      (when (= (length result1) (length result2))
 187        (every (lambda (r1 r2) (numerical-equal r1 r2 :test test))
 188               result1 result2)))
 189     ((and (arrayp result1) (arrayp result2))
 190      (array-equal result1 result2 :test test))
 191     (t (error "~A and ~A are not valid arguments." result1 result2))))
 192
 193 (defmacro assert-numerical-equal (expected form &rest extras)
 194   (expand-assert :equal form form expected extras :test #'numerical-equal))
 195
 196 ;;; Diagnostic functions
 197 ;;; Failing a unit test is only half the problem.
 198
 199 ;;; Relative floating point error.
 200 (defun float-error (float1 float2)
 201   "Return the relative error of the floating point numbers."
 202   (if (and (floatp float1) (floatp float2))
 203       (relative-error float1 float2)
 204       (error "~A and ~A are not float." float1 float2)))
 205
 206 ;;; Floating point error as a multiple of epsilon.
 207 (defun %float-error-epsilon (float1 float2 epsilon)
 208   "Return the relative error divided by epsilon."
 209   (cond
 210     (epsilon
 211      (/ (relative-error float1 float2) epsilon))
 212     ((or (typep float1 'short-float) (typep float2 'short-float))
 213      (/ (relative-error float1 float2) short-float-epsilon))
 214     ((or (typep float1 'single-float) (typep float2 'single-float))
 215      (/ (relative-error float1 float2) single-float-epsilon))
 216     ((or (typep float1 'double-float) (typep float2 'double-float))
 217      (/ (relative-error float1 float2) double-float-epsilon))
 218     (t (/ (relative-error float1 float2) long-float-epsilon))))
 219
 220 (defun float-error-epsilon (float1 float2 &optional (epsilon *epsilon*))
 221   "Return the relative error divided by epsilon."
 222   (if (and (floatp float1) (floatp float2))
 223       (%float-error-epsilon float1 float2 epsilon)
 224       (error "~A and ~A are not float." float1 float2)))
 225
 226 ;;; Relative complex floating point error.
 227 (defun %complex-error (complex1 complex2)
 228   "Internal version with no type checking."
 229   (complex
 230    (relative-error (realpart complex1) (realpart complex2))
 231    (relative-error (imagpart complex1) (imagpart complex2))))
 232
 233 (defun complex-error (complex1 complex2)
 234   "Return the relative error of each component of the complex numbers."
 235   (if (and (typep complex1 '(complex float))
 236            (typep complex2 '(complex float)))
 237       (%complex-error complex1 complex2)
 238       (error "~A and ~A are not complex float." complex1 complex2)))
 239
 240 ;;; Complex floating point error as a multiple of epsilon.
 241 (defun %complex-error-epsilon (complex1 complex2 epsilon)
 242   "Return the relative error of each component divided by epsilon."
 243   (complex
 244    (%float-error-epsilon (realpart complex1) (realpart complex2) epsilon)
 245    (%float-error-epsilon (imagpart complex1) (imagpart complex2) epsilon)))
 246
 247 (defun complex-error-epsilon (complex1 complex2 &optional (epsilon *epsilon*))
 248   "Return the relative error of each component divided by epsilon."
 249   (if (and (typep complex1 '(complex float))
 250            (typep complex2 '(complex float)))
 251       (%complex-error-epsilon complex1 complex2 epsilon)
 252       (error "~A and ~A are not complex float." complex1 complex2)))
 253
 254 ;;; Numeric relative error
 255 (defun number-error (number1 number2)
 256   "Return the relative error of the number."
 257   (cond
 258     ((and (floatp number1) (floatp number2))
 259      (relative-error number1 number2))
 260     ((and (typep number1 '(complex float)) (typep number2 '(complex float)))
 261      (%complex-error number1 number2))
 262     ((and (numberp number1) (numberp number2))
 263      (abs (- number1 number2)))
 264     (t (error "~A and ~A are not numbers." number1 number2))))
 265
 266 ;;; Numeric error as a multiple of epsilon.
 267 (defun number-error-epsilon (number1 number2 &optional (epsilon *epsilon*))
 268   "Return the relative error divided by epsilon."
 269   (cond
 270     ((and (floatp number1) (floatp number2))
 271      (%float-error-epsilon number1 number2 epsilon))
 272     ((and (typep number1 '(complex float)) (typep number2 '(complex float)))
 273      (%complex-error-epsilon number1 number2 epsilon))
 274     (t (error "~A and ~A are not float or complex float numbers." number1 number2))))
 275
 276 ;;; Sequence errors and the indices.
 277 (defun %sequence-error (sequence1 sequence2 test error-function)
 278   "Return a sequence of the indice and error between the sequences."
 279   (let ((n1 nil) (n2 nil)
 280         (errseq '()))
 281     (dotimes (index (length sequence1) errseq)
 282       (setf n1 (elt sequence1 index)
 283             n2 (elt sequence2 index))
 284       (unless (funcall test n1 n2)
 285         (push (list (1- index) n1 n2 (funcall error-function n1 n2))
 286               errseq)))))
 287
 288 (defun sequence-error (sequence1 sequence2 &key
 289                        (test #'number-equal)
 290                        (error-function #'number-error))
 291   "Return a sequence of the indice and error between the sequence elements."
 292   (cond
 293     ((not (typep sequence1 'sequence))
 294      (error "SEQUENCE1 is not a valid sequence."))
 295     ((not (typep sequence2 'sequence))
 296      (error "SEQUENCE2 is not a valid sequence."))
 297     ((not (= (length sequence1) (length sequence2)))
 298      (error "Lengths not equal. SEQUENCE1(~D) /= SEQUENCE2(~D)."
 299             (length sequence1) (length sequence2)))
 300     (t (%sequence-error sequence1 sequence2 test error-function))))
 301
 302 ;;; Array errors and the indices.
 303 (defun %array-indices (row-major-index dimensions)
 304   "Recursively calculate the indices from the row major index."
 305   (let ((remaining (rest dimensions)))
 306     (if remaining
 307         (multiple-value-bind (index remainder)
 308             (floor row-major-index (reduce #'* remaining))
 309           (cons index (%array-indices remainder remaining)))
 310         (cons row-major-index nil))))
 311
 312 (defun %array-error (array1 array2 test errfun)
 313   "Return a list of the indices, values and error of the elements that
 314 are not equal."
 315   (let ((dimensions (array-dimensions array1))
 316         (n1 nil) (n2 nil)
 317         (indices '())
 318         (errseq '()))
 319     (dotimes (index (array-total-size array1) errseq)
 320       (setf indices (%array-indices index dimensions)
 321             n1 (apply #'aref array1 indices)
 322             n2 (apply #'aref array2 indices))
 323       (unless (funcall test n1 n2)
 324         (push (list indices n1 n2 (funcall errfun n1 n2))
 325               errseq)))))
 326
 327 (defun array-error (array1 array2 &key
 328                     (test #'number-equal)
 329                     (error-function #'number-error))
 330   "Return a list of the indices and error between the array elements."
 331   (cond
 332     ((not (arrayp array1))
 333      (error "ARRAY1 is not an array."))
 334     ((not (arrayp array2))
 335      (error "ARRAY2 is not an array."))
 336     ((not (equal (array-dimensions array1) (array-dimensions array2)))
 337      (error "Arrays are not equal dimensions."))
 338     (t (%array-error array1 array2 test error-function))))
 339
 340 ;;; Floating point data functions
 341 (defun make-2d-list (rows columns &key (initial-element 0))
 342   "Return a nested list with INITIAL-ELEMENT."
 343   (mapcar (lambda (x) (make-list columns :initial-element x))
 344           (make-list rows :initial-element initial-element)))
 345
 346 (defun complex-random (limit &optional (state *random-state*))
 347   "Return a random complex number."
 348   (check-type limit complex)
 349   (complex
 350    (random (realpart limit) state)
 351    (random (imagpart limit) state)))
 352
 353 (defun make-random-list (size &optional (limit 1.0))
 354   "Return a list of random numbers."
 355   (mapcar (if (complexp limit) #'complex-random #'random)
 356           (make-list size :initial-element limit)))
 357
 358 (defun make-random-2d-list (rows columns &optional (limit 1.0))
 359   "Return a nested list of random numbers."
 360   (mapcar (lambda (x) (make-random-list columns x))
 361           (make-list rows :initial-element limit)))
 362
 363 (defun make-random-2d-array (rows columns &optional (limit 1.0))
 364   "Return a 2D array of random numbers."
 365   (let ((new-array (make-array (list rows columns)))
 366         (random-func (if (complexp limit)
 367                          #'complex-random
 368                          #'random)))
 369     (dotimes (i0 rows new-array)
 370       (dotimes (i1 columns)
 371         (setf (aref new-array i0 i1)
 372               (funcall random-func limit))))))