source: rt-tests.lisp @ 104efd

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

VALUE-OR-TINY returned a value that was too tiny.

  • qd-gamma.lisp::
    • Return sqrt(least-positive-normalized-double) instead of least-positive-normalized-double.
  • rt-tests.lisp::
    • Add test for this case.
  • Property mode set to 100644
File size: 41.5 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 47 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(rt:deftest oct.elliptic-pi.n0.q
1125    ;; Tests for random values for phi in [0, pi/2] and n in [0, 1]
1126    (loop for k from 0 below 100
1127       for phi = (random (/ +pi+ 2))
1128       for n = (random #q1)
1129       for epi = (elliptic-pi n phi 0)
1130       for true = (/ (atan (* (tan phi) (sqrt (- 1 n))))
1131                     (sqrt (- 1 n)))
1132       for result = (check-accuracy 204 epi true)
1133       unless (eq nil result)
1134       append (list (list (list k n phi) result)))
1135  nil)
1136
1137(rt:deftest oct.elliptic-pi.n1.q
1138    (loop for k from 0 below 100
1139       for phi = (random (/ +pi+ 2))
1140       for epi = (elliptic-pi 1 phi 0)
1141       for true = (tan phi)
1142       for result = (check-accuracy 194 epi true)
1143       unless (eq nil result)
1144       append (list (list (list k phi) result)))
1145  nil)
1146
1147(rt:deftest oct.elliptic-pi.n2.q
1148    (loop for k from 0 below 100
1149       for phi = (random (/ +pi+ 2))
1150       for n = (+ #q1 (random #q1))
1151       for epi = (elliptic-pi n phi 0)
1152       for true = (/ (atanh (* (tan phi) (sqrt (- n 1))))
1153                     (sqrt (- n 1)))
1154       for result = (check-accuracy 202 epi true)
1155       ;; Not sure if this formula holds when atanh gives a complex
1156       ;; result.  Wolfram doesn't say
1157       when (and (not (complexp true)) result)
1158       append (list (list (list k n phi) result)))
1159  nil)
1160
1161;; Tests for theta functions.
1162
1163(rt:deftest oct.theta3.1.d
1164    ;; A&S 16.38.5
1165    ;; sqrt(2*K/%pi) = theta3(0,q)
1166    (loop for k from 0 below 100
1167       for m = (random 1d0)
1168       for t3 = (elliptic-theta-3 0 (elliptic-nome m))
1169       for true = (sqrt (/ (* 2 (elliptic-k m)) (float-pi m)))
1170       for result = (check-accuracy 50.5 t3 true)
1171       when result
1172       append (list (list (list k m) result)))
1173  nil)
1174
1175(rt:deftest oct.theta3.1.q
1176    ;; A&S 16.38.5
1177    ;; sqrt(2*K/%pi) = theta3(0,q)
1178    (loop for k from 0 below 100
1179       for m = (random #q1)
1180       for t3 = (elliptic-theta-3 0 (elliptic-nome m))
1181       for true = (sqrt (/ (* 2 (elliptic-k m)) (float-pi m)))
1182       for result = (check-accuracy 205.7 t3 true)
1183       when result
1184       append (list (list (list k m) result)))
1185  nil)
1186
1187(rt:deftest oct.theta2.1.d
1188    ;; A&S 16.38.7
1189    ;; sqrt(2*sqrt(m)*K/%pi) = theta2(0,q)
1190    (loop for k from 0 below 100
1191       for m = (random 1d0)
1192       for t3 = (elliptic-theta-2 0 (elliptic-nome m))
1193       for true = (sqrt (/ (* 2 (sqrt m) (elliptic-k m)) (float-pi m)))
1194       for result = (check-accuracy 43.5 t3 true)
1195       when result
1196       append (list (list (list k m) result)))
1197  nil)
1198
1199(rt:deftest oct.theta2.1.q
1200    ;; A&S 16.38.7
1201    ;; sqrt(2*sqrt(m)*K/%pi) = theta2(0,q)
1202    (loop for k from 0 below 100
1203       for m = (random #q1)
1204       for t3 = (elliptic-theta-2 0 (elliptic-nome m))
1205       for true = (sqrt (/ (* 2 (sqrt m) (elliptic-k m)) (float-pi m)))
1206       for result = (check-accuracy 205 t3 true)
1207       when result
1208       append (list (list (list k m) result)))
1209  nil)
1210
1211(rt:deftest oct.theta4.1.d
1212    ;; A&S 16.38.8
1213    ;; sqrt(2*sqrt(1-m)*K/%pi) = theta2(0,q)
1214    (loop for k from 0 below 100
1215       for m = (random 1d0)
1216       for t3 = (elliptic-theta-4 0 (elliptic-nome m))
1217       for true = (sqrt (/ (* 2 (sqrt (- 1 m)) (elliptic-k m))
1218                           (float-pi m)))
1219       for result = (check-accuracy 49 t3 true)
1220       when result
1221       append (list (list (list k m) result)))
1222  nil)
1223
1224(rt:deftest oct.theta4.1.q
1225    ;; A&S 16.38.8
1226    ;; sqrt(2*sqrt(1-m)*K/%pi) = theta2(0,q)
1227    (loop for k from 0 below 100
1228       for m = (random #q1)
1229       for t3 = (elliptic-theta-4 0 (elliptic-nome m))
1230       for true = (sqrt (/ (* 2 (sqrt (- 1 m)) (elliptic-k m))
1231                           (float-pi m)))
1232       for result = (check-accuracy 204 t3 true)
1233       when result
1234       append (list (list (list k m) result)))
1235  nil)
1236
1237(rt:deftest lentz
1238    ;; This isn't really a test of cf-incomplete-gamma.  It's a test
1239    ;; that Lentz's algorithm works in this case.  For these args,
1240    ;; cf-incomplete-gamma used to generate an overflow or division by
1241    ;; zero because value-or-tiny was too tiny.
1242    (let ((g (cf-incomplete-gamma 3d0 5d0))
1243          (true (- 2 (* 37 (exp -5d0)))))
1244      (check-accuracy 53 g true))
1245  nil)
1246
1247(rt:deftest gamma.1.d
1248    (let ((g (gamma 0.5d0))
1249          (true (sqrt pi)))
1250      ;; This should give full accuracy but doesn't.
1251      (check-accuracy 51 g true))
1252  nil)
1253
1254(rt:deftest gamma.1.q
1255    (let ((g (gamma #q0.5))
1256          (true (sqrt +pi+)))
1257      ;; This should give full accuracy but doesn't.
1258      (check-accuracy 197 g true))
1259  nil)
1260
1261(rt:deftest gamma.2.d
1262    (loop for k from 0 below 100
1263       for y = (+ 1 (random 100d0))
1264       for g = (abs (gamma (complex 0 y)))
1265       for true = (sqrt (/ pi y (sinh (* pi y))))
1266       for result = (check-accuracy 44 g true)
1267       when result
1268       append (list (list (list k y) result)))
1269  nil)
1270
1271(rt:deftest gamma.2.q
1272    (loop for k from 0 below 100
1273       for y = (+ 1 (random #q100))
1274       for g = (abs (gamma (complex 0 y)))
1275       for true = (sqrt (/ +pi+ y (sinh (* +pi+ y))))
1276       for result = (check-accuracy 196 g true)
1277       when result
1278       append (list (list (list k y) result)))
1279  nil)
1280
1281(rt:deftest gamma.3.d
1282    (loop for k from 0 below 100
1283       for y = (+ 1 (random 100d0))
1284       for g = (abs (gamma (complex 1/2 y)))
1285       for true = (sqrt (/ pi (cosh (* pi y))))
1286       for result = (check-accuracy 44 g true)
1287       when result
1288       append (list (list (list k y) result)))
1289  nil)
1290
1291(rt:deftest gamma.3.q
1292    (loop for k from 0 below 100
1293       for y = (+ 1 (random #q100))
1294       for g = (abs (gamma (complex 1/2 y)))
1295       for true = (sqrt (/ +pi+ (cosh (* +pi+ y))))
1296       for result = (check-accuracy 196 g true)
1297       when result
1298       append (list (list (list k y) result)))
1299  nil)
1300
1301;; gamma_incomplete(2,z) = integrate(t*exp(-t), t, z, inf)
1302;;                       = (z+1)*exp(-z)
1303(rt:deftest gamma-incomplete-tail.1.d
1304    (let* ((z 5d0)
1305           (gi (incomplete-gamma-tail 2 z))
1306           (true (* (+ z 1) (exp (- z)))))
1307      (check-accuracy 52 gi true))
1308  nil)
1309
1310(rt:deftest gamma-incomplete-tail.2.d
1311    (let* ((z #c(1 5d0))
1312           (gi (incomplete-gamma-tail 2 z))
1313           (true (* (+ z 1) (exp (- z)))))
1314      (check-accuracy 50 gi true))
1315  nil)
1316
1317(rt:deftest gamma-incomplete-tail.1.q
1318    (let* ((z #q5)
1319           (gi (incomplete-gamma-tail 2 z))
1320           (true (* (+ z 1) (exp (- z)))))
1321      (check-accuracy 207 gi true))
1322  nil)
1323
1324(rt:deftest gamma-incomplete-tail.2.q
1325    (let* ((z #q(1 5))
1326           (gi (incomplete-gamma-tail 2 z))
1327           (true (* (+ z 1) (exp (- z)))))
1328      (check-accuracy 206 gi true))
1329  nil)
1330
1331(rt:deftest gamma-incomplete-tail.3.d
1332    (let* ((z -5d0)
1333           (gi (incomplete-gamma-tail 2 z))
1334           (true (* (+ z 1) (exp (- z)))))
1335      (check-accuracy 50 gi true))
1336  nil)
1337
1338(rt:deftest gamma-incomplete-tail.3.q
1339    (let* ((z #q-5)
1340           (gi (incomplete-gamma-tail 2 z))
1341           (true (* (+ z 1) (exp (- z)))))
1342      (check-accuracy 206 gi true))
1343  nil)
1344
1345;; See http://www.wolframalpha.com/input/?i=Gamma[1%2F2%2C-100%2Bi%2F%2810^10%29]
1346
1347(rt:deftest gamma-incomplete-tail.4.q
1348    (let* ((z #q(#q-100 #q1q-10))
1349           (gi (incomplete-gamma-tail 1/2 z))
1350           (true #q(#q-2.68811714181613544840818982228135651231579313476267430888499241497530341422025007816745898370049200133136q32
1351                    #q-2.70176456134384383878883307528351227886457379834795655467745609829086928772079968479767583764284583465328q42)))
1352      (check-accuracy 205 gi true))
1353  nil)
1354
1355
1356;; Fresnel integrals.
1357;;
1358;; For x small, Fresnel
1359;;
1360;;  S(z) = %pi/6*z^3*(1 - %pi^2*z^4/56 + %pi^4*z^8/2040 - ...)
1361;;
1362(defun fresnel-s-series (z)
1363  (let* ((fpi (float-pi z))
1364         (z^3 (expt z 3))
1365         (z^4 (* z^3 z)))
1366    (* fpi 1/6 z^3
1367       (+ 1 (/ (* fpi fpi z^4)
1368               -56)
1369          (/ (* (expt fpi 4) (expt z^4 2))
1370             7040)))))
1371
1372(rt:deftest fresnel-s.1d
1373    (let* ((z 1d-3)
1374           (s (fresnel-s z))
1375           (true (fresnel-s-series z)))
1376      (check-accuracy 52 s true))
1377  nil)
1378
1379(rt:deftest fresnel-s.2d
1380    (let* ((z #c(1d-3 1d-3))
1381           (s (fresnel-s z))
1382           (true (fresnel-s-series z)))
1383      (check-accuracy 52 s true))
1384  nil)
1385
1386(rt:deftest fresnel-s.1q
1387    (let* ((z #q1q-20)
1388           (s (fresnel-s z))
1389           (true (fresnel-s-series z)))
1390      (check-accuracy 212 s true))
1391  nil)
1392
1393(rt:deftest fresnel-s.2q
1394    (let* ((z #q(#q1q-3 #q1q-3))
1395           (s (fresnel-s z))
1396           (true (fresnel-s-series z)))
1397      (check-accuracy 212 s true))
1398  nil)
1399
1400(rt:deftest psi.1d
1401    (let* ((z 1d0)
1402           (p (psi z))
1403           (true (float (- +%gamma+) 1d0)))
1404      (check-accuracy 52 p true))
1405  nil)
1406
1407(rt:deftest psi.1q
1408    (let* ((z #q1)
1409           (p (psi z))
1410           (true (- +%gamma+)))
1411      (check-accuracy 208 p true))
1412  nil)
1413
1414(rt:deftest psi.2d
1415    (let* ((z (float 4/3 1d0))
1416           (p (psi z))
1417           (true (- 3
1418                    +%gamma+
1419                    (/ +pi+ (* 2 (sqrt #q3)))
1420                    (* 1.5 (log #q3)))))
1421      (check-accuracy 49.8 p true))
1422  nil)
1423
1424(rt:deftest psi.2q
1425    (let* ((z (float 4/3 #q1))
1426           (p (psi z))
1427           (true (- 3
1428                    +%gamma+
1429                    (/ +pi+ (* 2 (sqrt #q3)))
1430                    (* 1.5 (log #q3)))))
1431      (check-accuracy 205 p true))
1432  nil)
1433
1434(rt:deftest psi.3d
1435    (let* ((z (float -1/2 1d0))
1436           (p (psi z))
1437           (true (- 2
1438                    +%gamma+
1439                    (log #q4))))
1440      (check-accuracy 48 p true))
1441  nil)
1442
1443(rt:deftest psi.3q
1444    (let* ((z (float -1/2 #q1))
1445           (p (psi z))
1446           (true (- 2
1447                    +%gamma+
1448                    (log #q4))))
1449      (check-accuracy 204.1 p true))
1450  nil)
1451
1452(rt:deftest expintegral-e.1d
1453    (let* ((z 1d0)
1454           (e (exp-integral-e 0 z))
1455           (true (/ (exp (- z)) z)))
1456      (check-accuracy 53 e true))
1457  nil)
1458
1459(rt:deftest expintegral-e.1q
1460    (let* ((z #q1)
1461           (e (exp-integral-e 0 z))
1462           (true (/ (exp (- z)) z)))
1463      (check-accuracy 212 e true))
1464  nil)
1465
1466(rt:deftest expintegral-e.2d
1467    (let* ((z 15d0)
1468           (e (exp-integral-e 0 z))
1469           (true (/ (exp (- z)) z)))
1470      (check-accuracy 53 e true))
1471  nil)
1472
1473(rt:deftest expintegral-e.2q
1474    (let* ((z #q15)
1475           (e (exp-integral-e 0 z))
1476           (true (/ (exp (- z)) z)))
1477      (check-accuracy 212 e true))
1478  nil)
1479
1480(rt:deftest expintegral-e.3d
1481    (let* ((e (exp-integral-e 2 1d0))
1482           (true 0.14849550677592204791835999d0))
1483      (check-accuracy 47.5 e true))
1484  nil)
1485
1486(rt:deftest expintegral-e.4d
1487    (let* ((x .5d0)
1488           (e (exp-integral-e -2 x))
1489           (true (/ (* (exp (- x)) (+ (* x x x) (* 2 x x) (* 2 x)))
1490                    (expt x 4))))
1491      (check-accuracy 53 e true))
1492  nil)
1493
1494(rt:deftest expintegral-e.4q
1495    (let* ((x #q.5)
1496           (e (exp-integral-e -2 x))
1497           (true (/ (* (exp (- x)) (+ (* x x x) (* 2 x x) (* 2 x)))
1498                    (expt x 4))))
1499      (check-accuracy 210.8 e true))
1500  nil)
1501
1502(rt:deftest expintegral-e.5d
1503    (let* ((x .5d0)
1504           (e (exp-integral-e 2d0 x))
1505           (true #q0.3266438623245530177304015653336378358284946903290101))
1506      (check-accuracy 51.2 e true))
1507  nil)
1508
1509(rt:deftest expintegral-e.5q
1510    (let* ((x #q.5)
1511           (e (exp-integral-e #q2 x))
1512           (true #q0.326643862324553017730401565333637835828494690329010198058745549181386569998611289568))
1513      (check-accuracy 208.4 e true))
1514  nil)
1515
1516(rt:deftest expintegral-e.6d
1517    (let* ((x .5d0)
1518           (e (exp-integral-e 1d0 x))
1519           (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883))
1520      (check-accuracy 53.9 e true))
1521  nil)
1522
1523(rt:deftest expintegral-e.6q
1524    (let* ((x #q.5)
1525           (e (exp-integral-e #q1 x))
1526           (true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883))
1527      (check-accuracy 219.1 e true))
1528  nil)
1529
Note: See TracBrowser for help on using the repository browser.