1 ;;;; -*- Mode: lisp -*-
3 ;;;; Copyright (c) 2007 Raymond Toy
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
14 ;;;; The above copyright notice and this permission notice shall be
15 ;;;; included in all copies or substantial portions of the Software.
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.
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
)
40 (defun check-accuracy (limit est true
)
41 (let ((bits (bit-accuracy est true
)))
44 (list bits limit est true
)))))
46 (defvar *null
* (make-broadcast-stream))
48 ;;; Some simple tests from the Yozo Hida's qd package.
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
)))
55 (check-accuracy 213 val true
))
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
)))
63 (check-accuracy 202 val true
))
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
)))
71 (check-accuracy 211 val true
))
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
))
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
))
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
+))
95 (check-accuracy 212 y true
))
98 (rt:deftest oct.atan
.2
99 (let* ((arg (sqrt #q3
))
100 (y (/ (atan arg
) +pi
+))
102 (check-accuracy 212 y true
))
105 (rt:deftest oct.atan
.3
107 (y (/ (atan arg
) +pi
+))
109 (check-accuracy 212 y true
))
112 (rt:deftest oct.atan
.4
114 (y (/ (atan arg
) +pi
+))
116 (check-accuracy 212 y true
))
119 (rt:deftest oct.atan
.5
120 (let* ((arg #q-1q100
)
121 (y (/ (atan arg
) +pi
+))
123 (check-accuracy 212 y true
))
127 (defun atan-qd/duplication
(arg)
128 (make-instance 'qd-real
129 :value
(qdi::atan-qd
/duplication
(qd-value arg
))))
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
+))
137 (check-accuracy 212 y true
))
140 (rt:deftest oct.atan
/dup
.2
141 (let* ((arg (sqrt #q3
))
142 (y (/ (atan-qd/duplication arg
) +pi
+))
144 (check-accuracy 212 y true
))
147 (rt:deftest oct.atan
/dup
.3
149 (y (/ (atan-qd/duplication arg
) +pi
+))
151 (check-accuracy 212 y true
))
154 (rt:deftest oct.atan
/dup
.4
156 (y (/ (atan-qd/duplication arg
) +pi
+))
158 (check-accuracy 212 y true
))
161 (rt:deftest oct.atan
/dup
.5
162 (let* ((arg #q-1q100
)
163 (y (/ (atan-qd/duplication arg
) +pi
+))
165 (check-accuracy 212 y true
))
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
))))
174 (rt:deftest oct.atan
/cordic
.1
175 (let* ((arg (/ (sqrt #q3
)))
176 (y (/ (atan-qd/cordic arg
) +pi
+))
178 (check-accuracy 212 y true
))
181 (rt:deftest oct.atan
/cordic
.2
182 (let* ((arg (sqrt #q3
))
183 (y (/ (atan-qd/cordic arg
) +pi
+))
185 (check-accuracy 212 y true
))
188 (rt:deftest oct.atan
/cordic
.3
190 (y (/ (atan-qd/cordic arg
) +pi
+))
192 (check-accuracy 212 y true
))
195 (rt:deftest oct.atan
/cordic
.4
197 (y (/ (atan-qd/cordic arg
) +pi
+))
199 (check-accuracy 212 y true
))
202 (rt:deftest oct.atan
/cordic
.5
203 (let* ((arg #q-1q100
)
204 (y (/ (atan-qd/cordic arg
) +pi
+))
206 (check-accuracy 212 y true
))
210 ;;; Tests of sin where we know the analytical result.
211 (rt:deftest oct.sin
.1
212 (let* ((arg (/ +pi
+ 6))
215 (check-accuracy 212 y true
))
218 (rt:deftest oct.sin
.2
219 (let* ((arg (/ +pi
+ 4))
222 (check-accuracy 212 y true
))
225 (rt:deftest oct.sin
.3
226 (let* ((arg (/ +pi
+ 3))
228 (true (/ (sqrt #q3
) 2)))
229 (check-accuracy 212 y true
))
232 ;;; Tests of tan where we know the analytical result.
233 (rt:deftest oct.tan
.1
234 (let* ((arg (/ +pi
+ 6))
236 (true (/ (sqrt #q3
))))
237 (check-accuracy 212 y true
))
240 (rt:deftest oct.tan
.2
241 (let* ((arg (/ +pi
+ 4))
244 (check-accuracy 212 y true
))
247 (rt:deftest oct.tan
.3
248 (let* ((arg (/ +pi
+ 3))
251 (check-accuracy 212 y true
))
254 ;;; Tests of tan where we know the analytical result. Uses CORDIC
257 (defun tan/cordic
(arg)
258 (make-instance 'qd-real
259 :value
(qdi::tan-qd
/cordic
(qd-value arg
))))
261 (rt:deftest oct.tan
/cordic
.1
262 (let* ((arg (/ +pi
+ 6))
264 (true (/ (sqrt #q3
))))
265 (check-accuracy 211 y true
))
268 (rt:deftest oct.tan
/cordic
.2
269 (let* ((arg (/ +pi
+ 4))
272 (check-accuracy 211 y true
))
275 (rt:deftest oct.tan
/cordic
.3
276 (let* ((arg (/ +pi
+ 3))
279 (check-accuracy 210 y true
))
282 ;;; Tests of asin where we know the analytical result.
284 (rt:deftest oct.asin
.1
288 (check-accuracy 212 y true
))
291 (rt:deftest oct.asin
.2
292 (let* ((arg (sqrt #q
.5))
295 (check-accuracy 212 y true
))
298 (rt:deftest oct.asin
.3
299 (let* ((arg (/ (sqrt #q3
) 2))
302 (check-accuracy 212 y true
))
307 (rt:deftest oct.log
.1
310 (true (make-instance 'qd-real
:value qdi
::+qd-log2
+)))
311 (check-accuracy 212 y true
))
314 (rt:deftest oct.log
.2
317 (true (make-instance 'qd-real
:value qdi
::+qd-log10
+)))
318 (check-accuracy 207 y true
))
321 (rt:deftest oct.log
.3
322 (let* ((arg (+ 1 (scale-float #q1 -
80)))
324 (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25
))
325 (check-accuracy 212 y true
))
328 ;;; Tests of log using Newton iteration.
330 (defun log/newton
(arg)
331 (make-instance 'qd-real
332 :value
(qdi::log-qd
/newton
(qd-value arg
))))
334 (rt:deftest oct.log
/newton
.1
337 (true (make-instance 'qd-real
:value qdi
::+qd-log2
+)))
338 (check-accuracy 212 y true
))
341 (rt:deftest oct.log
/newton
.2
344 (true (make-instance 'qd-real
:value qdi
::+qd-log10
+)))
345 (check-accuracy 207 y true
))
348 (rt:deftest oct.log
/newton
.3
349 (let* ((arg (+ 1 (scale-float #q1 -
80)))
351 (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25
))
352 (check-accuracy 212 y true
))
355 ;;; Tests of log using AGM.
358 (make-instance 'qd-real
359 :value
(qdi::log-qd
/agm
(qd-value arg
))))
361 (rt:deftest oct.log
/agm
.1
364 (true (make-instance 'qd-real
:value qdi
::+qd-log2
+)))
365 (check-accuracy 203 y true
))
368 (rt:deftest oct.log
/agm
.2
371 (true (make-instance 'qd-real
:value qdi
::+qd-log10
+)))
372 (check-accuracy 205 y true
))
375 (rt:deftest oct.log
/agm
.3
376 (let* ((arg (+ 1 (scale-float #q1 -
80)))
378 (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25
))
379 (check-accuracy 123 y true
))
382 ;;; Tests of log using AGM2, a faster variaton of AGM.
384 (defun log/agm2
(arg)
385 (make-instance 'qd-real
386 :value
(qdi::log-qd
/agm2
(qd-value arg
))))
388 (rt:deftest oct.log
/agm2.1
391 (true (make-instance 'qd-real
:value qdi
::+qd-log2
+)))
392 (check-accuracy 203 y true
))
395 (rt:deftest oct.log
/agm2.2
398 (true (make-instance 'qd-real
:value qdi
::+qd-log10
+)))
399 (check-accuracy 205 y true
))
402 (rt:deftest oct.log
/agm2.3
403 (let* ((arg (+ 1 (scale-float #q1 -
80)))
405 (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25
))
406 (check-accuracy 123 y true
))
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
))))
414 (rt:deftest oct.log
/agm3.1
417 (true (make-instance 'qd-real
:value qdi
::+qd-log2
+)))
418 (check-accuracy 203 y true
))
421 (rt:deftest oct.log
/agm3.2
424 (true (make-instance 'qd-real
:value qdi
::+qd-log10
+)))
425 (check-accuracy 205 y true
))
428 (rt:deftest oct.log
/agm3.3
429 (let* ((arg (+ 1 (scale-float #q1 -
80)))
431 (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25
))
432 (check-accuracy 123 y true
))
435 ;;; Tests of sqrt to make sure we overflow or underflow where we
438 (rt:deftest oct.sqrt
.1
442 (check-accuracy 212 y true
))
445 (rt:deftest oct.sqrt
.2
449 (check-accuracy 212 y true
))
452 (rt:deftest oct.sqrt
.3
456 (check-accuracy 212 y true
))
459 (rt:deftest oct.sqrt
.4
460 (let* ((arg #q1q-200
)
463 (check-accuracy 212 y true
))
466 (rt:deftest oct.sqrt
.5
467 (let* ((arg #q1q-250
)
470 (check-accuracy 212 y true
))
473 ;;; Tests of log1p(x) = log(1+x), using the duplication formula.
475 (defun log1p/dup
(arg)
476 (make-instance 'qd-real
477 :value
(qdi::log1p-qd
/duplication
(qd-value arg
))))
479 (rt:deftest oct.log1p
.1
482 (true #q2.3025850929940456840179914546843642076011014886287729760333279009675726096773525q0
))
483 (check-accuracy 212 y true
))
486 (rt:deftest oct.log1p
.2
487 (let* ((arg (scale-float #q1 -
80))
489 (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25
))
490 (check-accuracy 212 y true
))
493 ;;; Tests of expm1(x) = exp(x) - 1, using a Taylor series with
494 ;;; argument reduction.
496 (defun expm1/series
(arg)
497 (make-instance 'qd-real
498 :value
(qdi::expm1-qd
/series
(qd-value arg
))))
500 (rt:deftest oct.expm1
/series
.1
502 (y (expm1/series arg
))
504 (check-accuracy 212 y true
))
507 (rt:deftest oct.expm1
/series
.2
509 (y (expm1/series arg
))
510 (true #q1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0
))
511 (check-accuracy 211 y true
))
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
))
521 ;;; Tests of expm1(x) = exp(x) - 1, using duplication formula.
523 (defun expm1/dup
(arg)
524 (make-instance 'qd-real
525 :value
(qdi::expm1-qd
/duplication
(qd-value arg
))))
528 (rt:deftest oct.expm1
/dup
.1
532 (check-accuracy 212 y true
))
535 (rt:deftest oct.expm1
/dup
.2
538 (true #q1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0
))
539 (check-accuracy 211 y true
))
542 (rt:deftest oct.expm1
/dup
.3
543 (let* ((arg (scale-float #q1 -
100))
545 (true #q7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31
))
546 (check-accuracy 211 y true
))