From 894142c36d25d0e906df98e246f9290970173f69 Mon Sep 17 00:00:00 2001 From: "Thomas M. Hermann" Date: Thu, 4 Jun 2009 14:03:07 -0500 Subject: [PATCH] Expanded the predicate and assertion interface. Expanded the testing interface with RATIONAL-EQUAL, FLOAT-EQUAL and NORM-EQUAL. These are implemented as generic functions that operate on atoms, sequences and arrays. In the case of sequences and arrays, RATIONAL-EQUAL and FLOAT-EQUAL perform and element-wise comparison for equality. NORM-EQUAL performs the test for equality on the relative error norm of the sequence or array. --- defpackage.lisp | 17 +- floating-point.lisp | 582 +++++++++++++++++++++++++++++++++++++++------------- lisp-unit.asd | 2 + rational.lisp | 88 ++++++++ 4 files changed, 541 insertions(+), 148 deletions(-) create mode 100644 rational.lisp diff --git a/defpackage.lisp b/defpackage.lisp index 396578a..42e8068 100644 --- a/defpackage.lisp +++ b/defpackage.lisp @@ -37,18 +37,21 @@ OTHER DEALINGS IN THE SOFTWARE. #:logically-equal #:set-equal #:use-debugger #:with-test-listener + ;; Rational predicates and assertions + #:rational-equal #:assert-rational-equal ;; Floating point parameters - #:*epsilon* #:*significant-figures* + #:*measure* #:*epsilon* #:*significant-figures* + ;; Floating point functions + #:default-epsilon #:relative-error + #:sumsq #:sump #:norm + #:relative-error-norm ;; Floating point predicates and assertions - #:float-equal #:assert-float-equal - #:complex-equal #:assert-complex-equal - #:number-equal #:assert-number-equal + #:float-equal #:assert-float-equal #:sigfig-equal #:assert-sigfig-equal + #:norm-equal #:assert-norm-equal + #:number-equal #:assert-number-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 ;; Floating point data functions #:make-2d-list diff --git a/floating-point.lisp b/floating-point.lisp index 8548c74..f81aed9 100644 --- a/floating-point.lisp +++ b/floating-point.lisp @@ -29,13 +29,89 @@ (common-lisp:in-package :lisp-unit) +(defparameter *measure* 1) + (defparameter *epsilon* nil "Set the error epsilon if the defaults are not acceptable.") (defparameter *significant-figures* 4 "Default to 4 significant figures.") -;;; (RELATIVE-ERROR x y) => number +(defgeneric default-epsilon (value) + (:documentation + "Return the default epsilon for the value.")) + +(defgeneric relative-error (exact approximate) + (:documentation + "Return the relative-error between the 2 quantities.")) + +(defgeneric float-equal (data1 data2 &optional epsilon) + (:documentation + "Return true if the floating point data is equal.")) + +(defgeneric sumsq (data) + (:documentation + "Return the scaling parameter and the sum of the squares of the ~ + data.")) + +(defgeneric sump (data p) + (:documentation + "Return the scaling parameter and the sum of the powers of p of the ~ + data.")) + +(defgeneric norm (data &optional measure) + (:documentation + "Return the element-wise norm of the data.")) + +(defgeneric relative-error-norm (exact approximate &optional measure) + (:documentation + "Return the relative error norm ")) + +(defgeneric norm-equal (data1 data2 &optional epsilon measure) + (:documentation + "Return true if the norm of the data is equal.")) + +(defgeneric numerical-equal (result1 result2 &key test) + (:documentation + "Return true if the results are numerically equal according to :TEST.")) + +;;; (DEFAULT-EPSILON value) => epsilon +(defmethod default-epsilon ((value float)) + "Return a default epsilon value based on the floating point type." + (typecase value + (short-float (* 2.0s0 short-float-epsilon)) + (single-float (* 2.0f0 single-float-epsilon)) + (double-float (* 2.0d0 double-float-epsilon)) + (long-float (* 2.0l0 long-float-epsilon)))) + +(defmethod default-epsilon ((value complex)) + "Return a default epsilon value based on the complex type." + (typecase value + ((complex short-float) (* 2.0s0 short-float-epsilon)) + ((complex single-float) (* 2.0f0 single-float-epsilon)) + ((complex double-float) (* 2.0d0 double-float-epsilon)) + ((complex long-float) (* 2.0l0 long-float-epsilon)) + (t 0))) + +(defmethod default-epsilon ((value list)) + "Return the default epsilon based on contents of the list." + (reduce (lambda (x y) (max x (default-epsilon y))) + value :initial-value 0)) + +(defmethod default-epsilon ((value vector)) + "Return the default epsilon based on the contents of the vector." + (reduce (lambda (x y) (max x (default-epsilon y))) + value :initial-value 0)) + +(defmethod default-epsilon ((value array)) + "Return the default epsilon based on the contents of the array." + (reduce (lambda (x y) (max x (default-epsilon y))) + (make-array (array-total-size value) + :element-type (array-element-type value) + :displaced-to value) + :initial-value 0)) + +;;; (RELATIVE-ERROR x y) => float ;;; [NumAlgoC] : Definition 1.3, pg. 2 ;;; modified with Definition 1.1, pg. 1 ;;; @@ -45,72 +121,363 @@ ;;; 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." +(defun %relative-error (exact approximate) + "Return the relative error of the numbers." (abs (if (or (zerop exact) (zerop approximate)) (- exact approximate) (/ (- exact approximate) exact)))) -;;; (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 - (< (relative-error float1 float2) epsilon)) - ((or (typep float1 'short-float) (typep float2 'short-float)) - (< (relative-error float1 float2) (* 2.0 short-float-epsilon))) - ((or (typep float1 'single-float) (typep float2 'single-float)) - (< (relative-error float1 float2) (* 2.0 single-float-epsilon))) - ((or (typep float1 'double-float) (typep float2 'double-float)) - (< (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 relative error between float1 and float2 is less -than some epsilon." - (when (and (floatp float1) (floatp float2)) - (%float-equal float1 float2 epsilon))) +(defmethod relative-error ((exact float) (approximate float)) + "Return the error delta between the exact and approximate floating +point value." + (%relative-error exact approximate)) + +(defmethod relative-error ((exact float) (approximate complex)) + "Return the relative error between the float and complex number." + (%relative-error exact approximate)) + +(defmethod relative-error ((exact complex) (approximate float)) + "Return the relative error between the float and complex number." + (%relative-error exact approximate)) + +(defmethod relative-error ((exact complex) (approximate complex)) + "Return the relative error of the complex numbers." + (if (or (typep exact '(complex float)) + (typep approximate '(complex float))) + (%relative-error exact approximate) + (error "Relative error is only applicable to complex values with ~ + floating point parts."))) + +;;; (FLOAT-EQUAL data1 data2 epsilon) => true or false +(defun %float-equal (data1 data2 epsilon) + "Return true if the relative error between the data is less than +epsilon." + (or + (and (zerop data1) (zerop data2)) + (< (%relative-error data1 data2) epsilon))) + +(defmethod float-equal ((data1 float) (data2 float) + &optional (epsilon *epsilon*)) + "Return true if the relative error between data1 and data2 is less +than epsilon." + (%float-equal data1 data2 + (or epsilon (max (default-epsilon data1) + (default-epsilon data2))))) + +(defmethod float-equal ((data1 float) (data2 rational) + &optional (epsilon *epsilon*)) + "Return true if the relative error between data1 and data2 is less +than epsilon." + (%float-equal data1 (float data2 data1) + (or epsilon (default-epsilon data1)))) + +(defmethod float-equal ((data1 rational) (data2 float) + &optional (epsilon *epsilon*)) + "Return true if the relative error between data1 and data2 is less +than epsilon." + (%float-equal (float data1 data2) data2 + (or epsilon (default-epsilon data2)))) + +(defmethod float-equal ((data1 float) (data2 complex) + &optional (epsilon *epsilon*)) + "Return true if the relative error between data1 and data2 is less +than epsilon." + (%float-equal data1 data2 + (or epsilon (max (default-epsilon data1) + (default-epsilon data2))))) + +(defmethod float-equal ((data1 complex) (data2 float) + &optional (epsilon *epsilon*)) + "Return true if the relative error between data1 and data2 is less +than epsilon." + (%float-equal data1 data2 + (or epsilon (max (default-epsilon data1) + (default-epsilon data2))))) + +(defmethod float-equal ((data1 complex) (data2 complex) + &optional (epsilon *epsilon*)) + "Return true if the relative error between data1 and data2 is less +than epsilon." + (< (relative-error data1 data2) + (or epsilon (max (default-epsilon data1) + (default-epsilon data2))))) + +(defun %seq-float-equal (seq1 seq2 epsilon) + "Return true if the element-wise comparison of relative error is +less than epsilon." + (or + (and (null seq1) (null seq2)) + (when (= (length seq1) (length seq2)) + (every + (lambda (d1 d2) (float-equal d1 d2 epsilon)) seq1 seq2)))) + +(defmethod float-equal ((data1 list) (data2 list) + &optional (epsilon *epsilon*)) + "Return true if the lists are equal in length and element-wise +comparison of the relative error is less than epsilon." + (%seq-float-equal data1 data2 epsilon)) + +(defmethod float-equal ((data1 list) (data2 vector) + &optional (epsilon *epsilon*)) + "Return true if the vector and the list are equal in length and +element-wise comparison of the relative error is less than epsilon." + (%seq-float-equal data1 data2 epsilon)) + +(defmethod float-equal ((data1 vector) (data2 list) + &optional (epsilon *epsilon*)) + "Return true if the vector and the list are equal in length and +element-wise comparison of the relative error is less than epsilon." + (%seq-float-equal data1 data2 epsilon)) + +(defmethod float-equal ((data1 vector) (data2 vector) + &optional (epsilon *epsilon*)) + "Return true if the vectors are equal in length and element-wise +comparison of the relative error is less than epsilon." + (%seq-float-equal data1 data2 epsilon)) + +(defmethod float-equal ((data1 array) (data2 array) + &optional (epsilon *epsilon*)) + "Return true if the arrays are equal in length and element-wise +comparison of the relative error is less than epsilon." + (when (equal (array-dimensions data1) + (array-dimensions data2)) + (%seq-float-equal + (make-array (array-total-size data1) + :element-type (array-element-type data1) + :displaced-to data1) + (make-array (array-total-size data2) + :element-type (array-element-type data2) + :displaced-to data2) + epsilon))) (defmacro assert-float-equal (expected form &rest extras) (expand-assert :equal form form expected extras :test #'float-equal)) -;;; (COMPLEX-EQUAL complex1 complex2 &optional epsilon) => true or false -(defun %complex-equal (complex1 complex2 epsilon) - "Internal version that does not verify arguments as (complex float)." - (and - (%float-equal (realpart complex1) (realpart complex2) epsilon) - (%float-equal (imagpart complex1) (imagpart complex2) epsilon))) - -(defun complex-equal (complex1 complex2 &optional (epsilon *epsilon*)) - "Return true if the relative error between Re(complex1), -Re(complex2) and between Im(complex1), Im(complex2) are each less than -epsilon." +;;; (SUMSQ data) => scale, sumsq +(defmethod sumsq ((data list)) + "Return the scaling parameter and the sum of the squares of the ~ + list." + (let ((scale 0) (sumsq 1) + (abs-val nil)) + (dolist (elm data (values scale sumsq)) + (when (< 0 (setf abs-val (abs elm))) + (if (< scale abs-val) + (setf sumsq (1+ (* sumsq (expt (/ scale abs-val) 2))) + scale abs-val) + (incf sumsq (expt (/ elm scale) 2))))))) + +(defmethod sumsq ((data vector)) + "Return the scaling parameter and the sum of the squares of the ~ + vector." + (let ((scale 0) (sumsq 1) + (size (length data)) + (abs-val nil)) + (dotimes (index size (values scale sumsq)) + (when (< 0 (setf abs-val (abs (svref data index)))) + (if (< scale abs-val) + (setf sumsq (1+ (* sumsq (expt (/ scale abs-val) 2))) + scale abs-val) + (incf sumsq (expt (/ (svref data index) scale) 2))))))) + +(defmethod sumsq ((data array)) + "Return the scaling parameter and the sum of the squares of the ~ + array." + (sumsq (make-array (array-total-size data) + :element-type (array-element-type data) + :displaced-to data))) + +;;; (SUMP data) => scale, sump +(defmethod sump ((data list) (p real)) + "Return the scaling parameter and the sum of the powers of p of the ~ + data." + (let ((scale 0) (sump 1) + (abs-val nil)) + (dolist (elm data (values scale sump)) + (when (< 0 (setf abs-val (abs elm))) + (if (< scale abs-val) + (setf sump (1+ (* sump (expt (/ scale abs-val) p))) + scale abs-val) + (incf sump (expt (/ elm scale) p))))))) + +(defmethod sump ((data vector) (p real)) + "Return the scaling parameter and the sum of the powers of p of the ~ + vector." + (let ((scale 0) (sump 1) + (size (length data)) + (abs-val nil)) + (dotimes (index size (values scale sump)) + (when (< 0 (setf abs-val (abs (svref data index)))) + (if (< scale abs-val) + (setf sump (1+ (* sump (expt (/ scale abs-val) p))) + scale abs-val) + (incf sump (expt (/ (svref data index) scale) p))))))) + +(defmethod sump ((data array) (p real)) + "Return the scaling parameter and the sum of the powers of p of the ~ + array." + (sump (make-array (array-total-size data) + :element-type (array-element-type data) + :displaced-to data) + p)) + +;;; (NORM data) => float +(defun %seq-1-norm (data) + "Return the Taxicab norm of the sequence." + (reduce (lambda (x y) (+ x (abs y))) + data :initial-value 0)) + +(defun %seq-2-norm (data) + "Return the Euclidean norm of the sequence." + (multiple-value-bind (scale sumsq) + (sumsq (map-into (make-array (length data)) #'abs data)) + (* scale (sqrt sumsq)))) + +(defun %seq-p-norm (data p) + "Return the p norm of the sequence." + (multiple-value-bind (scale sump) + (sump (map-into (make-array (length data)) #'abs data) p) + (* scale (expt sump (/ p))))) + +(defun %seq-inf-norm (data) + "Return the infinity, or maximum, norm of the sequence." + (reduce (lambda (x y) (max x (abs y))) + data :initial-value 0)) + +(defun %seq-norm (data measure) + "Return the norm of the sequence according to the measure." (cond - ((and (typep complex1 '(complex float)) - (typep complex2 '(complex float))) - (%complex-equal complex1 complex2 epsilon)) - ((and (complexp complex1) (complexp complex2)) - (= complex1 complex2)))) - -(defmacro assert-complex-equal (expected form &rest extras) - (expand-assert :equal form form expected extras :test #'complex-equal)) - -;;; (NUMBER-EQUAL number1 number2) => true or false -(defun number-equal (number1 number2 &optional - (epsilon *epsilon*) type-agreement) - "Return true if the numbers are equal using the appropriate -comparison." - (and - (or (not type-agreement) (eq (type-of number1) (type-of number2))) - (%complex-equal - (coerce number1 '(complex double-float)) - (coerce number2 '(complex double-float)) - epsilon))) - -(defmacro assert-number-equal (expected form &rest extras) - (expand-assert :equal form form expected extras :test #'number-equal)) + ((equalp measure 1) + (%seq-1-norm data)) + ((equalp measure 2) + (%seq-2-norm data)) + ((numberp measure) + (%seq-p-norm data measure)) + ((equalp measure :infinity) + (%seq-inf-norm data)) + (t (error "Unrecognized norm, ~A." measure)))) + +(defmethod norm ((data list) &optional (measure *measure*)) + "Return the norm of the list according to the measure." + (%seq-norm data measure)) + +(defmethod norm ((data vector) &optional (measure *measure*)) + "Return the norm of the vector according to the measure." + (%seq-norm data measure)) + +(defmethod norm ((data array) &optional (measure *measure*)) + "Return the entrywise norm of the array according to the measure." + (let ((flat-data (make-array (array-total-size data) + :element-type (array-element-type data) + :displaced-to data))) + (cond + ((and (numberp measure) (< 0 measure)) + (warn "Measure ~D results in an entrywise p-norm." measure) + (%seq-p-norm flat-data measure)) + ((equalp measure :frobenius) + (%seq-2-norm flat-data)) + ((equalp measure :max) + (%seq-inf-norm flat-data)) + (t (error "Unrecognized norm, ~A." measure))))) + +;;; (RELATIVE-ERROR-NORM exact approximate measure) => float +(defun %relative-error-norm (exact approximate measure) + "Return the relative error norm of the sequences." + (/ (norm (map-into (make-array (length exact)) + (lambda (x1 x2) (abs (- x1 x2))) + exact approximate) measure) + (norm exact measure))) + +(defmethod relative-error-norm ((exact list) (approximate list) + &optional (measure *measure*)) + "Return the relative error norm of the lists." + (if (= (length exact) (length approximate)) + (%relative-error-norm exact approximate measure) + (error "Lists are not equal in length."))) + +(defmethod relative-error-norm ((exact list) (approximate vector) + &optional (measure *measure*)) + "Return the relative error norm of the list and the vector." + (if (= (length exact) (length approximate)) + (%relative-error-norm exact approximate measure) + (error "The list and vector are not equal in length."))) + +(defmethod relative-error-norm ((exact vector) (approximate list) + &optional (measure *measure*)) + "Return the relative error norm of the list and the vector." + (if (= (length exact) (length approximate)) + (%relative-error-norm exact approximate measure) + (error "The list and vector are not equal in length."))) + +(defmethod relative-error-norm ((exact vector) (approximate vector) + &optional (measure *measure*)) + "Return the relative error norm of the vectors." + (if (= (length exact) (length approximate)) + (%relative-error-norm exact approximate measure) + (error "Vectors are not equal in length."))) + +(defmethod relative-error-norm ((exact array) (approximate vector) + &optional (measure *measure*)) + "Return the relative error norm of the arrays." + (if (equal (array-dimensions exact) + (array-dimensions approximate)) + (%relative-error-norm + (make-array (array-total-size exact) + :element-type (array-element-type exact) + :displaced-to exact) + (make-array (array-total-size approximate) + :element-type (array-element-type approximate) + :displaced-to approximate) + measure) + (error "Arrays are not equal dimensions."))) + +;;; (NORM-EQUAL data1 data2 epsilon measure) => boolean +(defun %norm-equal (seq1 seq2 epsilon measure) + "Return true if the relative error norm is less than epsilon." + (or + (and (null seq1) (null seq2)) + (< (%relative-error-norm seq1 seq2 measure) epsilon))) + +(defmethod norm-equal ((data1 list) (data2 list) &optional + (epsilon *epsilon*) (measure *measure*)) + "Return true if the lists are equal in length and the relative error +norm is less than epsilon." + (%norm-equal data1 data2 epsilon measure)) + +(defmethod norm-equal ((data1 list) (data2 vector) &optional + (epsilon *epsilon*) (measure *measure*)) + "Return true if the vector and the list are equal in length and the +relative error norm is less than epsilon." + (%norm-equal data1 data2 epsilon measure)) + +(defmethod norm-equal ((data1 vector) (data2 list) &optional + (epsilon *epsilon*) (measure *measure*)) + "Return true if the vector and the list are equal in length and the +relative error norm is less than epsilon." + (%norm-equal data1 data2 epsilon measure)) + +(defmethod norm-equal ((data1 vector) (data2 vector) &optional + (epsilon *epsilon*) (measure *measure*)) + "Return true if the vectors are equal in length and the relative +error norm is less than epsilon." + (%norm-equal data1 data2 epsilon measure)) + +(defmethod norm-equal ((data1 array) (data2 array) &optional + (epsilon *epsilon*) (measure *measure*)) + "Return true if the arrays are equal in length and the relative +error norm is less than epsilon." + (when (equal (array-dimensions data1) + (array-dimensions data2)) + (%norm-equal + (make-array (array-total-size data1) + :element-type (array-element-type data1) + :displaced-to data1) + (make-array (array-total-size data2) + :element-type (array-element-type data2) + :displaced-to data2) + epsilon measure))) + +(defmacro assert-norm-equal (expected form &rest extras) + (expand-assert :equal form form expected extras :test #'norm-equal)) ;;; (NORMALIZE-FLOAT significand &optional exponent) => significand,exponent ;;; [NumAlgoC] : Definition 1.7, pg. 4 @@ -129,6 +496,7 @@ comparison." (t (values (/ significand 10.0) (1+ exponent))))) ;;; (SIGFIG-EQUAL float1 float2 significant-figures) => true or false +;;; TODO : Convert to a generic function and specialize for sequences. (defun sigfig-equal (float1 float2 &optional (significant-figures *significant-figures*)) "Return true if the floating point numbers have equal significant figures." @@ -146,16 +514,25 @@ figures." (defmacro assert-sigfig-equal (expected form &rest extras) (expand-assert :equal form form expected extras :test #'sigfig-equal)) +;;; (NUMBER-EQUAL number1 number2) => true or false +(defun number-equal (number1 number2 &optional (epsilon *epsilon*) type-eq-p) + "Return true if the numbers are equal within some epsilon, +optionally requiring the types to be identical." + (and + (or (not type-eq-p) (eq (type-of number1) (type-of number2))) + (float-equal (coerce number1 '(complex double-float)) + (coerce number2 '(complex double-float)) + epsilon))) + +(defmacro assert-number-equal (expected form &rest extras) + (expand-assert :equal form form expected extras :test #'number-equal)) + ;;; (NUMERICAL-EQUAL result1 result2) => true or false ;;; ;;; This is a universal wrapper created by Liam Healy. It is ;;; implemented to support testing in GSLL. The interface is expanded, ;;; but backwards compatible with previous versions. ;;; -(defgeneric numerical-equal (result1 result2 &key test) - (:documentation - "Return true if the results are numerically equal according to :TEST.")) - (defmethod numerical-equal ((result1 number) (result2 number) &key (test #'number-equal)) "Return true if the the numbers are equal according to :TEST." @@ -211,83 +588,6 @@ is equal according to :TEST." ;;; Diagnostic functions ;;; Failing a unit test is only half the problem. -;;; Relative floating point error. -(defun float-error (float1 float2) - "Return the relative error of the floating point numbers." - (if (and (floatp float1) (floatp 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 relative error divided by epsilon." - (cond - (epsilon - (/ (relative-error float1 float2) epsilon)) - ((or (typep float1 'short-float) (typep float2 'short-float)) - (/ (relative-error float1 float2) short-float-epsilon)) - ((or (typep float1 'single-float) (typep float2 'single-float)) - (/ (relative-error float1 float2) single-float-epsilon)) - ((or (typep float1 'double-float) (typep float2 'double-float)) - (/ (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 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))) - -;;; Relative complex floating point error. -(defun %complex-error (complex1 complex2) - "Internal version with no type checking." - (complex - (relative-error (realpart complex1) (realpart complex2)) - (relative-error (imagpart complex1) (imagpart complex2)))) - -(defun complex-error (complex1 complex2) - "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) - (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 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 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 relative error -(defun number-error (number1 number2) - "Return the relative error of the number." - (cond - ((and (floatp number1) (floatp 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))) - (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 relative 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." @@ -302,7 +602,7 @@ is equal according to :TEST." (defun sequence-error (sequence1 sequence2 &key (test #'number-equal) - (error-function #'number-error)) + (error-function #'relative-error)) "Return a sequence of the indice and error between the sequence elements." (cond ((not (typep sequence1 'sequence)) @@ -341,7 +641,7 @@ are not equal." (defun array-error (array1 array2 &key (test #'number-equal) - (error-function #'number-error)) + (error-function #'relative-error)) "Return a list of the indices and error between the array elements." (cond ((not (arrayp array1)) diff --git a/lisp-unit.asd b/lisp-unit.asd index e7c5f4a..5c5b38d 100644 --- a/lisp-unit.asd +++ b/lisp-unit.asd @@ -34,6 +34,8 @@ :components ((:file "defpackage") (:file "lisp-unit" :depends-on ("defpackage")) + (:file "rational" + :depends-on ("defpackage")) (:file "floating-point" :depends-on ("defpackage" "lisp-unit")))) diff --git a/rational.lisp b/rational.lisp new file mode 100644 index 0000000..e31e74a --- /dev/null +++ b/rational.lisp @@ -0,0 +1,88 @@ +;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp -*- +;;; +;;; Floating tests and assertions for LISP-UNIT +;;; +;;; Copyright (c) 2009 Thomas M. Hermann +;;; +;;; Permission is hereby granted, free of charge, to any person obtaining +;;; a copy of this software and associated documentation files (the "Software"), +;;; to deal in the Software without restriction, including without limitation +;;; the rights to use, copy, modify, merge, publish, distribute, sublicense, +;;; and/or sell copies of the Software, and to permit persons to whom the +;;; Software is furnished to do so, subject to the following conditions: +;;; +;;; The above copyright notice and this permission notice shall be included +;;; in all copies or substantial portions of the Software. +;;; +;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +;;; OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +;;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +;;; THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR +;;; OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, +;;; ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +;;; OTHER DEALINGS IN THE SOFTWARE. +;;; + +(common-lisp:in-package :lisp-unit) + +(defgeneric rational-equal (data1 data2) + (:documentation + "Return true if the rational data is equal.")) + +;;; (RATIONAL-EQUAL data1 data2) => boolean +(defmethod rational-equal ((data1 rational) (data2 rational)) + "Return true if the rational numbers are equal." + (= data1 data2)) + +(defmethod rational-equal ((data1 complex) (data2 complex)) + "Return true if the complex parts are rational and equal." + (if (and (typep data1 '(complex rational)) + (typep data2 '(complex rational))) + (= data1 data2) + (error "Rational equality is not applicable to complex values ~ + with floating point parts."))) + +(defun %seq-rational-equal (seq1 seq2) + "Return true if the sequences are equal length and at each position +the corresponding elements are equal." + (or + (and (null seq1) (null seq2)) + (and + (= (length seq1) (length seq2)) + (every (lambda (d1 d2) (rational-equal d1 d2)) seq1 seq2)))) + +(defmethod rational-equal ((data1 list) (data2 list)) + "Return true if the lists are equal in length and element-wise +equal." + (%seq-rational-equal data1 data2)) + +(defmethod rational-equal ((data1 list) (data2 vector)) + "Return true if the vector and the list are equal in length and +element-wise equal." + (%seq-rational-equal data1 data2)) + +(defmethod rational-equal ((data1 vector) (data2 list)) + "Return true if the vector and the list are equal in length and +element-wise equal." + (%seq-rational-equal data1 data2)) + +(defmethod rational-equal ((data1 vector) (data2 vector)) + "Return true if the vectors are equal in length and element-wise +equal." + (%seq-rational-equal data1 data2)) + +(defmethod rational-equal ((data1 array) (data2 array)) + "Return true if the arrays are equal in dimension and element-wise +equal." + (when (equal (array-dimensions data1) + (array-dimensions data2)) + (%seq-rational-equal + (make-array (array-total-size data1) + :element-type (array-element-type data1) + :displaced-to data1) + (make-array (array-total-size data2) + :element-type (array-element-type data2) + :displaced-to data2)))) + +(defmacro assert-rational-equal (expected form &rest extras) + (expand-assert :equal form form expected extras :test #'rational-equal)) -- 2.11.4.GIT