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
  57 ;;; (FLOAT-EQUAL float1 float2 &optional epsilon) => true or false
  58 (defun %float-equal (float1 float2 epsilon)
  59   "Internal version that does not verify arguments as floats."
  60   (cond
  61     ((and (zerop float1) (zerop float2)))
  62     (epsilon
  63      (> epsilon (roundoff-error float1 float2)))
  64     ((and (typep float1 'double-float) (typep float2 'double-float))
  65      (> (* 2.0 double-float-epsilon) (roundoff-error float1 float2)))
  66     ((or (typep float1 'single-float) (typep float2 'single-float))
  67      (> (* 2.0 single-float-epsilon) (roundoff-error float1 float2)))
  68     (t nil)))
  69
  70 (defun float-equal (float1 float2 &optional (epsilon *epsilon*))
  71   "Return true if the absolute difference between float1 and float2 is
  72 less 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 absolute difference between Re(complex1),
  88 Re(complex2) and the absolute difference between Im(complex1),
  89 Im(complex2) is less than 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 (defun normalize-float (significand &optional (exponent 0))
 115   "Return the normalized floating point number and exponent."
 116   (cond
 117     ((zerop significand)
 118      (values significand 0))
 119     ((>= (abs significand) 10)
 120      (normalize-float (/ significand 10.0) (1+ exponent)))
 121     ((< (abs significand) 1)
 122      (normalize-float (* significand 10.0) (1- exponent)))
 123     (t (values significand exponent))))
 124
 125 ;;; (SIGFIG-EQUAL float1 float2 significant-figures) => true or false
 126 (defun sigfig-equal (float1 float2 &optional (significant-figures *significant-figures*))
 127   "Return true if the floating point numbers have equal significant
 128 figures."
 129   ;; Convert 5 to precision of FLOAT1 and 10 to precision of FLOAT2.
 130   ;; Then, rely on Rule of Float and Rational Contagion, CLHS 12.1.4.1,
 131   ;; to obtain a DELTA of the proper precision.
 132   (let ((delta (* (float 5 float1) (expt (float 10 float2) (- significant-figures)))))
 133     (if (or (zerop float1) (zerop float2))
 134         (< (abs (+ float1 float2)) delta)
 135         (multiple-value-bind (sig1 exp1) (normalize-float float1)
 136           (multiple-value-bind (sig2 exp2) (normalize-float float2)
 137             (and (= exp1 exp2)
 138                  (< (abs (- sig1 sig2)) delta)))))))
 139
 140 (defmacro assert-sigfig-equal (significant-figures expected form &rest extras)
 141   (expand-assert :equal form form expected extras
 142                  :test (lambda (f1 f2) (sigfig-equal f1 f2 significant-figures))))
 143
 144 ;;; (ARRAY-EQUAL array1 array2) => true or false
 145 (defun array-equal (array1 array2 &key (test #'number-equal))
 146   "Return true if the elements of the array are equal."
 147   (when (equal (array-dimensions array1) (array-dimensions array2))
 148     (every test
 149      (make-array (reduce #'* (array-dimensions array1))
 150                  :element-type (array-element-type array1)
 151                  :displaced-to array1)
 152      (make-array (reduce #'* (array-dimensions array2))
 153                  :element-type (array-element-type array2)
 154                  :displaced-to array2))))
 155
 156 (defmacro assert-array-equal (element-test expected form &rest extras)
 157   (expand-assert :equal form form expected extras
 158                  :test `(lambda (a1 a2) (array-equal a1 a2 :test ,element-test))))
 159
 160 ;;; (NUMERICAL-EQUAL result1 result2) => true or false
 161 ;;;
 162 ;;; This is a universal wrapper created by Liam Healy. It is
 163 ;;; implemented to support testing in GSLL. The interface is expanded,
 164 ;;; but backwards compatible with previous versions.
 165 ;;;
 166 (defun numerical-equal (result1 result2 &key (test #'number-equal))
 167   (cond
 168     ((and (numberp result1) (numberp result2))
 169      (funcall test result1 result2))
 170     ((and (typep result1 'sequence) (typep result2 'sequence))
 171      (when (= (length result1) (length result2))
 172        (every (lambda (r1 r2) (numerical-equal r1 r2 :test test))
 173               result1 result2)))
 174     ((and (arrayp result1) (arrayp result2))
 175      (array-equal result1 result2 :test test))
 176     (t (error "~A and/or ~A are not valid arguments." result1 result2))))
 177
 178 (defmacro assert-numerical-equal (expected form &rest extras)
 179   (expand-assert :equal form form expected extras :test #'numerical-equal))