external/oct/qd-test.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
  27 (in-package #:qdi)
  28
  29 ;; Compute to how many bits EST and TRUE are equal.  If they are
  30 ;; identical, return T.
  31 (defun bit-accuracy (est true)
  32   (let* ((diff (abs-qd (sub-qd est true)))
  33          (err (if (zerop-qd true)
  34                   (qd-0 diff)
  35                   (cl:/ (qd-0 diff) (abs (qd-0 true))))))
  36     (if (zerop (qd-0 diff))
  37         t
  38         (cl:- (log err 2d0)))))
  39
  40 (defun print-result (est true)
  41   (format t "est: ~/qdi::qd-format/~%" est)
  42   (format t "tru: ~/qdi::qd-format/~%" true)
  43   (format t "err: ~A~%" (qd-0 (sub-qd est true)))
  44   (format t "bits: ~,1f~%" (bit-accuracy est true)))
  45
  46 ;; Machin's formula for pi
  47 #+nil
  48 (defun atan-series (x)
  49   (let* ((d 1d0)
  50          (eps (make-qd-d (scale-float 1d0 -212)
  51                         0d0 0d0 0d0))
  52          (tmp x)
  53          (r (sqr-qd tmp))
  54          (s1 (make-qd-dd 0w0 0w0))
  55          (k 0)
  56          (sign 1))
  57     (loop while (qd-> tmp eps) do
  58           (incf k)
  59           (setf s1
  60                 (if (minusp sign)
  61                     (sub-qd s1 (div-qd tmp (make-qd-d d 0d0 0d0 0d0)))
  62                     (add-qd s1 (div-qd tmp (make-qd-d d 0d0 0d0 0d0)))))
  63           (incf d 2d0)
  64           (setf tmp (mul-qd tmp r))
  65           (setf sign (cl:- sign)))
  66     s1))
  67
  68 ;; pi =
  69 ;; 3.1415926535897932384626433832795028841971693993751058209749445923078L0
  70 (defun test2 ()
  71   ;; pi/4 = 4 * arctan(1/5) - arctan(1/239)
  72   ;;
  73   ;; Arctan is computed using the Taylor series
  74   ;;
  75   ;;   arctan(x) = x - x^3/3 + x^5/5 - x^7/7
  76   (flet ((atan-series (x)
  77            (let* ((d 1d0)
  78                   (eps (make-qd-d (scale-float 1d0 -212)))
  79                   (tmp x)
  80                   (r (sqr-qd tmp))
  81                   (s1 (make-qd-d 0d0))
  82                   (k 0)
  83                   (sign 1))
  84              (loop while (qd-> tmp eps) do
  85                    (incf k)
  86                    (setf s1
  87                          (if (minusp sign)
  88                              (sub-qd s1 (div-qd tmp (make-qd-d d)))
  89                              (add-qd s1 (div-qd tmp (make-qd-d d)))))
  90                    (incf d 2d0)
  91                    (setf tmp (mul-qd tmp r))
  92                    (setf sign (cl:- sign)))
  93              s1)))
  94     (let* ((x1 (div-qd +qd-one+
  95                        (make-qd-d 5d0)))
  96            (s1 (atan-series x1))
  97            (x2 (div-qd +qd-one+
  98                        (make-qd-d 239d0)))
  99            (s2 (atan-series x2))
 100            (p (mul-qd-d (sub-qd (mul-qd-d s1 4d0)
 101                                 s2)
 102                         4d0)))
 103       (format t "~2&pi via Machin's atan formula~%")
 104       (print-result p +qd-pi+)
 105       p)))
 106
 107 (defun test3 ()
 108   (declare (optimize (speed 3)))
 109   ;; Salamin-Brent Quadratic formula for pi
 110   (let* ((a +qd-one+)
 111          (b (sqrt-qd (make-qd-d 0.5d0)))
 112          (s (make-qd-d 0.5d0))
 113          (m 1d0)
 114          (p (div-qd (mul-qd-d (sqr-qd a) 2d0)
 115                     s)))
 116     (declare (double-float m))
 117     (dotimes (k 9)
 118       (setf m (cl:* 2 m))
 119       (let* ((a-new (mul-qd-d (add-qd a b) .5d0))
 120              (b-new (sqrt-qd (mul-qd a b)))
 121              (s-new (sub-qd s
 122                             (mul-qd-d (sub-qd (sqr-qd a-new)
 123                                               (sqr-qd b-new))
 124                                       m))))
 125         (setf a a-new)
 126         (setf b b-new)
 127         (setf s s-new)
 128         (setf p (div-qd (mul-qd-d (sqr-qd a) 2d0)
 129                         s))))
 130     (format t "~2&Salamin-Brent Quadratic formula for pi~%")
 131     (print-result p +qd-pi+)
 132     p))
 133
 134 (defun test4 ()
 135   (declare (optimize (speed 3)))
 136   ;; Borwein Quartic formula for pi
 137   (let* ((a (sub-qd (make-qd-d 6d0)
 138                     (mul-qd-d (sqrt-qd (make-qd-d 2d0))
 139                               4d0)))
 140          (y (sub-qd-d (sqrt-qd (make-qd-d 2d0))
 141                       1d0))
 142          (m 2d0)
 143          (p (div-qd +qd-one+
 144                     a)))
 145     (declare (double-float m))
 146     (dotimes (k 9)
 147       (setf m (cl:* 4 m))
 148       (let ((r (nroot-qd (sub-qd +qd-one+
 149                                  (sqr-qd (sqr-qd y)))
 150                          4)))
 151         (setf y (div-qd (sub-d-qd 1d0
 152                                   r)
 153                         (add-d-qd 1d0
 154                                 r)))
 155         (setf a (sub-qd (mul-qd a
 156                                 (sqr-qd (sqr-qd (add-qd-d y 1d0))))
 157                         (mul-qd-d (mul-qd y
 158                                           (add-qd-d (add-qd y (sqr-qd y))
 159                                                     1d0))
 160                                   m)))
 161         (setf p (div-qd +qd-one+
 162                         a))))
 163     (format t "~2&Borwein's Quartic formula for pi~%")
 164     (print-result p +qd-pi+)
 165     p))
 166
 167 ;; e =
 168 ;; 2.718281828459045235360287471352662497757247093699959574966967627724L0
 169 (defun test5 ()
 170   ;; Taylor series for e
 171   (let ((s (make-qd-d 2d0))
 172         (tmp +qd-one+)
 173         (n 1d0)
 174         (delta 0d0)
 175         (i 0))
 176     (loop while (qd-> tmp (make-qd-d 1d-100)) do
 177           (incf i)
 178           (incf n)
 179           (setf tmp (div-qd tmp
 180                             (make-qd-d (float n 1d0))))
 181           (setf s (add-qd s tmp)))
 182     (format t "~2&e via Taylor series~%")
 183     (print-result s +qd-e+)
 184     s))
 185
 186 ;; log(2) =
 187 ;; 0.6931471805599453094172321214581765680755001343602552541206800094934L0
 188 (defun test6 ()
 189   ;; Taylor series for log 2
 190   ;;
 191   ;; -log(1-x) = x + x^2/2 + x^3/3 + x^4/4 + ...
 192   ;;
 193   ;; with x = 1/2 to get log(1/2) = -log(2)
 194   (let ((s (make-qd-d .5d0))
 195         (tt (make-qd-d .5d0))
 196         (n 1d0)
 197         (i 0))
 198     (loop while (qd-> tt (make-qd-d 1d-100)) do
 199           (incf i)
 200           (incf n)
 201           (setf tt (mul-qd-d tt .5d0))
 202           (setf s (add-qd s
 203                           (div-qd tt (make-qd-d (float n 1d0))))))
 204     (format t "~2&log(2) via Taylor series~%")
 205     (print-result s +qd-log2+)
 206     s))
 207
 208 (defun test-atan (&optional (fun #'atan-qd))
 209   ;; Compute atan for known values
 210
 211   (format t "~2&atan via ~S~%" fun)
 212   ;; atan(1/sqrt(3)) = pi/6
 213   (let* ((arg (div-qd +qd-one+ (sqrt-qd (make-qd-d 3d0))))
 214          (y (div-qd (funcall fun arg) +qd-pi+))
 215          (true (div-qd +qd-one+ (make-qd-d 6d0))))
 216     (format t "atan(1/sqrt(3))/pi = ~/qdi::qd-format/~%" y)
 217     (format t "1/6                = ~/qdi::qd-format/~%" true)
 218     (format t "bits               = ~,1f~%"
 219             (bit-accuracy y true)))
 220   ;; atan(sqrt(3)) = %pi/3
 221   (let* ((arg (sqrt-qd (make-qd-d 3d0)))
 222          (y (div-qd (funcall fun arg) +qd-pi+))
 223          (true (div-qd +qd-one+ (make-qd-d 3d0))))
 224     (format t "atan(sqrt(3))/pi   = ~/qdi::qd-format/~%" y)
 225     (format t "1/3                = ~/qdi::qd-format/~%" true)
 226     (format t "bits               = ~,1f~%"
 227             (bit-accuracy y true)))
 228   ;; atan(1) = %pi/4
 229   (let* ((arg +qd-one+)
 230          (y (div-qd (funcall fun arg) +qd-pi+))
 231          (true (div-qd +qd-one+ (make-qd-d 4d0))))
 232     (format t "atan(1)/pi         = ~/qdi::qd-format/~%" y)
 233     (format t "1/4                = ~/qdi::qd-format/~%" true)
 234     (format t "bits               = ~,1f~%"
 235             (bit-accuracy y true))))
 236
 237 (defun test-sin ()
 238   (format t "~2&sin~%")
 239   (let* ((arg (div-qd +qd-pi+ (make-qd-d 6d0)))
 240          (y (sin-qd arg))
 241          (true (make-qd-d 0.5d0)))
 242     (format t "sin(pi/6)      = ~/qdi::qd-format/~%" y)
 243     (format t "1/2            = ~/qdi::qd-format/~%" true)
 244     (format t "bits           = ~,1f~%"
 245             (bit-accuracy y true)))
 246   (let* ((arg (div-qd +qd-pi+ (make-qd-d 4d0)))
 247          (y (sin-qd arg))
 248          (true (sqrt-qd (make-qd-d 0.5d0))))
 249     (format t "sin(pi/4)      = ~/qdi::qd-format/~%" y)
 250     (format t "1/sqrt(2)      = ~/qdi::qd-format/~%" true)
 251     (format t "bits           = ~,1f~%"
 252             (bit-accuracy y true)))
 253   (let* ((arg (div-qd +qd-pi+ (make-qd-d 3d0)))
 254          (y (sin-qd arg))
 255          (true (div-qd (sqrt-qd (make-qd-d 3d0)) (make-qd-d 2d0))))
 256     (format t "sin(pi/3)      = ~/qdi::qd-format/~%" y)
 257     (format t "sqrt(3)/2      = ~/qdi::qd-format/~%" true)
 258     (format t "bits           = ~,1f~%"
 259             (bit-accuracy y true)))
 260   )
 261
 262 (defun test-tan (&optional (f #'tan-qd))
 263   (format t "~2&tan via ~S~%" f)
 264   (let* ((arg (div-qd +qd-pi+ (make-qd-d 6d0)))
 265          (y (funcall f arg))
 266          (true (div-qd +qd-one+ (sqrt-qd (make-qd-d 3d0)))))
 267     (format t "tan(pi/6)      = ~/qdi::qd-format/~%" y)
 268     (format t "1/sqrt(3)      = ~/qdi::qd-format/~%" true)
 269     (format t "bits           = ~,1f~%"
 270             (bit-accuracy y true)))
 271   (let* ((arg (div-qd +qd-pi+ (make-qd-d 4d0)))
 272          (y (funcall f arg))
 273          (true +qd-one+))
 274     (format t "tan(pi/4)      = ~/qdi::qd-format/~%" y)
 275     (format t "1              = ~/qdi::qd-format/~%" true)
 276     (format t "bits           = ~,1f~%"
 277             (bit-accuracy y true)))
 278   (let* ((arg (div-qd +qd-pi+ (make-qd-d 3d0)))
 279          (y (funcall f arg))
 280          (true (sqrt-qd (make-qd-d 3d0))))
 281     (format t "tan(pi/3)      = ~/qdi::qd-format/~%" y)
 282     (format t "sqrt(3)        = ~/qdi::qd-format/~%" true)
 283     (format t "bits           = ~,1f~%"
 284             (bit-accuracy y true)))
 285   )
 286
 287 (defun test-asin ()
 288   (format t "~2&asin~%")
 289   (let* ((arg (make-qd-d 0.5d0))
 290          (y (asin-qd arg))
 291          (true (div-qd +qd-pi+ (make-qd-d 6d0))))
 292     (format t "asin(1/2)      = ~/qdi::qd-format/~%" y)
 293     (format t "pi/6           = ~/qdi::qd-format/~%" true)
 294     (format t "bits           = ~,1f~%"
 295             (bit-accuracy y true)))
 296   (let* ((arg (sqrt-qd (make-qd-d 0.5d0)))
 297          (y (asin-qd arg))
 298          (true (div-qd +qd-pi+ (make-qd-d 4d0))))
 299     (format t "asin(1/sqrt(2))= ~/qdi::qd-format/~%" y)
 300     (format t "pi/4           = ~/qdi::qd-format/~%" true)
 301     (format t "bits           = ~,1f~%"
 302             (bit-accuracy y true)))
 303   (let* ((arg (div-qd (sqrt-qd (make-qd-d 3d0)) (make-qd-d 2d0)))
 304          (y (asin-qd arg))
 305          (true (div-qd +qd-pi+ (make-qd-d 3d0))))
 306     (format t "asin(sqrt(3)/2)= ~/qdi::qd-format/~%" y)
 307     (format t "pi/3           = ~/qdi::qd-format/~%" true)
 308     (format t "bits           = ~,1f~%"
 309             (bit-accuracy y true)))
 310   )
 311
 312 (defun test-log (f)
 313   (format t "~2&Log via ~A~%" f)
 314   (let* ((arg (make-qd-d 2d0))
 315          (y (funcall f arg))
 316          (true +qd-log2+))
 317     (format t "log(2)      = ~/qdi::qd-format/~%" y)
 318     (format t "log(2)      = ~/qdi::qd-format/~%" true)
 319     (format t "bits        = ~,1f~%"
 320             (bit-accuracy y true)))
 321   (let* ((arg (make-qd-d 10d0))
 322          (y (funcall f arg))
 323          (true +qd-log10+))
 324     (format t "log(10)     = ~/qdi::qd-format/~%" y)
 325     (format t "log(10)     = ~/qdi::qd-format/~%" true)
 326     (format t "bits        = ~,1f~%"
 327             (bit-accuracy y true)))
 328   (let* ((arg (add-d-qd 1d0 (scale-float-qd (make-qd-d 1d0) -80)))
 329          (y (funcall f arg))
 330          (true (qd-from-string "8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25")))
 331     (format t "log(1+2^-80) = ~/qdi::qd-format/~%" y)
 332     (format t "log(1+2^-80) = ~/qdi::qd-format/~%" true)
 333     (format t "bits         = ~,1f~%"
 334             (bit-accuracy y true)))
 335   )
 336
 337 (defun test-log1p (f)
 338   (format t "~2&Log1p via ~A~%" f)
 339   (let* ((arg (make-qd-d 1d0))
 340          (y (funcall f arg))
 341          (true +qd-log2+))
 342     (format t "log1p(1)    = ~/qdi::qd-format/~%" y)
 343     (format t "log(2)      = ~/qdi::qd-format/~%" true)
 344     (format t "bits        = ~,1f~%"
 345             (bit-accuracy y true)))
 346   (let* ((arg (make-qd-d 9d0))
 347          (y (funcall f arg))
 348          (true (qd-from-string "2.3025850929940456840179914546843642076011014886287729760333279009675726096773525q0")))
 349     (format t "log1p(9)    = ~/qdi::qd-format/~%" y)
 350     (format t "log(10)     = ~/qdi::qd-format/~%" true)
 351     (format t "bits        = ~,1f~%"
 352             (bit-accuracy y true)))
 353   (let* ((arg (scale-float-qd (make-qd-d 1d0) -80))
 354          (y (funcall f arg))
 355          (true (qd-from-string "8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25")))
 356     (format t "log1p(2^-80) = ~/qdi::qd-format/~%" y)
 357     (format t "log(1+2^-80) = ~/qdi::qd-format/~%" true)
 358     (format t "bits         = ~,1f~%"
 359             (bit-accuracy y true)))
 360   )
 361
 362 (defun test-exp (f)
 363   (format t "~2&Exp via ~A~%" f)
 364   (let* ((arg +qd-zero+)
 365          (y (funcall f arg))
 366          (true +qd-zero+))
 367     (format t "exp(0)-1    = ~/qdi::qd-format/~%" y)
 368     (format t "0           = ~/qdi::qd-format/~%" true)
 369     (format t "bits        = ~,1f~%"
 370             (bit-accuracy y true)))
 371   (let* ((arg +qd-one+)
 372          (y (funcall f arg))
 373          (true (qd-from-string "1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0")))
 374     (format t "exp(1)-1    = ~/qdi::qd-format/~%" y)
 375     (format t "e-1         = ~/qdi::qd-format/~%" true)
 376     (format t "bits        = ~,1f~%"
 377             (bit-accuracy y true)))
 378
 379   (let* ((arg (scale-float-qd +qd-one+ -100))
 380          (y (funcall f arg))
 381          (true (qd-from-string "7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31")))
 382     (format t "exp(2^-80)-1 = ~/qdi::qd-format/~%" y)
 383     (format t "exp(2^-80)-1 = ~/qdi::qd-format/~%" true)
 384     (format t "bits         = ~,1f~%"
 385             (bit-accuracy y true)))
 386
 387   )
 388 (defun all-tests ()
 389   (test2)
 390   (test3)
 391   (test4)
 392   (test5)
 393   (test6)
 394   (test-sin)
 395   (test-asin)
 396   (dolist (f '(atan-qd/newton atan-qd/cordic atan-qd/duplication))
 397     (test-atan f))
 398   (dolist (f '(tan-qd/sincos tan-qd/cordic))
 399     (test-tan f))
 400   (dolist (f '(log-qd/newton log-qd/agm log-qd/agm2 log-qd/agm3 log-qd/halley))
 401     (test-log f))
 402   (dolist (f '(log1p-qd/duplication))
 403     (test-log1p f))
 404   (dolist (f (list #'(lambda (x)
 405                        (sub-qd-d (exp-qd/reduce x) 1d0))
 406                    #'expm1-qd/series
 407                    #'expm1-qd/duplication))
 408     (test-exp f))
 409   )