external/oct/rt-tests.lisp

   1 ;;;; -*- Mode: lisp -*-
   2 ;;;;
   3 ;;;; Copyright (c) 2007 Raymond Toy
   4 ;;;;
   5 ;;;; Permission is hereby granted, free of charge, to any person
   6 ;;;; obtaining a copy of this software and associated documentation
   7 ;;;; files (the "Software"), to deal in the Software without
   8 ;;;; restriction, including without limitation the rights to use,
   9 ;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
  10 ;;;; copies of the Software, and to permit persons to whom the
  11 ;;;; Software is furnished to do so, subject to the following
  12 ;;;; conditions:
  13 ;;;;
  14 ;;;; The above copyright notice and this permission notice shall be
  15 ;;;; included in all copies or substantial portions of the Software.
  16 ;;;;
  17 ;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  18 ;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  19 ;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  20 ;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  21 ;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  22 ;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  23 ;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  24 ;;;; OTHER DEALINGS IN THE SOFTWARE.
  25
  26 (in-package #:qd)
  27
  28 ;; Compute how many bits are the same for two numbers EST and TRUE.
  29 ;; Return T if they are identical.
  30 (defun bit-accuracy (est true)
  31   (let* ((diff (abs (- est true)))
  32          (err (float (if (zerop true)
  33                          diff
  34                          (/ diff (abs true)))
  35                      1d0)))
  36     (if (zerop diff)
  37         t
  38         (- (log err 2)))))
  39
  40 (defun check-accuracy (limit est true)
  41   (let ((bits (bit-accuracy est true)))
  42     (if (numberp bits)
  43         (if (< bits limit)
  44             (list bits limit est true)))))
  45
  46 (defvar *null* (make-broadcast-stream))
  47
  48 ;;; Some simple tests from the Yozo Hida's qd package.
  49
  50 ;; Pi via Machin's formula
  51 (rt:deftest oct.pi.machin
  52     (let* ((*standard-output* *null*)
  53            (val (make-instance 'qd-real :value (qdi::test2)))
  54            (true qd:+pi+))
  55       (check-accuracy 213 val true))
  56   nil)
  57
  58 ;; Pi via Salamin-Brent algorithm
  59 (rt:deftest oct.pi.salamin-brent
  60     (let* ((*standard-output* *null*)
  61            (val (make-instance 'qd-real :value (qdi::test3)))
  62            (true qd:+pi+))
  63       (check-accuracy 202 val true))
  64   nil)
  65
  66 ;; Pi via Borweign's Quartic formula
  67 (rt:deftest oct.pi.borweign
  68     (let* ((*standard-output* *null*)
  69            (val (make-instance 'qd-real :value (qdi::test4)))
  70            (true qd:+pi+))
  71       (check-accuracy 211 val true))
  72   nil)
  73
  74 ;; e via Taylor series
  75 (rt:deftest oct.e.taylor
  76     (let* ((*standard-output* *null*)
  77            (val (make-instance 'qd-real :value (qdi::test5)))
  78            (true (make-instance 'qd-real :value qdi::+qd-e+)))
  79       (check-accuracy 212 val true))
  80   nil)
  81
  82 ;; log(2) via Taylor series
  83 (rt:deftest oct.log2.taylor
  84     (let* ((*standard-output* *null*)
  85            (val (make-instance 'qd-real :value (qdi::test6)))
  86            (true (make-instance 'qd-real :value qdi::+qd-log2+)))
  87       (check-accuracy 212 val true))
  88   nil)
  89
  90 ;;; Tests of atan where we know the analytical result
  91 (rt:deftest oct.atan.1
  92     (let* ((arg (/ (sqrt #q3)))
  93            (y (/ (atan arg) +pi+))
  94            (true (/ #q6)))
  95       (check-accuracy 212 y true))
  96   nil)
  97
  98 (rt:deftest oct.atan.2
  99     (let* ((arg (sqrt #q3))
 100          (y (/ (atan arg) +pi+))
 101          (true (/ #q3)))
 102       (check-accuracy 212 y true))
 103   nil)
 104
 105 (rt:deftest oct.atan.3
 106     (let* ((arg #q1)
 107          (y (/ (atan arg) +pi+))
 108          (true (/ #q4)))
 109     (check-accuracy 212 y true))
 110   nil)
 111
 112 (rt:deftest oct.atan.4
 113     (let* ((arg #q1q100)
 114            (y (/ (atan arg) +pi+))
 115            (true #q.5))
 116       (check-accuracy 212 y true))
 117   nil)
 118
 119 (rt:deftest oct.atan.5
 120     (let* ((arg #q-1q100)
 121          (y (/ (atan arg) +pi+))
 122          (true #q-.5))
 123     (check-accuracy 212 y true))
 124   nil)
 125
 126
 127 (defun atan-qd/duplication (arg)
 128   (make-instance 'qd-real
 129                  :value (qdi::atan-qd/duplication (qd-value arg))))
 130
 131 ;;; Tests of atan where we know the analytical result.  Same tests,
 132 ;;; but using the atan duplication formula.
 133 (rt:deftest oct.atan/dup.1
 134     (let* ((arg (/ (sqrt #q3)))
 135            (y (/ (atan-qd/duplication arg) +pi+))
 136            (true (/ #q6)))
 137       (check-accuracy 212 y true))
 138   nil)
 139
 140 (rt:deftest oct.atan/dup.2
 141     (let* ((arg (sqrt #q3))
 142          (y (/ (atan-qd/duplication arg) +pi+))
 143          (true (/ #q3)))
 144       (check-accuracy 212 y true))
 145   nil)
 146
 147 (rt:deftest oct.atan/dup.3
 148     (let* ((arg #q1)
 149          (y (/ (atan-qd/duplication arg) +pi+))
 150          (true (/ #q4)))
 151     (check-accuracy 212 y true))
 152   nil)
 153
 154 (rt:deftest oct.atan/dup.4
 155     (let* ((arg #q1q100)
 156            (y (/ (atan-qd/duplication arg) +pi+))
 157            (true #q.5))
 158       (check-accuracy 212 y true))
 159   nil)
 160
 161 (rt:deftest oct.atan/dup.5
 162     (let* ((arg #q-1q100)
 163          (y (/ (atan-qd/duplication arg) +pi+))
 164          (true #q-.5))
 165     (check-accuracy 212 y true))
 166   nil)
 167
 168 ;;; Tests of atan where we know the analytical result.  Same tests,
 169 ;;; but using a CORDIC implementation.
 170 (defun atan-qd/cordic (arg)
 171   (make-instance 'qd-real
 172                  :value (qdi::atan-qd/cordic (qd-value arg))))
 173
 174 (rt:deftest oct.atan/cordic.1
 175     (let* ((arg (/ (sqrt #q3)))
 176            (y (/ (atan-qd/cordic arg) +pi+))
 177            (true (/ #q6)))
 178       (check-accuracy 212 y true))
 179   nil)
 180
 181 (rt:deftest oct.atan/cordic.2
 182     (let* ((arg (sqrt #q3))
 183          (y (/ (atan-qd/cordic arg) +pi+))
 184          (true (/ #q3)))
 185       (check-accuracy 212 y true))
 186   nil)
 187
 188 (rt:deftest oct.atan/cordic.3
 189     (let* ((arg #q1)
 190          (y (/ (atan-qd/cordic arg) +pi+))
 191          (true (/ #q4)))
 192     (check-accuracy 212 y true))
 193   nil)
 194
 195 (rt:deftest oct.atan/cordic.4
 196     (let* ((arg #q1q100)
 197            (y (/ (atan-qd/cordic arg) +pi+))
 198            (true #q.5))
 199       (check-accuracy 212 y true))
 200   nil)
 201
 202 (rt:deftest oct.atan/cordic.5
 203     (let* ((arg #q-1q100)
 204          (y (/ (atan-qd/cordic arg) +pi+))
 205          (true #q-.5))
 206     (check-accuracy 212 y true))
 207   nil)
 208
 209
 210 ;;; Tests of sin where we know the analytical result.
 211 (rt:deftest oct.sin.1
 212     (let* ((arg (/ +pi+ 6))
 213          (y (sin arg))
 214          (true #q.5))
 215     (check-accuracy 212 y true))
 216   nil)
 217
 218 (rt:deftest oct.sin.2
 219     (let* ((arg (/ +pi+ 4))
 220          (y (sin arg))
 221          (true (sqrt #q.5)))
 222     (check-accuracy 212 y true))
 223   nil)
 224
 225 (rt:deftest oct.sin.3
 226     (let* ((arg (/ +pi+ 3))
 227          (y (sin arg))
 228          (true (/ (sqrt #q3) 2)))
 229     (check-accuracy 212 y true))
 230   nil)
 231
 232 ;;; Tests of tan where we know the analytical result.
 233 (rt:deftest oct.tan.1
 234     (let* ((arg (/ +pi+ 6))
 235          (y (tan arg))
 236          (true (/ (sqrt #q3))))
 237     (check-accuracy 212 y true))
 238   nil)
 239
 240 (rt:deftest oct.tan.2
 241     (let* ((arg (/ +pi+ 4))
 242          (y (tan arg))
 243          (true #q1))
 244     (check-accuracy 212 y true))
 245   nil)
 246
 247 (rt:deftest oct.tan.3
 248   (let* ((arg (/ +pi+ 3))
 249          (y (tan arg))
 250          (true (sqrt #q3)))
 251     (check-accuracy 212 y true))
 252   nil)
 253
 254 ;;; Tests of tan where we know the analytical result.  Uses CORDIC
 255 ;;; algorithm.
 256
 257 (defun tan/cordic (arg)
 258   (make-instance 'qd-real
 259                  :value (qdi::tan-qd/cordic (qd-value arg))))
 260
 261 (rt:deftest oct.tan/cordic.1
 262     (let* ((arg (/ +pi+ 6))
 263          (y (tan/cordic arg))
 264          (true (/ (sqrt #q3))))
 265     (check-accuracy 211 y true))
 266   nil)
 267
 268 (rt:deftest oct.tan/cordic.2
 269     (let* ((arg (/ +pi+ 4))
 270          (y (tan/cordic arg))
 271          (true #q1))
 272     (check-accuracy 211 y true))
 273   nil)
 274
 275 (rt:deftest oct.tan/cordic.3
 276   (let* ((arg (/ +pi+ 3))
 277          (y (tan/cordic arg))
 278          (true (sqrt #q3)))
 279     (check-accuracy 210 y true))
 280   nil)
 281
 282 ;;; Tests of asin where we know the analytical result.
 283
 284 (rt:deftest oct.asin.1
 285     (let* ((arg #q.5)
 286          (y (asin arg))
 287          (true (/ +pi+ 6)))
 288     (check-accuracy 212 y true))
 289   nil)
 290
 291 (rt:deftest oct.asin.2
 292     (let* ((arg (sqrt #q.5))
 293            (y (asin arg))
 294            (true (/ +pi+ 4)))
 295       (check-accuracy 212 y true))
 296   nil)
 297
 298 (rt:deftest oct.asin.3
 299     (let* ((arg (/ (sqrt #q3) 2))
 300            (y (asin arg))
 301            (true (/ +pi+ 3)))
 302       (check-accuracy 212 y true))
 303   nil)
 304
 305 ;;; Tests of log.
 306
 307 (rt:deftest oct.log.1
 308     (let* ((arg #q2)
 309            (y (log arg))
 310            (true (make-instance 'qd-real :value qdi::+qd-log2+)))
 311       (check-accuracy 212 y true))
 312   nil)
 313
 314 (rt:deftest oct.log.2
 315     (let* ((arg #q10)
 316          (y (log arg))
 317          (true (make-instance 'qd-real :value qdi::+qd-log10+)))
 318     (check-accuracy 207 y true))
 319   nil)
 320
 321 (rt:deftest oct.log.3
 322     (let* ((arg (+ 1 (scale-float #q1 -80)))
 323            (y (log arg))
 324            (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25))
 325       (check-accuracy 212 y true))
 326   nil)
 327
 328 ;;; Tests of log using Newton iteration.
 329
 330 (defun log/newton (arg)
 331   (make-instance 'qd-real
 332                  :value (qdi::log-qd/newton (qd-value arg))))
 333
 334 (rt:deftest oct.log/newton.1
 335     (let* ((arg #q2)
 336            (y (log/newton arg))
 337            (true (make-instance 'qd-real :value qdi::+qd-log2+)))
 338       (check-accuracy 212 y true))
 339   nil)
 340
 341 (rt:deftest oct.log/newton.2
 342     (let* ((arg #q10)
 343          (y (log/newton arg))
 344          (true (make-instance 'qd-real :value qdi::+qd-log10+)))
 345     (check-accuracy 207 y true))
 346   nil)
 347
 348 (rt:deftest oct.log/newton.3
 349     (let* ((arg (+ 1 (scale-float #q1 -80)))
 350            (y (log/newton arg))
 351            (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25))
 352       (check-accuracy 212 y true))
 353   nil)
 354
 355 ;;; Tests of log using AGM.
 356
 357 (defun log/agm (arg)
 358   (make-instance 'qd-real
 359                  :value (qdi::log-qd/agm (qd-value arg))))
 360
 361 (rt:deftest oct.log/agm.1
 362     (let* ((arg #q2)
 363            (y (log/agm arg))
 364            (true (make-instance 'qd-real :value qdi::+qd-log2+)))
 365       (check-accuracy 203 y true))
 366   nil)
 367
 368 (rt:deftest oct.log/agm.2
 369     (let* ((arg #q10)
 370          (y (log/agm arg))
 371          (true (make-instance 'qd-real :value qdi::+qd-log10+)))
 372     (check-accuracy 205 y true))
 373   nil)
 374
 375 (rt:deftest oct.log/agm.3
 376     (let* ((arg (+ 1 (scale-float #q1 -80)))
 377            (y (log/agm arg))
 378            (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25))
 379       (check-accuracy 123 y true))
 380   nil)
 381
 382 ;;; Tests of log using AGM2, a faster variaton of AGM.
 383
 384 (defun log/agm2 (arg)
 385   (make-instance 'qd-real
 386                  :value (qdi::log-qd/agm2 (qd-value arg))))
 387
 388 (rt:deftest oct.log/agm2.1
 389     (let* ((arg #q2)
 390            (y (log/agm2 arg))
 391            (true (make-instance 'qd-real :value qdi::+qd-log2+)))
 392       (check-accuracy 203 y true))
 393   nil)
 394
 395 (rt:deftest oct.log/agm2.2
 396     (let* ((arg #q10)
 397          (y (log/agm2 arg))
 398          (true (make-instance 'qd-real :value qdi::+qd-log10+)))
 399     (check-accuracy 205 y true))
 400   nil)
 401
 402 (rt:deftest oct.log/agm2.3
 403     (let* ((arg (+ 1 (scale-float #q1 -80)))
 404            (y (log/agm2 arg))
 405            (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25))
 406       (check-accuracy 123 y true))
 407   nil)
 408
 409 ;;; Tests of log using AGM3, a faster variation of AGM2.
 410 (defun log/agm3 (arg)
 411   (make-instance 'qd-real
 412                  :value (qdi::log-qd/agm3 (qd-value arg))))
 413
 414 (rt:deftest oct.log/agm3.1
 415     (let* ((arg #q2)
 416            (y (log/agm3 arg))
 417            (true (make-instance 'qd-real :value qdi::+qd-log2+)))
 418       (check-accuracy 203 y true))
 419   nil)
 420
 421 (rt:deftest oct.log/agm3.2
 422     (let* ((arg #q10)
 423          (y (log/agm3 arg))
 424          (true (make-instance 'qd-real :value qdi::+qd-log10+)))
 425     (check-accuracy 205 y true))
 426   nil)
 427
 428 (rt:deftest oct.log/agm3.3
 429     (let* ((arg (+ 1 (scale-float #q1 -80)))
 430            (y (log/agm3 arg))
 431            (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25))
 432       (check-accuracy 123 y true))
 433   nil)
 434
 435 ;;; Tests of sqrt to make sure we overflow or underflow where we
 436 ;;; shouldn't.
 437
 438 (rt:deftest oct.sqrt.1
 439     (let* ((arg #q1q200)
 440            (y (sqrt arg))
 441            (true #q1q100))
 442       (check-accuracy 212 y true))
 443   nil)
 444
 445 (rt:deftest oct.sqrt.2
 446     (let* ((arg #q1q200)
 447            (y (sqrt arg))
 448            (true #q1q100))
 449       (check-accuracy 212 y true))
 450   nil)
 451
 452 (rt:deftest oct.sqrt.3
 453     (let* ((arg #q1q300)
 454            (y (sqrt arg))
 455            (true #q1q150))
 456       (check-accuracy 212 y true))
 457   nil)
 458
 459 (rt:deftest oct.sqrt.4
 460     (let* ((arg #q1q-200)
 461            (y (sqrt arg))
 462            (true #q1q-100))
 463       (check-accuracy 212 y true))
 464   nil)
 465
 466 (rt:deftest oct.sqrt.5
 467     (let* ((arg #q1q-250)
 468            (y (sqrt arg))
 469            (true #q1q-125))
 470       (check-accuracy 212 y true))
 471   nil)
 472
 473 ;;; Tests of log1p(x) = log(1+x), using the duplication formula.
 474
 475 (defun log1p/dup (arg)
 476   (make-instance 'qd-real
 477                  :value (qdi::log1p-qd/duplication (qd-value arg))))
 478
 479 (rt:deftest oct.log1p.1
 480     (let* ((arg #q9)
 481            (y (log1p/dup arg))
 482            (true #q2.3025850929940456840179914546843642076011014886287729760333279009675726096773525q0))
 483       (check-accuracy 212 y true))
 484   nil)
 485
 486 (rt:deftest oct.log1p.2
 487     (let* ((arg (scale-float #q1 -80))
 488            (y (log1p/dup arg))
 489            (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25))
 490       (check-accuracy 212 y true))
 491   nil)
 492
 493 ;;; Tests of expm1(x) = exp(x) - 1, using a Taylor series with
 494 ;;; argument reduction.
 495
 496 (defun expm1/series (arg)
 497   (make-instance 'qd-real
 498                  :value (qdi::expm1-qd/series (qd-value arg))))
 499
 500 (rt:deftest oct.expm1/series.1
 501   (let* ((arg #q0)
 502          (y (expm1/series arg))
 503          (true #q0))
 504     (check-accuracy 212 y true))
 505   nil)
 506
 507 (rt:deftest oct.expm1/series.2
 508   (let* ((arg #q1)
 509          (y (expm1/series arg))
 510          (true #q1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0))
 511     (check-accuracy 211 y true))
 512   nil)
 513
 514 (rt:deftest oct.expm1/series.3
 515     (let* ((arg (scale-float #q1 -100))
 516            (y (expm1/series arg))
 517            (true #q7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31))
 518       (check-accuracy 211 y true))
 519   nil)
 520
 521 ;;; Tests of expm1(x) = exp(x) - 1, using duplication formula.
 522
 523 (defun expm1/dup (arg)
 524   (make-instance 'qd-real
 525                  :value (qdi::expm1-qd/duplication (qd-value arg))))
 526
 527
 528 (rt:deftest oct.expm1/dup.1
 529   (let* ((arg #q0)
 530          (y (expm1/dup arg))
 531          (true #q0))
 532     (check-accuracy 212 y true))
 533   nil)
 534
 535 (rt:deftest oct.expm1/dup.2
 536   (let* ((arg #q1)
 537          (y (expm1/dup arg))
 538          (true #q1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0))
 539     (check-accuracy 211 y true))
 540   nil)
 541
 542 (rt:deftest oct.expm1/dup.3
 543     (let* ((arg (scale-float #q1 -100))
 544            (y (expm1/dup arg))
 545            (true #q7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31))
 546       (check-accuracy 211 y true))
 547   nil)