1 | ;;;; -*- Mode: lisp -*- |
---|
2 | ;;;; |
---|
3 | ;;;; Copyright (c) 2007,2011 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 #:oct) |
---|
27 | |
---|
28 | (eval-when (:compile-toplevel :load-toplevel :execute) |
---|
29 | (setf *readtable* *oct-readtable*)) |
---|
30 | |
---|
31 | ;; For the tests, we need to turn off underflow for clisp. |
---|
32 | #+clisp |
---|
33 | (ext:without-package-lock () |
---|
34 | (setq sys::*inhibit-floating-point-underflow* t)) |
---|
35 | |
---|
36 | ;; Compute how many bits are the same for two numbers EST and TRUE. |
---|
37 | ;; Return T if they are identical. |
---|
38 | (defun bit-accuracy (est true) |
---|
39 | (let* ((diff (abs (- est true))) |
---|
40 | (err (float (if (zerop true) |
---|
41 | diff |
---|
42 | (/ diff (abs true))) |
---|
43 | 1d0))) |
---|
44 | (if (zerop diff) |
---|
45 | t |
---|
46 | (- (log err 2))))) |
---|
47 | |
---|
48 | ;; Check actual value EST is with LIMIT bits of the true value TRUE. |
---|
49 | ;; If so, return NIL. Otherwise, return a list of the actual bits of |
---|
50 | ;; accuracy, the desired accuracy, and the values. This is mostly to |
---|
51 | ;; make it easy to see what the actual accuracy was and the arguments |
---|
52 | ;; for the test, which is important for the tests that use random |
---|
53 | ;; values. |
---|
54 | (defun check-accuracy (limit est true) |
---|
55 | (let ((bits (bit-accuracy est true))) |
---|
56 | (if (not (eq bits t)) |
---|
57 | (if (and (not (float-nan-p (realpart est))) |
---|
58 | (not (float-nan-p bits)) |
---|
59 | (< bits limit)) |
---|
60 | (list bits limit est true))))) |
---|
61 | |
---|
62 | (defvar *null* (make-broadcast-stream)) |
---|
63 | |
---|
64 | ;;; Some simple tests from the Yozo Hida's qd package. |
---|
65 | |
---|
66 | (rt:deftest float.1 |
---|
67 | (float 3/2) |
---|
68 | 1.5) |
---|
69 | |
---|
70 | (rt:deftest float.2 |
---|
71 | (float 3/2 1d0) |
---|
72 | 1.5d0) |
---|
73 | |
---|
74 | (rt:deftest float.3 |
---|
75 | (float 1.5d0) |
---|
76 | 1.5d0) |
---|
77 | |
---|
78 | (rt:deftest float.4 |
---|
79 | (= (float #q1.5) #q1.5) |
---|
80 | t) |
---|
81 | |
---|
82 | (rt:deftest ceiling-d.1 |
---|
83 | (multiple-value-list (ceiling -50d0)) |
---|
84 | (-50 0d0)) |
---|
85 | |
---|
86 | (rt:deftest ceiling-d.2 |
---|
87 | (let ((z -50.1d0)) |
---|
88 | (multiple-value-bind (res rem) |
---|
89 | (ceiling -50.1d0) |
---|
90 | (list res (= z (+ res rem))))) |
---|
91 | (-50 t)) |
---|
92 | |
---|
93 | (rt:deftest ceiling-q.1 |
---|
94 | (multiple-value-bind (res rem) |
---|
95 | (ceiling #q-50q0) |
---|
96 | (list res (zerop rem))) |
---|
97 | (-50 t)) |
---|
98 | |
---|
99 | (rt:deftest ceiling-q.2 |
---|
100 | (let ((z #q-50.1q0)) |
---|
101 | (multiple-value-bind (res rem) |
---|
102 | (ceiling z) |
---|
103 | (list res (= z (+ res rem))))) |
---|
104 | (-50 t)) |
---|
105 | |
---|
106 | (rt:deftest truncate-d.1 |
---|
107 | (multiple-value-list (truncate -50d0)) |
---|
108 | (-50 0d0)) |
---|
109 | |
---|
110 | (rt:deftest truncate-q.1 |
---|
111 | (multiple-value-bind (res rem) |
---|
112 | (truncate #q-50q0) |
---|
113 | (list res (zerop rem))) |
---|
114 | (-50 t)) |
---|
115 | |
---|
116 | (rt:deftest fceiling-d.1 |
---|
117 | (multiple-value-list (fceiling -50d0)) |
---|
118 | (-50d0 0d0)) |
---|
119 | |
---|
120 | (rt:deftest fceiling-d.2 |
---|
121 | (let ((z -50.1d0)) |
---|
122 | (multiple-value-bind (res rem) |
---|
123 | (fceiling -50.1d0) |
---|
124 | (list res (= z (+ res rem))))) |
---|
125 | (-50d0 t)) |
---|
126 | |
---|
127 | (rt:deftest fceiling-q.1 |
---|
128 | (multiple-value-bind (res rem) |
---|
129 | (fceiling #q-50q0) |
---|
130 | (list (= res -50) (zerop rem))) |
---|
131 | (t t)) |
---|
132 | |
---|
133 | (rt:deftest fceiling-q.2 |
---|
134 | (let ((z #q-50.1q0)) |
---|
135 | (multiple-value-bind (res rem) |
---|
136 | (fceiling z) |
---|
137 | (list (= res -50) (= z (+ res rem))))) |
---|
138 | (t t)) |
---|
139 | |
---|
140 | (rt:deftest ftruncate-d.1 |
---|
141 | (multiple-value-list (ftruncate -50d0)) |
---|
142 | (-50d0 0d0)) |
---|
143 | |
---|
144 | (rt:deftest ftruncate-q.1 |
---|
145 | (multiple-value-bind (res rem) |
---|
146 | (ftruncate #q-50q0) |
---|
147 | (list (= res -50) (zerop rem))) |
---|
148 | (t t)) |
---|
149 | |
---|
150 | ;; Pi via Machin's formula |
---|
151 | (rt:deftest oct.pi.machin |
---|
152 | (let* ((*standard-output* *null*) |
---|
153 | (val (make-instance 'qd-real :value (octi::test2 nil))) |
---|
154 | (true oct:+pi+)) |
---|
155 | (check-accuracy 213 val true)) |
---|
156 | nil) |
---|
157 | |
---|
158 | ;; Pi via Salamin-Brent algorithm |
---|
159 | (rt:deftest oct.pi.salamin-brent |
---|
160 | (let* ((*standard-output* *null*) |
---|
161 | (val (make-instance 'qd-real :value (octi::test3 nil))) |
---|
162 | (true oct:+pi+)) |
---|
163 | (check-accuracy 202 val true)) |
---|
164 | nil) |
---|
165 | |
---|
166 | ;; Pi via Borweign's Quartic formula |
---|
167 | (rt:deftest oct.pi.borweign |
---|
168 | (let* ((*standard-output* *null*) |
---|
169 | (val (make-instance 'qd-real :value (octi::test4 nil))) |
---|
170 | (true oct:+pi+)) |
---|
171 | (check-accuracy 211 val true)) |
---|
172 | nil) |
---|
173 | |
---|
174 | ;; e via Taylor series |
---|
175 | (rt:deftest oct.e.taylor |
---|
176 | (let* ((*standard-output* *null*) |
---|
177 | (val (make-instance 'qd-real :value (octi::test5 nil))) |
---|
178 | (true (make-instance 'qd-real :value octi::+qd-e+))) |
---|
179 | (check-accuracy 212 val true)) |
---|
180 | nil) |
---|
181 | |
---|
182 | ;; log(2) via Taylor series |
---|
183 | (rt:deftest oct.log2.taylor |
---|
184 | (let* ((*standard-output* *null*) |
---|
185 | (val (make-instance 'qd-real :value (octi::test6 nil))) |
---|
186 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
187 | (check-accuracy 212 val true)) |
---|
188 | nil) |
---|
189 | |
---|
190 | ;;; Tests of atan where we know the analytical result |
---|
191 | (rt:deftest oct.atan.1 |
---|
192 | (let* ((arg (/ (sqrt #q3))) |
---|
193 | (y (/ (atan arg) +pi+)) |
---|
194 | (true (/ #q6))) |
---|
195 | (check-accuracy 212 y true)) |
---|
196 | nil) |
---|
197 | |
---|
198 | (rt:deftest oct.atan.2 |
---|
199 | (let* ((arg (sqrt #q3)) |
---|
200 | (y (/ (atan arg) +pi+)) |
---|
201 | (true (/ #q3))) |
---|
202 | (check-accuracy 212 y true)) |
---|
203 | nil) |
---|
204 | |
---|
205 | (rt:deftest oct.atan.3 |
---|
206 | (let* ((arg #q1) |
---|
207 | (y (/ (atan arg) +pi+)) |
---|
208 | (true (/ #q4))) |
---|
209 | (check-accuracy 212 y true)) |
---|
210 | nil) |
---|
211 | |
---|
212 | (rt:deftest oct.atan.4 |
---|
213 | (let* ((arg #q1q100) |
---|
214 | (y (/ (atan arg) +pi+)) |
---|
215 | (true #q.5)) |
---|
216 | (check-accuracy 212 y true)) |
---|
217 | nil) |
---|
218 | |
---|
219 | (rt:deftest oct.atan.5 |
---|
220 | (let* ((arg #q-1q100) |
---|
221 | (y (/ (atan arg) +pi+)) |
---|
222 | (true #q-.5)) |
---|
223 | (check-accuracy 212 y true)) |
---|
224 | nil) |
---|
225 | |
---|
226 | |
---|
227 | (defun atan-qd/duplication (arg) |
---|
228 | (make-instance 'qd-real |
---|
229 | :value (octi::atan-qd/duplication (qd-value arg)))) |
---|
230 | |
---|
231 | ;;; Tests of atan where we know the analytical result. Same tests, |
---|
232 | ;;; but using the atan duplication formula. |
---|
233 | (rt:deftest oct.atan/dup.1 |
---|
234 | (let* ((arg (/ (sqrt #q3))) |
---|
235 | (y (/ (atan-qd/duplication arg) +pi+)) |
---|
236 | (true (/ #q6))) |
---|
237 | (check-accuracy 212 y true)) |
---|
238 | nil) |
---|
239 | |
---|
240 | (rt:deftest oct.atan/dup.2 |
---|
241 | (let* ((arg (sqrt #q3)) |
---|
242 | (y (/ (atan-qd/duplication arg) +pi+)) |
---|
243 | (true (/ #q3))) |
---|
244 | (check-accuracy 212 y true)) |
---|
245 | nil) |
---|
246 | |
---|
247 | (rt:deftest oct.atan/dup.3 |
---|
248 | (let* ((arg #q1) |
---|
249 | (y (/ (atan-qd/duplication arg) +pi+)) |
---|
250 | (true (/ #q4))) |
---|
251 | (check-accuracy 212 y true)) |
---|
252 | nil) |
---|
253 | |
---|
254 | (rt:deftest oct.atan/dup.4 |
---|
255 | (let* ((arg #q1q100) |
---|
256 | (y (/ (atan-qd/duplication arg) +pi+)) |
---|
257 | (true #q.5)) |
---|
258 | (check-accuracy 212 y true)) |
---|
259 | nil) |
---|
260 | |
---|
261 | (rt:deftest oct.atan/dup.5 |
---|
262 | (let* ((arg #q-1q100) |
---|
263 | (y (/ (atan-qd/duplication arg) +pi+)) |
---|
264 | (true #q-.5)) |
---|
265 | (check-accuracy 212 y true)) |
---|
266 | nil) |
---|
267 | |
---|
268 | ;;; Tests of atan where we know the analytical result. Same tests, |
---|
269 | ;;; but using a CORDIC implementation. |
---|
270 | (defun atan-qd/cordic (arg) |
---|
271 | (make-instance 'qd-real |
---|
272 | :value (octi::atan-qd/cordic (qd-value arg)))) |
---|
273 | |
---|
274 | (rt:deftest oct.atan/cordic.1 |
---|
275 | (let* ((arg (/ (sqrt #q3))) |
---|
276 | (y (/ (atan-qd/cordic arg) +pi+)) |
---|
277 | (true (/ #q6))) |
---|
278 | (check-accuracy 212 y true)) |
---|
279 | nil) |
---|
280 | |
---|
281 | (rt:deftest oct.atan/cordic.2 |
---|
282 | (let* ((arg (sqrt #q3)) |
---|
283 | (y (/ (atan-qd/cordic arg) +pi+)) |
---|
284 | (true (/ #q3))) |
---|
285 | (check-accuracy 212 y true)) |
---|
286 | nil) |
---|
287 | |
---|
288 | (rt:deftest oct.atan/cordic.3 |
---|
289 | (let* ((arg #q1) |
---|
290 | (y (/ (atan-qd/cordic arg) +pi+)) |
---|
291 | (true (/ #q4))) |
---|
292 | (check-accuracy 212 y true)) |
---|
293 | nil) |
---|
294 | |
---|
295 | (rt:deftest oct.atan/cordic.4 |
---|
296 | (let* ((arg #q1q100) |
---|
297 | (y (/ (atan-qd/cordic arg) +pi+)) |
---|
298 | (true #q.5)) |
---|
299 | (check-accuracy 212 y true)) |
---|
300 | nil) |
---|
301 | |
---|
302 | (rt:deftest oct.atan/cordic.5 |
---|
303 | (let* ((arg #q-1q100) |
---|
304 | (y (/ (atan-qd/cordic arg) +pi+)) |
---|
305 | (true #q-.5)) |
---|
306 | (check-accuracy 212 y true)) |
---|
307 | nil) |
---|
308 | |
---|
309 | |
---|
310 | ;;; Tests of sin where we know the analytical result. |
---|
311 | (rt:deftest oct.sin.1 |
---|
312 | (let* ((arg (/ +pi+ 6)) |
---|
313 | (y (sin arg)) |
---|
314 | (true #q.5)) |
---|
315 | (check-accuracy 212 y true)) |
---|
316 | nil) |
---|
317 | |
---|
318 | (rt:deftest oct.sin.2 |
---|
319 | (let* ((arg (/ +pi+ 4)) |
---|
320 | (y (sin arg)) |
---|
321 | (true (sqrt #q.5))) |
---|
322 | (check-accuracy 212 y true)) |
---|
323 | nil) |
---|
324 | |
---|
325 | (rt:deftest oct.sin.3 |
---|
326 | (let* ((arg (/ +pi+ 3)) |
---|
327 | (y (sin arg)) |
---|
328 | (true (/ (sqrt #q3) 2))) |
---|
329 | (check-accuracy 212 y true)) |
---|
330 | nil) |
---|
331 | |
---|
332 | (rt:deftest oct.big-sin.1 |
---|
333 | (let* ((arg (oct:make-qd (ash 1 120))) |
---|
334 | (y (sin arg)) |
---|
335 | (true #q3.778201093607520226555484700569229919605866976512306642257987199414885q-1)) |
---|
336 | (check-accuracy 205 y true)) |
---|
337 | nil) |
---|
338 | |
---|
339 | (rt:deftest oct.big-sin.2 |
---|
340 | (let* ((arg (oct:make-qd (ash 1 1023))) |
---|
341 | (y (sin arg)) |
---|
342 | (true #q5.631277798508840134529434079444683477103854907361251399182750155357133q-1)) |
---|
343 | (check-accuracy 205 y true)) |
---|
344 | nil) |
---|
345 | |
---|
346 | ;;; Tests of tan where we know the analytical result. |
---|
347 | (rt:deftest oct.tan.1 |
---|
348 | (let* ((arg (/ +pi+ 6)) |
---|
349 | (y (tan arg)) |
---|
350 | (true (/ (sqrt #q3)))) |
---|
351 | (check-accuracy 212 y true)) |
---|
352 | nil) |
---|
353 | |
---|
354 | (rt:deftest oct.tan.2 |
---|
355 | (let* ((arg (/ +pi+ 4)) |
---|
356 | (y (tan arg)) |
---|
357 | (true #q1)) |
---|
358 | (check-accuracy 212 y true)) |
---|
359 | nil) |
---|
360 | |
---|
361 | (rt:deftest oct.tan.3 |
---|
362 | (let* ((arg (/ +pi+ 3)) |
---|
363 | (y (tan arg)) |
---|
364 | (true (sqrt #q3))) |
---|
365 | (check-accuracy 212 y true)) |
---|
366 | nil) |
---|
367 | |
---|
368 | ;;; Tests of tan where we know the analytical result. Uses CORDIC |
---|
369 | ;;; algorithm. |
---|
370 | |
---|
371 | (defun tan/cordic (arg) |
---|
372 | (make-instance 'qd-real |
---|
373 | :value (octi::tan-qd/cordic (qd-value arg)))) |
---|
374 | |
---|
375 | (rt:deftest oct.tan/cordic.1 |
---|
376 | (let* ((arg (/ +pi+ 6)) |
---|
377 | (y (tan/cordic arg)) |
---|
378 | (true (/ (sqrt #q3)))) |
---|
379 | (check-accuracy 211 y true)) |
---|
380 | nil) |
---|
381 | |
---|
382 | (rt:deftest oct.tan/cordic.2 |
---|
383 | (let* ((arg (/ +pi+ 4)) |
---|
384 | (y (tan/cordic arg)) |
---|
385 | (true #q1)) |
---|
386 | (check-accuracy 211 y true)) |
---|
387 | nil) |
---|
388 | |
---|
389 | (rt:deftest oct.tan/cordic.3 |
---|
390 | (let* ((arg (/ +pi+ 3)) |
---|
391 | (y (tan/cordic arg)) |
---|
392 | (true (sqrt #q3))) |
---|
393 | (check-accuracy 210 y true)) |
---|
394 | nil) |
---|
395 | |
---|
396 | ;;; Tests of asin where we know the analytical result. |
---|
397 | |
---|
398 | (rt:deftest oct.asin.1 |
---|
399 | (let* ((arg #q.5) |
---|
400 | (y (asin arg)) |
---|
401 | (true (/ +pi+ 6))) |
---|
402 | (check-accuracy 212 y true)) |
---|
403 | nil) |
---|
404 | |
---|
405 | (rt:deftest oct.asin.2 |
---|
406 | (let* ((arg (sqrt #q.5)) |
---|
407 | (y (asin arg)) |
---|
408 | (true (/ +pi+ 4))) |
---|
409 | (check-accuracy 212 y true)) |
---|
410 | nil) |
---|
411 | |
---|
412 | (rt:deftest oct.asin.3 |
---|
413 | (let* ((arg (/ (sqrt #q3) 2)) |
---|
414 | (y (asin arg)) |
---|
415 | (true (/ +pi+ 3))) |
---|
416 | (check-accuracy 212 y true)) |
---|
417 | nil) |
---|
418 | |
---|
419 | ;;; Tests of log. |
---|
420 | |
---|
421 | (rt:deftest oct.log.1 |
---|
422 | (let* ((arg #q2) |
---|
423 | (y (log arg)) |
---|
424 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
425 | (check-accuracy 212 y true)) |
---|
426 | nil) |
---|
427 | |
---|
428 | (rt:deftest oct.log.2 |
---|
429 | (let* ((arg #q10) |
---|
430 | (y (log arg)) |
---|
431 | (true (make-instance 'qd-real :value octi::+qd-log10+))) |
---|
432 | (check-accuracy 207 y true)) |
---|
433 | nil) |
---|
434 | |
---|
435 | (rt:deftest oct.log.3 |
---|
436 | (let* ((arg (+ 1 (scale-float #q1 -80))) |
---|
437 | (y (log arg)) |
---|
438 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
439 | (check-accuracy 212 y true)) |
---|
440 | nil) |
---|
441 | |
---|
442 | ;;; Tests of log using Newton iteration. |
---|
443 | |
---|
444 | (defun log/newton (arg) |
---|
445 | (make-instance 'qd-real |
---|
446 | :value (octi::log-qd/newton (qd-value arg)))) |
---|
447 | |
---|
448 | (rt:deftest oct.log/newton.1 |
---|
449 | (let* ((arg #q2) |
---|
450 | (y (log/newton arg)) |
---|
451 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
452 | (check-accuracy 212 y true)) |
---|
453 | nil) |
---|
454 | |
---|
455 | (rt:deftest oct.log/newton.2 |
---|
456 | (let* ((arg #q10) |
---|
457 | (y (log/newton arg)) |
---|
458 | (true (make-instance 'qd-real :value octi::+qd-log10+))) |
---|
459 | (check-accuracy 207 y true)) |
---|
460 | nil) |
---|
461 | |
---|
462 | (rt:deftest oct.log/newton.3 |
---|
463 | (let* ((arg (+ 1 (scale-float #q1 -80))) |
---|
464 | (y (log/newton arg)) |
---|
465 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
466 | (check-accuracy 212 y true)) |
---|
467 | nil) |
---|
468 | |
---|
469 | ;;; Tests of log using AGM. |
---|
470 | |
---|
471 | (defun log/agm (arg) |
---|
472 | (make-instance 'qd-real |
---|
473 | :value (octi::log-qd/agm (qd-value arg)))) |
---|
474 | |
---|
475 | (rt:deftest oct.log/agm.1 |
---|
476 | (let* ((arg #q2) |
---|
477 | (y (log/agm arg)) |
---|
478 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
479 | (check-accuracy 203 y true)) |
---|
480 | nil) |
---|
481 | |
---|
482 | (rt:deftest oct.log/agm.2 |
---|
483 | (let* ((arg #q10) |
---|
484 | (y (log/agm arg)) |
---|
485 | (true (make-instance 'qd-real :value octi::+qd-log10+))) |
---|
486 | (check-accuracy 205 y true)) |
---|
487 | nil) |
---|
488 | |
---|
489 | (rt:deftest oct.log/agm.3 |
---|
490 | (let* ((arg (+ 1 (scale-float #q1 -80))) |
---|
491 | (y (log/agm arg)) |
---|
492 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
493 | (check-accuracy 123 y true)) |
---|
494 | nil) |
---|
495 | |
---|
496 | ;;; Tests of log using AGM2, a faster variaton of AGM. |
---|
497 | |
---|
498 | (defun log/agm2 (arg) |
---|
499 | (make-instance 'qd-real |
---|
500 | :value (octi::log-qd/agm2 (qd-value arg)))) |
---|
501 | |
---|
502 | (rt:deftest oct.log/agm2.1 |
---|
503 | (let* ((arg #q2) |
---|
504 | (y (log/agm2 arg)) |
---|
505 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
506 | (check-accuracy 203 y true)) |
---|
507 | nil) |
---|
508 | |
---|
509 | (rt:deftest oct.log/agm2.2 |
---|
510 | (let* ((arg #q10) |
---|
511 | (y (log/agm2 arg)) |
---|
512 | (true (make-instance 'qd-real :value octi::+qd-log10+))) |
---|
513 | (check-accuracy 205 y true)) |
---|
514 | nil) |
---|
515 | |
---|
516 | (rt:deftest oct.log/agm2.3 |
---|
517 | (let* ((arg (+ 1 (scale-float #q1 -80))) |
---|
518 | (y (log/agm2 arg)) |
---|
519 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
520 | (check-accuracy 123 y true)) |
---|
521 | nil) |
---|
522 | |
---|
523 | ;;; Tests of log using AGM3, a faster variation of AGM2. |
---|
524 | (defun log/agm3 (arg) |
---|
525 | (make-instance 'qd-real |
---|
526 | :value (octi::log-qd/agm3 (qd-value arg)))) |
---|
527 | |
---|
528 | (rt:deftest oct.log/agm3.1 |
---|
529 | (let* ((arg #q2) |
---|
530 | (y (log/agm3 arg)) |
---|
531 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
532 | (check-accuracy 203 y true)) |
---|
533 | nil) |
---|
534 | |
---|
535 | (rt:deftest oct.log/agm3.2 |
---|
536 | (let* ((arg #q10) |
---|
537 | (y (log/agm3 arg)) |
---|
538 | (true (make-instance 'qd-real :value octi::+qd-log10+))) |
---|
539 | (check-accuracy 205 y true)) |
---|
540 | nil) |
---|
541 | |
---|
542 | (rt:deftest oct.log/agm3.3 |
---|
543 | (let* ((arg (+ 1 (scale-float #q1 -80))) |
---|
544 | (y (log/agm3 arg)) |
---|
545 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
546 | (check-accuracy 123 y true)) |
---|
547 | nil) |
---|
548 | |
---|
549 | ;;; Tests of sqrt to make sure we overflow or underflow where we |
---|
550 | ;;; shouldn't. |
---|
551 | |
---|
552 | (rt:deftest oct.sqrt.1 |
---|
553 | (let* ((arg #q1q200) |
---|
554 | (y (sqrt arg)) |
---|
555 | (true #q1q100)) |
---|
556 | (check-accuracy 212 y true)) |
---|
557 | nil) |
---|
558 | |
---|
559 | (rt:deftest oct.sqrt.2 |
---|
560 | (let* ((arg #q1q200) |
---|
561 | (y (sqrt arg)) |
---|
562 | (true #q1q100)) |
---|
563 | (check-accuracy 212 y true)) |
---|
564 | nil) |
---|
565 | |
---|
566 | (rt:deftest oct.sqrt.3 |
---|
567 | (let* ((arg #q1q300) |
---|
568 | (y (sqrt arg)) |
---|
569 | (true #q1q150)) |
---|
570 | (check-accuracy 212 y true)) |
---|
571 | nil) |
---|
572 | |
---|
573 | (rt:deftest oct.sqrt.4 |
---|
574 | (let* ((arg #q1q-200) |
---|
575 | (y (sqrt arg)) |
---|
576 | (true #q1q-100)) |
---|
577 | (check-accuracy 212 y true)) |
---|
578 | nil) |
---|
579 | |
---|
580 | (rt:deftest oct.sqrt.5 |
---|
581 | (let* ((arg #q1q-250) |
---|
582 | (y (sqrt arg)) |
---|
583 | (true #q1q-125)) |
---|
584 | (check-accuracy 212 y true)) |
---|
585 | nil) |
---|
586 | |
---|
587 | ;;; Tests of log1p(x) = log(1+x), using the duplication formula. |
---|
588 | |
---|
589 | (defun log1p/dup (arg) |
---|
590 | (make-instance 'qd-real |
---|
591 | :value (octi::log1p-qd/duplication (qd-value arg)))) |
---|
592 | |
---|
593 | (rt:deftest oct.log1p.1 |
---|
594 | (let* ((arg #q9) |
---|
595 | (y (log1p/dup arg)) |
---|
596 | (true #q2.3025850929940456840179914546843642076011014886287729760333279009675726096773525q0)) |
---|
597 | (check-accuracy 212 y true)) |
---|
598 | nil) |
---|
599 | |
---|
600 | (rt:deftest oct.log1p.2 |
---|
601 | (let* ((arg (scale-float #q1 -80)) |
---|
602 | (y (log1p/dup arg)) |
---|
603 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
604 | (check-accuracy 212 y true)) |
---|
605 | nil) |
---|
606 | |
---|
607 | ;;; Tests of expm1(x) = exp(x) - 1, using a Taylor series with |
---|
608 | ;;; argument reduction. |
---|
609 | |
---|
610 | (defun expm1/series (arg) |
---|
611 | (make-instance 'qd-real |
---|
612 | :value (octi::expm1-qd/series (qd-value arg)))) |
---|
613 | |
---|
614 | (rt:deftest oct.expm1/series.1 |
---|
615 | (let* ((arg #q0) |
---|
616 | (y (expm1/series arg)) |
---|
617 | (true #q0)) |
---|
618 | (check-accuracy 212 y true)) |
---|
619 | nil) |
---|
620 | |
---|
621 | (rt:deftest oct.expm1/series.2 |
---|
622 | (let* ((arg #q1) |
---|
623 | (y (expm1/series arg)) |
---|
624 | (true #q1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0)) |
---|
625 | (check-accuracy 211 y true)) |
---|
626 | nil) |
---|
627 | |
---|
628 | (rt:deftest oct.expm1/series.3 |
---|
629 | (let* ((arg (scale-float #q1 -100)) |
---|
630 | (y (expm1/series arg)) |
---|
631 | (true #q7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31)) |
---|
632 | (check-accuracy 211 y true)) |
---|
633 | nil) |
---|
634 | |
---|
635 | ;;; Tests of expm1(x) = exp(x) - 1, using duplication formula. |
---|
636 | |
---|
637 | (defun expm1/dup (arg) |
---|
638 | (make-instance 'qd-real |
---|
639 | :value (octi::expm1-qd/duplication (qd-value arg)))) |
---|
640 | |
---|
641 | |
---|
642 | (rt:deftest oct.expm1/dup.1 |
---|
643 | (let* ((arg #q0) |
---|
644 | (y (expm1/dup arg)) |
---|
645 | (true #q0)) |
---|
646 | (check-accuracy 212 y true)) |
---|
647 | nil) |
---|
648 | |
---|
649 | (rt:deftest oct.expm1/dup.2 |
---|
650 | (let* ((arg #q1) |
---|
651 | (y (expm1/dup arg)) |
---|
652 | (true #q1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0)) |
---|
653 | (check-accuracy 211 y true)) |
---|
654 | nil) |
---|
655 | |
---|
656 | (rt:deftest oct.expm1/dup.3 |
---|
657 | (let* ((arg (scale-float #q1 -100)) |
---|
658 | (y (expm1/dup arg)) |
---|
659 | (true #q7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31)) |
---|
660 | (check-accuracy 211 y true)) |
---|
661 | nil) |
---|
662 | |
---|
663 | ;; If we screw up integer-decode-qd, printing is wrong. Here is one |
---|
664 | ;; case where integer-decode-qd was screwed up and printing the wrong |
---|
665 | ;; thing. |
---|
666 | (rt:deftest oct.integer-decode.1 |
---|
667 | (multiple-value-bind (frac exp s) |
---|
668 | (octi:integer-decode-qd (octi::%make-qd-d -0.03980126756814893d0 |
---|
669 | -2.7419792323327893d-18 |
---|
670 | 0d0 0d0)) |
---|
671 | (unless (and (eql frac 103329998279901916046530991816704) |
---|
672 | (eql exp -111) |
---|
673 | (eql s -1)) |
---|
674 | (list frac exp s))) |
---|
675 | nil) |
---|
676 | |
---|
677 | ;;; |
---|
678 | ;;; Add a few tests for the branch cuts. Many of these tests assume |
---|
679 | ;;; that Lisp has support for signed zeroes. If not, these tests are |
---|
680 | ;;; probably wrong. |
---|
681 | |
---|
682 | (defun check-signs (fun arg expected) |
---|
683 | (let* ((z (funcall fun arg)) |
---|
684 | (x (realpart z)) |
---|
685 | (y (imagpart z))) |
---|
686 | ;; If the Lisp doesn't support signed zeroes, then this test |
---|
687 | ;; should always pass. |
---|
688 | (if (or (eql -0d0 0d0) |
---|
689 | (and (= (float-sign x) (float-sign (realpart expected))) |
---|
690 | (= (float-sign y) (float-sign (imagpart expected))))) |
---|
691 | t |
---|
692 | (list z expected fun arg)))) |
---|
693 | |
---|
694 | ;; asin has a branch cut on the real axis |x|>1. For x < -1, it is |
---|
695 | ;; continuous with quadrant II; for x > 1, continuous with quadrant |
---|
696 | ;; IV. |
---|
697 | (rt:deftest oct.asin-branch-neg.1 |
---|
698 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
699 | (check-signs #'asin -2d0 true)) |
---|
700 | t) |
---|
701 | |
---|
702 | (rt:deftest oct.asin-branch-neg.2 |
---|
703 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
704 | (check-signs #'asin #q-2 true)) |
---|
705 | t) |
---|
706 | |
---|
707 | (rt:deftest oct.asin-branch-neg.3 |
---|
708 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
709 | (check-signs #'asin #c(-2d0 0d0) true)) |
---|
710 | t) |
---|
711 | |
---|
712 | (rt:deftest oct.asin-branch-neg.4 |
---|
713 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
714 | (check-signs #'asin #q(-2 0) true)) |
---|
715 | t) |
---|
716 | |
---|
717 | (rt:deftest oct.asin-branch-neg.5 |
---|
718 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
719 | (check-signs #'asin #c(-2d0 -0d0) (conjugate true))) |
---|
720 | t) |
---|
721 | |
---|
722 | (rt:deftest oct.asin-branch-neg.6 |
---|
723 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
724 | (check-signs #'asin #q(-2d0 -0d0) (conjugate true))) |
---|
725 | t) |
---|
726 | |
---|
727 | (rt:deftest oct.asin-branch-pos.1 |
---|
728 | (let ((true (cl:asin #c(2d0 -1d-20)))) |
---|
729 | (check-signs #'asin #c(2d0 0d0) (conjugate true))) |
---|
730 | t) |
---|
731 | |
---|
732 | (rt:deftest oct.asin-branch-pos.2 |
---|
733 | (let ((true (cl:asin #c(2d0 -1d-20)))) |
---|
734 | (check-signs #'asin #q(2 0d0) (conjugate true))) |
---|
735 | t) |
---|
736 | |
---|
737 | (rt:deftest oct.asin-branch-pos.3 |
---|
738 | (let ((true (cl:asin #c(2d0 -1d-20)))) |
---|
739 | (check-signs #'asin #c(2d0 -0d0) true)) |
---|
740 | t) |
---|
741 | |
---|
742 | (rt:deftest oct.asin-branch-pos.4 |
---|
743 | (let ((true (cl:asin #c(2d0 -1d-20)))) |
---|
744 | (check-signs #'asin #q(2d0 -0d0) true)) |
---|
745 | t) |
---|
746 | |
---|
747 | ;; acos branch cut is the real axis, |x| > 1. For x < -1, it is |
---|
748 | ;; continuous with quadrant II; for x > 1, quadrant IV. |
---|
749 | |
---|
750 | (rt:deftest oct.acos-branch-neg.1 |
---|
751 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
752 | (check-signs #'acos -2d0 true)) |
---|
753 | t) |
---|
754 | |
---|
755 | (rt:deftest oct.acos-branch-neg.2 |
---|
756 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
757 | (check-signs #'acos #q-2 true)) |
---|
758 | t) |
---|
759 | |
---|
760 | (rt:deftest oct.acos-branch-neg.3 |
---|
761 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
762 | (check-signs #'acos #c(-2d0 0d0) true)) |
---|
763 | t) |
---|
764 | |
---|
765 | (rt:deftest oct.acos-branch-neg.4 |
---|
766 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
767 | (check-signs #'acos #q(-2 0) true)) |
---|
768 | t) |
---|
769 | |
---|
770 | (rt:deftest oct.acos-branch-neg.5 |
---|
771 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
772 | (check-signs #'acos #c(-2d0 -0d0) (conjugate true))) |
---|
773 | t) |
---|
774 | |
---|
775 | (rt:deftest oct.acos-branch-neg.6 |
---|
776 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
777 | (check-signs #'acos #q(-2d0 -0d0) (conjugate true))) |
---|
778 | t) |
---|
779 | |
---|
780 | (rt:deftest oct.acos-branch-pos.1 |
---|
781 | (let ((true (cl:acos #c(2d0 -1d-20)))) |
---|
782 | (check-signs #'acos #c(2d0 0d0) (conjugate true))) |
---|
783 | t) |
---|
784 | |
---|
785 | (rt:deftest oct.acos-branch-pos.2 |
---|
786 | (let ((true (cl:acos #c(2d0 -1d-20)))) |
---|
787 | (check-signs #'acos #q(2 0d0) (conjugate true))) |
---|
788 | t) |
---|
789 | |
---|
790 | (rt:deftest oct.acos-branch-pos.3 |
---|
791 | (let ((true (cl:acos #c(2d0 -1d-20)))) |
---|
792 | (check-signs #'acos #c(2d0 -0d0) true)) |
---|
793 | t) |
---|
794 | |
---|
795 | (rt:deftest oct.acos-branch-pos.4 |
---|
796 | (let ((true (cl:acos #c(2d0 -1d-20)))) |
---|
797 | (check-signs #'acos #q(2d0 -0d0) true)) |
---|
798 | t) |
---|
799 | |
---|
800 | ;; atan branch cut is the imaginary axis, |y| > 1. For y < -1, it is |
---|
801 | ;; continuous with quadrant IV; for x > 1, quadrant II. |
---|
802 | |
---|
803 | (rt:deftest oct.atan-branch-neg.1 |
---|
804 | (let ((true (cl:atan #c(1d-20 -2d0)))) |
---|
805 | (check-signs #'atan #c(0d0 -2d0) true)) |
---|
806 | t) |
---|
807 | |
---|
808 | (rt:deftest oct.atan-branch-neg.2 |
---|
809 | (let ((true (cl:atan #c(1d-20 -2d0)))) |
---|
810 | (check-signs #'atan #q(0 -2) true)) |
---|
811 | t) |
---|
812 | |
---|
813 | (rt:deftest oct.atan-branch-neg.3 |
---|
814 | (let ((true (cl:atan #c(-1d-20 -2d0)))) |
---|
815 | (check-signs #'atan #c(-0d0 -2d0) true)) |
---|
816 | t) |
---|
817 | |
---|
818 | (rt:deftest oct.atan-branch-neg.4 |
---|
819 | (let ((true (cl:atan #c(-1d-20 -2d0)))) |
---|
820 | (check-signs #'atan #q(-0d0 -2d0) true)) |
---|
821 | t) |
---|
822 | |
---|
823 | (rt:deftest oct.atan-branch-pos.1 |
---|
824 | (let ((true (cl:atan #c(1d-20 2d0)))) |
---|
825 | (check-signs #'atan #c(0d0 2d0) true)) |
---|
826 | t) |
---|
827 | |
---|
828 | (rt:deftest oct.atan-branch-pos.2 |
---|
829 | (let ((true (cl:atan #c(1d-20 2d0)))) |
---|
830 | (check-signs #'atan #q(0d0 2d0) true)) |
---|
831 | t) |
---|
832 | |
---|
833 | (rt:deftest oct.atan-branch-pos.3 |
---|
834 | (let ((true (cl:atan #c(-1d-20 2d0)))) |
---|
835 | (check-signs #'atan #c(-0d0 2d0) true)) |
---|
836 | t) |
---|
837 | |
---|
838 | (rt:deftest oct.atan-branch-pos.4 |
---|
839 | (let ((true (cl:atan #c(-1d-20 2d0)))) |
---|
840 | (check-signs #'atan #q(-0d0 2d0) true)) |
---|
841 | t) |
---|
842 | |
---|
843 | ;; Test x < -1. CLHS says for x < -1, atanh is continuous with quadrant III. |
---|
844 | (rt:deftest oct.atanh-branch-neg.1 |
---|
845 | (let ((true (cl:atanh #c(-2d0 -1d-20)))) |
---|
846 | (check-signs #'atanh -2d0 true)) |
---|
847 | t) |
---|
848 | |
---|
849 | (rt:deftest oct.atanh-branch-neg.2 |
---|
850 | (let ((true (cl:atanh #c(-2d0 -1d-20)))) |
---|
851 | (check-signs #'atanh #q-2 true)) |
---|
852 | t) |
---|
853 | |
---|
854 | ;; Test x > 1. CLHS says for x > 1, atanh is continus with quadrant I. |
---|
855 | (rt:deftest oct.atanh-branch-pos.1 |
---|
856 | (let ((true (cl:atanh #c(2d0 1d-20)))) |
---|
857 | (check-signs #'atanh 2d0 true)) |
---|
858 | t) |
---|
859 | |
---|
860 | (rt:deftest oct.atanh-branch-pos.2 |
---|
861 | (let ((true (cl:atanh #c(2d0 1d-20)))) |
---|
862 | (check-signs #'atanh #q2 true)) |
---|
863 | t) |
---|
864 | |
---|
865 | ;; elliptic_k(-1) = gamma(1/4)^2/2^(5/2)/sqrt(%pi) |
---|
866 | (rt:deftest oct.elliptic-k.1d |
---|
867 | (let* ((val (elliptic-k -1d0)) |
---|
868 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395293535207125115147766480714547q0)) |
---|
869 | (check-accuracy 53 val true)) |
---|
870 | nil) |
---|
871 | |
---|
872 | (rt:deftest oct.elliptic-k.1q |
---|
873 | (let* ((val (elliptic-k #q-1q0)) |
---|
874 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395293535207125115147766480714547q0)) |
---|
875 | (check-accuracy 210 val true)) |
---|
876 | nil) |
---|
877 | |
---|
878 | ;; elliptic_k(1/2) = %pi^(3/2)/2/gamma(3/4)^2 |
---|
879 | (rt:deftest oct.elliptic-k.2d |
---|
880 | (let* ((val (elliptic-k 0.5d0)) |
---|
881 | (true #q1.854074677301371918433850347195260046217598823521766905585928045056021776838119978357271861650371897q0)) |
---|
882 | (check-accuracy 53 val true)) |
---|
883 | nil) |
---|
884 | |
---|
885 | (rt:deftest oct.elliptic-k.2q |
---|
886 | (let* ((val (elliptic-k #q.5)) |
---|
887 | (true #q1.854074677301371918433850347195260046217598823521766905585928045056021776838119978357271861650371897q0)) |
---|
888 | (check-accuracy 210 val true)) |
---|
889 | nil) |
---|
890 | |
---|
891 | ;; jacobi_sn(K,1/2) = 1, where K = elliptic_k(1/2) |
---|
892 | (rt:deftest oct.jacobi-sn.1d |
---|
893 | (let* ((ek (elliptic-k .5d0)) |
---|
894 | (val (jacobi-sn ek .5d0))) |
---|
895 | (check-accuracy 54 val 1d0)) |
---|
896 | nil) |
---|
897 | |
---|
898 | (rt:deftest oct.jacobi-sn.1q |
---|
899 | (let* ((ek (elliptic-k #q.5)) |
---|
900 | (val (jacobi-sn ek #q.5))) |
---|
901 | (check-accuracy 212 val #q1)) |
---|
902 | nil) |
---|
903 | |
---|
904 | ;; jacobi_cn(K,1/2) = 0 |
---|
905 | (rt:deftest oct.jacobi-cn.1d |
---|
906 | (let* ((ek (elliptic-k .5d0)) |
---|
907 | (val (jacobi-cn ek .5d0))) |
---|
908 | (check-accuracy 50 val 0d0)) |
---|
909 | nil) |
---|
910 | |
---|
911 | (rt:deftest oct.jacobi-cn.1q |
---|
912 | (let* ((ek (elliptic-k #q.5)) |
---|
913 | (val (jacobi-cn ek #q.5))) |
---|
914 | (check-accuracy 210 val #q0)) |
---|
915 | nil) |
---|
916 | |
---|
917 | ;; jacobi-dn(K, 1/2) = sqrt(1/2) |
---|
918 | (rt:deftest oct.jacobi-dn.1d |
---|
919 | (let* ((ek (elliptic-k .5d0)) |
---|
920 | (true (sqrt .5d0)) |
---|
921 | (val (jacobi-dn ek .5d0))) |
---|
922 | (check-accuracy 52 val true)) |
---|
923 | nil) |
---|
924 | |
---|
925 | (rt:deftest oct.jacobi-dn.1q |
---|
926 | (let* ((ek (elliptic-k #q.5)) |
---|
927 | (true (sqrt #q.5)) |
---|
928 | (val (jacobi-dn ek #q.5))) |
---|
929 | (check-accuracy 212 val true)) |
---|
930 | nil) |
---|
931 | |
---|
932 | (rt:deftest oct.carlson-rf.1d |
---|
933 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
934 | ;; = 1/4*beta(1/2,1/2) |
---|
935 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
936 | (let ((rf (carlson-rf 0d0 2d0 1d0)) |
---|
937 | (true 1.31102877714605990523241979494d0)) |
---|
938 | (check-accuracy 53 rf true)) |
---|
939 | nil) |
---|
940 | |
---|
941 | (rt:deftest oct.carlson-rf.1q |
---|
942 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
943 | (let ((rf (carlson-rf #q0 #q2 #q1)) |
---|
944 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0)) |
---|
945 | (check-accuracy 212 rf true)) |
---|
946 | nil) |
---|
947 | |
---|
948 | (rt:deftest oct.carlson-rd.1d |
---|
949 | ;; Rd(0,2,1) = 3*integrate(s^2/sqrt(1-s^4), s, 0 ,1) |
---|
950 | ;; = 3*beta(3/4,1/2)/4 |
---|
951 | ;; = 3*sqrt(%pi)*gamma(3/4)/gamma(1/4) |
---|
952 | (let ((rd (carlson-rd 0d0 2d0 1d0)) |
---|
953 | (true 1.7972103521033883d0)) |
---|
954 | (check-accuracy 51 rd true)) |
---|
955 | nil) |
---|
956 | |
---|
957 | (rt:deftest oct.carlson-rd.1q |
---|
958 | (let ((rd (carlson-rd #q0 #q2 #q1)) |
---|
959 | (true #q1.797210352103388311159883738420485817340818994823477337395512429419599q0)) |
---|
960 | (check-accuracy 212 rd true)) |
---|
961 | nil) |
---|
962 | |
---|
963 | ;; Test some of the contagion stuff. |
---|
964 | |
---|
965 | (rt:deftest oct.carlson-rf.contagion.1 |
---|
966 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
967 | ;; = 1/4*beta(1/2,1/2) |
---|
968 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
969 | (let ((rf (carlson-rf 0 2 1)) |
---|
970 | (true 1.31102877714605990523241979494d0)) |
---|
971 | (check-accuracy 23 rf true)) |
---|
972 | nil) |
---|
973 | |
---|
974 | (rt:deftest oct.carlson-rf.contagion.1d |
---|
975 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
976 | ;; = 1/4*beta(1/2,1/2) |
---|
977 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
978 | (let ((rf (carlson-rf 0d0 2 1)) |
---|
979 | (true 1.31102877714605990523241979494d0)) |
---|
980 | (check-accuracy 53 rf true)) |
---|
981 | nil) |
---|
982 | |
---|
983 | (rt:deftest oct.carlson-rf.contagion.2d |
---|
984 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
985 | ;; = 1/4*beta(1/2,1/2) |
---|
986 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
987 | (let ((rf (carlson-rf 0 2d0 1)) |
---|
988 | (true 1.31102877714605990523241979494d0)) |
---|
989 | (check-accuracy 53 rf true)) |
---|
990 | nil) |
---|
991 | |
---|
992 | (rt:deftest oct.carlson-rf.contagion.3d |
---|
993 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
994 | ;; = 1/4*beta(1/2,1/2) |
---|
995 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
996 | (let ((rf (carlson-rf 0 2 1d0)) |
---|
997 | (true 1.31102877714605990523241979494d0)) |
---|
998 | (check-accuracy 53 rf true)) |
---|
999 | nil) |
---|
1000 | |
---|
1001 | (rt:deftest oct.carlson-rf.contagion.1q |
---|
1002 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
1003 | ;; = 1/4*beta(1/2,1/2) |
---|
1004 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
1005 | (let ((rf (carlson-rf #q0q0 2 1)) |
---|
1006 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0)) |
---|
1007 | (check-accuracy 212 rf true)) |
---|
1008 | nil) |
---|
1009 | |
---|
1010 | (rt:deftest oct.carlson-rf.contagion.2q |
---|
1011 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
1012 | ;; = 1/4*beta(1/2,1/2) |
---|
1013 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
1014 | (let ((rf (carlson-rf 0 #q2q0 1)) |
---|
1015 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0)) |
---|
1016 | (check-accuracy 212 rf true)) |
---|
1017 | nil) |
---|
1018 | |
---|
1019 | (rt:deftest oct.carlson-rf.contagion.3q |
---|
1020 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
1021 | ;; = 1/4*beta(1/2,1/2) |
---|
1022 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
1023 | (let ((rf (carlson-rf 0 2 #q1q0)) |
---|
1024 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0)) |
---|
1025 | (check-accuracy 212 rf true)) |
---|
1026 | nil) |
---|
1027 | |
---|
1028 | ;; Elliptic integral of the third kind |
---|
1029 | |
---|
1030 | ;; elliptic-pi(0,phi,m) = elliptic-f(phi, m) |
---|
1031 | (rt:deftest oct.elliptic-pi.1d |
---|
1032 | (loop for k from 0 to 100 |
---|
1033 | for phi = (random (/ pi 2)) |
---|
1034 | for m = (random 1d0) |
---|
1035 | for epi = (elliptic-pi 0 phi m) |
---|
1036 | for ef = (elliptic-f phi m) |
---|
1037 | for result = (check-accuracy 48 epi ef) |
---|
1038 | unless (eq nil result) |
---|
1039 | append (list (list phi m) result)) |
---|
1040 | nil) |
---|
1041 | |
---|
1042 | (rt:deftest oct.elliptic-pi.1q |
---|
1043 | (loop for k from 0 below 100 |
---|
1044 | for phi = (random (/ +pi+ 2)) |
---|
1045 | for m = (random #q1) |
---|
1046 | for epi = (elliptic-pi 0 phi m) |
---|
1047 | for ef = (elliptic-f phi m) |
---|
1048 | for result = (check-accuracy 53 epi ef) |
---|
1049 | unless (eq nil result) |
---|
1050 | append (list (list phi m) result)) |
---|
1051 | nil) |
---|
1052 | |
---|
1053 | ;; DLMF 19.6.3 |
---|
1054 | ;; |
---|
1055 | ;; PI(n; pi/2 | 0) = pi/(2*sqrt(1-n)) |
---|
1056 | (rt:deftest oct.elliptic-pi.19.6.3.d |
---|
1057 | (loop for k from 0 below 100 |
---|
1058 | for n = (random 1d0) |
---|
1059 | for epi = (elliptic-pi n (/ pi 2) 0) |
---|
1060 | for true = (/ pi (* 2 (sqrt (- 1 n)))) |
---|
1061 | for result = (check-accuracy 47 epi true) |
---|
1062 | unless (eq nil result) |
---|
1063 | append (list (list (list k n) result))) |
---|
1064 | nil) |
---|
1065 | |
---|
1066 | (rt:deftest oct.elliptic-pi.19.6.3.q |
---|
1067 | (loop for k from 0 below 100 |
---|
1068 | for n = (random #q1) |
---|
1069 | for epi = (elliptic-pi n (/ (float-pi n) 2) 0) |
---|
1070 | for true = (/ (float-pi n) (* 2 (sqrt (- 1 n)))) |
---|
1071 | for result = (check-accuracy 208 epi true) |
---|
1072 | unless (eq nil result) |
---|
1073 | append (list (list (list k n) result))) |
---|
1074 | nil) |
---|
1075 | |
---|
1076 | ;; elliptic-pi(n, phi, 0) = |
---|
1077 | ;; atan(sqrt(1-n)*tan(phi))/sqrt(1-n) n < 1 |
---|
1078 | ;; atanh(sqrt(n-1)*tan(phi))/sqrt(n-1) n > 1 |
---|
1079 | ;; tan(phi) n = 1 |
---|
1080 | ;; |
---|
1081 | ;; These are easy to derive if you look at the integral: |
---|
1082 | ;; |
---|
1083 | ;; ellipti-pi(n, phi, 0) = integrate(1/(1-n*sin(t)^2), t, 0, phi) |
---|
1084 | ;; |
---|
1085 | ;; and this can be easily integrated to give the above expressions for |
---|
1086 | ;; the different values of n. |
---|
1087 | (rt:deftest oct.elliptic-pi.n0.d |
---|
1088 | ;; Tests for random values for phi in [0, pi/2] and n in [0, 1] |
---|
1089 | (loop for k from 0 below 100 |
---|
1090 | for phi = (random (/ pi 2)) |
---|
1091 | for n = (random 1d0) |
---|
1092 | for epi = (elliptic-pi n phi 0) |
---|
1093 | for true = (/ (atan (* (tan phi) (sqrt (- 1 n)))) |
---|
1094 | (sqrt (- 1 n))) |
---|
1095 | for result = (check-accuracy 46.5 epi true) |
---|
1096 | unless (eq nil result) |
---|
1097 | append (list (list (list k n phi) result))) |
---|
1098 | nil) |
---|
1099 | |
---|
1100 | (rt:deftest oct.elliptic-pi.n1.d |
---|
1101 | (loop for k from 0 below 100 |
---|
1102 | for phi = (random (/ pi 2)) |
---|
1103 | for epi = (elliptic-pi 1 phi 0) |
---|
1104 | for true = (tan phi) |
---|
1105 | for result = (check-accuracy 34.5 epi true) |
---|
1106 | unless (eq nil result) |
---|
1107 | append (list (list (list k phi) result))) |
---|
1108 | nil) |
---|
1109 | |
---|
1110 | (rt:deftest oct.elliptic-pi.n2.d |
---|
1111 | (loop for k from 0 below 100 |
---|
1112 | for phi = (random (/ pi 2)) |
---|
1113 | for n = (+ 1d0 (random 100d0)) |
---|
1114 | for epi = (elliptic-pi n phi 0) |
---|
1115 | for true = (/ (atanh (* (tan phi) (sqrt (- n 1)))) |
---|
1116 | (sqrt (- n 1))) |
---|
1117 | for result = (check-accuracy 45.85 epi true) |
---|
1118 | ;; Not sure if this formula holds when atanh gives a complex |
---|
1119 | ;; result. Wolfram doesn't say |
---|
1120 | when (and (not (complexp true)) result) |
---|
1121 | append (list (list (list k n phi) result))) |
---|
1122 | nil) |
---|
1123 | |
---|
1124 | ;; Failed test case: |
---|
1125 | ;; ((89 66.68551748022054d0 0.12266024127708153d0) |
---|
1126 | ;; (45.868614757480834d0 47 0.47787458521306514d0 |
---|
1127 | ;; 0.4778745852130726d0)) |
---|
1128 | ;; New threshold is 45.85 bits. |
---|
1129 | (rt:deftest oct.elliptic-pi.n2.d-1 |
---|
1130 | (let* ((n 66.68551748022054d0) |
---|
1131 | (phi 0.12266024127708153d0) |
---|
1132 | (epi (elliptic-pi n phi 0)) |
---|
1133 | (true (/ (atanh (* (tan phi) (sqrt (- n 1)))) |
---|
1134 | (sqrt (- n 1))))) |
---|
1135 | (check-accuracy 45.8686d0 epi true)) |
---|
1136 | nil) |
---|
1137 | |
---|
1138 | |
---|
1139 | (rt:deftest oct.elliptic-pi.n0.q |
---|
1140 | ;; Tests for random values for phi in [0, pi/2] and n in [0, 1] |
---|
1141 | (loop for k from 0 below 100 |
---|
1142 | for phi = (random (/ +pi+ 2)) |
---|
1143 | for n = (random #q1) |
---|
1144 | for epi = (elliptic-pi n phi 0) |
---|
1145 | for true = (/ (atan (* (tan phi) (sqrt (- 1 n)))) |
---|
1146 | (sqrt (- 1 n))) |
---|
1147 | for result = (check-accuracy 204 epi true) |
---|
1148 | unless (eq nil result) |
---|
1149 | append (list (list (list k n phi) result))) |
---|
1150 | nil) |
---|
1151 | |
---|
1152 | (rt:deftest oct.elliptic-pi.n1.q |
---|
1153 | (loop for k from 0 below 100 |
---|
1154 | for phi = (random (/ +pi+ 2)) |
---|
1155 | for epi = (elliptic-pi 1 phi 0) |
---|
1156 | for true = (tan phi) |
---|
1157 | for result = (check-accuracy 194 epi true) |
---|
1158 | unless (eq nil result) |
---|
1159 | append (list (list (list k phi) result))) |
---|
1160 | nil) |
---|
1161 | |
---|
1162 | (rt:deftest oct.elliptic-pi.n2.q |
---|
1163 | (loop for k from 0 below 100 |
---|
1164 | for phi = (random (/ +pi+ 2)) |
---|
1165 | for n = (+ #q1 (random #q1)) |
---|
1166 | for epi = (elliptic-pi n phi 0) |
---|
1167 | for true = (/ (atanh (* (tan phi) (sqrt (- n 1)))) |
---|
1168 | (sqrt (- n 1))) |
---|
1169 | for result = (check-accuracy 202 epi true) |
---|
1170 | ;; Not sure if this formula holds when atanh gives a complex |
---|
1171 | ;; result. Wolfram doesn't say |
---|
1172 | when (and (not (complexp true)) result) |
---|
1173 | append (list (list (list k n phi) result))) |
---|
1174 | nil) |
---|
1175 | |
---|
1176 | ;; Tests for theta functions. |
---|
1177 | |
---|
1178 | (rt:deftest oct.theta3.1.d |
---|
1179 | ;; A&S 16.38.5 |
---|
1180 | ;; sqrt(2*K/%pi) = theta3(0,q) |
---|
1181 | (loop for k from 0 below 100 |
---|
1182 | for m = (random 1d0) |
---|
1183 | for t3 = (elliptic-theta-3 0 (elliptic-nome m)) |
---|
1184 | for true = (sqrt (/ (* 2 (elliptic-k m)) (float-pi m))) |
---|
1185 | for result = (check-accuracy 50.5 t3 true) |
---|
1186 | when result |
---|
1187 | append (list (list (list k m) result))) |
---|
1188 | nil) |
---|
1189 | |
---|
1190 | (rt:deftest oct.theta3.1.q |
---|
1191 | ;; A&S 16.38.5 |
---|
1192 | ;; sqrt(2*K/%pi) = theta3(0,q) |
---|
1193 | (loop for k from 0 below 100 |
---|
1194 | for m = (random #q1) |
---|
1195 | for t3 = (elliptic-theta-3 0 (elliptic-nome m)) |
---|
1196 | for true = (sqrt (/ (* 2 (elliptic-k m)) (float-pi m))) |
---|
1197 | for result = (check-accuracy 205.7 t3 true) |
---|
1198 | when result |
---|
1199 | append (list (list (list k m) result))) |
---|
1200 | nil) |
---|
1201 | |
---|
1202 | (rt:deftest oct.theta2.1.d |
---|
1203 | ;; A&S 16.38.7 |
---|
1204 | ;; sqrt(2*sqrt(m)*K/%pi) = theta2(0,q) |
---|
1205 | (loop for k from 0 below 100 |
---|
1206 | for m = (random 1d0) |
---|
1207 | for t3 = (elliptic-theta-2 0 (elliptic-nome m)) |
---|
1208 | for true = (sqrt (/ (* 2 (sqrt m) (elliptic-k m)) (float-pi m))) |
---|
1209 | for result = (check-accuracy 43.5 t3 true) |
---|
1210 | when result |
---|
1211 | append (list (list (list k m) result))) |
---|
1212 | nil) |
---|
1213 | |
---|
1214 | (rt:deftest oct.theta2.1.q |
---|
1215 | ;; A&S 16.38.7 |
---|
1216 | ;; sqrt(2*sqrt(m)*K/%pi) = theta2(0,q) |
---|
1217 | (loop for k from 0 below 100 |
---|
1218 | for m = (random #q1) |
---|
1219 | for t3 = (elliptic-theta-2 0 (elliptic-nome m)) |
---|
1220 | for true = (sqrt (/ (* 2 (sqrt m) (elliptic-k m)) (float-pi m))) |
---|
1221 | for result = (check-accuracy 205 t3 true) |
---|
1222 | when result |
---|
1223 | append (list (list (list k m) result))) |
---|
1224 | nil) |
---|
1225 | |
---|
1226 | (rt:deftest oct.theta4.1.d |
---|
1227 | ;; A&S 16.38.8 |
---|
1228 | ;; sqrt(2*sqrt(1-m)*K/%pi) = theta2(0,q) |
---|
1229 | (loop for k from 0 below 100 |
---|
1230 | for m = (random 1d0) |
---|
1231 | for t3 = (elliptic-theta-4 0 (elliptic-nome m)) |
---|
1232 | for true = (sqrt (/ (* 2 (sqrt (- 1 m)) (elliptic-k m)) |
---|
1233 | (float-pi m))) |
---|
1234 | for result = (check-accuracy 49 t3 true) |
---|
1235 | when result |
---|
1236 | append (list (list (list k m) result))) |
---|
1237 | nil) |
---|
1238 | |
---|
1239 | (rt:deftest oct.theta4.1.q |
---|
1240 | ;; A&S 16.38.8 |
---|
1241 | ;; sqrt(2*sqrt(1-m)*K/%pi) = theta2(0,q) |
---|
1242 | (loop for k from 0 below 100 |
---|
1243 | for m = (random #q1) |
---|
1244 | for t3 = (elliptic-theta-4 0 (elliptic-nome m)) |
---|
1245 | for true = (sqrt (/ (* 2 (sqrt (- 1 m)) (elliptic-k m)) |
---|
1246 | (float-pi m))) |
---|
1247 | for result = (check-accuracy 204 t3 true) |
---|
1248 | when result |
---|
1249 | append (list (list (list k m) result))) |
---|
1250 | nil) |
---|
1251 | |
---|
1252 | (rt:deftest lentz |
---|
1253 | ;; This isn't really a test of cf-incomplete-gamma. It's a test |
---|
1254 | ;; that Lentz's algorithm works in this case. For these args, |
---|
1255 | ;; cf-incomplete-gamma used to generate an overflow or division by |
---|
1256 | ;; zero because value-or-tiny was too tiny. |
---|
1257 | (let ((g (cf-incomplete-gamma 3d0 5d0)) |
---|
1258 | (true (- 2 (* 37 (exp -5d0))))) |
---|
1259 | (check-accuracy 51.2 g true)) |
---|
1260 | nil) |
---|
1261 | |
---|
1262 | (rt:deftest gamma.1.d |
---|
1263 | (let ((g (gamma 0.5d0)) |
---|
1264 | (true (sqrt pi))) |
---|
1265 | ;; This should give full accuracy but doesn't. |
---|
1266 | (check-accuracy 51 g true)) |
---|
1267 | nil) |
---|
1268 | |
---|
1269 | (rt:deftest gamma.1.q |
---|
1270 | (let ((g (gamma #q0.5)) |
---|
1271 | (true (sqrt +pi+))) |
---|
1272 | ;; This should give full accuracy but doesn't. |
---|
1273 | (check-accuracy 197 g true)) |
---|
1274 | nil) |
---|
1275 | |
---|
1276 | (rt:deftest gamma.2.d |
---|
1277 | (loop for k from 0 below 100 |
---|
1278 | for y = (+ 1 (random 100d0)) |
---|
1279 | for g = (abs (gamma (complex 0 y))) |
---|
1280 | for true = (sqrt (/ pi y (sinh (* pi y)))) |
---|
1281 | for result = (check-accuracy 44 g true) |
---|
1282 | when result |
---|
1283 | append (list (list (list k y) result))) |
---|
1284 | nil) |
---|
1285 | |
---|
1286 | (rt:deftest gamma.2.q |
---|
1287 | (loop for k from 0 below 100 |
---|
1288 | for y = (+ 1 (random #q100)) |
---|
1289 | for g = (abs (gamma (complex 0 y))) |
---|
1290 | for true = (sqrt (/ +pi+ y (sinh (* +pi+ y)))) |
---|
1291 | for result = (check-accuracy 196 g true) |
---|
1292 | when result |
---|
1293 | append (list (list (list k y) result))) |
---|
1294 | nil) |
---|
1295 | |
---|
1296 | (rt:deftest gamma.3.d |
---|
1297 | (loop for k from 0 below 100 |
---|
1298 | for y = (+ 1 (random 100d0)) |
---|
1299 | for g = (abs (gamma (complex 1/2 y))) |
---|
1300 | for true = (sqrt (/ pi (cosh (* pi y)))) |
---|
1301 | for result = (check-accuracy 44 g true) |
---|
1302 | when result |
---|
1303 | append (list (list (list k y) result))) |
---|
1304 | nil) |
---|
1305 | |
---|
1306 | (rt:deftest gamma.3.q |
---|
1307 | (loop for k from 0 below 100 |
---|
1308 | for y = (+ 1 (random #q100)) |
---|
1309 | for g = (abs (gamma (complex 1/2 y))) |
---|
1310 | for true = (sqrt (/ +pi+ (cosh (* +pi+ y)))) |
---|
1311 | for result = (check-accuracy 196 g true) |
---|
1312 | when result |
---|
1313 | append (list (list (list k y) result))) |
---|
1314 | nil) |
---|
1315 | |
---|
1316 | ;; gamma_incomplete(2,z) = integrate(t*exp(-t), t, z, inf) |
---|
1317 | ;; = (z+1)*exp(-z) |
---|
1318 | (rt:deftest gamma-incomplete-tail.1.d |
---|
1319 | (let* ((z 5d0) |
---|
1320 | (gi (incomplete-gamma-tail 2 z)) |
---|
1321 | (true (* (+ z 1) (exp (- z))))) |
---|
1322 | (check-accuracy 52 gi true)) |
---|
1323 | nil) |
---|
1324 | |
---|
1325 | (rt:deftest gamma-incomplete-tail.2.d |
---|
1326 | (let* ((z #c(1 5d0)) |
---|
1327 | (gi (incomplete-gamma-tail 2 z)) |
---|
1328 | (true (* (+ z 1) (exp (- z))))) |
---|
1329 | (check-accuracy 50 gi true)) |
---|
1330 | nil) |
---|
1331 | |
---|
1332 | (rt:deftest gamma-incomplete-tail.1.q |
---|
1333 | (let* ((z #q5) |
---|
1334 | (gi (incomplete-gamma-tail 2 z)) |
---|
1335 | (true (* (+ z 1) (exp (- z))))) |
---|
1336 | (check-accuracy 207 gi true)) |
---|
1337 | nil) |
---|
1338 | |
---|
1339 | (rt:deftest gamma-incomplete-tail.2.q |
---|
1340 | (let* ((z #q(1 5)) |
---|
1341 | (gi (incomplete-gamma-tail 2 z)) |
---|
1342 | (true (* (+ z 1) (exp (- z))))) |
---|
1343 | (check-accuracy 206 gi true)) |
---|
1344 | nil) |
---|
1345 | |
---|
1346 | (rt:deftest gamma-incomplete-tail.3.d |
---|
1347 | (let* ((z -5d0) |
---|
1348 | (gi (incomplete-gamma-tail 2 z)) |
---|
1349 | (true (* (+ z 1) (exp (- z))))) |
---|
1350 | (check-accuracy 50 gi true)) |
---|
1351 | nil) |
---|
1352 | |
---|
1353 | (rt:deftest gamma-incomplete-tail.3.q |
---|
1354 | (let* ((z #q-5) |
---|
1355 | (gi (incomplete-gamma-tail 2 z)) |
---|
1356 | (true (* (+ z 1) (exp (- z))))) |
---|
1357 | (check-accuracy 206 gi true)) |
---|
1358 | nil) |
---|
1359 | |
---|
1360 | ;; See http://www.wolframalpha.com/input/?i=Gamma[1%2F2%2C-100%2Bi%2F%2810^10%29] |
---|
1361 | |
---|
1362 | (rt:deftest gamma-incomplete-tail.4.q |
---|
1363 | (let* ((z #q(#q-100 #q1q-10)) |
---|
1364 | (gi (incomplete-gamma-tail 1/2 z)) |
---|
1365 | (true #q(#q-2.68811714181613544840818982228135651231579313476267430888499241497530341422025007816745898370049200133136q32 |
---|
1366 | #q-2.70176456134384383878883307528351227886457379834795655467745609829086928772079968479767583764284583465328q42))) |
---|
1367 | (check-accuracy 205 gi true)) |
---|
1368 | nil) |
---|
1369 | |
---|
1370 | |
---|
1371 | ;; Fresnel integrals. |
---|
1372 | ;; |
---|
1373 | ;; For x small, Fresnel |
---|
1374 | ;; |
---|
1375 | ;; S(z) = %pi/6*z^3*(1 - %pi^2*z^4/56 + %pi^4*z^8/2040 - ...) |
---|
1376 | ;; |
---|
1377 | (defun fresnel-s-series (z) |
---|
1378 | (let* ((fpi (float-pi z)) |
---|
1379 | (z^3 (expt z 3)) |
---|
1380 | (z^4 (* z^3 z))) |
---|
1381 | (* fpi 1/6 z^3 |
---|
1382 | (+ 1 (/ (* fpi fpi z^4) |
---|
1383 | -56) |
---|
1384 | (/ (* (expt fpi 4) (expt z^4 2)) |
---|
1385 | 7040))))) |
---|
1386 | |
---|
1387 | (rt:deftest fresnel-s.1d |
---|
1388 | (let* ((z 1d-3) |
---|
1389 | (s (fresnel-s z)) |
---|
1390 | (true (fresnel-s-series z))) |
---|
1391 | (check-accuracy 52 s true)) |
---|
1392 | nil) |
---|
1393 | |
---|
1394 | (rt:deftest fresnel-s.2d |
---|
1395 | (let* ((z #c(1d-3 1d-3)) |
---|
1396 | (s (fresnel-s z)) |
---|
1397 | (true (fresnel-s-series z))) |
---|
1398 | (check-accuracy 52 s true)) |
---|
1399 | nil) |
---|
1400 | |
---|
1401 | (rt:deftest fresnel-s.1q |
---|
1402 | (let* ((z #q1q-20) |
---|
1403 | (s (fresnel-s z)) |
---|
1404 | (true (fresnel-s-series z))) |
---|
1405 | (check-accuracy 212 s true)) |
---|
1406 | nil) |
---|
1407 | |
---|
1408 | (rt:deftest fresnel-s.2q |
---|
1409 | (let* ((z #q(#q1q-3 #q1q-3)) |
---|
1410 | (s (fresnel-s z)) |
---|
1411 | (true (fresnel-s-series z))) |
---|
1412 | (check-accuracy 212 s true)) |
---|
1413 | nil) |
---|
1414 | |
---|
1415 | (rt:deftest psi.1d |
---|
1416 | (let* ((z 1d0) |
---|
1417 | (p (psi z)) |
---|
1418 | (true (float (- +%gamma+) 1d0))) |
---|
1419 | (check-accuracy 52 p true)) |
---|
1420 | nil) |
---|
1421 | |
---|
1422 | (rt:deftest psi.1q |
---|
1423 | (let* ((z #q1) |
---|
1424 | (p (psi z)) |
---|
1425 | (true (- +%gamma+))) |
---|
1426 | (check-accuracy 208 p true)) |
---|
1427 | nil) |
---|
1428 | |
---|
1429 | (rt:deftest psi.2d |
---|
1430 | (let* ((z (float 4/3 1d0)) |
---|
1431 | (p (psi z)) |
---|
1432 | (true (- 3 |
---|
1433 | +%gamma+ |
---|
1434 | (/ +pi+ (* 2 (sqrt #q3))) |
---|
1435 | (* 1.5 (log #q3))))) |
---|
1436 | (check-accuracy 49.8 p true)) |
---|
1437 | nil) |
---|
1438 | |
---|
1439 | (rt:deftest psi.2q |
---|
1440 | (let* ((z (float 4/3 #q1)) |
---|
1441 | (p (psi z)) |
---|
1442 | (true (- 3 |
---|
1443 | +%gamma+ |
---|
1444 | (/ +pi+ (* 2 (sqrt #q3))) |
---|
1445 | (* 1.5 (log #q3))))) |
---|
1446 | (check-accuracy 205 p true)) |
---|
1447 | nil) |
---|
1448 | |
---|
1449 | (rt:deftest psi.3d |
---|
1450 | (let* ((z (float -1/2 1d0)) |
---|
1451 | (p (psi z)) |
---|
1452 | (true (- 2 |
---|
1453 | +%gamma+ |
---|
1454 | (log #q4)))) |
---|
1455 | (check-accuracy 48 p true)) |
---|
1456 | nil) |
---|
1457 | |
---|
1458 | (rt:deftest psi.3q |
---|
1459 | (let* ((z (float -1/2 #q1)) |
---|
1460 | (p (psi z)) |
---|
1461 | (true (- 2 |
---|
1462 | +%gamma+ |
---|
1463 | (log #q4)))) |
---|
1464 | (check-accuracy 204.1 p true)) |
---|
1465 | nil) |
---|
1466 | |
---|
1467 | (rt:deftest expintegral-e.1d |
---|
1468 | (let* ((z 1d0) |
---|
1469 | (e (exp-integral-e 0 z)) |
---|
1470 | (true (/ (exp (- z)) z))) |
---|
1471 | (check-accuracy 53 e true)) |
---|
1472 | nil) |
---|
1473 | |
---|
1474 | (rt:deftest expintegral-e.1q |
---|
1475 | (let* ((z #q1) |
---|
1476 | (e (exp-integral-e 0 z)) |
---|
1477 | (true (/ (exp (- z)) z))) |
---|
1478 | (check-accuracy 212 e true)) |
---|
1479 | nil) |
---|
1480 | |
---|
1481 | (rt:deftest expintegral-e.2d |
---|
1482 | (let* ((z 15d0) |
---|
1483 | (e (exp-integral-e 0 z)) |
---|
1484 | (true (/ (exp (- z)) z))) |
---|
1485 | (check-accuracy 53 e true)) |
---|
1486 | nil) |
---|
1487 | |
---|
1488 | (rt:deftest expintegral-e.2q |
---|
1489 | (let* ((z #q15) |
---|
1490 | (e (exp-integral-e 0 z)) |
---|
1491 | (true (/ (exp (- z)) z))) |
---|
1492 | (check-accuracy 212 e true)) |
---|
1493 | nil) |
---|
1494 | |
---|
1495 | (rt:deftest expintegral-e.3d |
---|
1496 | (let* ((e (exp-integral-e 2 1d0)) |
---|
1497 | (true 0.14849550677592204791835999d0)) |
---|
1498 | (check-accuracy 47.5 e true)) |
---|
1499 | nil) |
---|
1500 | |
---|
1501 | (rt:deftest expintegral-e.4d |
---|
1502 | (let* ((x .5d0) |
---|
1503 | (e (exp-integral-e -2 x)) |
---|
1504 | (true (/ (* (exp (- x)) (+ (* x x x) (* 2 x x) (* 2 x))) |
---|
1505 | (expt x 4)))) |
---|
1506 | (check-accuracy 53 e true)) |
---|
1507 | nil) |
---|
1508 | |
---|
1509 | (rt:deftest expintegral-e.4q |
---|
1510 | (let* ((x #q.5) |
---|
1511 | (e (exp-integral-e -2 x)) |
---|
1512 | (true (/ (* (exp (- x)) (+ (* x x x) (* 2 x x) (* 2 x))) |
---|
1513 | (expt x 4)))) |
---|
1514 | (check-accuracy 210.8 e true)) |
---|
1515 | nil) |
---|
1516 | |
---|
1517 | (rt:deftest expintegral-e.5d |
---|
1518 | (let* ((x .5d0) |
---|
1519 | (e (exp-integral-e 2d0 x)) |
---|
1520 | (true #q0.3266438623245530177304015653336378358284946903290101)) |
---|
1521 | (check-accuracy 51.2 e true)) |
---|
1522 | nil) |
---|
1523 | |
---|
1524 | (rt:deftest expintegral-e.5q |
---|
1525 | (let* ((x #q.5) |
---|
1526 | (e (exp-integral-e #q2 x)) |
---|
1527 | (true #q0.326643862324553017730401565333637835828494690329010198058745549181386569998611289568)) |
---|
1528 | (check-accuracy 208.4 e true)) |
---|
1529 | nil) |
---|
1530 | |
---|
1531 | (rt:deftest expintegral-e.6d |
---|
1532 | (let* ((x .5d0) |
---|
1533 | (e (exp-integral-e 1d0 x)) |
---|
1534 | (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883)) |
---|
1535 | (check-accuracy 53.9 e true)) |
---|
1536 | nil) |
---|
1537 | |
---|
1538 | (rt:deftest expintegral-e.6q |
---|
1539 | (let* ((x #q.5) |
---|
1540 | (e (exp-integral-e #q1 x)) |
---|
1541 | (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883)) |
---|
1542 | (check-accuracy 219.1 e true)) |
---|
1543 | nil) |
---|
1544 | |
---|
1545 | |
---|
1546 | ;; Bessel J tests for negative order |
---|
1547 | (rt:deftest bessel-j.neg-order.d.1 |
---|
1548 | (let ((b (bessel-j -1d0 2d0)) |
---|
1549 | (true -0.5767248077568734d0)) |
---|
1550 | (check-accuracy 50.2 b true)) |
---|
1551 | nil) |
---|
1552 | |
---|
1553 | (rt:deftest bessel-j.neg-order.d.2 |
---|
1554 | (let ((b (bessel-j -1d0 1.5d0)) |
---|
1555 | (true -0.5579365079100996d0)) |
---|
1556 | (check-accuracy 50.5 b true)) |
---|
1557 | nil) |
---|
1558 | |
---|
1559 | (rt:deftest bessel-j.neg-order.d.3 |
---|
1560 | (let ((b (bessel-j -1.5d0 2d0)) |
---|
1561 | (true -0.3956232813587035d0)) |
---|
1562 | (check-accuracy 50.59 b true)) |
---|
1563 | nil) |
---|
1564 | |
---|
1565 | (rt:deftest bessel-j.neg-order.d.4 |
---|
1566 | (let ((b (bessel-j -1.8d0 1.5d0)) |
---|
1567 | (true -0.251327217627129314d0)) |
---|
1568 | (check-accuracy 49.98 b true)) |
---|
1569 | nil) |
---|
1570 | |
---|
1571 | (rt:deftest bessel-j.neg-order.d.5 |
---|
1572 | (let ((b (bessel-j -2d0 1.5d0)) |
---|
1573 | (true 0.2320876721442147d0)) |
---|
1574 | (check-accuracy 51.89 b true)) |
---|
1575 | nil) |
---|
1576 | |
---|
1577 | (rt:deftest bessel-j.neg-order.d.6 |
---|
1578 | (let ((b (bessel-j -2.5d0 1.5d0)) |
---|
1579 | (true 1.315037204805194d0)) |
---|
1580 | (check-accuracy 52.37 b true)) |
---|
1581 | nil) |
---|
1582 | |
---|
1583 | (rt:deftest bessel-j.neg-order.d.7 |
---|
1584 | (let ((b (bessel-j -2.3d0 1.5d0)) |
---|
1585 | (true 1.012178926325313d0)) |
---|
1586 | (check-accuracy 50.01 b true)) |
---|
1587 | nil) |
---|
1588 | |
---|
1589 | ;; Bessel-J tests for positive order |
---|
1590 | (rt:deftest bessel-j.pos-order.d.1 |
---|
1591 | (let ((b (bessel-j 1.5d0 1d0)) |
---|
1592 | (true 0.2402978391234270d0)) |
---|
1593 | (check-accuracy 51.83 b true)) |
---|
1594 | nil) |
---|
1595 | |
---|
1596 | (rt:deftest bessel-j.pos-order.d.2 |
---|
1597 | (let ((b (bessel-j 1.8d0 1d0)) |
---|
1598 | (true 0.1564953153109239d0)) |
---|
1599 | (check-accuracy 51.97 b true)) |
---|
1600 | nil) |
---|
1601 | |
---|
1602 | (rt:deftest bessel-j.pos-order.d.3 |
---|
1603 | (let ((b (bessel-j 2d0 1d0)) |
---|
1604 | (true 0.1149034849319005d0)) |
---|
1605 | (check-accuracy 51.87 b true)) |
---|
1606 | nil) |
---|
1607 | |
---|
1608 | (rt:deftest bessel-j.pos-order.d.4 |
---|
1609 | (let ((b (bessel-j 2.5d0 1d0)) |
---|
1610 | (true 0.04949681022847794d0)) |
---|
1611 | (check-accuracy 47.17 b true)) |
---|
1612 | nil) |
---|
1613 | |
---|
1614 | (rt:deftest bessel-j.pos-order.d.5 |
---|
1615 | (let ((b (bessel-j -2d0 1.5d0)) |
---|
1616 | (true 0.2320876721442147d0)) |
---|
1617 | (check-accuracy 51.89 b true)) |
---|
1618 | nil) |
---|
1619 | |
---|
1620 | ;; Bessel J for half integer order and real args |
---|
1621 | (rt:deftest bessel-j-1/2.d.1 |
---|
1622 | (loop for k from 0 below 100 |
---|
1623 | ;; x in [1,1+pi/2] because we don't want to test the Bessel |
---|
1624 | ;; series and we don't want to test near pi because sin(pi) |
---|
1625 | ;; = 0, where we will lose accuracy. |
---|
1626 | for x = (+ 1 (random (/ pi 2))) |
---|
1627 | for b = (bessel-j 0.5d0 x) |
---|
1628 | for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi))) |
---|
1629 | for result = (check-accuracy 48.42 b true) |
---|
1630 | when result |
---|
1631 | append (list (list (list k x) result))) |
---|
1632 | nil) |
---|
1633 | |
---|
1634 | (rt:deftest bessel-j-1/2.d.1.a |
---|
1635 | (let* ((x 2.3831631289164497d0) |
---|
1636 | (b (bessel-j 0.5d0 x)) |
---|
1637 | (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi))))) |
---|
1638 | (check-accuracy 48.42 b true)) |
---|
1639 | nil) |
---|
1640 | |
---|
1641 | (rt:deftest bessel-j-1/2.q.1 |
---|
1642 | (loop for k from 0 below 10 |
---|
1643 | ;; x in [1,1+pi/2] because we don't want to test the Bessel |
---|
1644 | ;; series and we don't want to test near pi because sin(pi) |
---|
1645 | ;; = 0, where we will lose accuracy. |
---|
1646 | for x = (+ 1 (random (/ (float-pi #q1) 2))) |
---|
1647 | for b = (bessel-j #q0.5 x) |
---|
1648 | for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))) |
---|
1649 | for result = (check-accuracy 173.28 b true) |
---|
1650 | when result |
---|
1651 | append (list (list (list k x) result))) |
---|
1652 | nil) |
---|
1653 | |
---|
1654 | (rt:deftest bessel-j-1/2.q.1.a |
---|
1655 | (let* ((x #q1.1288834862545916200627583005758663687705443417892789067029865493882q0) |
---|
1656 | (b (bessel-j #q0.5 x)) |
---|
1657 | (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))))) |
---|
1658 | (check-accuracy 182.92 b true)) |
---|
1659 | nil) |
---|
1660 | |
---|
1661 | (rt:deftest bessel-j-1/2.q.1.b |
---|
1662 | (let* ((x #q1.1288834862545916200627583005758663687705443417892789067029865493882q0) |
---|
1663 | (b (bessel-j #q0.5 x)) |
---|
1664 | (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1)))))) |
---|
1665 | (check-accuracy 173.28 b true)) |
---|
1666 | nil) |
---|
1667 | |
---|
1668 | ;; Bessel J for complex args |
---|
1669 | #+nil |
---|
1670 | (rt:deftest bessel-j-complex.pos-order.d.1 |
---|
1671 | (let ((b (bessel-j 0d0 #c(1d0 1))) |
---|
1672 | (true #c(0.9376084768060293d0 -0.4965299476091221d0))) |
---|
1673 | (check-accuracy 53 b true)) |
---|
1674 | nil) |
---|
1675 | |
---|