Context Navigation

root/rt-tests.lisp @ 5566bc397937c2e8979ee6847e8f59f279f1f643

Revision 5566bc397937c2e8979ee6847e8f59f279f1f643, 47.2 KB (checked in by Raymond Toy <toy.raymond@…>, 14 months ago)
Add more tests for Bessel J. Not all of them pass.
Property mode set to `100644`

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	;; = 3sqrt(%pi)gamma(3/4)/gamma(1/4)
952	(let ((rd (carlson-rd 0d0 2d0 1d0))
953	(true 1.7972103521033883d0))
954	(check-accuracy 51 rd true))
955	nil)
956
957	(rt:deftest oct.carlson-rd.1q
958	(let ((rd (carlson-rd #q0 #q2 #q1))
959	(true #q1.797210352103388311159883738420485817340818994823477337395512429419599q0))
960	(check-accuracy 212 rd true))
961	nil)
962
963	;; Test some of the contagion stuff.
964
965	(rt:deftest oct.carlson-rf.contagion.1
966	;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
967	;; = 1/4*beta(1/2,1/2)
968	;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
969	(let ((rf (carlson-rf 0 2 1))
970	(true 1.31102877714605990523241979494d0))
971	(check-accuracy 23 rf true))
972	nil)
973
974	(rt:deftest oct.carlson-rf.contagion.1d
975	;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
976	;; = 1/4*beta(1/2,1/2)
977	;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
978	(let ((rf (carlson-rf 0d0 2 1))
979	(true 1.31102877714605990523241979494d0))
980	(check-accuracy 53 rf true))
981	nil)
982
983	(rt:deftest oct.carlson-rf.contagion.2d
984	;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
985	;; = 1/4*beta(1/2,1/2)
986	;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
987	(let ((rf (carlson-rf 0 2d0 1))
988	(true 1.31102877714605990523241979494d0))
989	(check-accuracy 53 rf true))
990	nil)
991
992	(rt:deftest oct.carlson-rf.contagion.3d
993	;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
994	;; = 1/4*beta(1/2,1/2)
995	;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
996	(let ((rf (carlson-rf 0 2 1d0))
997	(true 1.31102877714605990523241979494d0))
998	(check-accuracy 53 rf true))
999	nil)
1000
1001	(rt:deftest oct.carlson-rf.contagion.1q
1002	;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
1003	;; = 1/4*beta(1/2,1/2)
1004	;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
1005	(let ((rf (carlson-rf #q0q0 2 1))
1006	(true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0))
1007	(check-accuracy 212 rf true))
1008	nil)
1009
1010	(rt:deftest oct.carlson-rf.contagion.2q
1011	;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
1012	;; = 1/4*beta(1/2,1/2)
1013	;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
1014	(let ((rf (carlson-rf 0 #q2q0 1))
1015	(true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0))
1016	(check-accuracy 212 rf true))
1017	nil)
1018
1019	(rt:deftest oct.carlson-rf.contagion.3q
1020	;; Rf(0,2,1) = integrate(1/sqrt(1-s^4), s, 0 ,1)
1021	;; = 1/4*beta(1/2,1/2)
1022	;; = sqrt(%pi)/4*gamma(1/4)/gamma(3/4)
1023	(let ((rf (carlson-rf 0 2 #q1q0))
1024	(true #q1.311028777146059905232419794945559706841377475715811581408410851900395q0))
1025	(check-accuracy 212 rf true))
1026	nil)
1027
1028	;; Elliptic integral of the third kind
1029
1030	;; elliptic-pi(0,phi,m) = elliptic-f(phi, m)
1031	(rt:deftest oct.elliptic-pi.1d
1032	(loop for k from 0 to 100
1033	for phi = (random (/ pi 2))
1034	for m = (random 1d0)
1035	for epi = (elliptic-pi 0 phi m)
1036	for ef = (elliptic-f phi m)
1037	for result = (check-accuracy 48 epi ef)
1038	unless (eq nil result)
1039	append (list (list phi m) result))
1040	nil)
1041
1042	(rt:deftest oct.elliptic-pi.1q
1043	(loop for k from 0 below 100
1044	for phi = (random (/ +pi+ 2))
1045	for m = (random #q1)
1046	for epi = (elliptic-pi 0 phi m)
1047	for ef = (elliptic-f phi m)
1048	for result = (check-accuracy 53 epi ef)
1049	unless (eq nil result)
1050	append (list (list phi m) result))
1051	nil)
1052
1053	;; DLMF 19.6.3
1054	;;
1055	;; PI(n; pi/2 \| 0) = pi/(2*sqrt(1-n))
1056	(rt:deftest oct.elliptic-pi.19.6.3.d
1057	(loop for k from 0 below 100
1058	for n = (random 1d0)
1059	for epi = (elliptic-pi n (/ pi 2) 0)
1060	for true = (/ pi (* 2 (sqrt (- 1 n))))
1061	for result = (check-accuracy 47 epi true)
1062	unless (eq nil result)
1063	append (list (list (list k n) result)))
1064	nil)
1065
1066	(rt:deftest oct.elliptic-pi.19.6.3.q
1067	(loop for k from 0 below 100
1068	for n = (random #q1)
1069	for epi = (elliptic-pi n (/ (float-pi n) 2) 0)
1070	for true = (/ (float-pi n) (* 2 (sqrt (- 1 n))))
1071	for result = (check-accuracy 208 epi true)
1072	unless (eq nil result)
1073	append (list (list (list k n) result)))
1074	nil)
1075
1076	;; elliptic-pi(n, phi, 0) =
1077	;; atan(sqrt(1-n)*tan(phi))/sqrt(1-n) n < 1
1078	;; atanh(sqrt(n-1)*tan(phi))/sqrt(n-1) n > 1
1079	;; tan(phi) n = 1
1080	;;
1081	;; These are easy to derive if you look at the integral:
1082	;;
1083	;; ellipti-pi(n, phi, 0) = integrate(1/(1-n*sin(t)^2), t, 0, phi)
1084	;;
1085	;; and this can be easily integrated to give the above expressions for
1086	;; the different values of n.
1087	(rt:deftest oct.elliptic-pi.n0.d
1088	;; Tests for random values for phi in [0, pi/2] and n in [0, 1]
1089	(loop for k from 0 below 100
1090	for phi = (random (/ pi 2))
1091	for n = (random 1d0)
1092	for epi = (elliptic-pi n phi 0)
1093	for true = (/ (atan (* (tan phi) (sqrt (- 1 n))))
1094	(sqrt (- 1 n)))
1095	for result = (check-accuracy 46.5 epi true)
1096	unless (eq nil result)
1097	append (list (list (list k n phi) result)))
1098	nil)
1099
1100	(rt:deftest oct.elliptic-pi.n1.d
1101	(loop for k from 0 below 100
1102	for phi = (random (/ pi 2))
1103	for epi = (elliptic-pi 1 phi 0)
1104	for true = (tan phi)
1105	for result = (check-accuracy 34.5 epi true)
1106	unless (eq nil result)
1107	append (list (list (list k phi) result)))
1108	nil)
1109
1110	(rt:deftest oct.elliptic-pi.n2.d
1111	(loop for k from 0 below 100
1112	for phi = (random (/ pi 2))
1113	for n = (+ 1d0 (random 100d0))
1114	for epi = (elliptic-pi n phi 0)
1115	for true = (/ (atanh (* (tan phi) (sqrt (- n 1))))
1116	(sqrt (- n 1)))
1117	for result = (check-accuracy 45.85 epi true)
1118	;; Not sure if this formula holds when atanh gives a complex
1119	;; result. Wolfram doesn't say
1120	when (and (not (complexp true)) result)
1121	append (list (list (list k n phi) result)))
1122	nil)
1123
1124	;; Failed test case:
1125	;; ((89 66.68551748022054d0 0.12266024127708153d0)
1126	;; (45.868614757480834d0 47 0.47787458521306514d0
1127	;; 0.4778745852130726d0))
1128	;; New threshold is 45.85 bits.
1129	(rt:deftest oct.elliptic-pi.n2.d-1
1130	(let* ((n 66.68551748022054d0)
1131	(phi 0.12266024127708153d0)
1132	(epi (elliptic-pi n phi 0))
1133	(true (/ (atanh (* (tan phi) (sqrt (- n 1))))
1134	(sqrt (- n 1)))))
1135	(check-accuracy 45.8686d0 epi true))
1136	nil)
1137
1138
1139	(rt:deftest oct.elliptic-pi.n0.q
1140	;; Tests for random values for phi in [0, pi/2] and n in [0, 1]
1141	(loop for k from 0 below 100
1142	for phi = (random (/ +pi+ 2))
1143	for n = (random #q1)
1144	for epi = (elliptic-pi n phi 0)
1145	for true = (/ (atan (* (tan phi) (sqrt (- 1 n))))
1146	(sqrt (- 1 n)))
1147	for result = (check-accuracy 204 epi true)
1148	unless (eq nil result)
1149	append (list (list (list k n phi) result)))
1150	nil)
1151
1152	(rt:deftest oct.elliptic-pi.n1.q
1153	(loop for k from 0 below 100
1154	for phi = (random (/ +pi+ 2))
1155	for epi = (elliptic-pi 1 phi 0)
1156	for true = (tan phi)
1157	for result = (check-accuracy 194 epi true)
1158	unless (eq nil result)
1159	append (list (list (list k phi) result)))
1160	nil)
1161
1162	(rt:deftest oct.elliptic-pi.n2.q
1163	(loop for k from 0 below 100
1164	for phi = (random (/ +pi+ 2))
1165	for n = (+ #q1 (random #q1))
1166	for epi = (elliptic-pi n phi 0)
1167	for true = (/ (atanh (* (tan phi) (sqrt (- n 1))))
1168	(sqrt (- n 1)))
1169	for result = (check-accuracy 202 epi true)
1170	;; Not sure if this formula holds when atanh gives a complex
1171	;; result. Wolfram doesn't say
1172	when (and (not (complexp true)) result)
1173	append (list (list (list k n phi) result)))
1174	nil)
1175
1176	;; Tests for theta functions.
1177
1178	(rt:deftest oct.theta3.1.d
1179	;; A&S 16.38.5
1180	;; sqrt(2*K/%pi) = theta3(0,q)
1181	(loop for k from 0 below 100
1182	for m = (random 1d0)
1183	for t3 = (elliptic-theta-3 0 (elliptic-nome m))
1184	for true = (sqrt (/ (* 2 (elliptic-k m)) (float-pi m)))
1185	for result = (check-accuracy 50.5 t3 true)
1186	when result
1187	append (list (list (list k m) result)))
1188	nil)
1189
1190	(rt:deftest oct.theta3.1.q
1191	;; A&S 16.38.5
1192	;; sqrt(2*K/%pi) = theta3(0,q)
1193	(loop for k from 0 below 100
1194	for m = (random #q1)
1195	for t3 = (elliptic-theta-3 0 (elliptic-nome m))
1196	for true = (sqrt (/ (* 2 (elliptic-k m)) (float-pi m)))
1197	for result = (check-accuracy 205.7 t3 true)
1198	when result
1199	append (list (list (list k m) result)))
1200	nil)
1201
1202	(rt:deftest oct.theta2.1.d
1203	;; A&S 16.38.7
1204	;; sqrt(2sqrt(m)K/%pi) = theta2(0,q)
1205	(loop for k from 0 below 100
1206	for m = (random 1d0)
1207	for t3 = (elliptic-theta-2 0 (elliptic-nome m))
1208	for true = (sqrt (/ (* 2 (sqrt m) (elliptic-k m)) (float-pi m)))
1209	for result = (check-accuracy 43.5 t3 true)
1210	when result
1211	append (list (list (list k m) result)))
1212	nil)
1213
1214	(rt:deftest oct.theta2.1.q
1215	;; A&S 16.38.7
1216	;; sqrt(2sqrt(m)K/%pi) = theta2(0,q)
1217	(loop for k from 0 below 100
1218	for m = (random #q1)
1219	for t3 = (elliptic-theta-2 0 (elliptic-nome m))
1220	for true = (sqrt (/ (* 2 (sqrt m) (elliptic-k m)) (float-pi m)))
1221	for result = (check-accuracy 205 t3 true)
1222	when result
1223	append (list (list (list k m) result)))
1224	nil)
1225
1226	(rt:deftest oct.theta4.1.d
1227	;; A&S 16.38.8
1228	;; sqrt(2sqrt(1-m)K/%pi) = theta2(0,q)
1229	(loop for k from 0 below 100
1230	for m = (random 1d0)
1231	for t3 = (elliptic-theta-4 0 (elliptic-nome m))
1232	for true = (sqrt (/ (* 2 (sqrt (- 1 m)) (elliptic-k m))
1233	(float-pi m)))
1234	for result = (check-accuracy 49 t3 true)
1235	when result
1236	append (list (list (list k m) result)))
1237	nil)
1238
1239	(rt:deftest oct.theta4.1.q
1240	;; A&S 16.38.8
1241	;; sqrt(2sqrt(1-m)K/%pi) = theta2(0,q)
1242	(loop for k from 0 below 100
1243	for m = (random #q1)
1244	for t3 = (elliptic-theta-4 0 (elliptic-nome m))
1245	for true = (sqrt (/ (* 2 (sqrt (- 1 m)) (elliptic-k m))
1246	(float-pi m)))
1247	for result = (check-accuracy 204 t3 true)
1248	when result
1249	append (list (list (list k m) result)))
1250	nil)
1251
1252	(rt:deftest lentz
1253	;; This isn't really a test of cf-incomplete-gamma. It's a test
1254	;; that Lentz's algorithm works in this case. For these args,
1255	;; cf-incomplete-gamma used to generate an overflow or division by
1256	;; zero because value-or-tiny was too tiny.
1257	(let ((g (cf-incomplete-gamma 3d0 5d0))
1258	(true (- 2 (* 37 (exp -5d0)))))
1259	(check-accuracy 51.2 g true))
1260	nil)
1261
1262	(rt:deftest gamma.1.d
1263	(let ((g (gamma 0.5d0))
1264	(true (sqrt pi)))
1265	;; This should give full accuracy but doesn't.
1266	(check-accuracy 51 g true))
1267	nil)
1268
1269	(rt:deftest gamma.1.q
1270	(let ((g (gamma #q0.5))
1271	(true (sqrt +pi+)))
1272	;; This should give full accuracy but doesn't.
1273	(check-accuracy 197 g true))
1274	nil)
1275
1276	(rt:deftest gamma.2.d
1277	(loop for k from 0 below 100
1278	for y = (+ 1 (random 100d0))
1279	for g = (abs (gamma (complex 0 y)))
1280	for true = (sqrt (/ pi y (sinh (* pi y))))
1281	for result = (check-accuracy 44 g true)
1282	when result
1283	append (list (list (list k y) result)))
1284	nil)
1285
1286	(rt:deftest gamma.2.q
1287	(loop for k from 0 below 100
1288	for y = (+ 1 (random #q100))
1289	for g = (abs (gamma (complex 0 y)))
1290	for true = (sqrt (/ +pi+ y (sinh (* +pi+ y))))
1291	for result = (check-accuracy 196 g true)
1292	when result
1293	append (list (list (list k y) result)))
1294	nil)
1295
1296	(rt:deftest gamma.3.d
1297	(loop for k from 0 below 100
1298	for y = (+ 1 (random 100d0))
1299	for g = (abs (gamma (complex 1/2 y)))
1300	for true = (sqrt (/ pi (cosh (* pi y))))
1301	for result = (check-accuracy 44 g true)
1302	when result
1303	append (list (list (list k y) result)))
1304	nil)
1305
1306	(rt:deftest gamma.3.q
1307	(loop for k from 0 below 100
1308	for y = (+ 1 (random #q100))
1309	for g = (abs (gamma (complex 1/2 y)))
1310	for true = (sqrt (/ +pi+ (cosh (* +pi+ y))))
1311	for result = (check-accuracy 196 g true)
1312	when result
1313	append (list (list (list k y) result)))
1314	nil)
1315
1316	;; gamma_incomplete(2,z) = integrate(t*exp(-t), t, z, inf)
1317	;; = (z+1)*exp(-z)
1318	(rt:deftest gamma-incomplete-tail.1.d
1319	(let* ((z 5d0)
1320	(gi (incomplete-gamma-tail 2 z))
1321	(true (* (+ z 1) (exp (- z)))))
1322	(check-accuracy 52 gi true))
1323	nil)
1324
1325	(rt:deftest gamma-incomplete-tail.2.d
1326	(let* ((z #c(1 5d0))
1327	(gi (incomplete-gamma-tail 2 z))
1328	(true (* (+ z 1) (exp (- z)))))
1329	(check-accuracy 50 gi true))
1330	nil)
1331
1332	(rt:deftest gamma-incomplete-tail.1.q
1333	(let* ((z #q5)
1334	(gi (incomplete-gamma-tail 2 z))
1335	(true (* (+ z 1) (exp (- z)))))
1336	(check-accuracy 207 gi true))
1337	nil)
1338
1339	(rt:deftest gamma-incomplete-tail.2.q
1340	(let* ((z #q(1 5))
1341	(gi (incomplete-gamma-tail 2 z))
1342	(true (* (+ z 1) (exp (- z)))))
1343	(check-accuracy 206 gi true))
1344	nil)
1345
1346	(rt:deftest gamma-incomplete-tail.3.d
1347	(let* ((z -5d0)
1348	(gi (incomplete-gamma-tail 2 z))
1349	(true (* (+ z 1) (exp (- z)))))
1350	(check-accuracy 50 gi true))
1351	nil)
1352
1353	(rt:deftest gamma-incomplete-tail.3.q
1354	(let* ((z #q-5)
1355	(gi (incomplete-gamma-tail 2 z))
1356	(true (* (+ z 1) (exp (- z)))))
1357	(check-accuracy 206 gi true))
1358	nil)
1359
1360	;; See http://www.wolframalpha.com/input/?i=Gamma[1%2F2%2C-100%2Bi%2F%2810^10%29]
1361
1362	(rt:deftest gamma-incomplete-tail.4.q
1363	(let* ((z #q(#q-100 #q1q-10))
1364	(gi (incomplete-gamma-tail 1/2 z))
1365	(true #q(#q-2.68811714181613544840818982228135651231579313476267430888499241497530341422025007816745898370049200133136q32
1366	#q-2.70176456134384383878883307528351227886457379834795655467745609829086928772079968479767583764284583465328q42)))
1367	(check-accuracy 205 gi true))
1368	nil)
1369
1370
1371	;; Fresnel integrals.
1372	;;
1373	;; For x small, Fresnel
1374	;;
1375	;; S(z) = %pi/6z^3(1 - %pi^2z^4/56 + %pi^4z^8/2040 - ...)
1376	;;
1377	(defun fresnel-s-series (z)
1378	(let* ((fpi (float-pi z))
1379	(z^3 (expt z 3))
1380	(z^4 (* z^3 z)))
1381	(* fpi 1/6 z^3
1382	(+ 1 (/ (* fpi fpi z^4)
1383	-56)
1384	(/ (* (expt fpi 4) (expt z^4 2))
1385	7040)))))
1386
1387	(rt:deftest fresnel-s.1d
1388	(let* ((z 1d-3)
1389	(s (fresnel-s z))
1390	(true (fresnel-s-series z)))
1391	(check-accuracy 52 s true))
1392	nil)
1393
1394	(rt:deftest fresnel-s.2d
1395	(let* ((z #c(1d-3 1d-3))
1396	(s (fresnel-s z))
1397	(true (fresnel-s-series z)))
1398	(check-accuracy 52 s true))
1399	nil)
1400
1401	(rt:deftest fresnel-s.1q
1402	(let* ((z #q1q-20)
1403	(s (fresnel-s z))
1404	(true (fresnel-s-series z)))
1405	(check-accuracy 212 s true))
1406	nil)
1407
1408	(rt:deftest fresnel-s.2q
1409	(let* ((z #q(#q1q-3 #q1q-3))
1410	(s (fresnel-s z))
1411	(true (fresnel-s-series z)))
1412	(check-accuracy 212 s true))
1413	nil)
1414
1415	(rt:deftest psi.1d
1416	(let* ((z 1d0)
1417	(p (psi z))
1418	(true (float (- +%gamma+) 1d0)))
1419	(check-accuracy 52 p true))
1420	nil)
1421
1422	(rt:deftest psi.1q
1423	(let* ((z #q1)
1424	(p (psi z))
1425	(true (- +%gamma+)))
1426	(check-accuracy 208 p true))
1427	nil)
1428
1429	(rt:deftest psi.2d
1430	(let* ((z (float 4/3 1d0))
1431	(p (psi z))
1432	(true (- 3
1433	+%gamma+
1434	(/ +pi+ (* 2 (sqrt #q3)))
1435	(* 1.5 (log #q3)))))
1436	(check-accuracy 49.8 p true))
1437	nil)
1438
1439	(rt:deftest psi.2q
1440	(let* ((z (float 4/3 #q1))
1441	(p (psi z))
1442	(true (- 3
1443	+%gamma+
1444	(/ +pi+ (* 2 (sqrt #q3)))
1445	(* 1.5 (log #q3)))))
1446	(check-accuracy 205 p true))
1447	nil)
1448
1449	(rt:deftest psi.3d
1450	(let* ((z (float -1/2 1d0))
1451	(p (psi z))
1452	(true (- 2
1453	+%gamma+
1454	(log #q4))))
1455	(check-accuracy 48 p true))
1456	nil)
1457
1458	(rt:deftest psi.3q
1459	(let* ((z (float -1/2 #q1))
1460	(p (psi z))
1461	(true (- 2
1462	+%gamma+
1463	(log #q4))))
1464	(check-accuracy 204.1 p true))
1465	nil)
1466
1467	(rt:deftest expintegral-e.1d
1468	(let* ((z 1d0)
1469	(e (exp-integral-e 0 z))
1470	(true (/ (exp (- z)) z)))
1471	(check-accuracy 53 e true))
1472	nil)
1473
1474	(rt:deftest expintegral-e.1q
1475	(let* ((z #q1)
1476	(e (exp-integral-e 0 z))
1477	(true (/ (exp (- z)) z)))
1478	(check-accuracy 212 e true))
1479	nil)
1480
1481	(rt:deftest expintegral-e.2d
1482	(let* ((z 15d0)
1483	(e (exp-integral-e 0 z))
1484	(true (/ (exp (- z)) z)))
1485	(check-accuracy 53 e true))
1486	nil)
1487
1488	(rt:deftest expintegral-e.2q
1489	(let* ((z #q15)
1490	(e (exp-integral-e 0 z))
1491	(true (/ (exp (- z)) z)))
1492	(check-accuracy 212 e true))
1493	nil)
1494
1495	(rt:deftest expintegral-e.3d
1496	(let* ((e (exp-integral-e 2 1d0))
1497	(true 0.14849550677592204791835999d0))
1498	(check-accuracy 47.5 e true))
1499	nil)
1500
1501	(rt:deftest expintegral-e.4d
1502	(let* ((x .5d0)
1503	(e (exp-integral-e -2 x))
1504	(true (/ (* (exp (- x)) (+ (* x x x) (* 2 x x) (* 2 x)))
1505	(expt x 4))))
1506	(check-accuracy 53 e true))
1507	nil)
1508
1509	(rt:deftest expintegral-e.4q
1510	(let* ((x #q.5)
1511	(e (exp-integral-e -2 x))
1512	(true (/ (* (exp (- x)) (+ (* x x x) (* 2 x x) (* 2 x)))
1513	(expt x 4))))
1514	(check-accuracy 210.8 e true))
1515	nil)
1516
1517	(rt:deftest expintegral-e.5d
1518	(let* ((x .5d0)
1519	(e (exp-integral-e 2d0 x))
1520	(true #q0.3266438623245530177304015653336378358284946903290101))
1521	(check-accuracy 51.2 e true))
1522	nil)
1523
1524	(rt:deftest expintegral-e.5q
1525	(let* ((x #q.5)
1526	(e (exp-integral-e #q2 x))
1527	(true #q0.326643862324553017730401565333637835828494690329010198058745549181386569998611289568))
1528	(check-accuracy 208.4 e true))
1529	nil)
1530
1531	(rt:deftest expintegral-e.6d
1532	(let* ((x .5d0)
1533	(e (exp-integral-e 1d0 x))
1534	(true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883))
1535	(check-accuracy 53.9 e true))
1536	nil)
1537
1538	(rt:deftest expintegral-e.6q
1539	(let* ((x #q.5)
1540	(e (exp-integral-e #q1 x))
1541	(true #q0.55977359477616081174679593931508523522684689031635351524829321910733989883))
1542	(check-accuracy 219.1 e true))
1543	nil)
1544
1545
1546	;; Bessel J tests for negative order
1547	(rt:deftest bessel-j.neg-order.d.1
1548	(let ((b (bessel-j -1d0 2d0))
1549	(true -0.5767248077568734d0))
1550	(check-accuracy 50.2 b true))
1551	nil)
1552
1553	(rt:deftest bessel-j.neg-order.d.2
1554	(let ((b (bessel-j -1d0 1.5d0))
1555	(true -0.5579365079100996d0))
1556	(check-accuracy 50.5 b true))
1557	nil)
1558
1559	(rt:deftest bessel-j.neg-order.d.3
1560	(let ((b (bessel-j -1.5d0 2d0))
1561	(true -0.3956232813587035d0))
1562	(check-accuracy 50.59 b true))
1563	nil)
1564
1565	(rt:deftest bessel-j.neg-order.d.4
1566	(let ((b (bessel-j -1.8d0 1.5d0))
1567	(true -0.251327217627129314d0))
1568	(check-accuracy 49.98 b true))
1569	nil)
1570
1571	(rt:deftest bessel-j.neg-order.d.5
1572	(let ((b (bessel-j -2d0 1.5d0))
1573	(true 0.2320876721442147d0))
1574	(check-accuracy 51.89 b true))
1575	nil)
1576
1577	(rt:deftest bessel-j.neg-order.d.6
1578	(let ((b (bessel-j -2.5d0 1.5d0))
1579	(true 1.315037204805194d0))
1580	(check-accuracy 52.37 b true))
1581	nil)
1582
1583	(rt:deftest bessel-j.neg-order.d.7
1584	(let ((b (bessel-j -2.3d0 1.5d0))
1585	(true 1.012178926325313d0))
1586	(check-accuracy 50.01 b true))
1587	nil)
1588
1589	;; Bessel-J tests for positive order
1590	(rt:deftest bessel-j.pos-order.d.1
1591	(let ((b (bessel-j 1.5d0 1d0))
1592	(true 0.2402978391234270d0))
1593	(check-accuracy 51.83 b true))
1594	nil)
1595
1596	(rt:deftest bessel-j.pos-order.d.2
1597	(let ((b (bessel-j 1.8d0 1d0))
1598	(true 0.1564953153109239d0))
1599	(check-accuracy 51.97 b true))
1600	nil)
1601
1602	(rt:deftest bessel-j.pos-order.d.3
1603	(let ((b (bessel-j 2d0 1d0))
1604	(true 0.1149034849319005d0))
1605	(check-accuracy 51.87 b true))
1606	nil)
1607
1608	(rt:deftest bessel-j.pos-order.d.4
1609	(let ((b (bessel-j 2.5d0 1d0))
1610	(true 0.04949681022847794d0))
1611	(check-accuracy 47.17 b true))
1612	nil)
1613
1614	(rt:deftest bessel-j.pos-order.d.5
1615	(let ((b (bessel-j -2d0 1.5d0))
1616	(true 0.2320876721442147d0))
1617	(check-accuracy 51.89 b true))
1618	nil)
1619
1620	;; Bessel J for half integer order and real args
1621	(rt:deftest bessel-j-1/2.d.1
1622	(loop for k from 0 below 100
1623	;; x in [1,1+pi/2] because we don't want to test the Bessel
1624	;; series and we don't want to test near pi because sin(pi)
1625	;; = 0, where we will lose accuracy.
1626	for x = (+ 1 (random (/ pi 2)))
1627	for b = (bessel-j 0.5d0 x)
1628	for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi)))
1629	for result = (check-accuracy 48.42 b true)
1630	when result
1631	append (list (list (list k x) result)))
1632	nil)
1633
1634	(rt:deftest bessel-j-1/2.d.1.a
1635	(let* ((x 2.3831631289164497d0)
1636	(b (bessel-j 0.5d0 x))
1637	(true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi)))))
1638	(check-accuracy 48.42 b true))
1639	nil)
1640
1641	(rt:deftest bessel-j-1/2.q.1
1642	(loop for k from 0 below 10
1643	;; x in [1,1+pi/2] because we don't want to test the Bessel
1644	;; series and we don't want to test near pi because sin(pi)
1645	;; = 0, where we will lose accuracy.
1646	for x = (+ 1 (random (/ (float-pi #q1) 2)))
1647	for b = (bessel-j #q0.5 x)
1648	for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1))))
1649	for result = (check-accuracy 169.45 b true)
1650	when result
1651	append (list (list (list k x) result)))
1652	nil)
1653
1654	(rt:deftest bessel-j-1/2.q.1.a
1655	(let* ((x #q1.1288834862545916200627583005758663687705443417892789067029865493882q0)
1656	(b (bessel-j #q0.5 x))
1657	(true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1))))))
1658	(check-accuracy 182.92 b true))
1659	nil)
1660
1661	(rt:deftest bessel-j-1/2.q.1.b
1662	(let* ((x #q1.1288834862545916200627583005758663687705443417892789067029865493882q0)
1663	(b (bessel-j #q0.5 x))
1664	(true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1))))))
1665	(check-accuracy 173.28 b true))
1666	nil)
1667
1668	(rt:deftest bessel-j-1/2.q.1.c
1669	(let* ((x #q1.0360263937639582798798376485114581552570020473846457752365459851056q0)
1670	(b (bessel-j #q0.5 x))
1671	(true (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1))))))
1672	(check-accuracy 169.45 b true))
1673	nil)
1674
1675	;; Bessel J for complex args
1676	(rt:deftest bessel-j-complex-arg.d.1
1677	(let ((b (bessel-j 0d0 #c(1d0 1)))
1678	(true #c(0.9376084768060293d0 -0.4965299476091221d0)))
1679	(check-accuracy 50.73 b true))
1680	nil)
1681
1682	(rt:deftest bessel-j-complex-arg.d.2
1683	(let ((b (bessel-j 1d0 #c(1d0 1)))
1684	(true #c(0.6141603349229036d0 0.3650280288270878d0)))
1685	(check-accuracy 52.51 b true))
1686	nil)
1687
1688	(rt:deftest bessel-j-complex-arg.d.3
1689	(let ((b (bessel-j 2d0 #c(1d0 1)))
1690	(true #c(0.0415798869439621d0 0.2473976415133063d0)))
1691	(check-accuracy 50.41 b true))
1692	nil)
1693
1694	(rt:deftest bessel-j-complex-arg.d.4
1695	(let ((b (bessel-j 2.3d0 #c(1d0 1)))
1696	(true #c(-0.0141615213034667d0 0.1677798241687935d0)))
1697	(check-accuracy 48.56 b true))
1698	nil)
1699
1700	(rt:deftest bessel-j-complex-arg.d.5
1701	(let ((b (bessel-j -2.3d0 #c(1d0 1)))
1702	(true #c(0.1920598664138632d0 -0.5158676904105332d0)))
1703	(check-accuracy 50.97 b true))
1704	nil)
1705
1706	(rt:deftest bessel-j-1/2-complex.d.1
1707	(loop for k from 0 below 10
1708	for x = (complex (random (/ pi 2))
1709	(random (/ pi 2)))
1710	for b = (bessel-j 0.5d0 x)
1711	for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 pi)))
1712	for result = (check-accuracy 49.8 b true)
1713	when result
1714	append (list (list (list k x) result)))
1715	nil)
1716
1717	(rt:deftest bessel-j-1/2-complex.q.1
1718	(loop for k from 0 below 10
1719	for x = (complex (random (/ (float-pi #q1) 2))
1720	(random (/ (float-pi #q1) 2)))
1721	for b = (bessel-j #q0.5 x)
1722	for true = (* (/ (sin x) (sqrt x)) (sqrt (/ 2 (float-pi #q1))))
1723	for result = (check-accuracy 212 b true)
1724	when result
1725	append (list (list (list k x) result)))
1726	nil)

Note: See TracBrowser for help on using the browser.

Download in other formats:

Original Format