external/oct/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
  27 (in-package #:qd)
  28
  29 (defun bit-accuracy (est true)
  30   (let* ((diff (abs (- est true)))
  31          (err (float (if (zerop true)
  32                          diff
  33                          (/ diff (abs true)))
  34                      1d0)))
  35     (if (zerop diff)
  36         t
  37         (- (log err 2)))))
  38
  39 (defun print-result (est true)
  40   (format t "est: ~A~%" est)
  41   (format t "tru: ~A~%" true)
  42   (format t "err: ~A~%" (float (- est true) 1d0))
  43   (format t "bits: ~,1f~%" (bit-accuracy est true)))
  44
  45 (defconstant +e+
  46   (make-instance 'qd-real :value qdi::+qd-e+))
  47
  48 (defconstant +log2+
  49   (make-instance 'qd-real :value qdi::+qd-log2+))
  50
  51 (defun test2 ()
  52   ;; pi/4 = 4 * arctan(1/5) - arctan(1/239)
  53   ;;
  54   ;; Arctan is computed using the Taylor series
  55   ;;
  56   ;;   arctan(x) = x - x^3/3 + x^5/5 - x^7/7
  57   (flet ((atan-series (x)
  58            (let* ((d 1d0)
  59                   (eps (float (scale-float 1d0 -212) #q1))
  60                   (tmp x)
  61                   (r (* tmp tmp))
  62                   (s1 #q0)
  63                   (k 0)
  64                   (sign 1))
  65              (loop while (> tmp eps) do
  66                    (incf k)
  67                    (setf s1
  68                          (if (minusp sign)
  69                              (- s1 (/ tmp d))
  70                              (+ s1 (/ tmp d))))
  71                    (incf d 2d0)
  72                    (setf tmp (* tmp r))
  73                    (setf sign (- sign)))
  74              s1)))
  75     (let* ((x1 (/ #q1 5))
  76            (s1 (atan-series x1))
  77            (x2 (/ #q1 239))
  78            (s2 (atan-series x2))
  79            (p (* (- (* s1 4)
  80                     s2)
  81                  4)))
  82       (format t "~2&pi via Machin's atan formula~%")
  83       (print-result p +pi+)
  84       p)))
  85
  86 (defun test3 ()
  87   (declare (optimize (speed 3)))
  88   ;; Salamin-Brent Quadratic formula for pi
  89   (let* ((a #q1)
  90          (b (sqrt #q.5))
  91          (s #q.5)
  92          (m 1d0)
  93          (p (/ (* (* a a)
  94                   2d0)
  95                s)))
  96     (declare (double-float m))
  97     (dotimes (k 9)
  98       (setf m (* 2 m))
  99       (let* ((a-new (* (+ a b) .5d0))
 100              (b-new (sqrt (* a b)))
 101              (s-new (- s
 102                        (* (- (* a-new a-new)
 103                              (* b-new b-new))
 104                           m))))
 105         (setf a a-new)
 106         (setf b b-new)
 107         (setf s s-new)
 108         (setf p (/ (* (* a a) 2d0)
 109                         s))))
 110     (format t "~2&Salamin-Brent Quadratic formula for pi~%")
 111     (print-result p +pi+)
 112     p))
 113
 114 (defun test4 ()
 115   (declare (optimize (speed 3)))
 116   ;; Borwein Quartic formula for pi
 117   (let* ((a (- 6
 118                (* (sqrt #q2)
 119                   4)))
 120          (y (- (sqrt #q2)
 121                1))
 122          (m 2d0)
 123          (p (/ a)))
 124     (declare (double-float m))
 125     (dotimes (k 9)
 126       (setf m (* 4 m))
 127       (let ((r (expt (- 1 (expt y 4))
 128                      1/4)))
 129         (setf y (/ (- 1d0 r)
 130                    (+ 1d0 r)))
 131         (setf a (- (* a
 132                       (expt (+ y 1d0) 4))
 133                    (* (* y
 134                          (+ (+ y (expt y 2))
 135                             1d0))
 136                       m)))
 137         (setf p (/ a))))
 138     (format t "~2&Borwein's Quartic formula for pi~%")
 139     (print-result p +pi+)
 140     p))
 141
 142 (defun test5 ()
 143   ;; Taylor series for e
 144   (let ((s #q2)
 145         (tmp #q1)
 146         (n 1d0)
 147         (delta 0d0)
 148         (i 0))
 149     (loop while (> tmp 1d-100) do
 150           (incf i)
 151           (incf n)
 152           (setf tmp (/ tmp n))
 153           (setf s (+ s tmp)))
 154     (format t "~2&e via Taylor series~%")
 155     (print-result s +e+)
 156     s))
 157
 158 (defun test6 ()
 159   ;; Taylor series for log 2
 160   ;;
 161   ;; -log(1-x) = x + x^2/2 + x^3/3 + x^4/4 + ...
 162   ;;
 163   ;; with x = 1/2 to get log(1/2) = -log(2)
 164   (let ((s #q.5)
 165         (tt #q.5)
 166         (n 1d0)
 167         (i 0))
 168     (loop while (> tt 1d-100) do
 169           (incf i)
 170           (incf n)
 171           (setf tt (* tt .5d0))
 172           (setf s (+ s
 173                      (/ tt n))))
 174     (format t "~2&log(2) via Taylor series~%")
 175     (print-result s +log2+)
 176     s))
 177
 178 (defun test-atan ()
 179   (let* ((arg (/ (sqrt #q3)))
 180          (y (/ (atan arg) +pi+))
 181          (true (/ #q6)))
 182     (format t "~2&atan for special args~%")
 183     (format t "atan(1/sqrt(3))/pi = 1/6~%")
 184     (print-result y true))
 185   ;; atan(sqrt(3)) = %pi/3
 186   (let* ((arg (sqrt #q3))
 187          (y (/ (atan arg) +pi+))
 188          (true (/ #q3)))
 189     (format t "atan(sqrt(3))/pi = 1/3~%")
 190     (print-result y true))
 191   ;; atan(1) = %pi/4
 192   (let* ((arg #q1)
 193          (y (/ (atan arg) +pi+))
 194          (true (/ #q4)))
 195     (format t "atan(1)/pi = 1/4~%")
 196     (print-result y true))
 197   (let* ((arg #q1q100)
 198          (y (/ (atan arg) +pi+))
 199          (true #q.5))
 200     (format t "atan(1q100)/pi = 1/2~%")
 201     (print-result y true))
 202   (let* ((arg #q-1q100)
 203          (y (/ (atan arg) +pi+))
 204          (true #q-.5))
 205     (format t "atan(-1q100)/pi = -1/2~%")
 206     (print-result y true)))
 207
 208 (defun test-sin ()
 209   (format t "~2&sin for special args~%")
 210   (let* ((arg (/ +pi+ 6))
 211          (y (sin arg))
 212          (true #q.5))
 213     (format t "sin(pi/6) = 1/2~%")
 214     (print-result y true))
 215   (let* ((arg (/ +pi+ 4))
 216          (y (sin arg))
 217          (true (sqrt #q.5)))
 218     (format t "sin(pi/4) = 1/sqrt(2)~%")
 219     (print-result y true))
 220   (let* ((arg (/ +pi+ 3))
 221          (y (sin arg))
 222          (true (/ (sqrt #q3) 2)))
 223     (format t "sin(pi/3) = sqrt(3)/2~%")
 224     (print-result y true)))
 225
 226 (defun test-tan ()
 227   (format t "~2&tan for special args~%")
 228   (let* ((arg (/ +pi+ 6))
 229          (y (tan arg))
 230          (true (/ (sqrt #q3))))
 231     (format t"tan(pi/6) = 1/sqrt(3)~%")
 232     (print-result y true))
 233   (let* ((arg (/ +pi+ 4))
 234          (y (tan arg))
 235          (true #q1))
 236     (format t "tan(pi/4) = 1~%")
 237     (print-result y true))
 238   (let* ((arg (/ +pi+ 3))
 239          (y (tan arg))
 240          (true (sqrt #q3)))
 241     (format t "tan(pi/3) = sqrt(3)~%")
 242     (print-result y true)))
 243
 244 (defun test-asin ()
 245   (format t "~2&asin for special args~%")
 246   (let* ((arg #q.5)
 247          (y (asin arg))
 248          (true (/ +pi+ 6)))
 249     (format t "asin(1/2) = pi/6~%")
 250     (print-result y true))
 251   (let* ((arg (sqrt #q.5))
 252          (y (asin arg))
 253          (true (/ +pi+ 4)))
 254     (format t "asin(1/sqrt(2) = pi/4~%")
 255     (print-result y true))
 256   (let* ((arg (/ (sqrt #q3) 2))
 257          (y (asin arg))
 258          (true (/ +pi+ 3)))
 259     (format t "asin(sqrt(3)/2) = pi/3~%")
 260     (print-result y true)))
 261
 262 (defun test-log ()
 263   (format t "~2&Log for special args~%")
 264   (let* ((arg #q2)
 265          (y (log arg))
 266          (true +log2+))
 267     (format t "log(2)~%")
 268     (print-result y true))
 269   (let* ((arg #q10)
 270          (y (log arg))
 271          (true (make-instance 'qd-real :value qdi::+qd-log10+)))
 272     (format t "log(10)~%")
 273     (print-result y true))
 274   (let* ((arg (+ 1 (scale-float #q1 -80)))
 275          (y (log arg))
 276          (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25))
 277     (format t "log(1+2^-80)~%")
 278     (print-result y true)))
 279
 280 (defun test-sqrt ()
 281   (format t "~2&Sqrt for special args~%")
 282   (dolist (f '((#q1q200 #q1q100)
 283                (#q1q300 #q1q150)
 284                (#q1q308 #q1q154)
 285                (#q1q-200 #q1q-100)
 286                (#q1q-250 #q1q-125)))
 287     (destructuring-bind (arg true)
 288         f
 289       (let ((y (sqrt arg)))
 290         (format t "sqrt(~/qdi::qd-format/)~%" (qd-value arg))
 291         (print-result y true)))))
 292
 293 (defun all-tests ()
 294   (test2)
 295   (test3)
 296   (test4)
 297   (test5)
 298   (test6)
 299   (test-atan)
 300   (test-sin)
 301   (test-tan)
 302   (test-asin)
 303   (test-log)
 304   (test-sqrt))