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 ;;; (NUMERICAL-EQUAL result1 result2) => true or false
 148 ;;;
 149 ;;; This is a universal wrapper created by Liam Healy. It is
 150 ;;; implemented to support testing in GSLL. The interface is expanded,
 151 ;;; but backwards compatible with previous versions.
 152 ;;;
 153 (defgeneric numerical-equal (result1 result2 &key test)
 154   (:documentation
 155    "Return true if the results are numerically equal according to :TEST."))
 156
 157 (defmethod numerical-equal ((result1 number) (result2 number)
 158                             &key (test #'number-equal))
 159   "Return true if the the numbers are equal according to :TEST."
 160   (funcall test result1 result2))
 161
 162 (defun %sequence-equal (seq1 seq2 test)
 163   "Return true if the sequences are equal in length and each element
 164 is equal according to :TEST."
 165   (when (= (length seq1) (length seq2))
 166     (every (lambda (s1 s2) (numerical-equal s1 s2 :test test))
 167            seq1 seq2)))
 168
 169 (defmethod numerical-equal ((result1 list) (result2 list)
 170                             &key (test #'number-equal))
 171   "Return true if the lists are equal in length and each element is
 172 equal according to :TEST."
 173   (%sequence-equal result1 result2 test))
 174
 175 (defmethod numerical-equal ((result1 vector) (result2 vector)
 176                             &key (test #'number-equal))
 177   "Return true if the vectors are equal in length and each element is
 178 equal according to :TEST."
 179   (%sequence-equal result1 result2 test))
 180
 181 (defmethod numerical-equal ((result1 list) (result2 vector)
 182                             &key (test #'number-equal))
 183   "Return true if every element of the list is equal to the
 184 corresponding element of the vector."
 185   (%sequence-equal result1 result2 test))
 186
 187 (defmethod numerical-equal ((result1 vector) (result2 list)
 188                             &key (test #'number-equal))
 189   "Return true if every element of the list is equla to the
 190 corresponding element of the vector."
 191   (%sequence-equal result1 result2 test))
 192
 193 (defmethod numerical-equal ((result1 array) (result2 array)
 194                             &key (test #'number-equal))
 195   "Return true if the arrays are equal in dimension and each element
 196 is equal according to :TEST."
 197   (when (equal (array-dimensions result1) (array-dimensions result2))
 198     (every test
 199            (make-array (array-total-size result1)
 200                        :element-type (array-element-type result1)
 201                        :displaced-to result1)
 202            (make-array (array-total-size result2)
 203                        :element-type (array-element-type result2)
 204                        :displaced-to result2))))
 205
 206 (defmacro assert-numerical-equal (expected form &rest extras)
 207   (expand-assert :equal form form expected extras :test #'numerical-equal))
 208
 209 ;;; Diagnostic functions
 210 ;;; Failing a unit test is only half the problem.
 211
 212 ;;; Relative floating point error.
 213 (defun float-error (float1 float2)
 214   "Return the relative error of the floating point numbers."
 215   (if (and (floatp float1) (floatp float2))
 216       (relative-error float1 float2)
 217       (error "~A and ~A are not float." float1 float2)))
 218
 219 ;;; Floating point error as a multiple of epsilon.
 220 (defun %float-error-epsilon (float1 float2 epsilon)
 221   "Return the relative error divided by epsilon."
 222   (cond
 223     (epsilon
 224      (/ (relative-error float1 float2) epsilon))
 225     ((or (typep float1 'short-float) (typep float2 'short-float))
 226      (/ (relative-error float1 float2) short-float-epsilon))
 227     ((or (typep float1 'single-float) (typep float2 'single-float))
 228      (/ (relative-error float1 float2) single-float-epsilon))
 229     ((or (typep float1 'double-float) (typep float2 'double-float))
 230      (/ (relative-error float1 float2) double-float-epsilon))
 231     (t (/ (relative-error float1 float2) long-float-epsilon))))
 232
 233 (defun float-error-epsilon (float1 float2 &optional (epsilon *epsilon*))
 234   "Return the relative error divided by epsilon."
 235   (if (and (floatp float1) (floatp float2))
 236       (%float-error-epsilon float1 float2 epsilon)
 237       (error "~A and ~A are not float." float1 float2)))
 238
 239 ;;; Relative complex floating point error.
 240 (defun %complex-error (complex1 complex2)
 241   "Internal version with no type checking."
 242   (complex
 243    (relative-error (realpart complex1) (realpart complex2))
 244    (relative-error (imagpart complex1) (imagpart complex2))))
 245
 246 (defun complex-error (complex1 complex2)
 247   "Return the relative error of each component of the complex numbers."
 248   (if (and (typep complex1 '(complex float))
 249            (typep complex2 '(complex float)))
 250       (%complex-error complex1 complex2)
 251       (error "~A and ~A are not complex float." complex1 complex2)))
 252
 253 ;;; Complex floating point error as a multiple of epsilon.
 254 (defun %complex-error-epsilon (complex1 complex2 epsilon)
 255   "Return the relative error of each component divided by epsilon."
 256   (complex
 257    (%float-error-epsilon (realpart complex1) (realpart complex2) epsilon)
 258    (%float-error-epsilon (imagpart complex1) (imagpart complex2) epsilon)))
 259
 260 (defun complex-error-epsilon (complex1 complex2 &optional (epsilon *epsilon*))
 261   "Return the relative error of each component divided by epsilon."
 262   (if (and (typep complex1 '(complex float))
 263            (typep complex2 '(complex float)))
 264       (%complex-error-epsilon complex1 complex2 epsilon)
 265       (error "~A and ~A are not complex float." complex1 complex2)))
 266
 267 ;;; Numeric relative error
 268 (defun number-error (number1 number2)
 269   "Return the relative error of the number."
 270   (cond
 271     ((and (floatp number1) (floatp number2))
 272      (relative-error number1 number2))
 273     ((and (typep number1 '(complex float)) (typep number2 '(complex float)))
 274      (%complex-error number1 number2))
 275     ((and (numberp number1) (numberp number2))
 276      (abs (- number1 number2)))
 277     (t (error "~A and ~A are not numbers." number1 number2))))
 278
 279 ;;; Numeric error as a multiple of epsilon.
 280 (defun number-error-epsilon (number1 number2 &optional (epsilon *epsilon*))
 281   "Return the relative error divided by epsilon."
 282   (cond
 283     ((and (floatp number1) (floatp number2))
 284      (%float-error-epsilon number1 number2 epsilon))
 285     ((and (typep number1 '(complex float)) (typep number2 '(complex float)))
 286      (%complex-error-epsilon number1 number2 epsilon))
 287     (t (error "~A and ~A are not float or complex float numbers." number1 number2))))
 288
 289 ;;; Sequence errors and the indices.
 290 (defun %sequence-error (sequence1 sequence2 test error-function)
 291   "Return a sequence of the indice and error between the sequences."
 292   (let ((n1 nil) (n2 nil)
 293         (errseq '()))
 294     (dotimes (index (length sequence1) errseq)
 295       (setf n1 (elt sequence1 index)
 296             n2 (elt sequence2 index))
 297       (unless (funcall test n1 n2)
 298         (push (list (1- index) n1 n2 (funcall error-function n1 n2))
 299               errseq)))))
 300
 301 (defun sequence-error (sequence1 sequence2 &key
 302                        (test #'number-equal)
 303                        (error-function #'number-error))
 304   "Return a sequence of the indice and error between the sequence elements."
 305   (cond
 306     ((not (typep sequence1 'sequence))
 307      (error "SEQUENCE1 is not a valid sequence."))
 308     ((not (typep sequence2 'sequence))
 309      (error "SEQUENCE2 is not a valid sequence."))
 310     ((not (= (length sequence1) (length sequence2)))
 311      (error "Lengths not equal. SEQUENCE1(~D) /= SEQUENCE2(~D)."
 312             (length sequence1) (length sequence2)))
 313     (t (%sequence-error sequence1 sequence2 test error-function))))
 314
 315 ;;; Array errors and the indices.
 316 (defun %array-indices (row-major-index dimensions)
 317   "Recursively calculate the indices from the row major index."
 318   (let ((remaining (rest dimensions)))
 319     (if remaining
 320         (multiple-value-bind (index remainder)
 321             (floor row-major-index (reduce #'* remaining))
 322           (cons index (%array-indices remainder remaining)))
 323         (cons row-major-index nil))))
 324
 325 (defun %array-error (array1 array2 test errfun)
 326   "Return a list of the indices, values and error of the elements that
 327 are not equal."
 328   (let ((dimensions (array-dimensions array1))
 329         (n1 nil) (n2 nil)
 330         (indices '())
 331         (errseq '()))
 332     (dotimes (index (array-total-size array1) errseq)
 333       (setf indices (%array-indices index dimensions)
 334             n1 (apply #'aref array1 indices)
 335             n2 (apply #'aref array2 indices))
 336       (unless (funcall test n1 n2)
 337         (push (list indices n1 n2 (funcall errfun n1 n2))
 338               errseq)))))
 339
 340 (defun array-error (array1 array2 &key
 341                     (test #'number-equal)
 342                     (error-function #'number-error))
 343   "Return a list of the indices and error between the array elements."
 344   (cond
 345     ((not (arrayp array1))
 346      (error "ARRAY1 is not an array."))
 347     ((not (arrayp array2))
 348      (error "ARRAY2 is not an array."))
 349     ((not (equal (array-dimensions array1) (array-dimensions array2)))
 350      (error "Arrays are not equal dimensions."))
 351     (t (%array-error array1 array2 test error-function))))
 352
 353 ;;; Floating point data functions
 354 (defun make-2d-list (rows columns &key (initial-element 0))
 355   "Return a nested list with INITIAL-ELEMENT."
 356   (mapcar (lambda (x) (make-list columns :initial-element x))
 357           (make-list rows :initial-element initial-element)))
 358
 359 (defun complex-random (limit &optional (state *random-state*))
 360   "Return a random complex number."
 361   (check-type limit complex)
 362   (complex
 363    (random (realpart limit) state)
 364    (random (imagpart limit) state)))
 365
 366 (defun make-random-list (size &optional (limit 1.0))
 367   "Return a list of random numbers."
 368   (mapcar (if (complexp limit) #'complex-random #'random)
 369           (make-list size :initial-element limit)))
 370
 371 (defun make-random-2d-list (rows columns &optional (limit 1.0))
 372   "Return a nested list of random numbers."
 373   (mapcar (lambda (x) (make-random-list columns x))
 374           (make-list rows :initial-element limit)))
 375
 376 (defun make-random-2d-array (rows columns &optional (limit 1.0))
 377   "Return a 2D array of random numbers."
 378   (let ((new-array (make-array (list rows columns)
 379                                :element-type (type-of limit)))
 380         (random-func (if (complexp limit)
 381                          #'complex-random
 382                          #'random)))
 383     (dotimes (i0 rows new-array)
 384       (dotimes (i1 columns)
 385         (setf (aref new-array i0 i1)
 386               (funcall random-func limit))))))