floating-point.lisp

   1 ;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp -*-
   2 ;;;
   3 ;;; Floating tests and assertions for LISP-UNIT
   4 ;;;
   5 ;;; Copyright (c) 2009 Thomas M. Hermann
   6 ;;;
   7 ;;; Permission is hereby granted, free of charge, to any person obtaining
   8 ;;; a copy of this software and associated documentation files (the "Software"),
   9 ;;; to deal in the Software without restriction, including without limitation
  10 ;;; the rights to use, copy, modify, merge, publish, distribute, sublicense,
  11 ;;; and/or sell copies of the Software, and to permit persons to whom the
  12 ;;; Software is furnished to do so, subject to the following conditions:
  13 ;;;
  14 ;;; The above copyright notice and this permission notice shall be included
  15 ;;; in all copies or substantial portions of the Software.
  16 ;;;
  17 ;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  18 ;;; OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19 ;;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  20 ;;; THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
  21 ;;; OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  22 ;;; ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  23 ;;; OTHER DEALINGS IN THE SOFTWARE.
  24 ;;;
  25 ;;; References
  26 ;;; [NumAlgoC] Gisela Engeln-Mullges and Frank Uhlig "Numerical
  27 ;;;            Algorithms with C", Springer, 1996
  28 ;;;            ISBN: 3-540-60530-4
  29
  30 (common-lisp:in-package :lisp-unit)
  31
  32 (defparameter *measure* 1)
  33
  34 (defparameter *epsilon* nil
  35   "Set the error epsilon if the defaults are not acceptable.")
  36
  37 (defparameter *significant-figures* 4
  38   "Default to 4 significant figures.")
  39
  40 (defgeneric default-epsilon (value)
  41   (:documentation
  42    "Return the default epsilon for the value."))
  43
  44 (defgeneric relative-error (exact approximate)
  45   (:documentation
  46    "Return the relative-error between the 2 quantities."))
  47
  48 (defgeneric float-equal (data1 data2 &optional epsilon)
  49   (:documentation
  50    "Return true if the floating point data is equal."))
  51
  52 (defgeneric sumsq (data)
  53   (:documentation
  54    "Return the scaling parameter and the sum of the squares of the ~
  55     data."))
  56
  57 (defgeneric sump (data p)
  58   (:documentation
  59    "Return the scaling parameter and the sum of the powers of p of the ~
  60     data."))
  61
  62 (defgeneric norm (data &optional measure)
  63   (:documentation
  64    "Return the element-wise norm of the data."))
  65
  66 (defgeneric relative-error-norm (exact approximate &optional measure)
  67   (:documentation
  68    "Return the relative error norm "))
  69
  70 (defgeneric norm-equal (data1 data2 &optional epsilon measure)
  71   (:documentation
  72    "Return true if the norm of the data is equal."))
  73
  74 (defgeneric sigfig-equal (data1 data2 &optional significant-figures)
  75   (:documentation
  76    "Return true if the data have equal significant figures."))
  77
  78 (defgeneric numerical-equal (result1 result2 &key test)
  79   (:documentation
  80    "Return true if the results are numerically equal according to :TEST."))
  81
  82 ;;; (DEFAULT-EPSILON value) => epsilon
  83 (defmethod default-epsilon ((value float))
  84   "Return a default epsilon value based on the floating point type."
  85   (typecase value
  86     (short-float  (* 2.0s0 short-float-epsilon))
  87     (single-float (* 2.0f0 single-float-epsilon))
  88     (double-float (* 2.0d0 double-float-epsilon))
  89     (long-float   (* 2.0l0 long-float-epsilon))))
  90
  91 (defmethod default-epsilon ((value complex))
  92   "Return a default epsilon value based on the complex type."
  93   (typecase value
  94     ((complex short-float)  (* 2.0s0 short-float-epsilon))
  95     ((complex single-float) (* 2.0f0 single-float-epsilon))
  96     ((complex double-float) (* 2.0d0 double-float-epsilon))
  97     ((complex long-float)   (* 2.0l0 long-float-epsilon))
  98     (t 0)))
  99
 100 (defmethod default-epsilon ((value list))
 101   "Return the default epsilon based on contents of the list."
 102   (reduce (lambda (x y) (max x (default-epsilon y)))
 103           value :initial-value 0))
 104
 105 (defmethod default-epsilon ((value vector))
 106   "Return the default epsilon based on the contents of the vector."
 107   (reduce (lambda (x y) (max x (default-epsilon y)))
 108           value :initial-value 0))
 109
 110 (defmethod default-epsilon ((value array))
 111   "Return the default epsilon based on the contents of the array."
 112   (reduce (lambda (x y) (max x (default-epsilon y)))
 113           (make-array (array-total-size value)
 114                       :element-type (array-element-type value)
 115                       :displaced-to value)
 116           :initial-value 0))
 117
 118 ;;; (RELATIVE-ERROR x y) => float
 119 ;;; [NumAlgoC] : Definition 1.3, pg. 2
 120 ;;;              modified with Definition 1.1, pg. 1
 121 ;;;
 122 ;;; The definition of relative error in this routine is modified from
 123 ;;; the Definition 1.3 in [NumAlgoC] for cases when either the exact
 124 ;;; or the approximate value equals zero. According to Definition 1.3,
 125 ;;; the relative error is identically equal to 1 in those cases. This
 126 ;;; function returns the absolue error in those cases. This is more
 127 ;;; useful for testing.
 128 (defun %relative-error (exact approximate)
 129   "Return the relative error of the numbers."
 130   (abs (if (or (zerop exact) (zerop approximate))
 131            (- exact approximate)
 132            (/ (- exact approximate) exact))))
 133
 134 (defmethod relative-error ((exact float) (approximate float))
 135   "Return the error delta between the exact and approximate floating
 136 point value."
 137   (%relative-error exact approximate))
 138
 139 (defmethod relative-error ((exact float) (approximate complex))
 140   "Return the relative error between the float and complex number."
 141   (%relative-error exact approximate))
 142
 143 (defmethod relative-error ((exact complex) (approximate float))
 144   "Return the relative error between the float and complex number."
 145   (%relative-error exact approximate))
 146
 147 (defmethod relative-error ((exact complex) (approximate complex))
 148   "Return the relative error of the complex numbers."
 149   (if (or (typep exact '(complex float))
 150           (typep approximate '(complex float)))
 151       (%relative-error exact approximate)
 152       (error "Relative error is only applicable to complex values with ~
 153               floating point parts.")))
 154
 155 ;;; (FLOAT-EQUAL data1 data2 epsilon) => true or false
 156 (defun %float-equal (data1 data2 epsilon)
 157   "Return true if the relative error between the data is less than
 158 epsilon."
 159   (or
 160    (and (zerop data1) (zerop data2))
 161    (< (%relative-error data1 data2) epsilon)))
 162
 163 (defmethod float-equal ((data1 float) (data2 float)
 164                         &optional (epsilon *epsilon*))
 165   "Return true if the relative error between data1 and data2 is less
 166 than epsilon."
 167   (%float-equal data1 data2
 168                 (or epsilon (max (default-epsilon data1)
 169                                  (default-epsilon data2)))))
 170
 171 (defmethod float-equal ((data1 float) (data2 rational)
 172                         &optional (epsilon *epsilon*))
 173   "Return true if the relative error between data1 and data2 is less
 174 than epsilon."
 175   (%float-equal data1 (float data2 data1)
 176                 (or epsilon (default-epsilon data1))))
 177
 178 (defmethod float-equal ((data1 rational) (data2 float)
 179                         &optional (epsilon *epsilon*))
 180   "Return true if the relative error between data1 and data2 is less
 181 than epsilon."
 182   (%float-equal (float data1 data2) data2
 183                 (or epsilon (default-epsilon data2))))
 184
 185 (defmethod float-equal ((data1 float) (data2 complex)
 186                         &optional (epsilon *epsilon*))
 187   "Return true if the relative error between data1 and data2 is less
 188 than epsilon."
 189   (%float-equal data1 data2
 190                 (or epsilon (max (default-epsilon data1)
 191                                  (default-epsilon data2)))))
 192
 193 (defmethod float-equal ((data1 complex) (data2 float)
 194                         &optional (epsilon *epsilon*))
 195   "Return true if the relative error between data1 and data2 is less
 196 than epsilon."
 197   (%float-equal data1 data2
 198                 (or epsilon (max (default-epsilon data1)
 199                                  (default-epsilon data2)))))
 200
 201 (defmethod float-equal ((data1 complex) (data2 complex)
 202                         &optional (epsilon *epsilon*))
 203   "Return true if the relative error between data1 and data2 is less
 204 than epsilon."
 205   (< (relative-error data1 data2)
 206      (or epsilon (max (default-epsilon data1)
 207                       (default-epsilon data2)))))
 208
 209 (defun %seq-float-equal (seq1 seq2 epsilon)
 210   "Return true if the element-wise comparison of relative error is
 211 less than epsilon."
 212   (or
 213    (and (null seq1) (null seq2))
 214    (when (= (length seq1) (length seq2))
 215      (every
 216       (lambda (d1 d2) (float-equal d1 d2 epsilon)) seq1 seq2))))
 217
 218 (defmethod float-equal ((data1 list) (data2 list)
 219                         &optional (epsilon *epsilon*))
 220   "Return true if the lists are equal in length and element-wise
 221 comparison of the relative error is less than epsilon."
 222   (%seq-float-equal data1 data2 epsilon))
 223
 224 (defmethod float-equal ((data1 list) (data2 vector)
 225                         &optional (epsilon *epsilon*))
 226   "Return true if the vector and the list are equal in length and
 227 element-wise comparison of the relative error is less than epsilon."
 228   (%seq-float-equal data1 data2 epsilon))
 229
 230 (defmethod float-equal ((data1 vector) (data2 list)
 231                         &optional (epsilon *epsilon*))
 232   "Return true if the vector and the list are equal in length and
 233 element-wise comparison of the relative error is less than epsilon."
 234   (%seq-float-equal data1 data2 epsilon))
 235
 236 (defmethod float-equal ((data1 vector) (data2 vector)
 237                         &optional (epsilon *epsilon*))
 238   "Return true if the vectors are equal in length and element-wise
 239 comparison of the relative error is less than epsilon."
 240   (%seq-float-equal data1 data2 epsilon))
 241
 242 (defmethod float-equal ((data1 array) (data2 array)
 243                         &optional (epsilon *epsilon*))
 244   "Return true if the arrays are equal in length and element-wise
 245 comparison of the relative error is less than epsilon."
 246   (when (equal (array-dimensions data1)
 247                (array-dimensions data2))
 248     (%seq-float-equal
 249      (make-array (array-total-size data1)
 250                  :element-type (array-element-type data1)
 251                  :displaced-to data1)
 252      (make-array (array-total-size data2)
 253                  :element-type (array-element-type data2)
 254                  :displaced-to data2)
 255      epsilon)))
 256
 257 (defmacro assert-float-equal (expected form &rest extras)
 258   (expand-assert :equal form form expected extras :test #'float-equal))
 259
 260 ;;; (SUMSQ data) => scale, sumsq
 261 (defmethod sumsq ((data list))
 262   "Return the scaling parameter and the sum of the squares of the ~
 263    list."
 264   (let ((scale 0) (sumsq 1)
 265         (abs-val nil))
 266     (dolist (elm data (values scale sumsq))
 267       (when (< 0 (setf abs-val (abs elm)))
 268         (if (< scale abs-val)
 269             (setf sumsq (1+ (* sumsq (expt (/ scale abs-val) 2)))
 270                   scale abs-val)
 271             (incf sumsq (expt (/ elm scale) 2)))))))
 272
 273 (defmethod sumsq ((data vector))
 274   "Return the scaling parameter and the sum of the squares of the ~
 275    vector."
 276   (let ((scale 0) (sumsq 1)
 277         (size (length data))
 278         (abs-val nil))
 279     (dotimes (index size (values scale sumsq))
 280       (when (< 0 (setf abs-val (abs (svref data index))))
 281         (if (< scale abs-val)
 282             (setf sumsq (1+ (* sumsq (expt (/ scale abs-val) 2)))
 283                   scale abs-val)
 284             (incf sumsq (expt (/ (svref data index) scale) 2)))))))
 285
 286 (defmethod sumsq ((data array))
 287   "Return the scaling parameter and the sum of the squares of the ~
 288    array."
 289   (sumsq (make-array (array-total-size data)
 290                      :element-type (array-element-type data)
 291                      :displaced-to data)))
 292
 293 ;;; (SUMP data) => scale, sump
 294 (defmethod sump ((data list) (p real))
 295   "Return the scaling parameter and the sum of the powers of p of the ~
 296    data."
 297   (let ((scale 0) (sump 1)
 298         (abs-val nil))
 299     (dolist (elm data (values scale sump))
 300       (when (< 0 (setf abs-val (abs elm)))
 301         (if (< scale abs-val)
 302             (setf sump  (1+ (* sump (expt (/ scale abs-val) p)))
 303                   scale abs-val)
 304             (incf sump (expt (/ elm scale) p)))))))
 305
 306 (defmethod sump ((data vector) (p real))
 307   "Return the scaling parameter and the sum of the powers of p of the ~
 308    vector."
 309   (let ((scale 0) (sump 1)
 310         (size (length data))
 311         (abs-val nil))
 312     (dotimes (index size (values scale sump))
 313       (when (< 0 (setf abs-val (abs (svref data index))))
 314         (if (< scale abs-val)
 315             (setf sump  (1+ (* sump (expt (/ scale abs-val) p)))
 316                   scale abs-val)
 317             (incf sump (expt (/ (svref data index) scale) p)))))))
 318
 319 (defmethod sump ((data array) (p real))
 320   "Return the scaling parameter and the sum of the powers of p of the ~
 321    array."
 322   (sump (make-array (array-total-size data)
 323                     :element-type (array-element-type data)
 324                     :displaced-to data)
 325         p))
 326
 327 ;;; (NORM data) => float
 328 (defun %seq-1-norm (data)
 329   "Return the Taxicab norm of the sequence."
 330   (reduce (lambda (x y) (+ x (abs y)))
 331           data :initial-value 0))
 332
 333 (defun %seq-2-norm (data)
 334   "Return the Euclidean norm of the sequence."
 335   (multiple-value-bind (scale sumsq)
 336       (sumsq (map-into (make-array (length data)) #'abs data))
 337     (* scale (sqrt sumsq))))
 338
 339 (defun %seq-p-norm (data p)
 340   "Return the p norm of the sequence."
 341   (multiple-value-bind (scale sump)
 342       (sump (map-into (make-array (length data)) #'abs data) p)
 343     (* scale (expt sump (/ p)))))
 344
 345 (defun %seq-inf-norm (data)
 346   "Return the infinity, or maximum, norm of the sequence."
 347   (reduce (lambda (x y) (max x (abs y)))
 348           data :initial-value 0))
 349
 350 (defun %seq-norm (data measure)
 351   "Return the norm of the sequence according to the measure."
 352   (cond
 353     ((equalp measure 1)
 354      (%seq-1-norm data))
 355     ((equalp measure 2)
 356      (%seq-2-norm data))
 357     ((numberp measure)
 358      (%seq-p-norm data measure))
 359     ((equalp measure :infinity)
 360      (%seq-inf-norm data))
 361     (t (error "Unrecognized norm, ~A." measure))))
 362
 363 (defmethod norm ((data list) &optional (measure *measure*))
 364   "Return the norm of the list according to the measure."
 365   (%seq-norm data measure))
 366
 367 (defmethod norm ((data vector) &optional (measure *measure*))
 368   "Return the norm of the vector according to the measure."
 369   (%seq-norm data measure))
 370
 371 (defmethod norm ((data array) &optional (measure *measure*))
 372   "Return the entrywise norm of the array according to the measure."
 373   (let ((flat-data (make-array (array-total-size data)
 374                                :element-type (array-element-type data)
 375                                :displaced-to data)))
 376     (cond
 377       ((and (numberp measure) (< 0 measure))
 378        (warn "Measure ~D results in an entrywise p-norm." measure)
 379        (%seq-p-norm flat-data measure))
 380       ((equalp measure :frobenius)
 381        (%seq-2-norm flat-data))
 382       ((equalp measure :max)
 383        (%seq-inf-norm flat-data))
 384       (t (error "Unrecognized norm, ~A." measure)))))
 385
 386 ;;; (RELATIVE-ERROR-NORM exact approximate measure) => float
 387 (defun %relative-error-norm (exact approximate measure)
 388   "Return the relative error norm of the sequences."
 389   (/ (norm (map-into (make-array (length exact))
 390                      (lambda (x1 x2) (abs (- x1 x2)))
 391                      exact approximate) measure)
 392      (norm exact measure)))
 393
 394 (defmethod relative-error-norm ((exact list) (approximate list)
 395                                 &optional (measure *measure*))
 396   "Return the relative error norm of the lists."
 397   (if (= (length exact) (length approximate))
 398       (%relative-error-norm exact approximate measure)
 399       (error "Lists are not equal in length.")))
 400
 401 (defmethod relative-error-norm ((exact list) (approximate vector)
 402                                 &optional (measure *measure*))
 403   "Return the relative error norm of the list and the vector."
 404   (if (= (length exact) (length approximate))
 405       (%relative-error-norm exact approximate measure)
 406       (error "The list and vector are not equal in length.")))
 407
 408 (defmethod relative-error-norm ((exact vector) (approximate list)
 409                                 &optional (measure *measure*))
 410   "Return the relative error norm of the list and the vector."
 411   (if (= (length exact) (length approximate))
 412       (%relative-error-norm exact approximate measure)
 413       (error "The list and vector are not equal in length.")))
 414
 415 (defmethod relative-error-norm ((exact vector) (approximate vector)
 416                                 &optional (measure *measure*))
 417   "Return the relative error norm of the vectors."
 418   (if (= (length exact) (length approximate))
 419       (%relative-error-norm exact approximate measure)
 420       (error "Vectors are not equal in length.")))
 421
 422 (defmethod relative-error-norm ((exact array) (approximate vector)
 423                                 &optional (measure *measure*))
 424   "Return the relative error norm of the arrays."
 425   (if (equal (array-dimensions exact)
 426              (array-dimensions approximate))
 427       (%relative-error-norm
 428        (make-array (array-total-size exact)
 429                    :element-type (array-element-type exact)
 430                    :displaced-to exact)
 431        (make-array (array-total-size approximate)
 432                    :element-type (array-element-type approximate)
 433                    :displaced-to approximate)
 434        measure)
 435       (error "Arrays are not equal dimensions.")))
 436
 437 ;;; (NORM-EQUAL data1 data2 epsilon measure) => boolean
 438 (defun %norm-equal (seq1 seq2 epsilon measure)
 439   "Return true if the relative error norm is less than epsilon."
 440   (or
 441    (and (null seq1) (null seq2))
 442    (< (%relative-error-norm seq1 seq2 measure) epsilon)))
 443
 444 (defmethod norm-equal ((data1 list) (data2 list) &optional
 445                        (epsilon *epsilon*) (measure *measure*))
 446   "Return true if the lists are equal in length and the relative error
 447 norm is less than epsilon."
 448   (%norm-equal data1 data2 epsilon measure))
 449
 450 (defmethod norm-equal ((data1 list) (data2 vector) &optional
 451                        (epsilon *epsilon*) (measure *measure*))
 452   "Return true if the vector and the list are equal in length and the
 453 relative error norm is less than epsilon."
 454   (%norm-equal data1 data2 epsilon measure))
 455
 456 (defmethod norm-equal ((data1 vector) (data2 list) &optional
 457                        (epsilon *epsilon*) (measure *measure*))
 458   "Return true if the vector and the list are equal in length and the
 459 relative error norm is less than epsilon."
 460   (%norm-equal data1 data2 epsilon measure))
 461
 462 (defmethod norm-equal ((data1 vector) (data2 vector) &optional
 463                        (epsilon *epsilon*) (measure *measure*))
 464   "Return true if the vectors are equal in length and the relative
 465 error norm is less than epsilon."
 466   (%norm-equal data1 data2 epsilon measure))
 467
 468 (defmethod norm-equal ((data1 array) (data2 array) &optional
 469                        (epsilon *epsilon*) (measure *measure*))
 470   "Return true if the arrays are equal in length and the relative
 471 error norm is less than epsilon."
 472   (when (equal (array-dimensions data1)
 473                (array-dimensions data2))
 474     (%norm-equal
 475      (make-array (array-total-size data1)
 476                  :element-type (array-element-type data1)
 477                  :displaced-to data1)
 478      (make-array (array-total-size data2)
 479                  :element-type (array-element-type data2)
 480                  :displaced-to data2)
 481      epsilon measure)))
 482
 483 (defmacro assert-norm-equal (expected form &rest extras)
 484   (expand-assert :equal form form expected extras :test #'norm-equal))
 485
 486 ;;; (NORMALIZE-FLOAT significand &optional exponent) => significand,exponent
 487 ;;; [NumAlgoC] : Definition 1.7, pg. 4
 488 ;;;
 489 ;;; To avoid using 0.1, first 1.0 <= significand < 10. On the final
 490 ;;; return, scale 0.1 <= significand < 1.
 491 (defun %normalize-float (significand &optional (exponent 0))
 492   "Return the normalized floating point number and exponent."
 493   (cond
 494     ((zerop significand)
 495      (values significand 0))
 496     ((>= (abs significand) 10)
 497      (%normalize-float (/ significand 10.0) (1+ exponent)))
 498     ((< (abs significand) 1)
 499      (%normalize-float (* significand 10.0) (1- exponent)))
 500     (t (values (/ significand 10.0) (1+ exponent)))))
 501
 502 ;;; (SIGFIG-EQUAL float1 float2 significant-figures) => true or false
 503 (defun %sigfig-equal (float1 float2 significant-figures)
 504   "Return true if the floating point numbers have equal significant
 505 figures."
 506   ;; Convert 0.5 to precision of FLOAT1 and 10 to precision of FLOAT2.
 507   ;; Then, rely on Rule of Float and Rational Contagion, CLHS 12.1.4.1,
 508   ;; to obtain a DELTA of the proper precision.
 509   (let ((delta (* (float 0.5 float1) (expt (float 10 float2) (- significant-figures)))))
 510     (if (or (zerop float1) (zerop float2))
 511         (< (abs (+ float1 float2)) delta)
 512         (multiple-value-bind (sig1 exp1) (%normalize-float float1)
 513           (multiple-value-bind (sig2 exp2) (%normalize-float float2)
 514             (and (= exp1 exp2)
 515                  (< (abs (- sig1 sig2)) delta)))))))
 516
 517 (defmethod sigfig-equal ((data1 float) (data2 float) &optional
 518                          (significant-figures *significant-figures*))
 519   "Return true if the floating point numbers have equal significant
 520 figures."
 521   (%sigfig-equal data1 data2 significant-figures))
 522
 523 (defun %seq-sigfig-equal (seq1 seq2 significant-figures)
 524   "Return true if the element-wise comparison is equal to the
 525 specified significant figures."
 526   (or
 527    (and (null seq1) (null seq2))
 528    (when (= (length seq1) (length seq2))
 529      (every
 530       (lambda (d1 d2) (sigfig-equal d1 d2 significant-figures))
 531       seq1 seq2))))
 532
 533 (defmethod sigfig-equal ((data1 list) (data2 list) &optional
 534                          (significant-figures *significant-figures*))
 535   "Return true if the lists are equal in length and the element-wise
 536 comparison is equal to significant figures."
 537   (%seq-sigfig-equal data1 data2 significant-figures))
 538
 539 (defmethod sigfig-equal ((data1 vector) (data2 list) &optional
 540                          (significant-figures *significant-figures*))
 541   "Return true if the vector and the list are equal in length and the
 542 element-wise comparison is equal to significant figures."
 543   (%seq-sigfig-equal data1 data2 significant-figures))
 544
 545 (defmethod sigfig-equal ((data1 list) (data2 vector) &optional
 546                          (significant-figures *significant-figures*))
 547   "Return true if the list and the vector are equal in length and the
 548 element-wise comparison is equal to significant figures."
 549   (%seq-sigfig-equal data1 data2 significant-figures))
 550
 551 (defmethod sigfig-equal ((data1 vector) (data2 vector) &optional
 552                          (significant-figures *significant-figures*))
 553   "Return true if the vectors are equal in length and the element-wise
 554 comparison is equal to significant figures."
 555   (%seq-sigfig-equal data1 data2 significant-figures))
 556
 557 (defmethod sigfig-equal ((data1 array) (data2 array) &optional
 558                          (significant-figures *significant-figures*))
 559   "Return true if the arrays are equal in length and the element-wise
 560 comparison is equal to significant figures."
 561   (when (equal (array-dimensions data1)
 562                (array-dimensions data2))
 563     (%seq-sigfig-equal
 564      (make-array (array-total-size data1)
 565                  :element-type (array-element-type data1)
 566                  :displaced-to data1)
 567      (make-array (array-total-size data2)
 568                  :element-type (array-element-type data2)
 569                  :displaced-to data2)
 570      significant-figures)))
 571
 572 (defmacro assert-sigfig-equal (expected form &rest extras)
 573   (expand-assert :equal form form expected extras :test #'sigfig-equal))
 574
 575 ;;; (NUMBER-EQUAL number1 number2) => true or false
 576 (defun number-equal (number1 number2 &optional (epsilon *epsilon*) type-eq-p)
 577   "Return true if the numbers are equal within some epsilon,
 578 optionally requiring the types to be identical."
 579   (and
 580    (or (not type-eq-p) (eq (type-of number1) (type-of number2)))
 581    (float-equal (coerce number1 '(complex double-float))
 582                 (coerce number2 '(complex double-float))
 583                 epsilon)))
 584
 585 (defmacro assert-number-equal (expected form &rest extras)
 586   (expand-assert :equal form form expected extras :test #'number-equal))
 587
 588 ;;; (NUMERICAL-EQUAL result1 result2) => true or false
 589 ;;;
 590 ;;; This is a universal wrapper created by Liam Healy. It is
 591 ;;; implemented to support testing in GSLL. The interface is expanded,
 592 ;;; but backwards compatible with previous versions.
 593 ;;;
 594 (defmethod numerical-equal ((result1 number) (result2 number)
 595                             &key (test #'number-equal))
 596   "Return true if the the numbers are equal according to :TEST."
 597   (funcall test result1 result2))
 598
 599 (defun %sequence-equal (seq1 seq2 test)
 600   "Return true if the sequences are equal in length and each element
 601 is equal according to :TEST."
 602   (when (= (length seq1) (length seq2))
 603     (every (lambda (s1 s2) (numerical-equal s1 s2 :test test))
 604            seq1 seq2)))
 605
 606 (defmethod numerical-equal ((result1 list) (result2 list)
 607                             &key (test #'number-equal))
 608   "Return true if the lists are equal in length and each element is
 609 equal according to :TEST."
 610   (%sequence-equal result1 result2 test))
 611
 612 (defmethod numerical-equal ((result1 vector) (result2 vector)
 613                             &key (test #'number-equal))
 614   "Return true if the vectors are equal in length and each element is
 615 equal according to :TEST."
 616   (%sequence-equal result1 result2 test))
 617
 618 (defmethod numerical-equal ((result1 list) (result2 vector)
 619                             &key (test #'number-equal))
 620   "Return true if every element of the list is equal to the
 621 corresponding element of the vector."
 622   (%sequence-equal result1 result2 test))
 623
 624 (defmethod numerical-equal ((result1 vector) (result2 list)
 625                             &key (test #'number-equal))
 626   "Return true if every element of the list is equla to the
 627 corresponding element of the vector."
 628   (%sequence-equal result1 result2 test))
 629
 630 (defmethod numerical-equal ((result1 array) (result2 array)
 631                             &key (test #'number-equal))
 632   "Return true if the arrays are equal in dimension and each element
 633 is equal according to :TEST."
 634   (when (equal (array-dimensions result1) (array-dimensions result2))
 635     (every test
 636            (make-array (array-total-size result1)
 637                        :element-type (array-element-type result1)
 638                        :displaced-to result1)
 639            (make-array (array-total-size result2)
 640                        :element-type (array-element-type result2)
 641                        :displaced-to result2))))
 642
 643 (defmacro assert-numerical-equal (expected form &rest extras)
 644   (expand-assert :equal form form expected extras :test #'numerical-equal))
 645
 646 ;;; Diagnostic functions
 647 ;;; Failing a unit test is only half the problem.
 648
 649 ;;; Sequence errors and the indices.
 650 (defun %sequence-error (sequence1 sequence2 test error-function)
 651   "Return a sequence of the indice and error between the sequences."
 652   (let ((n1 nil) (n2 nil)
 653         (errseq '()))
 654     (dotimes (index (length sequence1) errseq)
 655       (setf n1 (elt sequence1 index)
 656             n2 (elt sequence2 index))
 657       (unless (funcall test n1 n2)
 658         (push (list (1- index) n1 n2 (funcall error-function n1 n2))
 659               errseq)))))
 660
 661 (defun sequence-error (sequence1 sequence2 &key
 662                        (test #'number-equal)
 663                        (error-function #'relative-error))
 664   "Return a sequence of the indice and error between the sequence elements."
 665   (cond
 666     ((not (typep sequence1 'sequence))
 667      (error "SEQUENCE1 is not a valid sequence."))
 668     ((not (typep sequence2 'sequence))
 669      (error "SEQUENCE2 is not a valid sequence."))
 670     ((not (= (length sequence1) (length sequence2)))
 671      (error "Lengths not equal. SEQUENCE1(~D) /= SEQUENCE2(~D)."
 672             (length sequence1) (length sequence2)))
 673     (t (%sequence-error sequence1 sequence2 test error-function))))
 674
 675 ;;; Array errors and the indices.
 676 (defun %array-indices (row-major-index dimensions)
 677   "Recursively calculate the indices from the row major index."
 678   (let ((remaining (rest dimensions)))
 679     (if remaining
 680         (multiple-value-bind (index remainder)
 681             (floor row-major-index (reduce #'* remaining))
 682           (cons index (%array-indices remainder remaining)))
 683         (cons row-major-index nil))))
 684
 685 (defun %array-error (array1 array2 test errfun)
 686   "Return a list of the indices, values and error of the elements that
 687 are not equal."
 688   (let ((dimensions (array-dimensions array1))
 689         (n1 nil) (n2 nil)
 690         (indices '())
 691         (errseq '()))
 692     (dotimes (index (array-total-size array1) errseq)
 693       (setf indices (%array-indices index dimensions)
 694             n1 (apply #'aref array1 indices)
 695             n2 (apply #'aref array2 indices))
 696       (unless (funcall test n1 n2)
 697         (push (list indices n1 n2 (funcall errfun n1 n2))
 698               errseq)))))
 699
 700 (defun array-error (array1 array2 &key
 701                     (test #'number-equal)
 702                     (error-function #'relative-error))
 703   "Return a list of the indices and error between the array elements."
 704   (cond
 705     ((not (arrayp array1))
 706      (error "ARRAY1 is not an array."))
 707     ((not (arrayp array2))
 708      (error "ARRAY2 is not an array."))
 709     ((not (equal (array-dimensions array1) (array-dimensions array2)))
 710      (error "Arrays are not equal dimensions."))
 711     (t (%array-error array1 array2 test error-function))))
 712
 713 ;;; Floating point data functions
 714 (defun make-2d-list (rows columns &key (initial-element 0))
 715   "Return a nested list with INITIAL-ELEMENT."
 716   (mapcar (lambda (x) (make-list columns :initial-element x))
 717           (make-list rows :initial-element initial-element)))
 718
 719 (defun %complex-float-random (limit &optional (state *random-state*))
 720   "Return a random complex float number."
 721   (complex
 722    (random (realpart limit) state)
 723    (random (imagpart limit) state)))
 724
 725 (defun %complex-rational-random (limit &optional (state *random-state*))
 726   "Return a random complex rational number."
 727   (let ((imaglimit (imagpart limit)))
 728     (if (< 1 imaglimit)
 729         (complex
 730          (random (realpart limit) state)
 731          ;; Ensure that the imaginary part is not zero.
 732          (do ((im (random imaglimit state)
 733                   (random imaglimit state)))
 734              ((< 0 im) im)))
 735         (error "Imaginary part must be greater than 1."))))
 736
 737 (defun complex-random (limit &optional (state *random-state*))
 738   "Return a random complex number. "
 739   (check-type limit complex)
 740   (if (typep limit '(complex rational))
 741       (%complex-rational-random limit state)
 742       (%complex-float-random limit state)))
 743
 744 (defun make-random-list (size &optional (limit 1.0))
 745   "Return a list of random numbers."
 746   (mapcar (if (complexp limit) #'complex-random #'random)
 747           (make-list size :initial-element limit)))
 748
 749 (defun make-random-2d-list (rows columns &optional (limit 1.0))
 750   "Return a nested list of random numbers."
 751   (mapcar (lambda (x) (make-random-list columns x))
 752           (make-list rows :initial-element limit)))
 753
 754 (defun make-random-2d-array (rows columns &optional (limit 1.0))
 755   "Return a 2D array of random numbers."
 756   (let ((new-array (make-array (list rows columns)
 757                                :element-type (type-of limit)))
 758         (random-func (if (complexp limit)
 759                          #'complex-random
 760                          #'random)))
 761     (dotimes (i0 rows new-array)
 762       (dotimes (i1 columns)
 763         (setf (aref new-array i0 i1)
 764               (funcall random-func limit))))))