source: rt-tests.lisp @ 5566bc

Last change on this file since 5566bc was 5566bc, checked in by Raymond Toy <toy.raymond@…>, 3 years ago

Add more tests for Bessel J. Not all of them pass.

  • Property mode set to 100644
File size: 47.2 KB
Line 
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 169.45 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(rt:deftest bessel-j-1/2.q.1.c
1669    (let* ((x #q1.0360263937639582798798376485114581552570020473846457752365459851056q0)
1670           (b (bessel-j #q0.5 x))
1671           (true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1))))))
1672      (check-accuracy 169.45 b true))
1673  nil)
1674
1675;; Bessel J for complex args
1676(rt:deftest bessel-j-complex-arg.d.1
1677    (let ((b (bessel-j 0d0 #c(1d0 1)))
1678          (true #c(0.9376084768060293d0 -0.4965299476091221d0)))
1679      (check-accuracy 50.73 b true))
1680  nil)
1681
1682(rt:deftest bessel-j-complex-arg.d.2
1683    (let ((b (bessel-j 1d0 #c(1d0 1)))
1684          (true #c(0.6141603349229036d0 0.3650280288270878d0)))
1685      (check-accuracy 52.51 b true))
1686  nil)
1687
1688(rt:deftest bessel-j-complex-arg.d.3
1689    (let ((b (bessel-j 2d0 #c(1d0 1)))
1690          (true #c(0.0415798869439621d0 0.2473976415133063d0)))
1691      (check-accuracy 50.41 b true))
1692  nil)
1693
1694(rt:deftest bessel-j-complex-arg.d.4
1695    (let ((b (bessel-j 2.3d0 #c(1d0 1)))
1696          (true #c(-0.0141615213034667d0 0.1677798241687935d0)))
1697      (check-accuracy 48.56 b true))
1698  nil)
1699
1700(rt:deftest bessel-j-complex-arg.d.5
1701    (let ((b (bessel-j -2.3d0 #c(1d0 1)))
1702          (true #c(0.1920598664138632d0 -0.5158676904105332d0)))
1703      (check-accuracy 50.97 b true))
1704  nil)
1705
1706(rt:deftest bessel-j-1/2-complex.d.1
1707    (loop for k from 0 below 10
1708          for x = (complex (random (/ pi 2))
1709                           (random (/ pi 2)))
1710          for b = (bessel-j 0.5d0 x)
1711          for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi)))
1712          for result = (check-accuracy 49.8 b true)
1713          when result
1714            append (list (list (list k x) result)))
1715  nil)
1716
1717(rt:deftest bessel-j-1/2-complex.q.1
1718    (loop for k from 0 below 10
1719          for x = (complex (random (/ (float-pi #q1) 2))
1720                           (random (/ (float-pi #q1) 2)))
1721          for b = (bessel-j #q0.5 x)
1722          for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1))))
1723          for result = (check-accuracy 212 b true)
1724          when result
1725            append (list (list (list k x) result)))
1726  nil)
Note: See TracBrowser for help on using the repository browser.