From 8e3ea3e40aa33f5e7b7cd7dd948a7db9792aa3b6 Mon Sep 17 00:00:00 2001 From: "Thomas M. Hermann" Date: Fri, 3 Apr 2009 23:54:49 -0500 Subject: [PATCH] Use DO instead of DOTIMES in SEQUENCE-ERROR and ARRAY-ERROR. This commit contains the following trivial changes. - Renamed ROUNDOFF-ERROR to RELATIVE-ERROR. - Corrected NORMALIZE-FLOAT to return a significand in [0.1,1). Was returning a significand in (1,10]. - Added better references for the functions. - Added an export comment to better categorize exported symbols. --- defpackage.lisp | 3 +- floating-point.lisp | 141 +++++++++++++++++++++++++++++++--------------------- 2 files changed, 85 insertions(+), 59 deletions(-) diff --git a/defpackage.lisp b/defpackage.lisp index 08d599f..1205471 100644 --- a/defpackage.lisp +++ b/defpackage.lisp @@ -37,8 +37,9 @@ OTHER DEALINGS IN THE SOFTWARE. #:logically-equal #:set-equal #:use-debugger #:with-test-listener - ;; Floating point predicates and assertions + ;; Floating point parameters #:*epsilon* #:*significant-figures* + ;; Floating point predicates and assertions #:float-equal #:assert-float-equal #:complex-equal #:assert-complex-equal #:number-equal #:assert-number-equal diff --git a/floating-point.lisp b/floating-point.lisp index b6ad657..029c830 100644 --- a/floating-point.lisp +++ b/floating-point.lisp @@ -23,9 +23,9 @@ ;;; OTHER DEALINGS IN THE SOFTWARE. ;;; ;;; References -;;; [NumLinAlg] James W. Demmel "Applied Numerical Linear Algebra", -;;; Society for Industrial and Applied Mathematics, 1997 -;;; ISBN: 0-89871-389-7 +;;; [NumAlgoC] Gisela Engeln-Mullges and Frank Uhlig "Numerical +;;; Algorithms with C", Springer, 1996 +;;; ISBN: 3-540-60530-4 (common-lisp:in-package :lisp-unit) @@ -35,14 +35,22 @@ (defparameter *significant-figures* 4 "Default to 4 significant figures.") -;;; (ROUNDOFF-ERROR x y) => number -(defun roundoff-error (exact approximate) +;;; (RELATIVE-ERROR x y) => number +;;; [NumAlgoC] : Definition 1.3, pg. 2 +;;; modified with Definition 1.1, pg. 1 +;;; +;;; The definition of relative error in this routine is modified from +;;; the Definition 1.3 in [NumAlgoC] for cases when either the exact +;;; or the approximate value equals zero. According to Definition 1.3, +;;; the relative error is identically equal to 1 in those cases. This +;;; function returns the absolue error in those cases. This is more +;;; useful for testing. +(defun relative-error (exact approximate) "Return the error delta between the exact and approximate floating point value." - ;; [NumLinAlg] : Equation 1.1, pg. 12 (abs (if (or (zerop exact) (zerop approximate)) - (+ exact approximate) - (- (/ approximate exact) 1.0)))) + (- exact approximate) + (/ (- exact approximate) exact)))) ;;; (FLOAT-EQUAL float1 float2 &optional epsilon) => true or false (defun %float-equal (float1 float2 epsilon) @@ -50,18 +58,18 @@ point value." (cond ((and (zerop float1) (zerop float2))) (epsilon - (< (roundoff-error float1 float2) epsilon)) + (< (relative-error float1 float2) epsilon)) ((or (typep float1 'short-float) (typep float2 'short-float)) - (< (roundoff-error float1 float2) (* 2.0 short-float-epsilon))) + (< (relative-error float1 float2) (* 2.0 short-float-epsilon))) ((or (typep float1 'single-float) (typep float2 'single-float)) - (< (roundoff-error float1 float2) (* 2.0 single-float-epsilon))) + (< (relative-error float1 float2) (* 2.0 single-float-epsilon))) ((or (typep float1 'double-float) (typep float2 'double-float)) - (< (roundoff-error float1 float2) (* 2.0 double-float-epsilon))) - (t (< (roundoff-error float1 float2) (* 2.0 long-float-epsilon))))) + (< (relative-error float1 float2) (* 2.0 double-float-epsilon))) + (t (< (relative-error float1 float2) (* 2.0 long-float-epsilon))))) (defun float-equal (float1 float2 &optional (epsilon *epsilon*)) - "Return true if the absolute difference between float1 and float2 is -less than some epsilon." + "Return true if the relative error between float1 and float2 is less +than some epsilon." (when (and (floatp float1) (floatp float2)) (%float-equal float1 float2 epsilon))) @@ -76,9 +84,9 @@ less than some epsilon." (%float-equal (imagpart complex1) (imagpart complex2) epsilon))) (defun complex-equal (complex1 complex2 &optional (epsilon *epsilon*)) - "Return true if the absolute difference between Re(complex1), -Re(complex2) and the absolute difference between Im(complex1), -Im(complex2) is less than epsilon." + "Return true if the relative error between Re(complex1), +Re(complex2) and between Im(complex1), Im(complex2) are each less than +epsilon." (when (and (typep complex1 '(complex float)) (typep complex2 '(complex float))) (%complex-equal complex1 complex2 epsilon))) @@ -103,25 +111,29 @@ comparison." (expand-assert :equal form form expected extras :test #'number-equal)) ;;; (NORMALIZE-FLOAT significand &optional exponent) => significand,exponent +;;; [NumAlgoC] : Definition 1.7, pg. 4 +;;; +;;; To avoid using 0.1, first 1.0 <= significand < 10. On the final +;;; return, scale 0.1 <= significand < 1. (defun normalize-float (significand &optional (exponent 0)) "Return the normalized floating point number and exponent." (cond ((zerop significand) (values significand 0)) - ((>= (abs significand) 10) + ((> (abs significand) 10) (normalize-float (/ significand 10.0) (1+ exponent))) - ((< (abs significand) 1) + ((<= (abs significand) 1) (normalize-float (* significand 10.0) (1- exponent))) - (t (values significand exponent)))) + (t (values (/ significand 10.0) (1+ exponent))))) ;;; (SIGFIG-EQUAL float1 float2 significant-figures) => true or false (defun sigfig-equal (float1 float2 &optional (significant-figures *significant-figures*)) "Return true if the floating point numbers have equal significant figures." - ;; Convert 5 to precision of FLOAT1 and 10 to precision of FLOAT2. + ;; Convert 0.5 to precision of FLOAT1 and 10 to precision of FLOAT2. ;; Then, rely on Rule of Float and Rational Contagion, CLHS 12.1.4.1, ;; to obtain a DELTA of the proper precision. - (let ((delta (* (float 5 float1) (expt (float 10 float2) (- significant-figures))))) + (let ((delta (* (float 0.5 float1) (expt (float 10 float2) (- significant-figures))))) (if (or (zerop float1) (zerop float2)) (< (abs (+ float1 float2)) delta) (multiple-value-bind (sig1 exp1) (normalize-float float1) @@ -145,9 +157,9 @@ figures." :element-type (array-element-type array2) :displaced-to array2)))) -(defmacro assert-array-equal (element-test expected form &rest extras) +(defmacro assert-array-equal (test expected form &rest extras) (expand-assert :equal form form expected extras - :test `(lambda (a1 a2) (array-equal a1 a2 :test ,element-test)))) + :test `(lambda (a1 a2) (array-equal a1 a2 :test ,test)))) ;;; (NUMERICAL-EQUAL result1 result2) => true or false ;;; @@ -175,27 +187,27 @@ figures." ;;; Relative floating point error. (defun float-error (float1 float2) - "Return the roundoff error of the floating point numbers." + "Return the relative error of the floating point numbers." (if (and (floatp float1) (floatp float2)) - (roundoff-error float1 float2) + (relative-error float1 float2) (error "~A and ~A are not float." float1 float2))) ;;; Floating point error as a multiple of epsilon. (defun %float-error-epsilon (float1 float2 epsilon) - "Return the roundoff error divided by epsilon." + "Return the relative error divided by epsilon." (cond (epsilon - (/ (roundoff-error float1 float2) epsilon)) + (/ (relative-error float1 float2) epsilon)) ((or (typep float1 'short-float) (typep float2 'short-float)) - (/ (roundoff-error float1 float2) short-float-epsilon)) + (/ (relative-error float1 float2) short-float-epsilon)) ((or (typep float1 'single-float) (typep float2 'single-float)) - (/ (roundoff-error float1 float2) single-float-epsilon)) + (/ (relative-error float1 float2) single-float-epsilon)) ((or (typep float1 'double-float) (typep float2 'double-float)) - (/ (roundoff-error float1 float2) double-float-epsilon)) - (t (/ (roundoff-error float1 float2) long-float-epsilon)))) + (/ (relative-error float1 float2) double-float-epsilon)) + (t (/ (relative-error float1 float2) long-float-epsilon)))) (defun float-error-epsilon (float1 float2 &optional (epsilon *epsilon*)) - "Return the roundoff error divided by epsilon." + "Return the relative error divided by epsilon." (if (and (floatp float1) (floatp float2)) (%float-error-epsilon float1 float2 epsilon) (error "~A and ~A are not float." float1 float2))) @@ -204,11 +216,11 @@ figures." (defun %complex-error (complex1 complex2) "Internal version with no type checking." (complex - (roundoff-error (realpart complex1) (realpart complex2)) - (roundoff-error (imagpart complex1) (imagpart complex2)))) + (relative-error (realpart complex1) (realpart complex2)) + (relative-error (imagpart complex1) (imagpart complex2)))) (defun complex-error (complex1 complex2) - "Return the roundoff error of each component of the complex numbers." + "Return the relative error of each component of the complex numbers." (if (and (typep complex1 '(complex float)) (typep complex2 '(complex float))) (%complex-error complex1 complex2) @@ -216,33 +228,33 @@ figures." ;;; Complex floating point error as a multiple of epsilon. (defun %complex-error-epsilon (complex1 complex2 epsilon) - "Return the roundoff error of each component divided by epsilon." + "Return the relative error of each component divided by epsilon." (complex (%float-error-epsilon (realpart complex1) (realpart complex2) epsilon) (%float-error-epsilon (imagpart complex1) (imagpart complex2) epsilon))) (defun complex-error-epsilon (complex1 complex2 &optional (epsilon *epsilon*)) - "Return the roundoff error of each component divided by epsilon." + "Return the relative error of each component divided by epsilon." (if (and (typep complex1 '(complex float)) (typep complex2 '(complex float))) (%complex-error-epsilon complex1 complex2 epsilon) (error "~A and ~A are not complex float." complex1 complex2))) -;;; Numeric roundoff error +;;; Numeric relative error (defun number-error (number1 number2) - "Return the roundoff error of the number." + "Return the relative error of the number." (cond ((and (floatp number1) (floatp number2)) - (roundoff-error number1 number2)) + (relative-error number1 number2)) ((and (typep number1 '(complex float)) (typep number2 '(complex float))) (%complex-error number1 number2)) ((and (numberp number1) (numberp number2)) - (abs (- (/ number1 number2) 1))) + (abs (- number1 number2))) (t (error "~A and ~A are not numbers." number1 number2)))) ;;; Numeric error as a multiple of epsilon. (defun number-error-epsilon (number1 number2 &optional (epsilon *epsilon*)) - "Return the roundoff error divided by epsilon." + "Return the relative error divided by epsilon." (cond ((and (floatp number1) (floatp number2)) (%float-error-epsilon number1 number2 epsilon)) @@ -253,13 +265,18 @@ figures." ;;; Sequence errors and the indices. (defun %sequence-error (sequence1 sequence2 test error-function) "Return a sequence of the indice and error between the sequences." - (let ((errseq '())) - (dotimes (index (length sequence1) errseq) - (let ((n1 (elt sequence1 index)) - (n2 (elt sequence2 index))) - (unless (funcall test n1 n2) - (push (list index n1 n2 (funcall error-function n1 n2)) - errseq)))))) + (let ((end (1- (length sequence1))) + (errseq '())) + (do* ((index 0 (1+ index)) + (n1 (elt sequence1 index) (elt sequence1 index)) + (n2 (elt sequence2 index) (elt sequence2 index))) + ((>= index end) + (unless (funcall test n1 n2) + (push (list (1- index) n1 n2 (funcall error-function n1 n2)) + errseq))) + (unless (funcall test n1 n2) + (push (list (1- index) n1 n2 (funcall error-function n1 n2)) + errseq))))) (defun sequence-error (sequence1 sequence2 &key (test #'number-equal) @@ -289,14 +306,22 @@ figures." "Return a list of the indices, values and error of the elements that are not equal." (let ((dimensions (array-dimensions array1)) + (end (1- (array-total-size array1))) (errseq '())) - (dotimes (index (array-total-size array1) errseq) - (let* ((indices (%array-indices index dimensions)) - (n1 (apply #'aref array1 indices)) - (n2 (apply #'aref array2 indices))) - (unless (funcall test n1 n2) - (push (list indices n1 n2 (funcall errfun n1 n2)) - errseq)))))) + (do* ((index 0 (1+ index)) + (indices (%array-indices index dimensions) + (%array-indices index dimensions)) + (n1 (apply #'aref array1 indices) + (apply #'aref array1 indices)) + (n2 (apply #'aref array2 indices) + (apply #'aref array2 indices))) + ((>= index end) + (unless (funcall test n1 n2) + (push (list indices n1 n2 (funcall errfun n1 n2)) + errseq))) + (unless (funcall test n1 n2) + (push (list indices n1 n2 (funcall errfun n1 n2)) + errseq))))) (defun array-error (array1 array2 &key (test #'number-equal) -- 2.11.4.GIT