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 | (defun check-accuracy (limit est true) |
---|
49 | (let ((bits (bit-accuracy est true))) |
---|
50 | (if (numberp bits) |
---|
51 | (if (< bits limit) |
---|
52 | (list bits limit est true))))) |
---|
53 | |
---|
54 | (defvar *null* (make-broadcast-stream)) |
---|
55 | |
---|
56 | ;;; Some simple tests from the Yozo Hida's qd package. |
---|
57 | |
---|
58 | ;; Pi via Machin's formula |
---|
59 | (rt:deftest oct.pi.machin |
---|
60 | (let* ((*standard-output* *null*) |
---|
61 | (val (make-instance 'qd-real :value (octi::test2 nil))) |
---|
62 | (true oct:+pi+)) |
---|
63 | (check-accuracy 213 val true)) |
---|
64 | nil) |
---|
65 | |
---|
66 | ;; Pi via Salamin-Brent algorithm |
---|
67 | (rt:deftest oct.pi.salamin-brent |
---|
68 | (let* ((*standard-output* *null*) |
---|
69 | (val (make-instance 'qd-real :value (octi::test3 nil))) |
---|
70 | (true oct:+pi+)) |
---|
71 | (check-accuracy 202 val true)) |
---|
72 | nil) |
---|
73 | |
---|
74 | ;; Pi via Borweign's Quartic formula |
---|
75 | (rt:deftest oct.pi.borweign |
---|
76 | (let* ((*standard-output* *null*) |
---|
77 | (val (make-instance 'qd-real :value (octi::test4 nil))) |
---|
78 | (true oct:+pi+)) |
---|
79 | (check-accuracy 211 val true)) |
---|
80 | nil) |
---|
81 | |
---|
82 | ;; e via Taylor series |
---|
83 | (rt:deftest oct.e.taylor |
---|
84 | (let* ((*standard-output* *null*) |
---|
85 | (val (make-instance 'qd-real :value (octi::test5 nil))) |
---|
86 | (true (make-instance 'qd-real :value octi::+qd-e+))) |
---|
87 | (check-accuracy 212 val true)) |
---|
88 | nil) |
---|
89 | |
---|
90 | ;; log(2) via Taylor series |
---|
91 | (rt:deftest oct.log2.taylor |
---|
92 | (let* ((*standard-output* *null*) |
---|
93 | (val (make-instance 'qd-real :value (octi::test6 nil))) |
---|
94 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
95 | (check-accuracy 212 val true)) |
---|
96 | nil) |
---|
97 | |
---|
98 | ;;; Tests of atan where we know the analytical result |
---|
99 | (rt:deftest oct.atan.1 |
---|
100 | (let* ((arg (/ (sqrt #q3))) |
---|
101 | (y (/ (atan arg) +pi+)) |
---|
102 | (true (/ #q6))) |
---|
103 | (check-accuracy 212 y true)) |
---|
104 | nil) |
---|
105 | |
---|
106 | (rt:deftest oct.atan.2 |
---|
107 | (let* ((arg (sqrt #q3)) |
---|
108 | (y (/ (atan arg) +pi+)) |
---|
109 | (true (/ #q3))) |
---|
110 | (check-accuracy 212 y true)) |
---|
111 | nil) |
---|
112 | |
---|
113 | (rt:deftest oct.atan.3 |
---|
114 | (let* ((arg #q1) |
---|
115 | (y (/ (atan arg) +pi+)) |
---|
116 | (true (/ #q4))) |
---|
117 | (check-accuracy 212 y true)) |
---|
118 | nil) |
---|
119 | |
---|
120 | (rt:deftest oct.atan.4 |
---|
121 | (let* ((arg #q1q100) |
---|
122 | (y (/ (atan arg) +pi+)) |
---|
123 | (true #q.5)) |
---|
124 | (check-accuracy 212 y true)) |
---|
125 | nil) |
---|
126 | |
---|
127 | (rt:deftest oct.atan.5 |
---|
128 | (let* ((arg #q-1q100) |
---|
129 | (y (/ (atan arg) +pi+)) |
---|
130 | (true #q-.5)) |
---|
131 | (check-accuracy 212 y true)) |
---|
132 | nil) |
---|
133 | |
---|
134 | |
---|
135 | (defun atan-qd/duplication (arg) |
---|
136 | (make-instance 'qd-real |
---|
137 | :value (octi::atan-qd/duplication (qd-value arg)))) |
---|
138 | |
---|
139 | ;;; Tests of atan where we know the analytical result. Same tests, |
---|
140 | ;;; but using the atan duplication formula. |
---|
141 | (rt:deftest oct.atan/dup.1 |
---|
142 | (let* ((arg (/ (sqrt #q3))) |
---|
143 | (y (/ (atan-qd/duplication arg) +pi+)) |
---|
144 | (true (/ #q6))) |
---|
145 | (check-accuracy 212 y true)) |
---|
146 | nil) |
---|
147 | |
---|
148 | (rt:deftest oct.atan/dup.2 |
---|
149 | (let* ((arg (sqrt #q3)) |
---|
150 | (y (/ (atan-qd/duplication arg) +pi+)) |
---|
151 | (true (/ #q3))) |
---|
152 | (check-accuracy 212 y true)) |
---|
153 | nil) |
---|
154 | |
---|
155 | (rt:deftest oct.atan/dup.3 |
---|
156 | (let* ((arg #q1) |
---|
157 | (y (/ (atan-qd/duplication arg) +pi+)) |
---|
158 | (true (/ #q4))) |
---|
159 | (check-accuracy 212 y true)) |
---|
160 | nil) |
---|
161 | |
---|
162 | (rt:deftest oct.atan/dup.4 |
---|
163 | (let* ((arg #q1q100) |
---|
164 | (y (/ (atan-qd/duplication arg) +pi+)) |
---|
165 | (true #q.5)) |
---|
166 | (check-accuracy 212 y true)) |
---|
167 | nil) |
---|
168 | |
---|
169 | (rt:deftest oct.atan/dup.5 |
---|
170 | (let* ((arg #q-1q100) |
---|
171 | (y (/ (atan-qd/duplication arg) +pi+)) |
---|
172 | (true #q-.5)) |
---|
173 | (check-accuracy 212 y true)) |
---|
174 | nil) |
---|
175 | |
---|
176 | ;;; Tests of atan where we know the analytical result. Same tests, |
---|
177 | ;;; but using a CORDIC implementation. |
---|
178 | (defun atan-qd/cordic (arg) |
---|
179 | (make-instance 'qd-real |
---|
180 | :value (octi::atan-qd/cordic (qd-value arg)))) |
---|
181 | |
---|
182 | (rt:deftest oct.atan/cordic.1 |
---|
183 | (let* ((arg (/ (sqrt #q3))) |
---|
184 | (y (/ (atan-qd/cordic arg) +pi+)) |
---|
185 | (true (/ #q6))) |
---|
186 | (check-accuracy 212 y true)) |
---|
187 | nil) |
---|
188 | |
---|
189 | (rt:deftest oct.atan/cordic.2 |
---|
190 | (let* ((arg (sqrt #q3)) |
---|
191 | (y (/ (atan-qd/cordic arg) +pi+)) |
---|
192 | (true (/ #q3))) |
---|
193 | (check-accuracy 212 y true)) |
---|
194 | nil) |
---|
195 | |
---|
196 | (rt:deftest oct.atan/cordic.3 |
---|
197 | (let* ((arg #q1) |
---|
198 | (y (/ (atan-qd/cordic arg) +pi+)) |
---|
199 | (true (/ #q4))) |
---|
200 | (check-accuracy 212 y true)) |
---|
201 | nil) |
---|
202 | |
---|
203 | (rt:deftest oct.atan/cordic.4 |
---|
204 | (let* ((arg #q1q100) |
---|
205 | (y (/ (atan-qd/cordic arg) +pi+)) |
---|
206 | (true #q.5)) |
---|
207 | (check-accuracy 212 y true)) |
---|
208 | nil) |
---|
209 | |
---|
210 | (rt:deftest oct.atan/cordic.5 |
---|
211 | (let* ((arg #q-1q100) |
---|
212 | (y (/ (atan-qd/cordic arg) +pi+)) |
---|
213 | (true #q-.5)) |
---|
214 | (check-accuracy 212 y true)) |
---|
215 | nil) |
---|
216 | |
---|
217 | |
---|
218 | ;;; Tests of sin where we know the analytical result. |
---|
219 | (rt:deftest oct.sin.1 |
---|
220 | (let* ((arg (/ +pi+ 6)) |
---|
221 | (y (sin arg)) |
---|
222 | (true #q.5)) |
---|
223 | (check-accuracy 212 y true)) |
---|
224 | nil) |
---|
225 | |
---|
226 | (rt:deftest oct.sin.2 |
---|
227 | (let* ((arg (/ +pi+ 4)) |
---|
228 | (y (sin arg)) |
---|
229 | (true (sqrt #q.5))) |
---|
230 | (check-accuracy 212 y true)) |
---|
231 | nil) |
---|
232 | |
---|
233 | (rt:deftest oct.sin.3 |
---|
234 | (let* ((arg (/ +pi+ 3)) |
---|
235 | (y (sin arg)) |
---|
236 | (true (/ (sqrt #q3) 2))) |
---|
237 | (check-accuracy 212 y true)) |
---|
238 | nil) |
---|
239 | |
---|
240 | (rt:deftest oct.big-sin.1 |
---|
241 | (let* ((arg (oct:make-qd (ash 1 120))) |
---|
242 | (y (sin arg)) |
---|
243 | (true #q3.778201093607520226555484700569229919605866976512306642257987199414885q-1)) |
---|
244 | (check-accuracy 205 y true)) |
---|
245 | nil) |
---|
246 | |
---|
247 | (rt:deftest oct.big-sin.2 |
---|
248 | (let* ((arg (oct:make-qd (ash 1 1023))) |
---|
249 | (y (sin arg)) |
---|
250 | (true #q5.631277798508840134529434079444683477103854907361251399182750155357133q-1)) |
---|
251 | (check-accuracy 205 y true)) |
---|
252 | nil) |
---|
253 | |
---|
254 | ;;; Tests of tan where we know the analytical result. |
---|
255 | (rt:deftest oct.tan.1 |
---|
256 | (let* ((arg (/ +pi+ 6)) |
---|
257 | (y (tan arg)) |
---|
258 | (true (/ (sqrt #q3)))) |
---|
259 | (check-accuracy 212 y true)) |
---|
260 | nil) |
---|
261 | |
---|
262 | (rt:deftest oct.tan.2 |
---|
263 | (let* ((arg (/ +pi+ 4)) |
---|
264 | (y (tan arg)) |
---|
265 | (true #q1)) |
---|
266 | (check-accuracy 212 y true)) |
---|
267 | nil) |
---|
268 | |
---|
269 | (rt:deftest oct.tan.3 |
---|
270 | (let* ((arg (/ +pi+ 3)) |
---|
271 | (y (tan arg)) |
---|
272 | (true (sqrt #q3))) |
---|
273 | (check-accuracy 212 y true)) |
---|
274 | nil) |
---|
275 | |
---|
276 | ;;; Tests of tan where we know the analytical result. Uses CORDIC |
---|
277 | ;;; algorithm. |
---|
278 | |
---|
279 | (defun tan/cordic (arg) |
---|
280 | (make-instance 'qd-real |
---|
281 | :value (octi::tan-qd/cordic (qd-value arg)))) |
---|
282 | |
---|
283 | (rt:deftest oct.tan/cordic.1 |
---|
284 | (let* ((arg (/ +pi+ 6)) |
---|
285 | (y (tan/cordic arg)) |
---|
286 | (true (/ (sqrt #q3)))) |
---|
287 | (check-accuracy 211 y true)) |
---|
288 | nil) |
---|
289 | |
---|
290 | (rt:deftest oct.tan/cordic.2 |
---|
291 | (let* ((arg (/ +pi+ 4)) |
---|
292 | (y (tan/cordic arg)) |
---|
293 | (true #q1)) |
---|
294 | (check-accuracy 211 y true)) |
---|
295 | nil) |
---|
296 | |
---|
297 | (rt:deftest oct.tan/cordic.3 |
---|
298 | (let* ((arg (/ +pi+ 3)) |
---|
299 | (y (tan/cordic arg)) |
---|
300 | (true (sqrt #q3))) |
---|
301 | (check-accuracy 210 y true)) |
---|
302 | nil) |
---|
303 | |
---|
304 | ;;; Tests of asin where we know the analytical result. |
---|
305 | |
---|
306 | (rt:deftest oct.asin.1 |
---|
307 | (let* ((arg #q.5) |
---|
308 | (y (asin arg)) |
---|
309 | (true (/ +pi+ 6))) |
---|
310 | (check-accuracy 212 y true)) |
---|
311 | nil) |
---|
312 | |
---|
313 | (rt:deftest oct.asin.2 |
---|
314 | (let* ((arg (sqrt #q.5)) |
---|
315 | (y (asin arg)) |
---|
316 | (true (/ +pi+ 4))) |
---|
317 | (check-accuracy 212 y true)) |
---|
318 | nil) |
---|
319 | |
---|
320 | (rt:deftest oct.asin.3 |
---|
321 | (let* ((arg (/ (sqrt #q3) 2)) |
---|
322 | (y (asin arg)) |
---|
323 | (true (/ +pi+ 3))) |
---|
324 | (check-accuracy 212 y true)) |
---|
325 | nil) |
---|
326 | |
---|
327 | ;;; Tests of log. |
---|
328 | |
---|
329 | (rt:deftest oct.log.1 |
---|
330 | (let* ((arg #q2) |
---|
331 | (y (log arg)) |
---|
332 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
333 | (check-accuracy 212 y true)) |
---|
334 | nil) |
---|
335 | |
---|
336 | (rt:deftest oct.log.2 |
---|
337 | (let* ((arg #q10) |
---|
338 | (y (log arg)) |
---|
339 | (true (make-instance 'qd-real :value octi::+qd-log10+))) |
---|
340 | (check-accuracy 207 y true)) |
---|
341 | nil) |
---|
342 | |
---|
343 | (rt:deftest oct.log.3 |
---|
344 | (let* ((arg (+ 1 (scale-float #q1 -80))) |
---|
345 | (y (log arg)) |
---|
346 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
347 | (check-accuracy 212 y true)) |
---|
348 | nil) |
---|
349 | |
---|
350 | ;;; Tests of log using Newton iteration. |
---|
351 | |
---|
352 | (defun log/newton (arg) |
---|
353 | (make-instance 'qd-real |
---|
354 | :value (octi::log-qd/newton (qd-value arg)))) |
---|
355 | |
---|
356 | (rt:deftest oct.log/newton.1 |
---|
357 | (let* ((arg #q2) |
---|
358 | (y (log/newton arg)) |
---|
359 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
360 | (check-accuracy 212 y true)) |
---|
361 | nil) |
---|
362 | |
---|
363 | (rt:deftest oct.log/newton.2 |
---|
364 | (let* ((arg #q10) |
---|
365 | (y (log/newton arg)) |
---|
366 | (true (make-instance 'qd-real :value octi::+qd-log10+))) |
---|
367 | (check-accuracy 207 y true)) |
---|
368 | nil) |
---|
369 | |
---|
370 | (rt:deftest oct.log/newton.3 |
---|
371 | (let* ((arg (+ 1 (scale-float #q1 -80))) |
---|
372 | (y (log/newton arg)) |
---|
373 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
374 | (check-accuracy 212 y true)) |
---|
375 | nil) |
---|
376 | |
---|
377 | ;;; Tests of log using AGM. |
---|
378 | |
---|
379 | (defun log/agm (arg) |
---|
380 | (make-instance 'qd-real |
---|
381 | :value (octi::log-qd/agm (qd-value arg)))) |
---|
382 | |
---|
383 | (rt:deftest oct.log/agm.1 |
---|
384 | (let* ((arg #q2) |
---|
385 | (y (log/agm arg)) |
---|
386 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
387 | (check-accuracy 203 y true)) |
---|
388 | nil) |
---|
389 | |
---|
390 | (rt:deftest oct.log/agm.2 |
---|
391 | (let* ((arg #q10) |
---|
392 | (y (log/agm arg)) |
---|
393 | (true (make-instance 'qd-real :value octi::+qd-log10+))) |
---|
394 | (check-accuracy 205 y true)) |
---|
395 | nil) |
---|
396 | |
---|
397 | (rt:deftest oct.log/agm.3 |
---|
398 | (let* ((arg (+ 1 (scale-float #q1 -80))) |
---|
399 | (y (log/agm arg)) |
---|
400 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
401 | (check-accuracy 123 y true)) |
---|
402 | nil) |
---|
403 | |
---|
404 | ;;; Tests of log using AGM2, a faster variaton of AGM. |
---|
405 | |
---|
406 | (defun log/agm2 (arg) |
---|
407 | (make-instance 'qd-real |
---|
408 | :value (octi::log-qd/agm2 (qd-value arg)))) |
---|
409 | |
---|
410 | (rt:deftest oct.log/agm2.1 |
---|
411 | (let* ((arg #q2) |
---|
412 | (y (log/agm2 arg)) |
---|
413 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
414 | (check-accuracy 203 y true)) |
---|
415 | nil) |
---|
416 | |
---|
417 | (rt:deftest oct.log/agm2.2 |
---|
418 | (let* ((arg #q10) |
---|
419 | (y (log/agm2 arg)) |
---|
420 | (true (make-instance 'qd-real :value octi::+qd-log10+))) |
---|
421 | (check-accuracy 205 y true)) |
---|
422 | nil) |
---|
423 | |
---|
424 | (rt:deftest oct.log/agm2.3 |
---|
425 | (let* ((arg (+ 1 (scale-float #q1 -80))) |
---|
426 | (y (log/agm2 arg)) |
---|
427 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
428 | (check-accuracy 123 y true)) |
---|
429 | nil) |
---|
430 | |
---|
431 | ;;; Tests of log using AGM3, a faster variation of AGM2. |
---|
432 | (defun log/agm3 (arg) |
---|
433 | (make-instance 'qd-real |
---|
434 | :value (octi::log-qd/agm3 (qd-value arg)))) |
---|
435 | |
---|
436 | (rt:deftest oct.log/agm3.1 |
---|
437 | (let* ((arg #q2) |
---|
438 | (y (log/agm3 arg)) |
---|
439 | (true (make-instance 'qd-real :value octi::+qd-log2+))) |
---|
440 | (check-accuracy 203 y true)) |
---|
441 | nil) |
---|
442 | |
---|
443 | (rt:deftest oct.log/agm3.2 |
---|
444 | (let* ((arg #q10) |
---|
445 | (y (log/agm3 arg)) |
---|
446 | (true (make-instance 'qd-real :value octi::+qd-log10+))) |
---|
447 | (check-accuracy 205 y true)) |
---|
448 | nil) |
---|
449 | |
---|
450 | (rt:deftest oct.log/agm3.3 |
---|
451 | (let* ((arg (+ 1 (scale-float #q1 -80))) |
---|
452 | (y (log/agm3 arg)) |
---|
453 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
454 | (check-accuracy 123 y true)) |
---|
455 | nil) |
---|
456 | |
---|
457 | ;;; Tests of sqrt to make sure we overflow or underflow where we |
---|
458 | ;;; shouldn't. |
---|
459 | |
---|
460 | (rt:deftest oct.sqrt.1 |
---|
461 | (let* ((arg #q1q200) |
---|
462 | (y (sqrt arg)) |
---|
463 | (true #q1q100)) |
---|
464 | (check-accuracy 212 y true)) |
---|
465 | nil) |
---|
466 | |
---|
467 | (rt:deftest oct.sqrt.2 |
---|
468 | (let* ((arg #q1q200) |
---|
469 | (y (sqrt arg)) |
---|
470 | (true #q1q100)) |
---|
471 | (check-accuracy 212 y true)) |
---|
472 | nil) |
---|
473 | |
---|
474 | (rt:deftest oct.sqrt.3 |
---|
475 | (let* ((arg #q1q300) |
---|
476 | (y (sqrt arg)) |
---|
477 | (true #q1q150)) |
---|
478 | (check-accuracy 212 y true)) |
---|
479 | nil) |
---|
480 | |
---|
481 | (rt:deftest oct.sqrt.4 |
---|
482 | (let* ((arg #q1q-200) |
---|
483 | (y (sqrt arg)) |
---|
484 | (true #q1q-100)) |
---|
485 | (check-accuracy 212 y true)) |
---|
486 | nil) |
---|
487 | |
---|
488 | (rt:deftest oct.sqrt.5 |
---|
489 | (let* ((arg #q1q-250) |
---|
490 | (y (sqrt arg)) |
---|
491 | (true #q1q-125)) |
---|
492 | (check-accuracy 212 y true)) |
---|
493 | nil) |
---|
494 | |
---|
495 | ;;; Tests of log1p(x) = log(1+x), using the duplication formula. |
---|
496 | |
---|
497 | (defun log1p/dup (arg) |
---|
498 | (make-instance 'qd-real |
---|
499 | :value (octi::log1p-qd/duplication (qd-value arg)))) |
---|
500 | |
---|
501 | (rt:deftest oct.log1p.1 |
---|
502 | (let* ((arg #q9) |
---|
503 | (y (log1p/dup arg)) |
---|
504 | (true #q2.3025850929940456840179914546843642076011014886287729760333279009675726096773525q0)) |
---|
505 | (check-accuracy 212 y true)) |
---|
506 | nil) |
---|
507 | |
---|
508 | (rt:deftest oct.log1p.2 |
---|
509 | (let* ((arg (scale-float #q1 -80)) |
---|
510 | (y (log1p/dup arg)) |
---|
511 | (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25)) |
---|
512 | (check-accuracy 212 y true)) |
---|
513 | nil) |
---|
514 | |
---|
515 | ;;; Tests of expm1(x) = exp(x) - 1, using a Taylor series with |
---|
516 | ;;; argument reduction. |
---|
517 | |
---|
518 | (defun expm1/series (arg) |
---|
519 | (make-instance 'qd-real |
---|
520 | :value (octi::expm1-qd/series (qd-value arg)))) |
---|
521 | |
---|
522 | (rt:deftest oct.expm1/series.1 |
---|
523 | (let* ((arg #q0) |
---|
524 | (y (expm1/series arg)) |
---|
525 | (true #q0)) |
---|
526 | (check-accuracy 212 y true)) |
---|
527 | nil) |
---|
528 | |
---|
529 | (rt:deftest oct.expm1/series.2 |
---|
530 | (let* ((arg #q1) |
---|
531 | (y (expm1/series arg)) |
---|
532 | (true #q1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0)) |
---|
533 | (check-accuracy 211 y true)) |
---|
534 | nil) |
---|
535 | |
---|
536 | (rt:deftest oct.expm1/series.3 |
---|
537 | (let* ((arg (scale-float #q1 -100)) |
---|
538 | (y (expm1/series arg)) |
---|
539 | (true #q7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31)) |
---|
540 | (check-accuracy 211 y true)) |
---|
541 | nil) |
---|
542 | |
---|
543 | ;;; Tests of expm1(x) = exp(x) - 1, using duplication formula. |
---|
544 | |
---|
545 | (defun expm1/dup (arg) |
---|
546 | (make-instance 'qd-real |
---|
547 | :value (octi::expm1-qd/duplication (qd-value arg)))) |
---|
548 | |
---|
549 | |
---|
550 | (rt:deftest oct.expm1/dup.1 |
---|
551 | (let* ((arg #q0) |
---|
552 | (y (expm1/dup arg)) |
---|
553 | (true #q0)) |
---|
554 | (check-accuracy 212 y true)) |
---|
555 | nil) |
---|
556 | |
---|
557 | (rt:deftest oct.expm1/dup.2 |
---|
558 | (let* ((arg #q1) |
---|
559 | (y (expm1/dup arg)) |
---|
560 | (true #q1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0)) |
---|
561 | (check-accuracy 211 y true)) |
---|
562 | nil) |
---|
563 | |
---|
564 | (rt:deftest oct.expm1/dup.3 |
---|
565 | (let* ((arg (scale-float #q1 -100)) |
---|
566 | (y (expm1/dup arg)) |
---|
567 | (true #q7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31)) |
---|
568 | (check-accuracy 211 y true)) |
---|
569 | nil) |
---|
570 | |
---|
571 | ;; If we screw up integer-decode-qd, printing is wrong. Here is one |
---|
572 | ;; case where integer-decode-qd was screwed up and printing the wrong |
---|
573 | ;; thing. |
---|
574 | (rt:deftest oct.integer-decode.1 |
---|
575 | (multiple-value-bind (frac exp s) |
---|
576 | (octi:integer-decode-qd (octi::%make-qd-d -0.03980126756814893d0 |
---|
577 | -2.7419792323327893d-18 |
---|
578 | 0d0 0d0)) |
---|
579 | (unless (and (eql frac 103329998279901916046530991816704) |
---|
580 | (eql exp -111) |
---|
581 | (eql s -1)) |
---|
582 | (list frac exp s))) |
---|
583 | nil) |
---|
584 | |
---|
585 | ;;; |
---|
586 | ;;; Add a few tests for the branch cuts. Many of these tests assume |
---|
587 | ;;; that Lisp has support for signed zeroes. If not, these tests are |
---|
588 | ;;; probably wrong. |
---|
589 | |
---|
590 | (defun check-signs (fun arg expected) |
---|
591 | (let* ((z (funcall fun arg)) |
---|
592 | (x (realpart z)) |
---|
593 | (y (imagpart z))) |
---|
594 | ;; If the Lisp doesn't support signed zeroes, then this test |
---|
595 | ;; should always pass. |
---|
596 | (if (or (eql -0d0 0d0) |
---|
597 | (and (= (float-sign x) (float-sign (realpart expected))) |
---|
598 | (= (float-sign y) (float-sign (imagpart expected))))) |
---|
599 | t |
---|
600 | (list z expected fun arg)))) |
---|
601 | |
---|
602 | ;; asin has a branch cut on the real axis |x|>1. For x < -1, it is |
---|
603 | ;; continuous with quadrant II; for x > 1, continuous with quadrant |
---|
604 | ;; IV. |
---|
605 | (rt:deftest oct.asin-branch-neg.1 |
---|
606 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
607 | (check-signs #'asin -2d0 true)) |
---|
608 | t) |
---|
609 | |
---|
610 | (rt:deftest oct.asin-branch-neg.2 |
---|
611 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
612 | (check-signs #'asin #q-2 true)) |
---|
613 | t) |
---|
614 | |
---|
615 | (rt:deftest oct.asin-branch-neg.3 |
---|
616 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
617 | (check-signs #'asin #c(-2d0 0d0) true)) |
---|
618 | t) |
---|
619 | |
---|
620 | (rt:deftest oct.asin-branch-neg.4 |
---|
621 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
622 | (check-signs #'asin #q(-2 0) true)) |
---|
623 | t) |
---|
624 | |
---|
625 | (rt:deftest oct.asin-branch-neg.5 |
---|
626 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
627 | (check-signs #'asin #c(-2d0 -0d0) (conjugate true))) |
---|
628 | t) |
---|
629 | |
---|
630 | (rt:deftest oct.asin-branch-neg.6 |
---|
631 | (let ((true (cl:asin #c(-2d0 1d-20)))) |
---|
632 | (check-signs #'asin #q(-2d0 -0d0) (conjugate true))) |
---|
633 | t) |
---|
634 | |
---|
635 | (rt:deftest oct.asin-branch-pos.1 |
---|
636 | (let ((true (cl:asin #c(2d0 -1d-20)))) |
---|
637 | (check-signs #'asin #c(2d0 0d0) (conjugate true))) |
---|
638 | t) |
---|
639 | |
---|
640 | (rt:deftest oct.asin-branch-pos.2 |
---|
641 | (let ((true (cl:asin #c(2d0 -1d-20)))) |
---|
642 | (check-signs #'asin #q(2 0d0) (conjugate true))) |
---|
643 | t) |
---|
644 | |
---|
645 | (rt:deftest oct.asin-branch-pos.3 |
---|
646 | (let ((true (cl:asin #c(2d0 -1d-20)))) |
---|
647 | (check-signs #'asin #c(2d0 -0d0) true)) |
---|
648 | t) |
---|
649 | |
---|
650 | (rt:deftest oct.asin-branch-pos.4 |
---|
651 | (let ((true (cl:asin #c(2d0 -1d-20)))) |
---|
652 | (check-signs #'asin #q(2d0 -0d0) true)) |
---|
653 | t) |
---|
654 | |
---|
655 | ;; acos branch cut is the real axis, |x| > 1. For x < -1, it is |
---|
656 | ;; continuous with quadrant II; for x > 1, quadrant IV. |
---|
657 | |
---|
658 | (rt:deftest oct.acos-branch-neg.1 |
---|
659 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
660 | (check-signs #'acos -2d0 true)) |
---|
661 | t) |
---|
662 | |
---|
663 | (rt:deftest oct.acos-branch-neg.2 |
---|
664 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
665 | (check-signs #'acos #q-2 true)) |
---|
666 | t) |
---|
667 | |
---|
668 | (rt:deftest oct.acos-branch-neg.3 |
---|
669 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
670 | (check-signs #'acos #c(-2d0 0d0) true)) |
---|
671 | t) |
---|
672 | |
---|
673 | (rt:deftest oct.acos-branch-neg.4 |
---|
674 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
675 | (check-signs #'acos #q(-2 0) true)) |
---|
676 | t) |
---|
677 | |
---|
678 | (rt:deftest oct.acos-branch-neg.5 |
---|
679 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
680 | (check-signs #'acos #c(-2d0 -0d0) (conjugate true))) |
---|
681 | t) |
---|
682 | |
---|
683 | (rt:deftest oct.acos-branch-neg.6 |
---|
684 | (let ((true (cl:acos #c(-2d0 1d-20)))) |
---|
685 | (check-signs #'acos #q(-2d0 -0d0) (conjugate true))) |
---|
686 | t) |
---|
687 | |
---|
688 | (rt:deftest oct.acos-branch-pos.1 |
---|
689 | (let ((true (cl:acos #c(2d0 -1d-20)))) |
---|
690 | (check-signs #'acos #c(2d0 0d0) (conjugate true))) |
---|
691 | t) |
---|
692 | |
---|
693 | (rt:deftest oct.acos-branch-pos.2 |
---|
694 | (let ((true (cl:acos #c(2d0 -1d-20)))) |
---|
695 | (check-signs #'acos #q(2 0d0) (conjugate true))) |
---|
696 | t) |
---|
697 | |
---|
698 | (rt:deftest oct.acos-branch-pos.3 |
---|
699 | (let ((true (cl:acos #c(2d0 -1d-20)))) |
---|
700 | (check-signs #'acos #c(2d0 -0d0) true)) |
---|
701 | t) |
---|
702 | |
---|
703 | (rt:deftest oct.acos-branch-pos.4 |
---|
704 | (let ((true (cl:acos #c(2d0 -1d-20)))) |
---|
705 | (check-signs #'acos #q(2d0 -0d0) true)) |
---|
706 | t) |
---|
707 | |
---|
708 | ;; atan branch cut is the imaginary axis, |y| > 1. For y < -1, it is |
---|
709 | ;; continuous with quadrant IV; for x > 1, quadrant II. |
---|
710 | |
---|
711 | (rt:deftest oct.atan-branch-neg.1 |
---|
712 | (let ((true (cl:atan #c(1d-20 -2d0)))) |
---|
713 | (check-signs #'atan #c(0d0 -2d0) true)) |
---|
714 | t) |
---|
715 | |
---|
716 | (rt:deftest oct.atan-branch-neg.2 |
---|
717 | (let ((true (cl:atan #c(1d-20 -2d0)))) |
---|
718 | (check-signs #'atan #q(0 -2) true)) |
---|
719 | t) |
---|
720 | |
---|
721 | (rt:deftest oct.atan-branch-neg.3 |
---|
722 | (let ((true (cl:atan #c(-1d-20 -2d0)))) |
---|
723 | (check-signs #'atan #c(-0d0 -2d0) true)) |
---|
724 | t) |
---|
725 | |
---|
726 | (rt:deftest oct.atan-branch-neg.4 |
---|
727 | (let ((true (cl:atan #c(-1d-20 -2d0)))) |
---|
728 | (check-signs #'atan #q(-0d0 -2d0) true)) |
---|
729 | t) |
---|
730 | |
---|
731 | (rt:deftest oct.atan-branch-pos.1 |
---|
732 | (let ((true (cl:atan #c(1d-20 2d0)))) |
---|
733 | (check-signs #'atan #c(0d0 2d0) true)) |
---|
734 | t) |
---|
735 | |
---|
736 | (rt:deftest oct.atan-branch-pos.2 |
---|
737 | (let ((true (cl:atan #c(1d-20 2d0)))) |
---|
738 | (check-signs #'atan #q(0d0 2d0) true)) |
---|
739 | t) |
---|
740 | |
---|
741 | (rt:deftest oct.atan-branch-pos.3 |
---|
742 | (let ((true (cl:atan #c(-1d-20 2d0)))) |
---|
743 | (check-signs #'atan #c(-0d0 2d0) true)) |
---|
744 | t) |
---|
745 | |
---|
746 | (rt:deftest oct.atan-branch-pos.4 |
---|
747 | (let ((true (cl:atan #c(-1d-20 2d0)))) |
---|
748 | (check-signs #'atan #q(-0d0 2d0) true)) |
---|
749 | t) |
---|
750 | |
---|
751 | ;; Test x < -1 |
---|
752 | (rt:deftest oct.atanh-branch-neg.1 |
---|
753 | (let ((true (cl:atanh #c(-2d0 -1d-20)))) |
---|
754 | (check-signs #'atanh -2d0 true)) |
---|
755 | t) |
---|
756 | |
---|
757 | (rt:deftest oct.atanh-branch-neg.2 |
---|
758 | (let ((true (cl:atanh #c(-2d0 -1d-20)))) |
---|
759 | (check-signs #'atanh #q-2 true)) |
---|
760 | t) |
---|
761 | |
---|
762 | ;; Test x > 1 |
---|
763 | (rt:deftest oct.atanh-branch-pos.1 |
---|
764 | (let ((true (cl:atanh #c(2d0 1d-20)))) |
---|
765 | (check-signs #'atanh 2d0 true)) |
---|
766 | t) |
---|
767 | |
---|
768 | (rt:deftest oct.atanh-branch-pos.2 |
---|
769 | (let ((true (cl:atanh #c(2d0 1d-20)))) |
---|
770 | (check-signs #'atanh #q2 true)) |
---|
771 | t) |
---|
772 | |
---|
773 | ;; elliptic_k(-1) = gamma(1/4)^2/2^(5/2)/sqrt(%pi) |
---|
774 | (rt:deftest oct.elliptic-k.1d |
---|
775 | (let* ((val (elliptic-k -1d0)) |
---|
776 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395293535207125115147766480714547q0)) |
---|
777 | (check-accuracy 53 val true)) |
---|
778 | nil) |
---|
779 | |
---|
780 | (rt:deftest oct.elliptic-k.1q |
---|
781 | (let* ((val (elliptic-k #q-1q0)) |
---|
782 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395293535207125115147766480714547q0)) |
---|
783 | (check-accuracy 210 val true)) |
---|
784 | nil) |
---|
785 | |
---|
786 | ;; elliptic_k(1/2) = %pi^(3/2)/2/gamma(3/4)^2 |
---|
787 | (rt:deftest oct.elliptic-k.2d |
---|
788 | (let* ((val (elliptic-k 0.5d0)) |
---|
789 | (true #q1.854074677301371918433850347195260046217598823521766905585928045056021776838119978357271861650371897q0)) |
---|
790 | (check-accuracy 53 val true)) |
---|
791 | nil) |
---|
792 | |
---|
793 | (rt:deftest oct.elliptic-k.2q |
---|
794 | (let* ((val (elliptic-k #q.5)) |
---|
795 | (true #q1.854074677301371918433850347195260046217598823521766905585928045056021776838119978357271861650371897q0)) |
---|
796 | (check-accuracy 210 val true)) |
---|
797 | nil) |
---|
798 | |
---|
799 | ;; jacobi_sn(K,1/2) = 1, where K = elliptic_k(1/2) |
---|
800 | (rt:deftest oct.jacobi-sn.1d |
---|
801 | (let* ((ek (elliptic-k .5d0)) |
---|
802 | (val (jacobi-sn ek .5d0))) |
---|
803 | (check-accuracy 54 val 1d0)) |
---|
804 | nil) |
---|
805 | |
---|
806 | (rt:deftest oct.jacobi-sn.1q |
---|
807 | (let* ((ek (elliptic-k #q.5)) |
---|
808 | (val (jacobi-sn ek #q.5))) |
---|
809 | (check-accuracy 212 val #q1)) |
---|
810 | nil) |
---|
811 | |
---|
812 | ;; jacobi_cn(K,1/2) = 0 |
---|
813 | (rt:deftest oct.jacobi-cn.1d |
---|
814 | (let* ((ek (elliptic-k .5d0)) |
---|
815 | (val (jacobi-cn ek .5d0))) |
---|
816 | (check-accuracy 50 val 0d0)) |
---|
817 | nil) |
---|
818 | |
---|
819 | (rt:deftest oct.jacobi-sn.1q |
---|
820 | (let* ((ek (elliptic-k #q.5)) |
---|
821 | (val (jacobi-cn ek #q.5))) |
---|
822 | (check-accuracy 210 val #q0)) |
---|
823 | nil) |
---|
824 | |
---|
825 | ;; jacobi-dn(K, 1/2) = sqrt(1/2) |
---|
826 | (rt:deftest oct.jacobi-dn.1d |
---|
827 | (let* ((ek (elliptic-k .5d0)) |
---|
828 | (true (sqrt .5d0)) |
---|
829 | (val (jacobi-dn ek .5d0))) |
---|
830 | (check-accuracy 52 val true)) |
---|
831 | nil) |
---|
832 | |
---|
833 | (rt:deftest oct.jacobi-dn.1q |
---|
834 | (let* ((ek (elliptic-k #q.5)) |
---|
835 | (true (sqrt #q.5)) |
---|
836 | (val (jacobi-dn ek #q.5))) |
---|
837 | (check-accuracy 212 val true)) |
---|
838 | nil) |
---|
839 | |
---|
840 | (rt:deftest oct.carlson-rf.1d |
---|
841 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
842 | ;; = 1/4*beta(1/2,1/2) |
---|
843 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
844 | (let ((rf (carlson-rf 0d0 2d0 1d0)) |
---|
845 | (true 1.31102877714605990523241979494d0)) |
---|
846 | (check-accuracy 53 rf true)) |
---|
847 | nil) |
---|
848 | |
---|
849 | (rt:deftest oct.carlson-rf.1q |
---|
850 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
851 | (let ((rf (carlson-rf #q0 #q2 #q1)) |
---|
852 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0)) |
---|
853 | (check-accuracy 212 rf true)) |
---|
854 | nil) |
---|
855 | |
---|
856 | (rt:deftest oct.carlson-rd.1d |
---|
857 | ;; Rd(0,2,1) = 3*integrate(s^2/sqrt(1-s^4), s, 0 ,1) |
---|
858 | ;; = 3*beta(3/4,1/2)/4 |
---|
859 | ;; = 3*sqrt(%pi)*gamma(3/4)/gamma(1/4) |
---|
860 | (let ((rd (carlson-rd 0d0 2d0 1d0)) |
---|
861 | (true 1.7972103521033883d0)) |
---|
862 | (check-accuracy 51 rd true)) |
---|
863 | nil) |
---|
864 | |
---|
865 | (rt:deftest oct.carlson-rd.1q |
---|
866 | (let ((rd (carlson-rd #q0 #q2 #q1)) |
---|
867 | (true #q1.797210352103388311159883738420485817340818994823477337395512429419599q0)) |
---|
868 | (check-accuracy 212 rd true)) |
---|
869 | nil) |
---|
870 | |
---|
871 | ;; Test some of the contagion stuff. |
---|
872 | |
---|
873 | (rt:deftest oct.carlson-rf.contagion.1 |
---|
874 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
875 | ;; = 1/4*beta(1/2,1/2) |
---|
876 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
877 | (let ((rf (carlson-rf 0 2 1)) |
---|
878 | (true 1.31102877714605990523241979494d0)) |
---|
879 | (check-accuracy 23 rf true)) |
---|
880 | nil) |
---|
881 | |
---|
882 | (rt:deftest oct.carlson-rf.contagion.1d |
---|
883 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
884 | ;; = 1/4*beta(1/2,1/2) |
---|
885 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
886 | (let ((rf (carlson-rf 0d0 2 1)) |
---|
887 | (true 1.31102877714605990523241979494d0)) |
---|
888 | (check-accuracy 53 rf true)) |
---|
889 | nil) |
---|
890 | |
---|
891 | (rt:deftest oct.carlson-rf.contagion.2d |
---|
892 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
893 | ;; = 1/4*beta(1/2,1/2) |
---|
894 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
895 | (let ((rf (carlson-rf 0 2d0 1)) |
---|
896 | (true 1.31102877714605990523241979494d0)) |
---|
897 | (check-accuracy 53 rf true)) |
---|
898 | nil) |
---|
899 | |
---|
900 | (rt:deftest oct.carlson-rf.contagion.3d |
---|
901 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
902 | ;; = 1/4*beta(1/2,1/2) |
---|
903 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
904 | (let ((rf (carlson-rf 0 2 1d0)) |
---|
905 | (true 1.31102877714605990523241979494d0)) |
---|
906 | (check-accuracy 53 rf true)) |
---|
907 | nil) |
---|
908 | |
---|
909 | (rt:deftest oct.carlson-rf.contagion.1q |
---|
910 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
911 | ;; = 1/4*beta(1/2,1/2) |
---|
912 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
913 | (let ((rf (carlson-rf #q0q0 2 1)) |
---|
914 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0)) |
---|
915 | (check-accuracy 212 rf true)) |
---|
916 | nil) |
---|
917 | |
---|
918 | (rt:deftest oct.carlson-rf.contagion.2q |
---|
919 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
920 | ;; = 1/4*beta(1/2,1/2) |
---|
921 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
922 | (let ((rf (carlson-rf 0 #q2q0 1)) |
---|
923 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0)) |
---|
924 | (check-accuracy 212 rf true)) |
---|
925 | nil) |
---|
926 | |
---|
927 | (rt:deftest oct.carlson-rf.contagion.3q |
---|
928 | ;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1) |
---|
929 | ;; = 1/4*beta(1/2,1/2) |
---|
930 | ;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4) |
---|
931 | (let ((rf (carlson-rf 0 2 #q1q0)) |
---|
932 | (true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0)) |
---|
933 | (check-accuracy 212 rf true)) |
---|
934 | nil) |
---|
935 | |
---|
936 | ;; Elliptic integral of the third kind |
---|
937 | |
---|
938 | ;; elliptic-pi(0,phi,m) = elliptic-f(phi, m) |
---|
939 | (rt:deftest oct.elliptic-pi.1d |
---|
940 | (loop for k from 0 to 100 |
---|
941 | for phi = (random (/ pi 2)) |
---|
942 | for m = (random 1d0) |
---|
943 | for epi = (elliptic-pi 0 phi m) |
---|
944 | for ef = (elliptic-f phi m) |
---|
945 | for result = (check-accuracy 48 epi ef) |
---|
946 | unless (eq nil result) |
---|
947 | append (list (list phi m) result)) |
---|
948 | nil) |
---|
949 | |
---|
950 | (rt:deftest oct.elliptic-pi.1q |
---|
951 | (loop for k from 0 below 100 |
---|
952 | for phi = (random (/ +pi+ 2)) |
---|
953 | for m = (random #q1) |
---|
954 | for epi = (elliptic-pi 0 phi m) |
---|
955 | for ef = (elliptic-f phi m) |
---|
956 | for result = (check-accuracy 53 epi ef) |
---|
957 | unless (eq nil result) |
---|
958 | append (list (list phi m) result)) |
---|
959 | nil) |
---|
960 | |
---|
961 | ;; DLMF 19.6.3 |
---|
962 | ;; |
---|
963 | ;; PI(n; pi/2 | 0) = pi/(2*sqrt(1-n)) |
---|
964 | (rt:deftest oct.elliptic-pi.19.6.3.d |
---|
965 | (loop for k from 0 below 100 |
---|
966 | for n = (random 1d0) |
---|
967 | for epi = (elliptic-pi n (/ pi 2) 0) |
---|
968 | for true = (/ pi (* 2 (sqrt (- 1 n)))) |
---|
969 | for result = (check-accuracy 49 epi true) |
---|
970 | unless (eq nil result) |
---|
971 | append (list (list (list k n) result))) |
---|
972 | nil) |
---|
973 | |
---|
974 | (rt:deftest oct.elliptic-pi.19.6.3.q |
---|
975 | (loop for k from 0 below 100 |
---|
976 | for n = (random #q1) |
---|
977 | for epi = (elliptic-pi n (/ (float-pi n) 2) 0) |
---|
978 | for true = (/ (float-pi n) (* 2 (sqrt (- 1 n)))) |
---|
979 | for result = (check-accuracy 209 epi true) |
---|
980 | unless (eq nil result) |
---|
981 | append (list (list (list k n) result))) |
---|
982 | nil) |
---|
983 | |
---|
984 | ;; elliptic-pi(n, phi, 0) = |
---|
985 | ;; atan(sqrt(1-n)*tan(phi))/sqrt(1-n) n < 1 |
---|
986 | ;; atanh(sqrt(n-1)*tan(phi))/sqrt(n-1) n > 1 |
---|
987 | ;; tan(phi) n = 1 |
---|
988 | ;; |
---|
989 | ;; These are easy to derive if you look at the integral: |
---|
990 | ;; |
---|
991 | ;; ellipti-pi(n, phi, 0) = integrate(1/(1-n*sin(t)^2), t, 0, phi) |
---|
992 | ;; |
---|
993 | ;; and this can be easily integrated to give the above expressions for |
---|
994 | ;; the different values of n. |
---|
995 | (rt:deftest oct.elliptic-pi.n0.d |
---|
996 | ;; Tests for random values for phi in [0, pi/2] and n in [0, 1] |
---|
997 | (loop for k from 0 below 100 |
---|
998 | for phi = (random (/ pi 2)) |
---|
999 | for n = (random 1d0) |
---|
1000 | for epi = (elliptic-pi n phi 0) |
---|
1001 | for true = (/ (atan (* (tan phi) (sqrt (- 1 n)))) |
---|
1002 | (sqrt (- 1 n))) |
---|
1003 | for result = (check-accuracy 47.5 epi true) |
---|
1004 | unless (eq nil result) |
---|
1005 | append (list (list (list k n phi) result))) |
---|
1006 | nil) |
---|
1007 | |
---|
1008 | (rt:deftest oct.elliptic-pi.n1.d |
---|
1009 | (loop for k from 0 below 100 |
---|
1010 | for phi = (random (/ pi 2)) |
---|
1011 | for epi = (elliptic-pi 1 phi 0) |
---|
1012 | for true = (tan phi) |
---|
1013 | for result = (check-accuracy 36 epi true) |
---|
1014 | unless (eq nil result) |
---|
1015 | append (list (list (list k phi) result))) |
---|
1016 | nil) |
---|
1017 | |
---|
1018 | (rt:deftest oct.elliptic-pi.n2.d |
---|
1019 | (loop for k from 0 below 100 |
---|
1020 | for phi = (random (/ pi 2)) |
---|
1021 | for n = (+ 1d0 (random 100d0)) |
---|
1022 | for epi = (elliptic-pi n phi 0) |
---|
1023 | for true = (/ (atanh (* (tan phi) (sqrt (- n 1)))) |
---|
1024 | (sqrt (- n 1))) |
---|
1025 | for result = (check-accuracy 47 epi true) |
---|
1026 | ;; Not sure if this formula holds when atanh gives a complex |
---|
1027 | ;; result. Wolfram doesn't say |
---|
1028 | when (and (not (complexp true)) result) |
---|
1029 | append (list (list (list k n phi) result))) |
---|
1030 | nil) |
---|
1031 | |
---|
1032 | (rt:deftest oct.elliptic-pi.n0.q |
---|
1033 | ;; Tests for random values for phi in [0, pi/2] and n in [0, 1] |
---|
1034 | (loop for k from 0 below 100 |
---|
1035 | for phi = (random (/ +pi+ 2)) |
---|
1036 | for n = (random #q1) |
---|
1037 | for epi = (elliptic-pi n phi 0) |
---|
1038 | for true = (/ (atan (* (tan phi) (sqrt (- 1 n)))) |
---|
1039 | (sqrt (- 1 n))) |
---|
1040 | for result = (check-accuracy 208 epi true) |
---|
1041 | unless (eq nil result) |
---|
1042 | append (list (list (list k n phi) result))) |
---|
1043 | nil) |
---|
1044 | |
---|
1045 | (rt:deftest oct.elliptic-pi.n1.q |
---|
1046 | (loop for k from 0 below 100 |
---|
1047 | for phi = (random (/ +pi+ 2)) |
---|
1048 | for epi = (elliptic-pi 1 phi 0) |
---|
1049 | for true = (tan phi) |
---|
1050 | for result = (check-accuracy 194 epi true) |
---|
1051 | unless (eq nil result) |
---|
1052 | append (list (list (list k phi) result))) |
---|
1053 | nil) |
---|
1054 | |
---|
1055 | (rt:deftest oct.elliptic-pi.n2.q |
---|
1056 | (loop for k from 0 below 100 |
---|
1057 | for phi = (random (/ +pi+ 2)) |
---|
1058 | for n = (+ #q1 (random #q1)) |
---|
1059 | for epi = (elliptic-pi n phi 0) |
---|
1060 | for true = (/ (atanh (* (tan phi) (sqrt (- n 1)))) |
---|
1061 | (sqrt (- n 1))) |
---|
1062 | for result = (check-accuracy 206 epi true) |
---|
1063 | ;; Not sure if this formula holds when atanh gives a complex |
---|
1064 | ;; result. Wolfram doesn't say |
---|
1065 | when (and (not (complexp true)) result) |
---|
1066 | append (list (list (list k n phi) result))) |
---|
1067 | nil) |
---|
1068 | |
---|
1069 | ;; Tests for theta functions. |
---|
1070 | |
---|
1071 | (rt:deftest oct.theta3.1.d |
---|
1072 | ;; A&S 16.38.5 |
---|
1073 | ;; sqrt(2*K/%pi) = theta3(0,q) |
---|
1074 | (loop for k from 0 below 100 |
---|
1075 | for m = (random 1d0) |
---|
1076 | for t3 = (theta3 0 (elliptic-nome m)) |
---|
1077 | for true = (sqrt (/ (* 2 (elliptic-k m)) (float-pi m))) |
---|
1078 | for result = (check-accuracy 51 t3 true) |
---|
1079 | when result |
---|
1080 | append (list (list (list k m) result))) |
---|
1081 | nil) |
---|
1082 | |
---|
1083 | (rt:deftest oct.theta3.1.q |
---|
1084 | ;; A&S 16.38.5 |
---|
1085 | ;; sqrt(2*K/%pi) = theta3(0,q) |
---|
1086 | (loop for k from 0 below 100 |
---|
1087 | for m = (random #q1) |
---|
1088 | for t3 = (theta3 0 (elliptic-nome m)) |
---|
1089 | for true = (sqrt (/ (* 2 (elliptic-k m)) (float-pi m))) |
---|
1090 | for result = (check-accuracy 206 t3 true) |
---|
1091 | when result |
---|
1092 | append (list (list (list k m) result))) |
---|
1093 | nil) |
---|
1094 | |
---|
1095 | (rt:deftest oct.theta2.1.d |
---|
1096 | ;; A&S 16.38.7 |
---|
1097 | ;; sqrt(2*sqrt(m)*K/%pi) = theta2(0,q) |
---|
1098 | (loop for k from 0 below 100 |
---|
1099 | for m = (random 1d0) |
---|
1100 | for t3 = (theta2 0 (elliptic-nome m)) |
---|
1101 | for true = (sqrt (/ (* 2 (sqrt m) (elliptic-k m)) (float-pi m))) |
---|
1102 | for result = (check-accuracy 49 t3 true) |
---|
1103 | when result |
---|
1104 | append (list (list (list k m) result))) |
---|
1105 | nil) |
---|
1106 | |
---|
1107 | (rt:deftest oct.theta2.1.q |
---|
1108 | ;; A&S 16.38.7 |
---|
1109 | ;; sqrt(2*sqrt(m)*K/%pi) = theta2(0,q) |
---|
1110 | (loop for k from 0 below 100 |
---|
1111 | for m = (random #q1) |
---|
1112 | for t3 = (theta2 0 (elliptic-nome m)) |
---|
1113 | for true = (sqrt (/ (* 2 (sqrt m) (elliptic-k m)) (float-pi m))) |
---|
1114 | for result = (check-accuracy 206 t3 true) |
---|
1115 | when result |
---|
1116 | append (list (list (list k m) result))) |
---|
1117 | nil) |
---|
1118 | |
---|
1119 | (rt:deftest oct.theta4.1.d |
---|
1120 | ;; A&S 16.38.8 |
---|
1121 | ;; sqrt(2*sqrt(1-m)*K/%pi) = theta2(0,q) |
---|
1122 | (loop for k from 0 below 100 |
---|
1123 | for m = (random 1d0) |
---|
1124 | for t3 = (theta4 0 (elliptic-nome m)) |
---|
1125 | for true = (sqrt (/ (* 2 (sqrt (- 1 m)) (elliptic-k m)) |
---|
1126 | (float-pi m))) |
---|
1127 | for result = (check-accuracy 49 t3 true) |
---|
1128 | when result |
---|
1129 | append (list (list (list k m) result))) |
---|
1130 | nil) |
---|
1131 | |
---|
1132 | (rt:deftest oct.theta4.1.q |
---|
1133 | ;; A&S 16.38.8 |
---|
1134 | ;; sqrt(2*sqrt(1-m)*K/%pi) = theta2(0,q) |
---|
1135 | (loop for k from 0 below 100 |
---|
1136 | for m = (random #q1) |
---|
1137 | for t3 = (theta4 0 (elliptic-nome m)) |
---|
1138 | for true = (sqrt (/ (* 2 (sqrt (- 1 m)) (elliptic-k m)) |
---|
1139 | (float-pi m))) |
---|
1140 | for result = (check-accuracy 204 t3 true) |
---|
1141 | when result |
---|
1142 | append (list (list (list k m) result))) |
---|
1143 | nil) |
---|