root/rt-tests.lisp @ 8ade177a0ce9bbb89b963ff29e46f38f377e9530

Revision 8ade177a0ce9bbb89b963ff29e46f38f377e9530, 31.9 KB (checked in by Raymond Toy <toy.raymond@…>, 3 years ago)

Add Elliptic theta functions and tests.

oct.asd:
o Add qd-theta.

qd-theta.lisp:
o New file with Elliptic theta functions and elliptic nome function.

rt-tests.lisp:
o Tests for theta functions.
o Relax accuracy requirements for some of the tests os that they can

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