floating-point.lisp

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