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 ;;; (ROUNDOFF-ERROR x y) => number
  11 ;;; Return the error delta between the exact and approximate floating
  12 ;;; point value.
  13 ;;; [NumLinAlg] : Equation 1.1, pg. 12
  14 (defun roundoff-error (exact approximate)
  15   "Returned the error delta between the exact and approximate floating
  16 point value."
  17   (abs (if (or (zerop exact) (zerop approximate))
  18            (+ exact approximate)
  19            (- (/ approximate exact) 1.0))))
  20
  21 ;;; (FLOAT-EQUAL float1 float2 &optional epsilon) => true or false
  22 ;;; Return true if the absolute difference between float1 and float2
  23 ;;; is less than epsilon. If an epsilon is not specified and either
  24 ;;; float1 or float2 is single precision, the single-float-epsilon is
  25 ;;; used.
  26 (defun float-equal (float1 float2 &optional epsilon)
  27   "Return true if the absolute difference between float1 and float2 is
  28 less than some epsilon."
  29   (and
  30    (floatp float1)
  31    (floatp float2)
  32    (cond
  33      ((and (zerop float1) (zerop float2)))
  34      (epsilon
  35       (> epsilon (roundoff-error float1 float2)))
  36      ((and (typep float1 'double-float) (typep float2 'double-float))
  37       (> (* 2.0 double-float-epsilon) (roundoff-error float1 float2)))
  38      ((or (typep float1 'single-float) (typep float2 'single-float))
  39       (> (* 2.0 single-float-epsilon) (roundoff-error float1 float2)))
  40      (t nil))))
  41
  42 ;;; (COMPLEX-EQUAL complex1 complex2 &optional epsilon) => true or false
  43 ;;; Return true if the absolute difference of the real components and
  44 ;;; the absolute difference of the imaginary components is less then
  45 ;;; epsilon. If an epsilon is not specified and either complex1 or
  46 ;;; complex2 is (complex single-float), the single-float-epsilon is
  47 ;;; used.
  48 (defun complex-equal (complex1 complex2 &optional epsilon)
  49   "Return true if the absolute difference between Re(complex1),
  50 Re(complex2) and the absolute difference between Im(complex1),
  51 Im(complex2) is less than epsilon."
  52   (and
  53    (typep complex1 '(complex float))
  54    (typep complex2 '(complex float))
  55    (float-equal (realpart complex1) (realpart complex2) epsilon)
  56    (float-equal (imagpart complex1) (imagpart complex2) epsilon)))
  57
  58 ;;; (NUMBER-EQUAL number1 number2) => true or false
  59 ;;; Return true if the numbers are equal using the appropriate
  60 ;;; comparison.
  61 (defun number-equal (number1 number2 &optional epsilon)
  62   "Return true if the numbers are equal using the appropriate
  63 comparison."
  64   (cond
  65     ((and (floatp number1) (floatp number2))
  66      (float-equal number1 number2 epsilon))
  67     ((and (typep number1 '(complex float)) (typep number2 '(complex float)))
  68      (complex-equal number1 number2 epsilon))
  69     ((and (numberp number1) (numberp number2))
  70      (= number1 number2))
  71     (t (error "~A and ~A are not numbers." number1 number2))))
  72
  73 ;;; (ELEMENT-EQUAL array1 array2 indice dimensions) => true or false
  74 ;;; A utility function for ARRAY-EQUAL.
  75 (defun element-equal (array1 array2 indices dimensions &key (test #'number-equal))
  76   "Return true if the index of array1 equals array2."
  77   (let* ((rank (first dimensions))
  78          (remaining (rest dimensions))
  79          (update-result
  80           (if remaining
  81               (lambda (index)
  82                 (element-equal array1 array2
  83                                (cons index indices) remaining :test test))
  84               (lambda (index)
  85                 (funcall test
  86                          (apply #'aref array1 index (reverse indices))
  87                          (apply #'aref array2 index (reverse indices)))))))
  88     (do ((index 0 (1+ index))
  89          (result t (funcall update-result index)))
  90         ((or (not result) (>= index rank)) result))))
  91
  92 ;;; (ARRAY-EQUAL array1 array2) => true or false
  93 ;;; Return true of the elements of the array are equal.
  94 (defun array-equal (array1 array2 &key (test #'number-equal))
  95   "Return true if the elements of the array are equal."
  96   (when (equal (array-dimensions array1) (array-dimensions array2))
  97     (element-equal array1 array2 nil (array-dimensions array1) :test test)))
  98
  99 ;;; (NORMALIZE-FLOAT significand &optional exponent) => significand,exponent
 100 (defun normalize-float (significand &optional (exponent 0))
 101   "Return the normalized floating point number and exponent."
 102   (cond
 103     ((zerop significand)
 104      (values significand 0))
 105     ((>= (abs significand) 10)
 106      (normalize-float (/ significand 10.0) (1+ exponent)))
 107     ((< (abs significand) 1)
 108      (normalize-float (* significand 10.0) (1- exponent)))
 109     (t (values significand exponent))))
 110
 111 ;;; (SIGNIFICANT-FIGURES-EQUAL float1 float2 significant-figures) => true or false
 112 (defun significant-figures-equal (float1 float2 significant-figures)
 113   "Return true if the floating point numbers have equal significant
 114 figures."
 115   ;; Convert 5 to precision of FLOAT1 and 10 to precision of
 116   ;; FLOAT2. Then, rely on Rule of Float and Rational Contagion, CLHS
 117   ;; 12.1.4.1, to obtain a DELTA of the proper precision.
 118   (let ((delta (* (float 5 float1) (expt (float 10 float2) (- significant-figures)))))
 119     (if (or (zerop float1) (zerop float2))
 120         (< (abs (+ float1 float2)) delta)
 121         (multiple-value-bind (sig1 exp1) (normalize-float float1)
 122           (multiple-value-bind (sig2 exp2) (normalize-float float2)
 123             (and (= exp1 exp2)
 124                  (< (abs (- sig1 sig2)) delta)))))))
 125
 126 (defun 2-sigfig-equal (float1 float2)
 127   "Return true if the floats are equal to 2 significant figures."
 128   (significant-figures-equal float1 float2 2))
 129
 130 (defun 3-sigfig-equal (float1 float2)
 131   "Return true if the floats are equal to 3 significant figures."
 132   (significant-figures-equal float1 float2 3))
 133
 134 (defun 4-sigfig-equal (float1 float2)
 135   "Return true if the floats are equal to 4 significant figures."
 136   (significant-figures-equal float1 float2 4))
 137
 138 (defun 5-sigfig-equal (float1 float2)
 139   "Return true if the floats are equal to 5 significant figures."
 140   (significant-figures-equal float1 float2 5))
 141
 142 (defun 6-sigfig-equal (float1 float2)
 143   "Return true if the floats are equal to 6 significant figures."
 144   (significant-figures-equal float1 float2 6))