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)
 300         (sumsq 1)
 301         (abs-val nil))
 302     (dolist (elm data (values scale sumsq))
 303       (when (plusp (setq abs-val (abs elm)))
 304         (if (< scale abs-val)
 305             (setq
 306              sumsq (1+ (* sumsq (expt (/ scale abs-val) 2)))
 307              scale abs-val)
 308             (setq sumsq (+ sumsq (expt (/ elm scale) 2))))))))
 309
 310 (defmethod sumsq ((data vector))
 311   "Return the scaling parameter and the sum of the squares of the ~
 312    vector."
 313   (let ((scale 0)
 314         (sumsq 1)
 315         (abs-val nil)
 316         (size (length data)))
 317     (dotimes (index size (values scale sumsq))
 318       (when (plusp (setq abs-val (abs (svref data index))))
 319         (if (< scale abs-val)
 320             (setq
 321              sumsq (1+ (* sumsq (expt (/ scale abs-val) 2)))
 322              scale abs-val)
 323             (setq
 324              sumsq
 325              (+ sumsq (expt (/ (svref data index) scale) 2))))))))
 326
 327 (defmethod sumsq ((data array))
 328   "Return the scaling parameter and the sum of the squares of the ~
 329    array."
 330   (sumsq (make-array (array-total-size data)
 331                      :element-type (array-element-type data)
 332                      :displaced-to data)))
 333
 334 ;;; (SUMP data) => scale, sump
 335 (defmethod sump ((data list) (p real))
 336   "Return the scaling parameter and the sum of the powers of p of the ~
 337    data."
 338   (let ((scale 0)
 339         (sump 1)
 340         (abs-val nil))
 341     (dolist (elm data (values scale sump))
 342       (when (plusp (setq abs-val (abs elm)))
 343         (if (< scale abs-val)
 344             (setq
 345              sump (1+ (* sump (expt (/ scale abs-val) p)))
 346              scale abs-val)
 347             (setq sump (+ sump (expt (/ elm scale) p))))))))
 348
 349 (defmethod sump ((data vector) (p real))
 350   "Return the scaling parameter and the sum of the powers of p of the ~
 351    vector."
 352   (let ((scale 0)
 353         (sump 1)
 354         (abs-val nil)
 355         (size (length data)))
 356     (dotimes (index size (values scale sump))
 357       (when (plusp (setq abs-val (abs (svref data index))))
 358         (if (< scale abs-val)
 359             (setq
 360              sump (1+ (* sump (expt (/ scale abs-val) p)))
 361              scale abs-val)
 362             (setq
 363              sump
 364              (+ sump (expt (/ (svref data index) scale) p))))))))
 365
 366 (defmethod sump ((data array) (p real))
 367   "Return the scaling parameter and the sum of the powers of p of the ~
 368    array."
 369   (sump (make-array (array-total-size data)
 370                     :element-type (array-element-type data)
 371                     :displaced-to data)
 372         p))
 373
 374 ;;; (NORM data) => float
 375
 376 (defgeneric %norm (data measure)
 377   (:documentation
 378    "Return the norm of the data according to measure."))
 379
 380 (defmethod %norm ((data list) (measure (eql 1)))
 381   "Return the Taxicab norm of the list."
 382   (loop for item in data sum (abs item)))
 383
 384 (defmethod %norm ((data vector) (measure (eql 1)))
 385   "Return the Taxicab norm of the vector."
 386   (loop for item across data sum (abs item)))
 387
 388 (defmethod %norm ((data list) (measure (eql 2)))
 389   "Return the Euclidean norm of the list."
 390   (multiple-value-bind (scale sumsq)
 391       (sumsq (map-into (make-array (length data)) #'abs data))
 392     (* scale (sqrt sumsq))))
 393
 394 (defmethod %norm ((data vector) (measure (eql 2)))
 395   "Return the Euclidean norm of the vector."
 396   (multiple-value-bind (scale sumsq)
 397       (sumsq (map-into (make-array (length data)) #'abs data))
 398     (* scale (sqrt sumsq))))
 399
 400 (defmethod %norm ((data list) (measure integer))
 401   "Return the Euclidean norm of the list."
 402   (multiple-value-bind (scale sump)
 403       (sump (map-into (make-array (length data)) #'abs data)
 404             measure)
 405     (* scale (expt sump (/ measure)))))
 406
 407 (defmethod %norm ((data vector) (measure integer))
 408   "Return the Euclidean norm of the vector."
 409   (multiple-value-bind (scale sump)
 410       (sump (map-into (make-array (length data)) #'abs data)
 411             measure)
 412     (* scale (expt sump (/ measure)))))
 413
 414 (defmethod %norm ((data list) (measure (eql :infinity)))
 415   "Return the infinity, or maximum, norm of the list."
 416   (loop for item in data maximize (abs item)))
 417
 418 (defmethod %norm ((data vector) (measure (eql :infinity)))
 419   "Return the infinity, or maximum, norm of the vector."
 420   (loop for item across data maximize (abs item)))
 421
 422 (defmethod norm ((data list) &optional (measure *measure*))
 423   "Return the norm of the list according to the measure."
 424   (%norm data measure))
 425
 426 (defmethod norm ((data vector) &optional (measure *measure*))
 427   "Return the norm of the vector according to the measure."
 428   (%norm data measure))
 429
 430 ;;; FIXME : Is the entrywise norm of an array useful or confusing?
 431 (defmethod norm ((data array) &optional (measure *measure*))
 432   "Return the entrywise norm of the array according to the measure."
 433   (%norm
 434    (make-array
 435     (array-total-size data)
 436     :element-type (array-element-type data)
 437     :displaced-to data)
 438    measure))
 439
 440 ;;; (RELATIVE-ERROR-NORM exact approximate measure) => float
 441 (defun %relative-error-norm (exact approximate measure)
 442   "Return the relative error norm of the sequences."
 443   (/ (norm (map-into (make-array (length exact))
 444                      (lambda (x1 x2) (abs (- x1 x2)))
 445                      exact approximate) measure)
 446      (norm exact measure)))
 447
 448 (defmethod relative-error-norm ((exact list) (approximate list)
 449                                 &optional (measure *measure*))
 450   "Return the relative error norm of the lists."
 451   (if (= (length exact) (length approximate))
 452       (%relative-error-norm exact approximate measure)
 453       (error "Lists are not equal in length.")))
 454
 455 (defmethod relative-error-norm ((exact list) (approximate vector)
 456                                 &optional (measure *measure*))
 457   "Return the relative error norm of the list and the vector."
 458   (if (= (length exact) (length approximate))
 459       (%relative-error-norm exact approximate measure)
 460       (error "The list and vector are not equal in length.")))
 461
 462 (defmethod relative-error-norm ((exact vector) (approximate list)
 463                                 &optional (measure *measure*))
 464   "Return the relative error norm of the list and the vector."
 465   (if (= (length exact) (length approximate))
 466       (%relative-error-norm exact approximate measure)
 467       (error "The list and vector are not equal in length.")))
 468
 469 (defmethod relative-error-norm ((exact vector) (approximate vector)
 470                                 &optional (measure *measure*))
 471   "Return the relative error norm of the vectors."
 472   (if (= (length exact) (length approximate))
 473       (%relative-error-norm exact approximate measure)
 474       (error "Vectors are not equal in length.")))
 475
 476 (defmethod relative-error-norm ((exact array) (approximate vector)
 477                                 &optional (measure *measure*))
 478   "Return the relative error norm of the arrays."
 479   (if (equal (array-dimensions exact)
 480              (array-dimensions approximate))
 481       (%relative-error-norm
 482        (make-array (array-total-size exact)
 483                    :element-type (array-element-type exact)
 484                    :displaced-to exact)
 485        (make-array (array-total-size approximate)
 486                    :element-type (array-element-type approximate)
 487                    :displaced-to approximate)
 488        measure)
 489       (error "Arrays are not equal dimensions.")))
 490
 491 ;;; (NORM-EQUAL data1 data2 epsilon measure) => boolean
 492 (defun %norm-equal (seq1 seq2 epsilon measure)
 493   "Return true if the relative error norm is less than epsilon."
 494   (or
 495    (and (null seq1) (null seq2))
 496    (< (%relative-error-norm seq1 seq2 measure) epsilon)))
 497
 498 (defmethod norm-equal ((data1 list) (data2 list) &optional
 499                        (epsilon *epsilon*) (measure *measure*))
 500   "Return true if the lists are equal in length and the relative error
 501 norm is less than epsilon."
 502   (%norm-equal data1 data2 epsilon measure))
 503
 504 (defmethod norm-equal ((data1 list) (data2 vector) &optional
 505                        (epsilon *epsilon*) (measure *measure*))
 506   "Return true if the vector and the list are equal in length and the
 507 relative error norm is less than epsilon."
 508   (%norm-equal data1 data2 epsilon measure))
 509
 510 (defmethod norm-equal ((data1 vector) (data2 list) &optional
 511                        (epsilon *epsilon*) (measure *measure*))
 512   "Return true if the vector and the list are equal in length and the
 513 relative error norm is less than epsilon."
 514   (%norm-equal data1 data2 epsilon measure))
 515
 516 (defmethod norm-equal ((data1 vector) (data2 vector) &optional
 517                        (epsilon *epsilon*) (measure *measure*))
 518   "Return true if the vectors are equal in length and the relative
 519 error norm is less than epsilon."
 520   (%norm-equal data1 data2 epsilon measure))
 521
 522 (defmethod norm-equal ((data1 array) (data2 array) &optional
 523                        (epsilon *epsilon*) (measure *measure*))
 524   "Return true if the arrays are equal in length and the relative
 525 error norm is less than epsilon."
 526   (when (equal (array-dimensions data1)
 527                (array-dimensions data2))
 528     (%norm-equal
 529      (make-array (array-total-size data1)
 530                  :element-type (array-element-type data1)
 531                  :displaced-to data1)
 532      (make-array (array-total-size data2)
 533                  :element-type (array-element-type data2)
 534                  :displaced-to data2)
 535      epsilon measure)))
 536
 537 (defmacro assert-norm-equal (expected form &rest extras)
 538   `(expand-assert :equal ,form ,form ,expected ,extras :test #'norm-equal))
 539
 540 ;;; (NORMALIZE-FLOAT significand &optional exponent) => significand,exponent
 541 ;;; [NumAlgoC] : Definition 1.7, pg. 4
 542 ;;;
 543 ;;; To avoid using 0.1, first 1.0 <= significand < 10. On the final
 544 ;;; return, scale 0.1 <= significand < 1.
 545 (defun %normalize-float (significand &optional (exponent 0))
 546   "Return the normalized floating point number and exponent."
 547   ;;; FIXME : Use the LOOP.
 548   (cond
 549     ((zerop significand)
 550      (values significand 0))
 551     ((>= (abs significand) 10)
 552      (%normalize-float (/ significand 10.0) (1+ exponent)))
 553     ((< (abs significand) 1)
 554      (%normalize-float (* significand 10.0) (1- exponent)))
 555     (t (values (/ significand 10.0) (1+ exponent)))))
 556
 557 ;;; (SIGFIG-EQUAL float1 float2 significant-figures) => true or false
 558 (defun %sigfig-equal (float1 float2 significant-figures)
 559   "Return true if the floating point numbers have equal significant
 560 figures."
 561   (if (or (zerop float1) (zerop float2))
 562       (< (abs (+ float1 float2)) (* 5D-1 (expt 1D1 (- significant-figures))))
 563       (multiple-value-bind (sig1 exp1) (%normalize-float float1)
 564         (multiple-value-bind (sig2 exp2) (%normalize-float float2)
 565           (= (round (* sig1 (expt 1D1 significant-figures)))
 566              (round (* sig2 (expt 1D1 (- significant-figures (- exp1 exp2))))))))))
 567
 568 (defmethod sigfig-equal ((data1 float) (data2 float) &optional
 569                          (significant-figures *significant-figures*))
 570   "Return true if the floating point numbers have equal significant
 571 figures."
 572   (%sigfig-equal data1 data2 significant-figures))
 573
 574 (defun %seq-sigfig-equal (seq1 seq2 significant-figures)
 575   "Return true if the element-wise comparison is equal to the
 576 specified significant figures."
 577   (or
 578    (and (null seq1) (null seq2))
 579    (when (= (length seq1) (length seq2))
 580      (every
 581       (lambda (d1 d2) (sigfig-equal d1 d2 significant-figures))
 582       seq1 seq2))))
 583
 584 (defmethod sigfig-equal ((data1 list) (data2 list) &optional
 585                          (significant-figures *significant-figures*))
 586   "Return true if the lists are equal in length and the element-wise
 587 comparison is equal to significant figures."
 588   (%seq-sigfig-equal data1 data2 significant-figures))
 589
 590 (defmethod sigfig-equal ((data1 vector) (data2 list) &optional
 591                          (significant-figures *significant-figures*))
 592   "Return true if the vector and the list are equal in length and the
 593 element-wise comparison is equal to significant figures."
 594   (%seq-sigfig-equal data1 data2 significant-figures))
 595
 596 (defmethod sigfig-equal ((data1 list) (data2 vector) &optional
 597                          (significant-figures *significant-figures*))
 598   "Return true if the list and the vector are equal in length and the
 599 element-wise comparison is equal to significant figures."
 600   (%seq-sigfig-equal data1 data2 significant-figures))
 601
 602 (defmethod sigfig-equal ((data1 vector) (data2 vector) &optional
 603                          (significant-figures *significant-figures*))
 604   "Return true if the vectors are equal in length and the element-wise
 605 comparison is equal to significant figures."
 606   (%seq-sigfig-equal data1 data2 significant-figures))
 607
 608 (defmethod sigfig-equal ((data1 array) (data2 array) &optional
 609                          (significant-figures *significant-figures*))
 610   "Return true if the arrays are equal in length and the element-wise
 611 comparison is equal to significant figures."
 612   (when (equal (array-dimensions data1)
 613                (array-dimensions data2))
 614     (%seq-sigfig-equal
 615      (make-array (array-total-size data1)
 616                  :element-type (array-element-type data1)
 617                  :displaced-to data1)
 618      (make-array (array-total-size data2)
 619                  :element-type (array-element-type data2)
 620                  :displaced-to data2)
 621      significant-figures)))
 622
 623 (defmacro assert-sigfig-equal (expected form &rest extras)
 624   `(expand-assert :equal ,form ,form ,expected ,extras :test #'sigfig-equal))
 625
 626 ;;; (NUMBER-EQUAL number1 number2) => true or false
 627 (defun number-equal (number1 number2 &optional (epsilon *epsilon*) type-eq-p)
 628   "Return true if the numbers are equal within some epsilon,
 629 optionally requiring the types to be identical."
 630   (and
 631    (or (not type-eq-p) (eq (type-of number1) (type-of number2)))
 632    (float-equal (coerce number1 '(complex double-float))
 633                 (coerce number2 '(complex double-float))
 634                 epsilon)))
 635
 636 (defmacro assert-number-equal (expected form &rest extras)
 637   `(expand-assert :equal ,form ,form ,expected ,extras :test #'number-equal))
 638
 639 ;;; (NUMERICAL-EQUAL result1 result2) => true or false
 640 ;;;
 641 ;;; This is a universal wrapper created by Liam Healy. It is
 642 ;;; implemented to support testing in GSLL. The interface is expanded,
 643 ;;; but backwards compatible with previous versions.
 644 ;;;
 645 (defmethod numerical-equal ((result1 number) (result2 number)
 646                             &key (test #'number-equal))
 647   "Return true if the the numbers are equal according to :TEST."
 648   (funcall test result1 result2))
 649
 650 (defun %sequence-equal (seq1 seq2 test)
 651   "Return true if the sequences are equal in length and each element
 652 is equal according to :TEST."
 653   (when (= (length seq1) (length seq2))
 654     (every (lambda (s1 s2) (numerical-equal s1 s2 :test test))
 655            seq1 seq2)))
 656
 657 (defmethod numerical-equal ((result1 list) (result2 list)
 658                             &key (test #'number-equal))
 659   "Return true if the lists are equal in length and each element is
 660 equal according to :TEST."
 661   (%sequence-equal result1 result2 test))
 662
 663 (defmethod numerical-equal ((result1 vector) (result2 vector)
 664                             &key (test #'number-equal))
 665   "Return true if the vectors are equal in length and each element is
 666 equal according to :TEST."
 667   (%sequence-equal result1 result2 test))
 668
 669 (defmethod numerical-equal ((result1 list) (result2 vector)
 670                             &key (test #'number-equal))
 671   "Return true if every element of the list is equal to the
 672 corresponding element of the vector."
 673   (%sequence-equal result1 result2 test))
 674
 675 (defmethod numerical-equal ((result1 vector) (result2 list)
 676                             &key (test #'number-equal))
 677   "Return true if every element of the list is equla to the
 678 corresponding element of the vector."
 679   (%sequence-equal result1 result2 test))
 680
 681 (defmethod numerical-equal ((result1 array) (result2 array)
 682                             &key (test #'number-equal))
 683   "Return true if the arrays are equal in dimension and each element
 684 is equal according to :TEST."
 685   (when (equal (array-dimensions result1) (array-dimensions result2))
 686     (every test
 687            (make-array (array-total-size result1)
 688                        :element-type (array-element-type result1)
 689                        :displaced-to result1)
 690            (make-array (array-total-size result2)
 691                        :element-type (array-element-type result2)
 692                        :displaced-to result2))))
 693
 694 (defmacro assert-numerical-equal (expected form &rest extras)
 695   `(expand-assert :equal ,form ,form ,expected ,extras :test #'numerical-equal))
 696
 697 ;;; FIXME : Audit and move the diagnostic functions to a separate
 698 ;;; file.
 699
 700 ;;; Diagnostic functions
 701 ;;; Failing a unit test is only half the problem.
 702
 703 ;;; Sequence errors and the indices.
 704 (defun %sequence-error (sequence1 sequence2 test error-function)
 705   "Return a sequence of the indice and error between the sequences."
 706   (let ((n1 nil) (n2 nil)
 707         (errseq '()))
 708     (dotimes (index (length sequence1) errseq)
 709       (setf n1 (elt sequence1 index)
 710             n2 (elt sequence2 index))
 711       (unless (funcall test n1 n2)
 712         (push (list (1- index) n1 n2 (funcall error-function n1 n2))
 713               errseq)))))
 714
 715 (defun sequence-error (sequence1 sequence2 &key
 716                        (test #'number-equal)
 717                        (error-function #'relative-error))
 718   "Return a sequence of the indice and error between the sequence elements."
 719   (cond
 720     ((not (typep sequence1 'sequence))
 721      (error "SEQUENCE1 is not a valid sequence."))
 722     ((not (typep sequence2 'sequence))
 723      (error "SEQUENCE2 is not a valid sequence."))
 724     ((not (= (length sequence1) (length sequence2)))
 725      (error "Lengths not equal. SEQUENCE1(~D) /= SEQUENCE2(~D)."
 726             (length sequence1) (length sequence2)))
 727     (t (%sequence-error sequence1 sequence2 test error-function))))
 728
 729 ;;; Array errors and the indices.
 730 (defun %array-indices (row-major-index dimensions)
 731   "Recursively calculate the indices from the row major index."
 732   (let ((remaining (rest dimensions)))
 733     (if remaining
 734         (multiple-value-bind (index remainder)
 735             (floor row-major-index (reduce #'* remaining))
 736           (cons index (%array-indices remainder remaining)))
 737         (cons row-major-index nil))))
 738
 739 (defun %array-error (array1 array2 test errfun)
 740   "Return a list of the indices, values and error of the elements that
 741 are not equal."
 742   (let ((dimensions (array-dimensions array1))
 743         (n1 nil) (n2 nil)
 744         (indices '())
 745         (errseq '()))
 746     (dotimes (index (array-total-size array1) errseq)
 747       (setf indices (%array-indices index dimensions)
 748             n1 (apply #'aref array1 indices)
 749             n2 (apply #'aref array2 indices))
 750       (unless (funcall test n1 n2)
 751         (push (list indices n1 n2 (funcall errfun n1 n2))
 752               errseq)))))
 753
 754 (defun array-error (array1 array2 &key
 755                     (test #'number-equal)
 756                     (error-function #'relative-error))
 757   "Return a list of the indices and error between the array elements."
 758   (cond
 759     ((not (arrayp array1))
 760      (error "ARRAY1 is not an array."))
 761     ((not (arrayp array2))
 762      (error "ARRAY2 is not an array."))
 763     ((not (equal (array-dimensions array1) (array-dimensions array2)))
 764      (error "Arrays are not equal dimensions."))
 765     (t (%array-error array1 array2 test error-function))))
 766
 767 ;;; Floating point data functions
 768 (defun make-2d-list (rows columns &key (initial-element 0))
 769   "Return a nested list with INITIAL-ELEMENT."
 770   (mapcar (lambda (x) (make-list columns :initial-element x))
 771           (make-list rows :initial-element initial-element)))
 772
 773 (defun %complex-float-random (limit &optional (state *random-state*))
 774   "Return a random complex float number."
 775   (complex
 776    (random (realpart limit) state)
 777    (random (imagpart limit) state)))
 778
 779 (defun %complex-rational-random (limit &optional (state *random-state*))
 780   "Return a random complex rational number."
 781   (let ((imaglimit (imagpart limit)))
 782     (if (< 1 imaglimit)
 783         (complex
 784          (random (realpart limit) state)
 785          ;; Ensure that the imaginary part is not zero.
 786          (do ((im (random imaglimit state)
 787                   (random imaglimit state)))
 788              ((< 0 im) im)))
 789         (error "Imaginary part must be greater than 1."))))
 790
 791 (defun complex-random (limit &optional (state *random-state*))
 792   "Return a random complex number. "
 793   (check-type limit complex)
 794   (if (typep limit '(complex rational))
 795       (%complex-rational-random limit state)
 796       (%complex-float-random limit state)))
 797
 798 (defun make-random-list (size &optional (limit 1.0))
 799   "Return a list of random numbers."
 800   (mapcar (if (complexp limit) #'complex-random #'random)
 801           (make-list size :initial-element limit)))
 802
 803 (defun make-random-2d-list (rows columns &optional (limit 1.0))
 804   "Return a nested list of random numbers."
 805   (mapcar (lambda (x) (make-random-list columns x))
 806           (make-list rows :initial-element limit)))
 807
 808 (defun make-random-2d-array (rows columns &optional (limit 1.0))
 809   "Return a 2D array of random numbers."
 810   (let ((new-array (make-array (list rows columns)
 811                                :element-type (type-of limit)))
 812         (random-func (if (complexp limit)
 813                          #'complex-random
 814                          #'random)))
 815     (dotimes (i0 rows new-array)
 816       (dotimes (i1 columns)
 817         (setf (aref new-array i0 i1)
 818               (funcall random-func limit))))))