floating-point.lisp

   1 ;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp -*-
   2 ;;;;
   3 ;;;; Additional LISP-UNIT tests and assertions
   4 ;;;;
   5 ;;;; Copyright (c) 2009, Thomas M. Hermann
   6 ;;;; All rights reserved.
   7 ;;;;
   8 ;;;; Redistribution and  use  in  source  and  binary  forms, with or without
   9 ;;;; modification, are permitted  provided  that the following conditions are
  10 ;;;; met:
  11 ;;;;
  12 ;;;;   o  Redistributions of  source  code  must  retain  the above copyright
  13 ;;;;      notice, this list of conditions and the following disclaimer.
  14 ;;;;   o  Redistributions in binary  form  must reproduce the above copyright
  15 ;;;;      notice, this list of  conditions  and  the  following disclaimer in
  16 ;;;;      the  documentation  and/or   other   materials  provided  with  the
  17 ;;;;      distribution.
  18 ;;;;   o  The names of the contributors may not be used to endorse or promote
  19 ;;;;      products derived from this software without  specific prior written
  20 ;;;;      permission.
  21 ;;;;
  22 ;;;; THIS SOFTWARE IS  PROVIDED  BY  THE  COPYRIGHT  HOLDERS AND CONTRIBUTORS
  23 ;;;; "AS IS"  AND  ANY  EXPRESS  OR  IMPLIED  WARRANTIES, INCLUDING,  BUT NOT
  24 ;;;; LIMITED TO, THE IMPLIED WARRANTIES  OF MERCHANTABILITY AND FITNESS FOR A
  25 ;;;; PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
  26 ;;;; OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  27 ;;;; EXEMPLARY, OR  CONSEQUENTIAL  DAMAGES  (INCLUDING,  BUT  NOT LIMITED TO,
  28 ;;;; PROCUREMENT OF  SUBSTITUTE  GOODS  OR  SERVICES;  LOSS  OF USE, DATA, OR
  29 ;;;; PROFITS; OR BUSINESS INTERRUPTION)  HOWEVER  CAUSED AND ON ANY THEORY OF
  30 ;;;; LIABILITY, WHETHER  IN  CONTRACT,  STRICT  LIABILITY, OR TORT (INCLUDING
  31 ;;;; NEGLIGENCE OR  OTHERWISE)  ARISING  IN  ANY  WAY  OUT OF THE USE OF THIS
  32 ;;;; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  33 ;;;;
  34 ;;;; References
  35 ;;;; [NumLinAlg] James W. Demmel "Applied Numerical Linear Algebra",
  36 ;;;;             Society for Industrial and Applied Mathematics, 1997
  37 ;;;;             ISBN: 0-89871-389-7
  38
  39 (common-lisp:in-package :lisp-unit)
  40
  41 (defparameter *epsilon* nil
  42   "Set the error epsilon if the defaults are not acceptable.")
  43
  44 (defparameter *significant-figures* 4
  45   "Default to 4 significant figures.")
  46
  47 ;;; (ROUNDOFF-ERROR x y) => number
  48 (defun roundoff-error (exact approximate)
  49   "Return the error delta between the exact and approximate floating
  50 point value."
  51   ;; [NumLinAlg] : Equation 1.1, pg. 12
  52   (abs (if (or (zerop exact) (zerop approximate))
  53            (+ exact approximate)
  54            (- (/ approximate exact) 1.0))))
  55
  56 ;;; (FLOAT-EQUAL float1 float2 &optional epsilon) => true or false
  57 (defun %float-equal (float1 float2 epsilon)
  58   "Internal version that does not verify arguments as floats."
  59   (cond
  60     ((and (zerop float1) (zerop float2)))
  61     (epsilon
  62      (< (roundoff-error float1 float2) epsilon))
  63     ((or (typep float1 'short-float) (typep float2 'short-float))
  64      (< (roundoff-error float1 float2) (* 2.0 short-float-epsilon)))
  65     ((or (typep float1 'single-float) (typep float2 'single-float))
  66      (< (roundoff-error float1 float2) (* 2.0 single-float-epsilon)))
  67     ((or (typep float1 'double-float) (typep float2 'double-float))
  68      (< (roundoff-error float1 float2) (* 2.0 double-float-epsilon)))
  69     (t (< (roundoff-error float1 float2) (* 2.0 long-float-epsilon)))))
  70
  71 (defun float-equal (float1 float2 &optional (epsilon *epsilon*))
  72   "Return true if the absolute difference between float1 and float2 is
  73 less than some epsilon."
  74   (when (and (floatp float1) (floatp float2))
  75     (%float-equal float1 float2 epsilon)))
  76
  77 (defmacro assert-float-equal (expected form &rest extras)
  78   (expand-assert :equal form form expected extras :test #'float-equal))
  79
  80 ;;; (COMPLEX-EQUAL complex1 complex2 &optional epsilon) => true or false
  81 (defun %complex-equal (complex1 complex2 epsilon)
  82   "Internal version that does not verify arguments as (complex float)."
  83   (and
  84    (%float-equal (realpart complex1) (realpart complex2) epsilon)
  85    (%float-equal (imagpart complex1) (imagpart complex2) epsilon)))
  86
  87 (defun complex-equal (complex1 complex2 &optional (epsilon *epsilon*))
  88   "Return true if the absolute difference between Re(complex1),
  89 Re(complex2) and the absolute difference between Im(complex1),
  90 Im(complex2) is less than epsilon."
  91   (when (and (typep complex1 '(complex float))
  92              (typep complex2 '(complex float)))
  93     (%complex-equal complex1 complex2 epsilon)))
  94
  95 (defmacro assert-complex-equal (expected form &rest extras)
  96   (expand-assert :equal form form expected extras :test #'complex-equal))
  97
  98 ;;; (NUMBER-EQUAL number1 number2) => true or false
  99 (defun number-equal (number1 number2 &optional (epsilon *epsilon*))
 100   "Return true if the numbers are equal using the appropriate
 101 comparison."
 102   (cond
 103     ((and (floatp number1) (floatp number2))
 104      (%float-equal number1 number2 epsilon))
 105     ((and (typep number1 '(complex float)) (typep number2 '(complex float)))
 106      (%complex-equal number1 number2 epsilon))
 107     ((and (numberp number1) (numberp number2))
 108      (= number1 number2))
 109     (t (error "~A and ~A are not numbers." number1 number2))))
 110
 111 (defmacro assert-number-equal (expected form &rest extras)
 112   (expand-assert :equal form form expected extras :test #'number-equal))
 113
 114 ;;; (NORMALIZE-FLOAT significand &optional exponent) => significand,exponent
 115 (defun normalize-float (significand &optional (exponent 0))
 116   "Return the normalized floating point number and exponent."
 117   (cond
 118     ((zerop significand)
 119      (values significand 0))
 120     ((>= (abs significand) 10)
 121      (normalize-float (/ significand 10.0) (1+ exponent)))
 122     ((< (abs significand) 1)
 123      (normalize-float (* significand 10.0) (1- exponent)))
 124     (t (values significand exponent))))
 125
 126 ;;; (SIGFIG-EQUAL float1 float2 significant-figures) => true or false
 127 (defun sigfig-equal (float1 float2 &optional (significant-figures *significant-figures*))
 128   "Return true if the floating point numbers have equal significant
 129 figures."
 130   ;; Convert 5 to precision of FLOAT1 and 10 to precision of FLOAT2.
 131   ;; Then, rely on Rule of Float and Rational Contagion, CLHS 12.1.4.1,
 132   ;; to obtain a DELTA of the proper precision.
 133   (let ((delta (* (float 5 float1) (expt (float 10 float2) (- significant-figures)))))
 134     (if (or (zerop float1) (zerop float2))
 135         (< (abs (+ float1 float2)) delta)
 136         (multiple-value-bind (sig1 exp1) (normalize-float float1)
 137           (multiple-value-bind (sig2 exp2) (normalize-float float2)
 138             (and (= exp1 exp2)
 139                  (< (abs (- sig1 sig2)) delta)))))))
 140
 141 (defmacro assert-sigfig-equal (significant-figures expected form &rest extras)
 142   (expand-assert :equal form form expected extras
 143                  :test (lambda (f1 f2) (sigfig-equal f1 f2 significant-figures))))
 144
 145 ;;; (ARRAY-EQUAL array1 array2) => true or false
 146 (defun array-equal (array1 array2 &key (test #'number-equal))
 147   "Return true if the elements of the array are equal."
 148   (when (equal (array-dimensions array1) (array-dimensions array2))
 149     (every test
 150      (make-array (array-total-size array1)
 151                  :element-type (array-element-type array1)
 152                  :displaced-to array1)
 153      (make-array (array-total-size array2)
 154                  :element-type (array-element-type array2)
 155                  :displaced-to array2))))
 156
 157 (defmacro assert-array-equal (element-test expected form &rest extras)
 158   (expand-assert :equal form form expected extras
 159                  :test `(lambda (a1 a2) (array-equal a1 a2 :test ,element-test))))
 160
 161 ;;; (NUMERICAL-EQUAL result1 result2) => true or false
 162 ;;;
 163 ;;; This is a universal wrapper created by Liam Healy. It is
 164 ;;; implemented to support testing in GSLL. The interface is expanded,
 165 ;;; but backwards compatible with previous versions.
 166 ;;;
 167 (defun numerical-equal (result1 result2 &key (test #'number-equal))
 168   (cond
 169     ((and (numberp result1) (numberp result2))
 170      (funcall test result1 result2))
 171     ((and (typep result1 'sequence) (typep result2 'sequence))
 172      (when (= (length result1) (length result2))
 173        (every (lambda (r1 r2) (numerical-equal r1 r2 :test test))
 174               result1 result2)))
 175     ((and (arrayp result1) (arrayp result2))
 176      (array-equal result1 result2 :test test))
 177     (t (error "~A and ~A are not valid arguments." result1 result2))))
 178
 179 (defmacro assert-numerical-equal (expected form &rest extras)
 180   (expand-assert :equal form form expected extras :test #'numerical-equal))
 181
 182 ;;; Diagnostic functions
 183 ;;; Failing a unit test is only half the problem.
 184
 185 ;;; Relative floating point error.
 186 (defun float-error (float1 float2)
 187   "Return the roundoff error of the floating point numbers."
 188   (if (and (floatp float1) (floatp float2))
 189       (roundoff-error float1 float2)
 190       (error "~A and ~A are not float." float1 float2)))
 191
 192 ;;; Floating point error as a multiple of epsilon.
 193 (defun %float-error-epsilon (float1 float2 epsilon)
 194   "Return the roundoff error divided by epsilon."
 195   (cond
 196     (epsilon
 197      (/ (roundoff-error float1 float2) epsilon))
 198     ((or (typep float1 'short-float) (typep float2 'short-float))
 199      (/ (roundoff-error float1 float2) short-float-epsilon))
 200     ((or (typep float1 'single-float) (typep float2 'single-float))
 201      (/ (roundoff-error float1 float2) single-float-epsilon))
 202     ((or (typep float1 'double-float) (typep float2 'double-float))
 203      (/ (roundoff-error float1 float2) double-float-epsilon))
 204     (t (/ (roundoff-error float1 float2) long-float-epsilon))))
 205
 206 (defun float-error-epsilon (float1 float2 &optional (epsilon *epsilon*))
 207   "Return the roundoff error divided by epsilon."
 208   (if (and (floatp float1) (floatp float2))
 209       (%float-error-epsilon float1 float2 epsilon)
 210       (error "~A and ~A are not float." float1 float2)))
 211
 212 ;;; Relative complex floating point error.
 213 (defun %complex-error (complex1 complex2)
 214   "Internal version with no type checking."
 215   (complex
 216    (roundoff-error (realpart complex1) (realpart complex2))
 217    (roundoff-error (imagpart complex1) (imagpart complex2))))
 218
 219 (defun complex-error (complex1 complex2)
 220   "Return the roundoff error of each component of the complex numbers."
 221   (if (and (typep complex1 '(complex float))
 222            (typep complex2 '(complex float)))
 223       (%complex-error complex1 complex2)
 224       (error "~A and ~A are not complex float." complex1 complex2)))
 225
 226 ;;; Complex floating point error as a multiple of epsilon.
 227 (defun %complex-error-epsilon (complex1 complex2 epsilon)
 228   "Return the roundoff error of each component divided by epsilon."
 229   (complex
 230    (%float-error-epsilon (realpart complex1) (realpart complex2) epsilon)
 231    (%float-error-epsilon (imagpart complex1) (imagpart complex2) epsilon)))
 232
 233 (defun complex-error-epsilon (complex1 complex2 &optional (epsilon *epsilon*))
 234   "Return the roundoff error of each component divided by epsilon."
 235   (if (and (typep complex1 '(complex float))
 236            (typep complex2 '(complex float)))
 237       (%complex-error-epsilon complex1 complex2 epsilon)
 238       (error "~A and ~A are not complex float." complex1 complex2)))
 239
 240 ;;; Numeric roundoff error
 241 (defun number-error (number1 number2)
 242   "Return the roundoff error of the number."
 243   (cond
 244     ((and (floatp number1) (floatp number2))
 245      (roundoff-error number1 number2))
 246     ((and (typep number1 '(complex float)) (typep number2 '(complex float)))
 247      (%complex-error number1 number2))
 248     ((and (numberp number1) (numberp number2))
 249      (abs (- (/ number1 number2) 1)))
 250     (t (error "~A and ~A are not numbers." number1 number2))))
 251
 252 ;;; Numeric error as a multiple of epsilon.
 253 (defun number-error-epsilon (number1 number2 &optional (epsilon *epsilon*))
 254   "Return the roundoff error divided by epsilon."
 255   (cond
 256     ((and (floatp number1) (floatp number2))
 257      (%float-error-epsilon number1 number2 epsilon))
 258     ((and (typep number1 '(complex float)) (typep number2 '(complex float)))
 259      (%complex-error-epsilon number1 number2 epsilon))
 260     (t (error "~A and ~A are not float or complex float numbers." number1 number2))))
 261
 262 ;;; Sequence errors and the indices.
 263 (defun %sequence-error (sequence1 sequence2 test error-function)
 264   "Return a sequence of the indice and error between the sequences."
 265   (let ((errseq '()))
 266     (dotimes (index (length sequence1) errseq)
 267       (let ((n1 (elt sequence1 index))
 268             (n2 (elt sequence2 index)))
 269         (unless (funcall test n1 n2)
 270           (push (list index n1 n2 (funcall error-function n1 n2))
 271                 errseq))))))
 272
 273 (defun sequence-error (sequence1 sequence2 &key
 274                        (test #'number-equal)
 275                        (error-function #'number-error))
 276   "Return a sequence of the indice and error between the sequence elements."
 277   (cond
 278     ((not (typep sequence1 'sequence))
 279      (error "SEQUENCE1 is not a valid sequence."))
 280     ((not (typep sequence2 'sequence))
 281      (error "SEQUENCE2 is not a valid sequence."))
 282     ((not (= (length sequence1) (length sequence2)))
 283      (error "Lengths not equal. SEQUENCE1(~D) /= SEQUENCE2(~D)."
 284             (length sequence1) (length sequence2)))
 285     (t (%sequence-error sequence1 sequence2 test error-function))))
 286
 287 ;;; Array errors and the indices.
 288 (defun %array-indices (row-major-index dimensions)
 289   "Recursively calculate the indices from the row major index."
 290   (let ((remaining (rest dimensions)))
 291     (if remaining
 292         (multiple-value-bind (index remainder)
 293             (floor row-major-index (reduce #'* remaining))
 294           (cons index (%array-indices remainder remaining)))
 295         (cons row-major-index nil))))
 296
 297 (defun %array-error (array1 array2 test errfun)
 298   "Return a list of the indices, values and error of the elements that
 299 are not equal."
 300   (let ((dimensions (array-dimensions array1))
 301         (errseq '()))
 302     (dotimes (index (array-total-size array1) errseq)
 303       (let* ((indices (%array-indices index dimensions))
 304              (n1 (apply #'aref array1 indices))
 305              (n2 (apply #'aref array2 indices)))
 306         (unless (funcall test n1 n2)
 307           (push (list indices n1 n2 (funcall errfun n1 n2))
 308                 errseq))))))
 309
 310 (defun array-error (array1 array2 &key
 311                     (test #'number-equal)
 312                     (error-function #'number-error))
 313   "Return a list of the indices and error between the array elements."
 314   (cond
 315     ((not (arrayp array1))
 316      (error "ARRAY1 is not an array."))
 317     ((not (arrayp array2))
 318      (error "ARRAY2 is not an array."))
 319     ((not (equal (array-dimensions array1) (array-dimensions array2)))
 320      (error "Arrays are not equal dimensions."))
 321     (t (%array-error array1 array2 test error-function))))