floating-point.lisp

   1 ;;;-*- Mode: Lisp; Syntax: ANSI-Common-Lisp -*-
   2
   3 ;;;; References
   4 ;;;; [NumLinAlg] James W. Demmel "Applied Numerical Linear Algebra",
   5 ;;;;             Society for Industrial and Applied Mathematics, 1997
   6 ;;;;             ISBN: 0-89871-389-7
   7
   8 (common-lisp:in-package :lisp-unit)
   9
  10 (defparameter *epsilon* nil
  11   "Set the error epsilon if the defaults are not acceptable.")
  12
  13 (defparameter *significant-figures* 4
  14   "Default to 4 significant figures.")
  15
  16 ;;; (ROUNDOFF-ERROR x y) => number
  17 ;;; Return the error delta between the exact and approximate floating
  18 ;;; point value.
  19 ;;; [NumLinAlg] : Equation 1.1, pg. 12
  20 (defun roundoff-error (exact approximate)
  21   "Return the error delta between the exact and approximate floating
  22 point value."
  23   (abs (if (or (zerop exact) (zerop approximate))
  24            (+ exact approximate)
  25            (- (/ approximate exact) 1.0))))
  26
  27 ;;; (FLOAT-EQUAL float1 float2 &optional epsilon) => true or false
  28 ;;; Return true if the absolute difference between float1 and float2
  29 ;;; is less than epsilon. If an epsilon is not specified and either
  30 ;;; float1 or float2 is single precision, the single-float-epsilon is
  31 ;;; used.
  32 (defun float-equal (float1 float2 &optional (epsilon *epsilon*))
  33   "Return true if the absolute difference between float1 and float2 is
  34 less than some epsilon."
  35   (and
  36    (floatp float1)
  37    (floatp float2)
  38    (cond
  39      ((and (zerop float1) (zerop float2)))
  40      (epsilon
  41       (> epsilon (roundoff-error float1 float2)))
  42      ((and (typep float1 'double-float) (typep float2 'double-float))
  43       (> (* 2.0 double-float-epsilon) (roundoff-error float1 float2)))
  44      ((or (typep float1 'single-float) (typep float2 'single-float))
  45       (> (* 2.0 single-float-epsilon) (roundoff-error float1 float2)))
  46      (t nil))))
  47
  48 (defmacro assert-float-equal (expected form &rest extras)
  49   (expand-assert :equal form form expected extras :test #'float-equal))
  50
  51 ;;; (COMPLEX-EQUAL complex1 complex2 &optional epsilon) => true or false
  52 ;;; Return true if the absolute difference of the real components and
  53 ;;; the absolute difference of the imaginary components is less then
  54 ;;; epsilon. If an epsilon is not specified and either complex1 or
  55 ;;; complex2 is (complex single-float), the single-float-epsilon is
  56 ;;; used.
  57 (defun complex-equal (complex1 complex2 &optional (epsilon *epsilon*))
  58   "Return true if the absolute difference between Re(complex1),
  59 Re(complex2) and the absolute difference between Im(complex1),
  60 Im(complex2) is less than epsilon."
  61   (and
  62    (typep complex1 '(complex float))
  63    (typep complex2 '(complex float))
  64    (float-equal (realpart complex1) (realpart complex2) epsilon)
  65    (float-equal (imagpart complex1) (imagpart complex2) epsilon)))
  66
  67 (defmacro assert-complex-equal (expected form &rest extras)
  68   (expand-assert :equal form form expected extras :test #'complex-equal))
  69
  70 ;;; (NUMBER-EQUAL number1 number2) => true or false
  71 ;;; Return true if the numbers are equal using the appropriate
  72 ;;; comparison.
  73 (defun number-equal (number1 number2 &optional (epsilon *epsilon*))
  74   "Return true if the numbers are equal using the appropriate
  75 comparison."
  76   (cond
  77     ((and (floatp number1) (floatp number2))
  78      (float-equal number1 number2 epsilon))
  79     ((and (typep number1 '(complex float)) (typep number2 '(complex float)))
  80      (complex-equal number1 number2 epsilon))
  81     ((and (numberp number1) (numberp number2))
  82      (= number1 number2))
  83     (t (error "~A and ~A are not numbers." number1 number2))))
  84
  85 (defmacro assert-number-equal (expected form &rest extras)
  86   (expand-assert :equal form form expected extras :test #'number-equal))
  87
  88 ;;; (NORMALIZE-FLOAT significand &optional exponent) => significand,exponent
  89 (defun normalize-float (significand &optional (exponent 0))
  90   "Return the normalized floating point number and exponent."
  91   (cond
  92     ((zerop significand)
  93      (values significand 0))
  94     ((>= (abs significand) 10)
  95      (normalize-float (/ significand 10.0) (1+ exponent)))
  96     ((< (abs significand) 1)
  97      (normalize-float (* significand 10.0) (1- exponent)))
  98     (t (values significand exponent))))
  99
 100 ;;; (SIGFIG-EQUAL float1 float2 significant-figures) => true or false
 101 (defun sigfig-equal (float1 float2 &optional (significant-figures *significant-figures*))
 102   "Return true if the floating point numbers have equal significant
 103 figures."
 104   ;; Convert 5 to precision of FLOAT1 and 10 to precision of FLOAT2.
 105   ;; Then, rely on Rule of Float and Rational Contagion, CLHS 12.1.4.1,
 106   ;; to obtain a DELTA of the proper precision.
 107   (let ((delta (* (float 5 float1) (expt (float 10 float2) (- significant-figures)))))
 108     (if (or (zerop float1) (zerop float2))
 109         (< (abs (+ float1 float2)) delta)
 110         (multiple-value-bind (sig1 exp1) (normalize-float float1)
 111           (multiple-value-bind (sig2 exp2) (normalize-float float2)
 112             (and (= exp1 exp2)
 113                  (< (abs (- sig1 sig2)) delta)))))))
 114
 115 (defmacro assert-sigfig-equal (significant-figures expected form &rest extras)
 116   (expand-assert :equal form form expected extras
 117                  :test (lambda (f1 f2) (sigfig-equal f1 f2 significant-figures))))
 118
 119 ;;; (ARRAY-EQUAL array1 array2) => true or false
 120 ;;; Return true of the elements of the array are equal.
 121 (defun array-equal (array1 array2 &key (test #'number-equal))
 122   "Return true if the elements of the array are equal."
 123   (when (equal (array-dimensions array1) (array-dimensions array2))
 124     (every test
 125      (make-array (reduce #'* (array-dimensions array1)) :displaced-to array1)
 126      (make-array (reduce #'* (array-dimensions array2)) :displaced-to array2))))
 127
 128 (defmacro assert-array-equal (element-test expected form &rest extras)
 129   (expand-assert :equal form form expected extras
 130                  :test `(lambda (a1 a2) (array-equal a1 a2 :test ,element-test))))
 131
 132 ;;; (NUMERICAL-EQUAL result1 result2) => true or false
 133 ;;; This is a universal wrapper function used by Liam Healy. It is
 134 ;;; implemented to support testing in GSLL. While the interface is
 135 ;;; identical to previous versions, the implementation details are
 136 ;;; slightly different to use the routines here.
 137 (defun numerical-equal (result1 result2 &key (test #'number-equal))
 138   (cond
 139     ((and (numberp result1) (numberp result2))
 140      (funcall test result1 result2))
 141     ((and (typep result1 'sequence) (typep result2 'sequence))
 142      (when (= (length result1) (length result2))
 143        (every (lambda (r1 r2) (numerical-equal r1 r2 :test test))
 144               result1 result2)))
 145     ((and (arrayp result1) (arrayp result2))
 146      (array-equal result1 result2 :test test))
 147     (t (error "~A and/or ~A are not valid arguments." result1 result2))))
 148
 149 (defmacro assert-numerical-equal (expected form &rest extras)
 150   (expand-assert :equal form form expected extras :test #'numerical-equal))