From c91614523fc460da6e09ff226eb3ee2a0f272145 Mon Sep 17 00:00:00 2001 From: "Thomas M. Hermann" Date: Sun, 15 Mar 2009 23:25:51 -0500 Subject: [PATCH] Implemented diagnostic functions for numbers, sequences and arrays. Two types of diagnostic functions are implemented for numbers. The first type returns the roundoff error between 2 numbers. The second type returns the roundoff error as a multiple of the appropriate precision float-epsilon value. For sequences and arrays, the diagnostic functions return a list of the elements that are not equal, noting the element indices, values and numeric error. Each function has 2 keywords. First, a :TEST keyword to specify the equality test, default to NUMBER-EQUAL. Second, an :ERROR-FUNCTION keyword specifying which type of error should be reported, defaulting to #'NUMBER-ERROR. NOTE: This is an initial implementation and subject to change. --- defpackage.lisp | 7 ++- floating-point.lisp | 160 +++++++++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 157 insertions(+), 10 deletions(-) diff --git a/defpackage.lisp b/defpackage.lisp index b5fcb2d..9f6ff44 100644 --- a/defpackage.lisp +++ b/defpackage.lisp @@ -22,6 +22,11 @@ #:number-equal #:assert-number-equal #:sigfig-equal #:assert-sigfig-equal #:array-equal #:assert-array-equal - #:numerical-equal #:assert-numerical-equal)) + #:numerical-equal #:assert-numerical-equal + ;; Floating point diagnostic functions + #:float-error #:float-error-epsilon + #:complex-error #:complex-error-epsilon + #:number-error #:number-error-epsilon + #:sequence-error #:array-error)) (pushnew :lisp-unit common-lisp:*features*) diff --git a/floating-point.lisp b/floating-point.lisp index afcc033..557d4db 100644 --- a/floating-point.lisp +++ b/floating-point.lisp @@ -53,19 +53,20 @@ point value." (+ exact approximate) (- (/ approximate exact) 1.0)))) - ;;; (FLOAT-EQUAL float1 float2 &optional epsilon) => true or false (defun %float-equal (float1 float2 epsilon) "Internal version that does not verify arguments as floats." (cond ((and (zerop float1) (zerop float2))) (epsilon - (> epsilon (roundoff-error float1 float2))) - ((and (typep float1 'double-float) (typep float2 'double-float)) - (> (* 2.0 double-float-epsilon) (roundoff-error float1 float2))) + (< (roundoff-error float1 float2) epsilon)) + ((or (typep float1 'short-float) (typep float2 'short-float)) + (< (roundoff-error float1 float2) (* 2.0 short-float-epsilon))) ((or (typep float1 'single-float) (typep float2 'single-float)) - (> (* 2.0 single-float-epsilon) (roundoff-error float1 float2))) - (t nil))) + (< (roundoff-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))))) (defun float-equal (float1 float2 &optional (epsilon *epsilon*)) "Return true if the absolute difference between float1 and float2 is @@ -146,10 +147,10 @@ figures." "Return true if the elements of the array are equal." (when (equal (array-dimensions array1) (array-dimensions array2)) (every test - (make-array (reduce #'* (array-dimensions array1)) + (make-array (array-total-size array1) :element-type (array-element-type array1) :displaced-to array1) - (make-array (reduce #'* (array-dimensions array2)) + (make-array (array-total-size array2) :element-type (array-element-type array2) :displaced-to array2)))) @@ -173,7 +174,148 @@ figures." result1 result2))) ((and (arrayp result1) (arrayp result2)) (array-equal result1 result2 :test test)) - (t (error "~A and/or ~A are not valid arguments." result1 result2)))) + (t (error "~A and ~A are not valid arguments." result1 result2)))) (defmacro assert-numerical-equal (expected form &rest extras) (expand-assert :equal form form expected extras :test #'numerical-equal)) + +;;; Diagnostic functions +;;; Failing a unit test is only half the problem. + +;;; Relative floating point error. +(defun float-error (float1 float2) + "Return the roundoff error of the floating point numbers." + (if (and (floatp float1) (floatp float2)) + (roundoff-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." + (cond + (epsilon + (/ (roundoff-error float1 float2) epsilon)) + ((or (typep float1 'short-float) (typep float2 'short-float)) + (/ (roundoff-error float1 float2) short-float-epsilon)) + ((or (typep float1 'single-float) (typep float2 'single-float)) + (/ (roundoff-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)))) + +(defun float-error-epsilon (float1 float2 &optional (epsilon *epsilon*)) + "Return the roundoff 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))) + +;;; Relative complex floating point error. +(defun %complex-error (complex1 complex2) + "Internal version with no type checking." + (complex + (roundoff-error (realpart complex1) (realpart complex2)) + (roundoff-error (imagpart complex1) (imagpart complex2)))) + +(defun complex-error (complex1 complex2) + "Return the roundoff error of each component of the complex numbers." + (if (and (typep complex1 '(complex float)) + (typep complex2 '(complex float))) + (%complex-error complex1 complex2) + (error "~A and ~A are not complex float." complex1 complex2))) + +;;; 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." + (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." + (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 +(defun number-error (number1 number2) + "Return the roundoff error of the number." + (cond + ((and (floatp number1) (floatp number2)) + (roundoff-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))) + (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." + (cond + ((and (floatp number1) (floatp number2)) + (%float-error-epsilon number1 number2 epsilon)) + ((and (typep number1 '(complex float)) (typep number2 '(complex float))) + (%complex-error-epsilon number1 number2 epsilon)) + (t (error "~A and ~A are not float or complex float numbers." number1 number2)))) + +;;; 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)))))) + +(defun sequence-error (sequence1 sequence2 &key + (test #'number-equal) + (error-function #'number-error)) + "Return a sequence of the indice and error between the sequence elements." + (cond + ((not (typep sequence1 'sequence)) + (error "SEQUENCE1 is not a valid sequence.")) + ((not (typep sequence2 'sequence)) + (error "SEQUENCE2 is not a valid sequence.")) + ((not (= (length sequence1) (length sequence2))) + (error "Lengths not equal. SEQUENCE1(~D) /= SEQUENCE2(~D)." + (length sequence1) (length sequence2))) + (t (%sequence-error sequence1 sequence2 test error-function)))) + +;;; Array errors and the indices. +(defun %array-indices (row-major-index dimensions) + "Recursively calculate the indices from the row major index." + (let ((remaining (rest dimensions))) + (if remaining + (multiple-value-bind (index remainder) + (floor row-major-index (reduce #'* remaining)) + (cons index (%array-indices remainder remaining))) + (cons row-major-index nil)))) + +(defun %array-error (array1 array2 test errfun) + "Return a list of the indices, values and error of the elements that +are not equal." + (let ((dimensions (array-dimensions 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)))))) + +(defun array-error (array1 array2 &key + (test #'number-equal) + (error-function #'number-error)) + "Return a list of the indices and error between the array elements." + (cond + ((not (arrayp array1)) + (error "ARRAY1 is not an array.")) + ((not (arrayp array2)) + (error "ARRAY2 is not an array.")) + ((not (equal (array-dimensions array1) (array-dimensions array2))) + (error "Arrays are not equal dimensions.")) + (t (%array-error array1 array2 test error-function)))) -- 2.11.4.GIT