source: rt-tests.lisp @ 4b332e

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

Fix bug in psi for -n/2 for n odd which was causing an overflow. Add
tests too.

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