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 &optional (epsilon *epsilon*))
  58   "Return true if the absolute difference between float1 and float2 is
  59 less than some epsilon."
  60   (when (and (floatp float1) (floatp float2))
  61     (cond
  62       ((and (zerop float1) (zerop float2)))
  63       (epsilon
  64        (> epsilon (roundoff-error float1 float2)))
  65       ((and (typep float1 'double-float) (typep float2 'double-float))
  66        (> (* 2.0 double-float-epsilon) (roundoff-error float1 float2)))
  67       ((or (typep float1 'single-float) (typep float2 'single-float))
  68        (> (* 2.0 single-float-epsilon) (roundoff-error float1 float2)))
  69       (t nil))))
  70
  71 (defmacro assert-float-equal (expected form &rest extras)
  72   (expand-assert :equal form form expected extras :test #'float-equal))
  73
  74 ;;; (COMPLEX-EQUAL complex1 complex2 &optional epsilon) => true or false
  75 (defun complex-equal (complex1 complex2 &optional (epsilon *epsilon*))
  76   "Return true if the absolute difference between Re(complex1),
  77 Re(complex2) and the absolute difference between Im(complex1),
  78 Im(complex2) is less than epsilon."
  79   (and
  80    (typep complex1 '(complex float))
  81    (typep complex2 '(complex float))
  82    (float-equal (realpart complex1) (realpart complex2) epsilon)
  83    (float-equal (imagpart complex1) (imagpart complex2) epsilon)))
  84
  85 (defmacro assert-complex-equal (expected form &rest extras)
  86   (expand-assert :equal form form expected extras :test #'complex-equal))
  87
  88 ;;; (NUMBER-EQUAL number1 number2) => true or false
  89 (defun number-equal (number1 number2 &optional (epsilon *epsilon*))
  90   "Return true if the numbers are equal using the appropriate
  91 comparison."
  92   (cond
  93     ((and (floatp number1) (floatp number2))
  94      (float-equal number1 number2 epsilon))
  95     ((and (typep number1 '(complex float)) (typep number2 '(complex float)))
  96      (complex-equal number1 number2 epsilon))
  97     ((and (numberp number1) (numberp number2))
  98      (= number1 number2))
  99     (t (error "~A and ~A are not numbers." number1 number2))))
 100
 101 (defmacro assert-number-equal (expected form &rest extras)
 102   (expand-assert :equal form form expected extras :test #'number-equal))
 103
 104 ;;; (NORMALIZE-FLOAT significand &optional exponent) => significand,exponent
 105 (defun normalize-float (significand &optional (exponent 0))
 106   "Return the normalized floating point number and exponent."
 107   (cond
 108     ((zerop significand)
 109      (values significand 0))
 110     ((>= (abs significand) 10)
 111      (normalize-float (/ significand 10.0) (1+ exponent)))
 112     ((< (abs significand) 1)
 113      (normalize-float (* significand 10.0) (1- exponent)))
 114     (t (values significand exponent))))
 115
 116 ;;; (SIGFIG-EQUAL float1 float2 significant-figures) => true or false
 117 (defun sigfig-equal (float1 float2 &optional (significant-figures *significant-figures*))
 118   "Return true if the floating point numbers have equal significant
 119 figures."
 120   ;; Convert 5 to precision of FLOAT1 and 10 to precision of FLOAT2.
 121   ;; Then, rely on Rule of Float and Rational Contagion, CLHS 12.1.4.1,
 122   ;; to obtain a DELTA of the proper precision.
 123   (let ((delta (* (float 5 float1) (expt (float 10 float2) (- significant-figures)))))
 124     (if (or (zerop float1) (zerop float2))
 125         (< (abs (+ float1 float2)) delta)
 126         (multiple-value-bind (sig1 exp1) (normalize-float float1)
 127           (multiple-value-bind (sig2 exp2) (normalize-float float2)
 128             (and (= exp1 exp2)
 129                  (< (abs (- sig1 sig2)) delta)))))))
 130
 131 (defmacro assert-sigfig-equal (significant-figures expected form &rest extras)
 132   (expand-assert :equal form form expected extras
 133                  :test (lambda (f1 f2) (sigfig-equal f1 f2 significant-figures))))
 134
 135 ;;; (ARRAY-EQUAL array1 array2) => true or false
 136 (defun array-equal (array1 array2 &key (test #'number-equal))
 137   "Return true if the elements of the array are equal."
 138   (when (equal (array-dimensions array1) (array-dimensions array2))
 139     (every test
 140      (make-array (reduce #'* (array-dimensions array1))
 141                  :element-type (array-element-type array1)
 142                  :displaced-to array1)
 143      (make-array (reduce #'* (array-dimensions array2))
 144                  :element-type (array-element-type array2)
 145                  :displaced-to array2))))
 146
 147 (defmacro assert-array-equal (element-test expected form &rest extras)
 148   (expand-assert :equal form form expected extras
 149                  :test `(lambda (a1 a2) (array-equal a1 a2 :test ,element-test))))
 150
 151 ;;; (NUMERICAL-EQUAL result1 result2) => true or false
 152 ;;;
 153 ;;; This is a universal wrapper created by Liam Healy. It is
 154 ;;; implemented to support testing in GSLL. The interface is expanded,
 155 ;;; but backwards compatible with previous versions.
 156 ;;;
 157 (defun numerical-equal (result1 result2 &key (test #'number-equal))
 158   (cond
 159     ((and (numberp result1) (numberp result2))
 160      (funcall test result1 result2))
 161     ((and (typep result1 'sequence) (typep result2 'sequence))
 162      (when (= (length result1) (length result2))
 163        (every (lambda (r1 r2) (numerical-equal r1 r2 :test test))
 164               result1 result2)))
 165     ((and (arrayp result1) (arrayp result2))
 166      (array-equal result1 result2 :test test))
 167     (t (error "~A and/or ~A are not valid arguments." result1 result2))))
 168
 169 (defmacro assert-numerical-equal (expected form &rest extras)
 170   (expand-assert :equal form form expected extras :test #'numerical-equal))