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