| blob | 8d4b3e48674ce3f5e27a8060232ec6156a145683 |
1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10-*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
4 ;;; ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;; ;;;;;
8 ;;; Copyright (c) 2001 by Raymond Toy. Replaced everything and added ;;;;;
9 ;;; support for symbolic manipulation of all 12 Jacobian elliptic ;;;;;
10 ;;; functions and the complete and incomplete elliptic integral ;;;;;
11 ;;; of the first, second and third kinds. ;;;;;
12 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
15 ;;(macsyma-module ellipt)
17 ;;;
18 ;;; Jacobian elliptic functions and elliptic integrals.
19 ;;;
20 ;;; References:
21 ;;;
22 ;;; [1] Abramowitz and Stegun
23 ;;; [2] Lawden, Elliptic Functions and Applications, Springer-Verlag, 1989
24 ;;; [3] Whittaker & Watson, A Course of Modern Analysis
25 ;;;
26 ;;; We use the definitions from Abramowitz and Stegun where our
27 ;;; sn(u,m) is sn(u|m). That is, the second arg is the parameter,
28 ;;; instead of the modulus k or modular angle alpha.
29 ;;;
30 ;;; Note that m = k^2 and k = sin(alpha).
31 ;;;
33 ;;
34 ;; Routines for computing the basic elliptic functions sn, cn, and dn.
35 ;;
36 ;;
37 ;; A&S gives several methods for computing elliptic functions
38 ;; including the AGM method (16.4) and ascending and descending Landen
39 ;; transformations (16.12 and 16.14). The latter are actually quite
40 ;; fast, only requiring simple arithmetic and square roots for the
41 ;; transformation until the last step. The AGM requires evaluation of
42 ;; several trigonometric functions at each stage.
43 ;;
44 ;; However, the Landen transformations appear to have some round-off
45 ;; issues. For example, using the ascending transform to compute cn,
46 ;; cn(100,.7) > 1e10. This is clearly not right since |cn| <= 1.
47 ;;
54 ;; A&S 16.14.1
55 ;;
56 ;; Take care in computing this transform. For the case where
57 ;; m is complex, we should compute sqrt(mu1) first as
58 ;; (1-sqrt(m))/(1+sqrt(m)), and then square this to get mu1.
59 ;; If not, we may choose the wrong branch when computing
60 ;; sqrt(mu1).
69 ;; Note: Don't calculate mu first, as given in 16.12.1. We
70 ;; should calculate sqrt(mu) = (1-sqrt(m1)/(1+sqrt(m1)), and
71 ;; then compute mu = sqrt(mu)^2. If we calculate mu first,
72 ;; sqrt(mu) loses information when m or m1 is complex.
79 ;;
80 ;; This appears to work quite well for both real and complex values
81 ;; of u.
84 ;; A&S 16.6.1
87 ;; A&S 16.6.1
98 ;; jacobi_sn(u,0) = sin(u). Should we use A&S 16.13.1 if m
99 ;; is small enough?
100 ;;
101 ;; sn(u,m) = sin(u) - 1/4*m(u-sin(u)*cos(u))*cos(u)
104 ;; jacobi_sn(u,1) = tanh(u). Should we use A&S 16.15.1 if m
105 ;; is close enough to 1?
106 ;;
107 ;; sn(u,m) = tanh(u) + 1/4*(1-m)*(sinh(u)*cosh(u)-u)*sech(u)^2
110 ;; Use the ascending Landen transformation to compute sn.
114 s)))))
118 ;; jacobi_dn(u,0) = 1. Should we use A&S 16.13.3 for small m?
119 ;;
120 ;; dn(u,m) = 1 - 1/2*m*sin(u)^2
123 ;; jacobi_dn(u,1) = sech(u). Should we use A&S 16.15.3 if m
124 ;; is close enough to 1?
125 ;;
126 ;; dn(u,m) = sech(u) + 1/4*(1-m)*(sinh(u)*cosh(u)+u)*tanh(u)*sech(u)
129 ;; Use the Gauss transformation from
130 ;; http://functions.wolfram.com/09.29.16.0013.01:
131 ;;
132 ;;
133 ;; dn((1+sqrt(m))*z, 4*sqrt(m)/(1+sqrt(m))^2)
134 ;; = (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2)
135 ;;
136 ;; So
137 ;;
138 ;; dn(y, mu) = (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2)
139 ;;
140 ;; where z = y/(1+sqrt(m)) and mu=4*sqrt(m)/(1+sqrt(m))^2.
141 ;;
142 ;; Solve for m, and we get
143 ;;
144 ;; sqrt(m) = -(mu+2*sqrt(1-mu)-2)/mu or (-mu+2*sqrt(1-mu)+2)/mu.
145 ;;
146 ;; I don't think it matters which sqrt we use, so I (rtoy)
147 ;; arbitrarily choose the first one above.
148 ;;
149 ;; Note that (1-sqrt(1-mu))/(1+sqrt(1-mu)) is the same as
150 ;; -(mu+2*sqrt(1-mu)-2)/mu. Also, the former is more
151 ;; accurate for small mu.
163 ;; jacobi_cn(u,0) = cos(u). Should we use A&S 16.13.2 for
164 ;; small m?
165 ;;
166 ;; cn(u,m) = cos(u) + 1/4*m*(u-sin(u)*cos(u))*sin(u)
169 ;; jacobi_cn(u,1) = sech(u). Should we use A&S 16.15.3 if m
170 ;; is close enough to 1?
171 ;;
172 ;; cn(u,m) = sech(u) - 1/4*(1-m)*(sinh(u)*cosh(u)-u)*tanh(u)*sech(u)
175 ;; Use the ascending Landen transformation, A&S 16.14.3.
181 d)))))))
185 ;; Tell maxima what the derivatives are.
186 ;;
187 ;; Lawden says the derivative wrt to k but that's not what we want.
188 ;;
189 ;; Here's the derivation we used, based on how Lawden gets his results.
190 ;;
191 ;; Let
192 ;;
193 ;; diff(sn(u,m),m) = s
194 ;; diff(cn(u,m),m) = p
195 ;; diff(dn(u,m),m) = q
196 ;;
197 ;; From the derivatives of sn, cn, dn wrt to u, we have
198 ;;
199 ;; diff(sn(u,m),u) = cn(u)*dn(u)
200 ;; diff(cn(u,m),u) = -cn(u)*dn(u)
201 ;; diff(dn(u,m),u) = -m*sn(u)*cn(u)
202 ;;
204 ;; Differentiate these wrt to m:
205 ;;
206 ;; diff(s,u) = p*dn + cn*q
207 ;; diff(p,u) = -p*dn - q*dn
208 ;; diff(q,u) = -sn*cn - m*s*cn - m*sn*q
209 ;;
210 ;; Also recall that
211 ;;
212 ;; sn(u)^2 + cn(u)^2 = 1
213 ;; dn(u)^2 + m*sn(u)^2 = 1
214 ;;
215 ;; Differentiate these wrt to m:
216 ;;
217 ;; sn*s + cn*p = 0
218 ;; 2*dn*q + sn^2 + 2*m*sn*s = 0
219 ;;
220 ;; Thus,
221 ;;
222 ;; p = -s*sn/cn
223 ;; q = -m*s*sn/dn - sn^2/dn/2
224 ;;
225 ;; So
226 ;; diff(s,u) = -s*sn*dn/cn - m*s*sn*cn/dn - sn^2*cn/dn/2
227 ;;
228 ;; or
229 ;;
230 ;; diff(s,u) + s*(sn*dn/cn + m*sn*cn/dn) = -1/2*sn^2*cn/dn
231 ;;
232 ;; diff(s,u) + s*sn/cn/dn*(dn^2 + m*cn^2) = -1/2*sn^2*cn/dn
233 ;;
234 ;; Multiply through by the integrating factor 1/cn/dn:
235 ;;
236 ;; diff(s/cn/dn, u) = -1/2*sn^2/dn^2 = -1/2*sd^2.
237 ;;
238 ;; Integrate this to get
239 ;;
240 ;; s/cn/dn = C + -1/2*int sd^2
241 ;;
242 ;; It can be shown that C is zero.
243 ;;
244 ;; We know that (by differentiating this expression)
245 ;;
246 ;; int dn^2 = (1-m)*u+m*sn*cd + m*(1-m)*int sd^2
247 ;;
248 ;; or
249 ;;
250 ;; int sd^2 = 1/m/(1-m)*int dn^2 - u/m - sn*cd/(1-m)
251 ;;
252 ;; Thus, we get
253 ;;
254 ;; s/cn/dn = u/(2*m) + sn*cd/(2*(1-m)) - 1/2/m/(1-m)*int dn^2
255 ;;
256 ;; or
257 ;;
258 ;; s = 1/(2*m)*u*cn*dn + 1/(2*(1-m))*sn*cn^2 - 1/2/(m*(1-m))*cn*dn*E(u)
259 ;;
260 ;; where E(u) = int dn^2 = elliptic_e(am(u)) = elliptic_e(asin(sn(u)))
261 ;;
262 ;; This is our desired result:
263 ;;
264 ;; s = 1/(2*m)*cn*dn*[u - elliptic_e(asin(sn(u)),m)/(1-m)] + sn*cn^2/2/(1-m)
265 ;;
266 ;;
267 ;; Since diff(cn(u,m),m) = p = -s*sn/cn, we have
268 ;;
269 ;; p = -1/(2*m)*sn*dn[u - elliptic_e(asin(sn(u)),m)/(1-m)] - sn^2*cn/2/(1-m)
270 ;;
271 ;; diff(dn(u,m),m) = q = -m*s*sn/dn - sn^2/dn/2
272 ;;
273 ;; q = -1/2*sn*cn*[u-elliptic_e(asin(sn),m)/(1-m)] - m*sn^2*cn^2/dn/2/(1-m)
274 ;;
275 ;; - sn^2/dn/2
276 ;;
277 ;; = -1/2*sn*cn*[u-elliptic_e(asin(sn),m)/(1-m)] + dn*sn^2/2/(m-1)
278 ;;
291 grad)
305 grad)
320 grad)
322 ;; The inverse elliptic functions.
323 ;;
324 ;; F(phi|m) = asn(sin(phi),m)
325 ;;
326 ;; so asn(u,m) = F(asin(u)|m)
329 ;; Lawden 3.1.2:
330 ;; inverse_jacobi_sn(x) = integrate(1/sqrt(1-t^2)/sqrt(1-m*t^2),t,0,x)
331 ;; -> 1/sqrt(1-x^2)/sqrt(1-m*x^2)
337 ;; diff(F(asin(u)|m),m)
349 grad)
351 ;; Let u = inverse_jacobi_cn(x). Then jacobi_cn(u) = x or
352 ;; sqrt(1-jacobi_sn(u)^2) = x. Or
353 ;;
354 ;; jacobi_sn(u) = sqrt(1-x^2)
355 ;;
356 ;; So u = inverse_jacobi_sn(sqrt(1-x^2),m) = inverse_jacob_cn(x,m)
357 ;;
360 ;; Whittaker and Watson, 22.121
361 ;; inverse_jacobi_cn(u,m) = integrate(1/sqrt(1-t^2)/sqrt(1-m+m*t^2), t, u, 1)
362 ;; -> -1/sqrt(1-x^2)/sqrt(1-m+m*x^2)
384 m)
389 m)))))))
390 grad)
392 ;; Let u = inverse_jacobi_dn(x). Then
393 ;;
394 ;; jacobi_dn(u) = x or
395 ;;
396 ;; x^2 = jacobi_dn(u)^2 = 1 - m*jacobi_sn(u)^2
397 ;;
398 ;; so jacobi_sn(u) = sqrt(1-x^2)/sqrt(m)
399 ;;
400 ;; or u = inverse_jacobi_sn(sqrt(1-x^2)/sqrt(m))
403 ;; Whittaker and Watson, 22.121
404 ;; inverse_jacobi_dn(u,m) = integrate(1/sqrt(1-t^2)/sqrt(t^2-(1-m)), t, u, 1)
405 ;; -> -1/sqrt(1-x^2)/sqrt(x^2+m-1)
439 m)
447 m))))))))
448 grad)
451 ;; Possible forms of a complex number:
452 ;;
453 ;; 2.3
454 ;; $%i
455 ;; ((mplus simp) 2.3 ((mtimes simp) 2.3 $%i))
456 ;; ((mplus simp) 2.3 $%i))
457 ;; ((mtimes simp) 2.3 $%i)
458 ;;
461 ;; Is argument u a complex number with real and imagpart satisfying predicate ntypep?
484 ;; Convert a Lisp number to a maxima number
487 (t (merror (intl:gettext "COMPLEXIFY: argument must be a Lisp real or complex number.~%COMPLEXIFY: found: ~:M") x))))
490 ;; Replace elliptic_kc(m) in the expression with sym. Check to see
491 ;; if the resulting expression is linear in sym and the constant
492 ;; term is zero. If so, return the coefficient of sym, i.e, the
493 ;; coefficient of elliptic_kc(m).
500 nil)))
503 ;; Replace elliptic_kc(m) in the expression with sym. Check to see
504 ;; if the resulting expression is linear in sym and the constant
505 ;; term is zero. If so, return the coefficient of sym, i.e, the
506 ;; coefficient of elliptic_kc(m), and the constant term. Otherwise,
507 ;; return NIL.
514 nil)))
516 ;; Tell maxima how to simplify the functions
525 args
531 args
534 ;; A&S 16.5.1
537 ;; A&S 16.6.1
540 ;; A&S 16.6.1
547 %inverse_jacobi_ns
548 %inverse_jacobi_cn
549 %inverse_jacobi_nc
550 %inverse_jacobi_dn
551 %inverse_jacobi_nd
552 %inverse_jacobi_sc
553 %inverse_jacobi_cs
554 %inverse_jacobi_sd
555 %inverse_jacobi_ds
556 %inverse_jacobi_cd
557 %inverse_jacobi_dc))
561 (%inverse_jacobi_sn
562 ;; jacobi_sn(inverse_jacobi_sn(u,m), m) = u
563 inv-arg)
564 (%inverse_jacobi_ns
565 ;; inverse_jacobi_ns(u,m) = inverse_jacobi_sn(1/u,m)
567 (%inverse_jacobi_cn
568 ;; sn(x)^2 + cn(x)^2 = 1 so sn(x) = sqrt(1-cn(x)^2)
570 (%inverse_jacobi_nc
571 ;; inverse_jacobi_nc(u) = inverse_jacobi_cn(1/u)
573 m))
574 (%inverse_jacobi_dn
575 ;; dn(x)^2 + m*sn(x)^2 = 1 so
576 ;; sn(x) = 1/sqrt(m)*sqrt(1-dn(x)^2)
579 (%inverse_jacobi_nd
580 ;; inverse_jacobi_nd(u) = inverse_jacobi_dn(1/u)
582 m))
583 (%inverse_jacobi_sc
584 ;; See below for inverse_jacobi_sc.
586 (%inverse_jacobi_cs
587 ;; inverse_jacobi_cs(u) = inverse_jacobi_sc(1/u)
589 m))
590 (%inverse_jacobi_sd
591 ;; See below for inverse_jacobi_sd
593 (%inverse_jacobi_ds
594 ;; inverse_jacobi_ds(u) = inverse_jacobi_sd(1/u)
596 m))
597 (%inverse_jacobi_cd
598 ;; See below
601 (%inverse_jacobi_dc
603 ;; A&S 16.20.1 (Jacobi's Imaginary transformation)
608 ;; sn(m*K+u)
609 ;;
610 ;; A&S 16.8.1
612 coef
616 ;; sn(4*m*K + u) = sn(u), sn(0) = 0
618 0
621 ;; sn(4*m*K + K + u) = sn(K+u) = cd(u)
622 ;; sn(K) = 1
624 1
627 ;; sn(4*m*K+2*K + u) = sn(2*K+u) = -sn(u)
628 ;; sn(2*K) = 0
630 0
633 ;; sn(4*m*K+3*K+u) = sn(2*K + K + u) = -sn(K+u) = -cd(u)
634 ;; sn(3*K) = -1
636 -1
640 ;; A&S 16.5.2
641 ;;
642 ;; sn(1/2*K) = 1/sqrt(1+sqrt(1-m))
648 ;; A&S 16.5.2
649 ;;
650 ;; sn(1/2*K + K) = cd(1/2*K,m)
653 m))
657 ;; Nothing to do
667 args
673 args
676 ;; A&S 16.5.1
679 ;; A&S 16.6.2
682 ;; A&S 16.6.2
689 %inverse_jacobi_ns
690 %inverse_jacobi_cn
691 %inverse_jacobi_nc
692 %inverse_jacobi_dn
693 %inverse_jacobi_nd
694 %inverse_jacobi_sc
695 %inverse_jacobi_cs
696 %inverse_jacobi_sd
697 %inverse_jacobi_ds
698 %inverse_jacobi_cd
699 %inverse_jacobi_dc))
704 ;; I'm lazy. Use cn(x) = sqrt(1-sn(x)^2). Hope
705 ;; this is right.
708 ;; A&S 16.20.2 (Jacobi's Imaginary transformation)
712 ;; cn(m*K+u)
713 ;;
714 ;; A&S 16.8.2
716 coef
720 ;; cn(4*m*K + u) = cn(u),
721 ;; cn(0) = 1
723 1
726 ;; cn(4*m*K + K + u) = cn(K+u) = -sqrt(m1)*sd(u)
727 ;; cn(K) = 0
729 0
733 ;; cn(4*m*K + 2*K + u) = cn(2*K+u) = -cn(u)
734 ;; cn(2*K) = -1
736 -1
739 ;; cn(4*m*K + 3*K + u) = cn(2*K + K + u) =
740 ;; -cn(K+u) = sqrt(m1)*sd(u)
741 ;;
742 ;; cn(3*K) = 0
744 0
749 ;; A&S 16.5.2
750 ;; cn(1/2*K) = (1-m)^(1/4)/sqrt(1+sqrt(1-m))
768 args
774 args
777 ;; A&S 16.5.1
780 ;; A&S 16.6.3
783 ;; A&S 16.6.3
790 %inverse_jacobi_ns
791 %inverse_jacobi_cn
792 %inverse_jacobi_nc
793 %inverse_jacobi_dn
794 %inverse_jacobi_nd
795 %inverse_jacobi_sc
796 %inverse_jacobi_cs
797 %inverse_jacobi_sd
798 %inverse_jacobi_ds
799 %inverse_jacobi_cd
800 %inverse_jacobi_dc))
803 ;; jacobi_dn(inverse_jacobi_dn(u,m), m) = u
806 ;; Express in terms of sn:
807 ;; dn(x) = sqrt(1-m*sn(x)^2)
812 ;; A&S 16.5.3
813 ;; dn(sqrt(1-m),m) = K(m)
815 ;; A&S 16.20.2 (Jacobi's Imaginary transformation)
820 ;; A&S 16.8.3
821 ;;
822 ;; dn(m*K+u) has period 2K
823 ;;
825 coef
829 ;; dn(2*m*K + u) = dn(u)
830 ;; dn(0) = 1
832 1
833 ;; dn(4*m*K+2*K + u) = dn(2*K+u) = dn(u)
836 ;; dn(2*m*K + K + u) = dn(K + u) = sqrt(1-m)*nd(u)
837 ;; dn(K) = sqrt(1-m)
844 ;; A&S 16.5.2
845 ;; dn(1/2*K) = (1-m)^(1/4)
852 ;; Should we simplify the inverse elliptic functions into the
853 ;; appropriate incomplete elliptic integral? I think we should leave
854 ;; it, but perhaps allow some way to do that transformation if
855 ;; desired.
859 ;; To numerically evaluate inverse_jacobi_sn (asn), use
860 ;;
861 ;; asn(x,m) = F(asin(x),m)
862 ;;
863 ;; But F(phi,m) = sin(phi)*rf(cos(phi)^2, 1-m*sin(phi)^2,1). Thus
864 ;;
865 ;; asn(x,m) = F(asin(x),m)
866 ;; = x*rf(1-x^2,1-m*x^2,1)
867 ;;
868 ;; I (rtoy) am not 100% about the first identity above for all
869 ;; complex values of x and m, but tests seem to indicate that it
870 ;; produces the correct value as verified by verifying
871 ;; jacobi_sn(inverse_jacobi_sn(x,m),m) = x.
876 (* uu
882 args
897 args
905 ;; asn(0,m) = 0
908 ;; asn(1,m) = elliptic_kc(m)
911 ;; asn(-1,m) = -elliptic_kc(m)
914 ;; asn(x,0) = F(asin(x),0) = asin(x)
917 ;; asn(x,1) = F(asin(x),1) = log(tan(pi/4+asin(x)/2))
923 ;; inverse_jacobi_sn(sn(u)) = u
926 ;; Nothing to do
932 ;; Numerically evaluate acn
933 ;;
934 ;; acn(x,m) = F(acos(x),m)
938 args
946 args
950 ;; asn(x,0) = F(acos(x),0) = acos(x)
953 ;; asn(x,1) = F(asin(x),1) = log(tan(pi/2+asin(x)/2))
963 ;; inverse_jacobi_cn(cn(u)) = u
966 ;; Nothing to do
976 args
986 args
989 ;; x = dn(u,1) = sech(u). so u = asech(x)
992 ;; jacobi_dn(0,m) = 1
995 ;; jacobi_dn(K(m),m) = sqrt(1-m) so
996 ;; inverse_jacobi_dn(sqrt(1-m),m) = K(m)
1002 ;; inverse_jacobi_dn(dn(u)) = u
1005 ;; Nothing to do
1008 ;;;; Elliptic integrals
1016 "Compute Carlson's incomplete or complete elliptic integral of the
1017 first kind:
1019 INF
1020 /
1021 [ 1
1022 RF(x, y, z) = I ----------------------------------- dt
1023 ] SQRT(x + t) SQRT(y + t) SQRT(z + t)
1024 /
1025 0
1027 x, y, and z may be complex.
1028 "
1056 ;; Elliptic integral of the first kind (Legendre's form):
1057 ;;
1058 ;;
1059 ;; phi
1060 ;; /
1061 ;; [ 1
1062 ;; I ------------------- ds
1063 ;; ] 2
1064 ;; / SQRT(1 - m SIN (s))
1065 ;; 0
1073 ;; A&S 17.4.15
1074 ;;
1075 ;; F(phi|m) = 1/sqrt(m)*F(theta|1/m)
1076 ;;
1077 ;; with sin(theta) = sqrt(m)*sin(phi)
1081 ;; A&S 17.4.17
1087 root)
1089 root))))
1091 ;; A&S 17.4.19
1092 phi)
1094 ;; A&S 17.4.21
1095 ;;
1097 ;; = log(tan(pi/4+pi/2))
1102 ;; A&S 17.4.3
1111 (* sin-phi
1121 phi-arg m-arg)))))
1128 (* sin-phi
1133 ;; Elliptic F is quasi-periodic wrt to z:
1134 ;;
1135 ;; F(z|m) = F(z - pi*round(Re(z)/pi)|m) + 2*round(Re(z)/pi)*K(m)
1139 0
1143 ;; Complete elliptic integral of the first kind
1147 ;; A&S 17.4.17
1153 root)
1155 root))))
1157 ;; A&S 17.4.19
1169 ;; Elliptic integral of the second kind (Legendre's form):
1170 ;;
1171 ;;
1172 ;; phi
1173 ;; /
1174 ;; [ 2
1175 ;; I SQRT(1 - m SIN (s)) ds
1176 ;; ]
1177 ;; /
1178 ;; 0
1184 ;; A&S 17.4.23
1185 phi)
1187 ;; A&S 17.4.25
1200 ;; Elliptic E is quasi-periodic wrt to phi:
1201 ;;
1202 ;; E(z|m) = E(z - %pi*round(Re(z)/%pi)|m) + 2*round(Re(z)/%pi)*E(m)
1207 ;; Complete version
1211 ;; A&S 17.4.23
1214 ;; A&S 17.4.25
1223 ;; Define the elliptic integrals for maxima
1224 ;;
1225 ;; We use the definitions given in A&S 17.2.6 and 17.2.8. In particular:
1226 ;;
1227 ;; phi
1228 ;; /
1229 ;; [ 1
1230 ;; F(phi|m) = I ------------------- ds
1231 ;; ] 2
1232 ;; / SQRT(1 - m SIN (s))
1233 ;; 0
1234 ;;
1235 ;; and
1236 ;;
1237 ;; phi
1238 ;; /
1239 ;; [ 2
1240 ;; E(phi|m) = I SQRT(1 - m SIN (s)) ds
1241 ;; ]
1242 ;; /
1243 ;; 0
1244 ;;
1245 ;; That is, we do not use the modular angle, alpha, as the second arg;
1246 ;; the parameter m = sin(alpha)^2 is used.
1247 ;;
1250 ;; The derivative of F(phi|m) wrt to phi is easy. The derivative wrt
1251 ;; to m is harder. Here is a derivation. Hope I got it right.
1252 ;;
1253 ;; diff(integrate(1/sqrt(1-m*sin(x)^2),x,0,phi), m);
1254 ;;
1255 ;; PHI
1256 ;; / 2
1257 ;; [ SIN (x)
1258 ;; I ------------------ dx
1259 ;; ] 2 3/2
1260 ;; / (1 - m SIN (x))
1261 ;; 0
1262 ;; --------------------------
1263 ;; 2
1264 ;;
1265 ;;
1266 ;; Now use the following relationship that is easily verified:
1267 ;;
1268 ;; 2 2
1269 ;; (1 - m) SIN (x) COS (x) COS(x) SIN(x)
1270 ;; ------------------- = ------------------- - DIFF(-------------------, x)
1271 ;; 2 2 2
1272 ;; SQRT(1 - m SIN (x)) SQRT(1 - m SIN (x)) SQRT(1 - m SIN (x))
1273 ;;
1274 ;;
1275 ;; Now integrate this to get:
1276 ;;
1277 ;;
1278 ;; PHI
1279 ;; / 2
1280 ;; [ SIN (x)
1281 ;; (1 - m) I ------------------- dx =
1282 ;; ] 2
1283 ;; / SQRT(1 - m SIN (x))
1284 ;; 0
1286 ;;
1287 ;; PHI
1288 ;; / 2
1289 ;; [ COS (x)
1290 ;; + I ------------------- dx
1291 ;; ] 2
1292 ;; / SQRT(1 - m SIN (x))
1293 ;; 0
1294 ;; COS(PHI) SIN(PHI)
1295 ;; - ---------------------
1296 ;; 2
1297 ;; SQRT(1 - m SIN (PHI))
1298 ;;
1299 ;; Use the fact that cos(x)^2 = 1 - sin(x)^2 to show that this
1300 ;; integral on the RHS is:
1301 ;;
1302 ;;
1303 ;; (1 - m) elliptic_F(PHI, m) + elliptic_E(PHI, m)
1304 ;; -------------------------------------------
1305 ;; m
1306 ;; So, finally, we have
1307 ;;
1308 ;;
1309 ;;
1310 ;; d
1311 ;; 2 -- (elliptic_F(PHI, m)) =
1312 ;; dm
1313 ;;
1314 ;; elliptic_E(PHI, m) - (1 - m) elliptic_F(PHI, m) COS(PHI) SIN(PHI)
1315 ;; ---------------------------------------------- - ---------------------
1316 ;; m 2
1317 ;; SQRT(1 - m SIN (PHI))
1318 ;; ----------------------------------------------------------------------
1319 ;; 1 - m
1323 ;; diff wrt phi
1324 ;; 1/sqrt(1-m*sin(phi)^2)
1328 ;; diff wrt m
1341 grad)
1343 ;;
1344 ;; The derivative of E(phi|m) wrt to m is much simpler to derive than F(phi|m).
1345 ;;
1346 ;; Take the derivative of the definition to get
1347 ;;
1348 ;; PHI
1349 ;; / 2
1350 ;; [ SIN (x)
1351 ;; I ------------------- dx
1352 ;; ] 2
1353 ;; / SQRT(1 - m SIN (x))
1354 ;; 0
1355 ;; - ---------------------------
1356 ;; 2
1357 ;;
1358 ;; It is easy to see that
1359 ;;
1360 ;; PHI
1361 ;; / 2
1362 ;; [ SIN (x)
1363 ;; elliptic_F(PHI, m) - m I ------------------- dx = elliptic_E(PHI, m)
1364 ;; ] 2
1365 ;; / SQRT(1 - m SIN (x))
1366 ;; 0
1367 ;;
1368 ;; So we finally have
1369 ;;
1370 ;; d elliptic_E(PHI, m) - elliptic_F(PHI, m)
1371 ;; -- (elliptic_E(PHI, m)) = ---------------------------------------
1372 ;; dm 2 m
1376 ;; sqrt(1-m*sin(phi)^2)
1380 ;; diff wrt m
1384 grad)
1389 ;; Numerically evaluate it
1393 args
1401 args
1407 ;; A&S 17.4.19
1408 phi)
1410 ;; A&S 17.4.21. Let's pick the log tan form. But this
1411 ;; isn't right if we know that abs(phi) > %pi/2, where
1412 ;; elliptic_f is undefined (or infinity).
1420 phi m))))
1422 ;; Complete elliptic integral
1425 ;; Nothing to do
1431 ;; Numerically evaluate it
1434 (complexify (bigfloat::bf-elliptic-e (complex ($float ($realpart phi)) ($float ($imagpart phi)))
1441 args
1447 ;; A&S 17.4.23
1448 phi)
1450 ;; A&S 17.4.25, but handle periodicity:
1451 ;; elliptic_e(x,m) = elliptic_e(x-%pi*round(x/%pi), m)
1452 ;; + 2*round(x/%pi)*elliptic_ec(m)
1453 ;;
1454 ;; Or
1455 ;;
1456 ;; elliptic_e(x,1) = sin(x-%pi*round(x/%pi)) + 2*round(x/%pi)*elliptic_ec(m)
1457 ;;
1461 mult-pi)))
1466 ;; Complete elliptic integral
1472 ;; Handle the case where phi is a number where we can apply
1473 ;; the periodicity property without blowing up the
1474 ;; expression.
1479 m)
1484 ;; Nothing to do
1487 ;; Complete elliptic integrals
1488 ;;
1489 ;; elliptic_kc(m) = elliptic_f(%pi/2, m)
1490 ;;
1491 ;; elliptic_ec(m) = elliptic_e(%pi/2, m)
1492 ;;
1494 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1496 ;;; We support a simplim%function. The function is looked up in simplimit and
1497 ;;; handles specific values of the function.
1502 ;; Look for the limit of the argument
1505 ;; For an argument 1 return $infinity.
1508 ;; All other cases are handled by the simplifier of the function.
1511 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1516 ;; elliptic_kc(1) is complex infinity. Maxima can not handle
1517 ;; infinities correctly, throw a Maxima error.
1520 m))
1522 ;; Numerically evaluate it
1525 (complexify (bigfloat::bf-elliptic-k (complex ($float ($realpart m)) ($float ($imagpart m))))))
1528 args
1533 ;; http://functions.wolfram.com/EllipticIntegrals/EllipticK/03/01/
1534 ;;
1535 ;; elliptic_kc(1/2) = 8*%pi^(3/2)/gamma(-1/4)^2
1539 ;; elliptic_kc(-1) = gamma(1/4)^2/(4*sqrt(2*%pi))
1543 ;; elliptic_kc(17-12*sqrt(2)) = 2*(2+sqrt(2))*%pi^(3/2)/gamma(-1/4)^2
1552 ;; elliptic_kc((sqrt(3)+2)/4) = sqrt(%pi)*gamma(1/3)/gamma(5/6).
1553 ;;
1554 ;; First evaluate this integral, where y = sqrt(1+t^3).
1555 ;;
1556 ;; integrate(1/y,t,-1,inf) = integrate(1/y,t,-1,0) + integrate(1/y,t,0,inf).
1557 ;;
1558 ;; The second integral, maxima gives beta(1/6,1/3)/3.
1559 ;;
1560 ;; For the first, we can use the change of variable x=-u^(1/3) to get
1561 ;;
1562 ;; integrate(1/sqrt(1-u)/u^(2/3),u,0,1)
1563 ;;
1564 ;; which is a beta integral that maxima can evaluate to
1565 ;; beta(1/3,1/2)/3. Then we see the value of the initial
1566 ;; integral is
1567 ;;
1568 ;; beta(1/6,1/3)/3 + beta(1/3,1/2)/3
1569 ;;
1570 ;; (Thanks to Guilherme Namen for this derivation on the mailing list, 2023-03-09.)
1571 ;;
1572 ;; We can simplify this expression by converting to gamma functions:
1573 ;;
1574 ;; beta(1/6,1/3)/3 + beta(1/3,1/2)/3 =
1575 ;; (gamma(1/3)*(gamma(1/6)*gamma(5/6)+%pi))/(3*sqrt(%pi)*gamma(5/6));
1576 ;;
1577 ;; Using the reflection formula gamma(1-z)*gamma(z) =
1578 ;; %pi/sin(%pi*z), we can write gamma(1/6)*gamma(5/6) =
1579 ;; %pi/sin(%pi*1/6) = 2*%pi. Finally, we have
1580 ;;
1581 ;; sqrt(%pi)*gamma(1/3)/gamma(5/6);
1582 ;;
1583 ;; All that remains is to show that integrate(1/y,t) can be
1584 ;; written as an inverse_jacobi_cn function with modulus
1585 ;; (sqrt(3)+2)/4.
1586 ;;
1587 ;; First apply the substitution
1588 ;;
1589 ;; s = (t+sqrt(3)+1)/(t-sqrt(3)+1). We then have the integral
1590 ;;
1591 ;; C*integrate(1/sqrt(s^2-1)/sqrt(s^2+4*sqrt(3)+7),s)
1592 ;;
1593 ;; where C is some constant. From A&S 14.4.49, we can see
1594 ;; this integral is the inverse_jacobi_nc function with
1595 ;; modulus of (4*sqrt(3)+7)/(4*sqrt(3)+7+1) =
1596 ;; (sqrt(3)+2)/4.
1600 ($hypergeometric_representation
1601 ;; See http://functions.wolfram.com/08.02.26.0001.01
1602 ;;
1603 ;; elliptic_kc(z) = %pi/2*%f[2,1]([1/2,1/2],[1], z)
1604 ;;
1609 m)))
1611 ;; Nothing to do
1616 ;; diff wrt m
1625 grad)
1630 ;; Numerically evaluate it
1634 args
1638 args
1640 ;; Some special cases we know about.
1646 ;; elliptic_ec(1/2). Use the identity
1647 ;;
1648 ;; elliptic_ec(z)*elliptic_kc(1-z) - elliptic_kc(z)*elliptic_kc(1-z)
1649 ;; + elliptic_ec(1-z)*elliptic_kc(z) = %pi/2;
1650 ;;
1651 ;; Let z = 1/2 to get
1652 ;;
1653 ;; %pi^(3/2)*'elliptic_ec(1/2)/gamma(3/4)^2-%pi^3/(4*gamma(3/4)^4) = %pi/2
1654 ;;
1655 ;; since we know that elliptic_kc(1/2) = %pi^(3/2)/(2*gamma(3/4)^2). Hence
1656 ;;
1657 ;; elliptic_ec(1/2)
1658 ;; = (2*%pi*gamma(3/4)^4+%pi^3)/(4*%pi^(3/2)*gamma(3/4)^2)
1659 ;; = gamma(3/4)^2/(2*sqrt(%pi))+%pi^(3/2)/(4*gamma(3/4)^2)
1660 ;;
1666 ;; elliptic_ec(-1). Use the identity
1667 ;; http://functions.wolfram.com/08.01.17.0002.01
1668 ;;
1669 ;;
1670 ;; elliptic_ec(z) = sqrt(1 - z)*elliptic_ec(z/(z-1))
1671 ;;
1672 ;; Let z = -1 to get
1673 ;;
1674 ;; elliptic_ec(-1) = sqrt(2)*elliptic_ec(1/2)
1675 ;;
1676 ;; Should we expand out elliptic_ec(1/2) using the above result?
1679 ($hypergeometric_representation
1680 ;; See http://functions.wolfram.com/08.01.26.0001.01
1681 ;;
1682 ;; elliptic_ec(z) = %pi/2*%f[2,1]([-1/2,1/2],[1], z)
1683 ;;
1688 m)))
1690 ;; Nothing to do
1698 m)))
1700 grad)
1702 ;;
1703 ;; Elliptic integral of the third kind:
1704 ;;
1705 ;; (A&S 17.2.14)
1706 ;;
1707 ;; phi
1708 ;; /
1709 ;; [ 1
1710 ;; PI(n;phi|m) = I ----------------------------------- ds
1711 ;; ] 2 2
1712 ;; / SQRT(1 - m SIN (s)) (1 - n SIN (s))
1713 ;; 0
1714 ;;
1715 ;; As with E and F, we do not use the modular angle alpha but the
1716 ;; parameter m = sin(alpha)^2.
1717 ;;
1722 ;; Numerically evaluate it
1726 args
1736 args
1743 ;; 3 cases depending on n < 1, n > 1, or n = 1.
1746 ($positive
1750 ($negative
1754 ($zero
1757 ;; Nothing to do
1760 ;; Complete elliptic-pi. That is phi = %pi/2. Then
1761 ;; elliptic_pi(n,m)
1762 ;; = Rf(0, 1-m,1) + Rj(0,1-m,1-n)*n/3;
1767 ;; To compute elliptic_pi for all z, we use the property
1768 ;; (http://functions.wolfram.com/08.06.16.0002.01)
1769 ;;
1770 ;; elliptic_pi(n, z + %pi*k, m)
1771 ;; = 2*k*elliptic_pi(n, %pi/2, m) + elliptic_pi(n, z, m)
1772 ;;
1773 ;; So we are left with computing the integral for 0 <= z < %pi. Using
1774 ;; Carlson's formulation produces the wrong values for %pi/2 < z <
1775 ;; %pi. How to do that?
1776 ;;
1777 ;; Let
1778 ;;
1779 ;; I(a,b) = integrate(1/(1-n*sin(x)^2)/sqrt(1 - m*sin(x)^2), x, a, b)
1780 ;;
1781 ;; That is, I(a,b) is the integral for the elliptic_pi function but
1782 ;; with a lower limit of a and an upper limit of b.
1783 ;;
1784 ;; Then, we want to compute I(0, z), with %pi <= z < %pi. Let w = z +
1785 ;; %pi/2, 0 <= w < %pi/2. Then
1786 ;;
1787 ;; I(0, w+%pi/2) = I(0, %pi/2) + I(%pi/2, w+%pi/2)
1788 ;;
1789 ;; To evaluate I(%pi/2, w+%pi/2), use a change of variables:
1790 ;;
1791 ;; changevar('integrate(1/(1-n*sin(x)^2)/sqrt(1 - m*sin(x)^2), x, %pi/2, w + %pi/2),
1792 ;; x-%pi+u,u,x)
1793 ;;
1794 ;; = integrate(-1/(sqrt(1-m*sin(u)^2)*(1-n*sin(u)^2)),u,%pi/2-w,%pi/2)
1795 ;; = I(%pi/2-w,%pi/2)
1796 ;; = I(0,%pi/2) - I(0,%pi/2-w)
1797 ;;
1798 ;; Thus,
1799 ;;
1800 ;; I(0,%pi/2+w) = 2*I(0,%pi/2) - I(0,%pi/2-w)
1801 ;;
1802 ;; This allows us to compute the general result with 0 <= z < %pi
1803 ;;
1804 ;; I(0, k*%pi + z) = 2*k*I(0,%pi/2) + I(0,z);
1805 ;;
1806 ;; If 0 <= z < %pi/2, then the we are done. If %pi/2 <= z < %pi, let
1807 ;; z = w+%pi/2. Then
1808 ;;
1809 ;; I(0,z) = 2*I(0,%pi/2) - I(0,%pi/2-w)
1810 ;;
1811 ;; Or, since w = z-%pi/2:
1812 ;;
1813 ;; I(0,z) = 2*I(0,%pi/2) - I(0,%pi-z)
1816 ;; elliptic_pi(n, -phi, m) = -elliptic_pi(n, phi, m). That is, it
1817 ;; is an odd function of phi.
1821 ;; Note: Carlson's DRJ has n defined as the negative of the n given
1822 ;; in A&S.
1824 ;; elliptic_pi(n,phi,m) =
1825 ;; sin(phi)*Rf(cos(phi)^2, 1-m*sin(phi)^2, 1)
1826 ;; - (-n / 3) * sin(phi)^3
1827 ;; * Rj(cos(phi)^2, 1-m*sin(phi)^2, 1, 1-n*sin(phi)^2)
1838 ;; FIXME: Reducing the arg by pi has significant round-off.
1839 ;; Consider doing something better.
1846 ;;; Deriviatives from functions.wolfram.com
1847 ;;; http://functions.wolfram.com/EllipticIntegrals/EllipticPi3/20/
1850 ;Derivative wrt first argument
1869 ;derivative wrt second argument
1876 ;Derivative wrt third argument
1887 grad)
1890 ;; Translation of Jim FitzSimons' bigfloat implementation of elliptic
1891 ;; integrals from http://www.getnet.com/~cherry/elliptbf3.mac.
1892 ;;
1893 ;; The algorithms are based on B.C. Carlson's "Numerical Computation
1894 ;; of Real or Complex Elliptic Integrals". These are updated to the
1895 ;; algorithms in Journal of Computational and Applied Mathematics 118
1896 ;; (2000) 71-85 "Reduction Theorems for Elliptic Integrands with the
1897 ;; Square Root of two quadritic factors"
1898 ;;
1900 ;; NOTE: Despite the names indicating these are for bigfloat numbers,
1901 ;; the algorithms and routines are generic and will work with floats
1902 ;; and bigfloats.
1905 ;; Compute error tolerance as sqrt(2^(-fpprec)). Not sure this is
1906 ;; quite right, but it makes the routines more accurate as fpprec
1907 ;; increases.
1912 args))))
1914 ;; rc(x,y) = integrate(1/2*(t+x)^(-1/2)/(t+y), t, 0, inf)
1915 ;;
1916 ;; log(x) = (x-1)*rc(((1+x)/2)^2, x), x > 0
1917 ;; asin(x) = x * rc(1-x^2, 1), |x|<= 1
1918 ;; acos(x) = sqrt(1-x^2)*rc(x^2,1), 0 <= x <=1
1919 ;; atan(x) = x * rc(1,1+x^2)
1920 ;; asinh(x) = x * rc(1+x^2,1)
1921 ;; acosh(x) = sqrt(x^2-1) * rc(x^2,1), x >= 1
1922 ;; atanh(x) = x * rc(1,1-x^2), |x|<=1
1926 xn z w a an pwr4 n epslon lambda sn s)
1943 do
1950 ;; c2=3/10,c3=1/7,c4=3/8,c5=9/22,c6=159/208,c7=9/8
1963 ;; See https://dlmf.nist.gov/19.16.E5:
1964 ;;
1965 ;; rd(x,y,z) = integrate(3/2/sqrt(t+x)/sqrt(t+y)/sqrt(t+z)/(t+z), t, 0, inf)
1966 ;;
1967 ;; rd(1,1,1) = 1
1968 ;; E(K) = rf(0, 1-K^2, 1) - (K^2/3)*rd(0,1-K^2,1)
1969 ;;
1970 ;; B = integrate(s^2/sqrt(1-s^4), s, 0 ,1)
1971 ;; = beta(3/4,1/2)/4
1972 ;; = sqrt(%pi)*gamma(3/4)/gamma(1/4)
1973 ;; = 1/3*rd(0,2,1)
1988 xnroot ynroot znroot lam)
1990 do
2005 ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26
2012 zndev))
2031 ;; See https://dlmf.nist.gov/19.16.E1
2032 ;;
2033 ;; rf(x,y,z) = 1/2*integrate(1/(sqrt(t+x)*sqrt(t+y)*sqrt(t+z)), t, 0, inf);
2034 ;;
2035 ;; rf(1,1,1) = 1
2048 xnroot ynroot znroot lam)
2050 do
2063 ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26
2094 xnroot ynroot znroot pnroot lam dn)
2096 do
2108 (/ (* power4
2110 dn)))
2133 pndev))
2181 ;; elliptic_f(phi,m) = sin(phi)*rf(cos(phi)^2, 1-m*sin(phi)^2,1)
2185 ;; F(z|1) = log(tan(z/2+%pi/4))
2191 ;; Handle periodicity (see elliptic-f)
2196 0
2199 ;; elliptic_kc(k) = rf(0, 1-k^2,1)
2200 ;;
2201 ;; or
2202 ;; elliptic_kc(m) = rf(0, 1-m,1)
2215 ;; elliptic_e(phi, k) = sin(phi)*rf(cos(phi)^2,1-k^2*sin(phi)^2,1)
2216 ;; - (k^2/3)*sin(phi)^3*rd(cos(phi)^2, 1-k^2*sin(phi)^2,1)
2217 ;;
2218 ;;
2219 ;; or
2220 ;; elliptic_e(phi, m) = sin(phi)*rf(cos(phi)^2,1-m*sin(phi)^2,1)
2221 ;; - (m/3)*sin(phi)^3*rd(cos(phi)^2, 1-m*sin(phi)^2,1)
2222 ;;
2231 ;; Elliptic E is quasi-periodic wrt to phi:
2232 ;;
2233 ;; E(z|m) = E(z - %pi*round(Re(z)/%pi)|m) + 2*round(Re(z)/%pi)*E(m)
2240 ;; elliptic_ec(k) = rf(0,1-k^2,1) - (k^2/3)*rd(0,1-k^2,1);
2241 ;;
2242 ;; or
2243 ;; elliptic_ec(m) = rf(0,1-m,1) - (m/3)*rd(0,1-m,1);
2264 ;; Note: Carlson's DRJ has n defined as the negative of the n given
2265 ;; in A&S.
2277 ;; FIXME: Reducing the arg by pi has significant round-off.
2278 ;; Consider doing something better.
2286 ;; Compute inverse_jacobi_sn, for float or bigfloat args.
2292 ;; Compute inverse_jacobi_dn. We use the following identity
2293 ;; from Gradshteyn & Ryzhik, 8.153.6
2294 ;;
2295 ;; w = dn(z|m) = cn(sqrt(m)*z, 1/m)
2296 ;;
2297 ;; Solve for z to get
2298 ;;
2299 ;; z = inverse_jacobi_dn(w,m)
2300 ;; = 1/sqrt(m) * inverse_jacobi_cn(w, 1/m)
2305 ;; jacobi_dn(x,1) = sech(x) so the inverse is asech(x)
2308 ;; We should do something better to make sure that things
2309 ;; that should be real are real.
2317 ;; Define Carlson's elliptic integrals.
2323 ;; If z is a complex rational, convert to a
2324 ;; complex double-float. Otherwise, leave it as
2325 ;; is. If we don't do this, %i is handled as
2326 ;; #c(0 1), which makes bf-rc use single-float
2327 ;; arithmetic instead of the desired
2328 ;; double-float.
2332 z)))
2335 ;; See comments from bf-rc
2342 args
2346 args
2350 ;; rc(0, 1) = %pi/2
2354 ;; rc(0,1/4) = %pi
2358 ;; rc(2,1) = 1/2*integrate(1/sqrt(t+2)/(t+1), t, 0, inf)
2359 ;; = (log(sqrt(2)+1)-log(sqrt(2)-1))/2
2360 ;; ratsimp(logcontract(%)),algebraic:
2361 ;; = -log(3-2^(3/2))/2
2362 ;; = -log(sqrt(3-2^(3/2)))
2363 ;; = -log(sqrt(2)-1)
2364 ;; = log(1/(sqrt(2)-1))
2365 ;; ratsimp(%),algebraic;
2366 ;; = log(sqrt(2)+1)
2370 ;; rc(%i, %i+1) = 1/2*integrate(1/sqrt(t+%i)/(t+%i+1), t, 0, inf)
2371 ;; = %pi/2-atan((-1)^(1/4))
2372 ;; ratsimp(logcontract(ratsimp(rectform(%o42)))),algebraic;
2373 ;; = (%i*log(3-2^(3/2))+%pi)/4
2374 ;; = (%i*log(3-2^(3/2)))/4+%pi/4
2375 ;; = %i*log(sqrt(3-2^(3/2)))/2+%pi/4
2376 ;; sqrtdenest(%);
2377 ;; = %pi/4 + %i*log(sqrt(2)-1)/2
2384 ;; rc(0,%i) = 1/2*integrate(1/(sqrt(t)*(t+%i)), t, 0, inf)
2385 ;; = -((sqrt(2)*%i-sqrt(2))*%pi)/4
2386 ;; = ((1-%i)*%pi)/2^(3/2)
2392 ;; carlson_rc(x,x) = 1/2*integrate(1/sqrt(t+x)/(t+x), t, 0, inf)
2393 ;; = 1/sqrt(x)
2397 ;; Rc(((1+x)/2)^2,x) = log(x)/(x-1) for x > 0.
2398 ;;
2399 ;; This is done by looking at Rc(x,y) and seeing if
2400 ;; ((1+y)/2)^2 is the same as x.
2412 ;; See https://dlmf.nist.gov/19.20.E18
2416 ;; Rd(1,1,1) = 1
2420 ;; Rd(x,x,x) = x^(-3/2)
2424 ;; Rd(0,y,y) = 3/4*%pi*y^(-3/2)
2426 '$%pi
2429 ;; Rd(x,y,y) = 3/(2*(y-x))*(Rc(x, y) - sqrt(x)/y)
2433 y))))
2435 ;; Rd(x,x,z) = 3/(z-x)*(Rc(z,x) - 1/sqrt(z))
2444 ;; Rd(0,2,1) = 3*(gamma(3/4)^2)/sqrt(2*%pi)
2445 ;; See https://dlmf.nist.gov/19.20.E22.
2446 ;;
2447 ;; But that's the same as
2448 ;; 3*sqrt(%pi)*gamma(3/4)/gamma(1/4). We can see this by
2449 ;; taking the ratio to get
2450 ;; gamma(1/4)*gamma(3/4)/sqrt(2)*%pi. But
2451 ;; gamma(1/4)*gamma(3/4) = beta(1/4,3/4) = sqrt(2)*%pi.
2452 ;; Hence, the ratio is 1.
2453 ;;
2454 ;; Note also that Rd(x,y,z) = Rd(y,x,z)
2461 ;; 1/3*m*Rd(0,1-m,1) = K(m) - E(m).
2462 ;; See https://dlmf.nist.gov/19.25.E1
2463 ;;
2464 ;; Thus, Rd(0,y,1) = 3/(1-y)*(K(1-y) - E(1-y))
2465 ;;
2466 ;; Note that Rd(x,y,z) = Rd(y,x,z).
2475 ;; 1/3*m*(1-m)*Rd(0,1,1-m) = E(m) - (1-m)*K(m)
2476 ;; See https://dlmf.nist.gov/19.25.E1
2477 ;;
2478 ;; Thus
2479 ;; Rd(0,1,z) = 3/(z*(1-z))*(E(1-z) - z*K(1-z))
2480 ;; Recall that Rd(x,y,z) = Rd(y,x,z).
2491 args
2495 args
2506 ;; See https://dlmf.nist.gov/19.20.i
2509 ;; Rf(x,x,x) = x^(-1/2)
2513 ;; Rf(0,y,y) = 1/2*%pi*y^(-1/2)
2520 args
2530 ;; Rf(0,1,2) = (gamma(1/4))^2/(4*sqrt(2*%pi))
2531 ;;
2532 ;; And Rf is symmetric in all the args, so check every
2533 ;; permutation too. This could probably be simplified
2534 ;; without consing all the lists, but I'm lazy.
2539 args
2549 ;; rf(%i, -%i, 0)
2550 ;; = 1/2*integrate(1/sqrt(t^2+1)/sqrt(t),t,0,inf)
2551 ;; = beta(1/4,1/4)/4;
2552 ;; makegamma(%)
2553 ;; = gamma(1/4)^2/(4*sqrt(%pi))
2554 ;;
2555 ;; Rf is symmetric, so check all the permutations too.
2561 args
2562 ;; Check that x = 0 and z = 1, and
2563 ;; return y.
2566 y)))
2573 ;; Rf(0,1-m,1) = elliptic_kc(m).
2574 ;; See https://dlmf.nist.gov/19.25.E1
2578 args
2588 ;; rf(%i, -%i, 0)
2589 ;; = 1/2*integrate(1/sqrt(t^2+1)/sqrt(t),t,0,inf)
2590 ;; = beta(1/4,1/4)/4;
2591 ;; makegamma(%)
2592 ;; = gamma(1/4)^2/(4*sqrt(%pi))
2593 ;;
2594 ;; Rf is symmetric, so check all the permutations too.
2603 args
2607 args
2619 ;; See https://dlmf.nist.gov/19.20.iii
2623 ;; Rj(x,x,x,x) = x^(-3/2)
2626 ;; Rj(x,y,z,z) = Rd(x,y,z)
2630 ;; Rj(0,y,y,p) = 3*%pi/(2*(y*sqrt(p)+p*sqrt(y)))
2636 ;; Rj(x,y,y,p) = 3/(p-y)*(Rc(x,y) - Rc(x,p))
2642 ;; Rj(x,y,y,y) = Rd(x,y,y)
2650 args
2654 args
2659 ;;; Other Jacobian elliptic functions
2661 ;; jacobi_ns(u,m) = 1/jacobi_sn(u,m)
2664 ;; diff wrt u
2667 ;; diff wrt m
2680 m)))))))
2681 grad)
2691 args
2700 args
2705 ;; A&S 16.6.10
2708 ;; A&S 16.6.10
2713 ;; ns is odd
2718 %inverse_jacobi_ns
2719 %inverse_jacobi_cn
2720 %inverse_jacobi_nc
2721 %inverse_jacobi_dn
2722 %inverse_jacobi_nd
2723 %inverse_jacobi_sc
2724 %inverse_jacobi_cs
2725 %inverse_jacobi_sd
2726 %inverse_jacobi_ds
2727 %inverse_jacobi_cd
2728 %inverse_jacobi_dc))
2733 ;; Express in terms of sn:
2734 ;; ns(x) = 1/sn(x)
2736 ;; A&S 16.20 (Jacobi's Imaginary transformation)
2738 ;; ns(i*u) = 1/sn(i*u) = -i/sc(u,m1) = -i*cs(u,m1)
2742 ;; A&S 16.8.10
2743 ;;
2744 ;; ns(m*K+u) = 1/sn(m*K+u)
2745 ;;
2747 coef
2751 ;; ns(4*m*K+u) = ns(u)
2752 ;; ns(0) = infinity
2757 ;; ns(4*m*K + K + u) = ns(K+u) = dc(u)
2758 ;; ns(K) = 1
2760 1
2763 ;; ns(4*m*K+2*K + u) = ns(2*K+u) = -ns(u)
2764 ;; ns(2*K) = infinity
2769 ;; ns(4*m*K+3*K+u) = ns(2*K + K + u) = -ns(K+u) = -dc(u)
2770 ;; ns(3*K) = -1
2772 -1
2780 ;; Nothing to do
2783 ;; jacobi_nc(u,m) = 1/jacobi_cn(u,m)
2786 ;; wrt u
2789 ;; wrt m
2800 grad)
2810 args
2819 args
2826 ;; A&S 16.6.8
2829 ;; A&S 16.6.8
2832 ;; nc is even
2837 %inverse_jacobi_ns
2838 %inverse_jacobi_cn
2839 %inverse_jacobi_nc
2840 %inverse_jacobi_dn
2841 %inverse_jacobi_nd
2842 %inverse_jacobi_sc
2843 %inverse_jacobi_cs
2844 %inverse_jacobi_sd
2845 %inverse_jacobi_ds
2846 %inverse_jacobi_cd
2847 %inverse_jacobi_dc))
2852 ;; Express in terms of cn:
2853 ;; nc(x) = 1/cn(x)
2855 ;; A&S 16.20 (Jacobi's Imaginary transformation)
2857 ;; nc(i*u) = 1/cn(i*u) = 1/nc(u,1-m) = cn(u,1-m)
2860 ;; A&S 16.8.8
2861 ;;
2862 ;; nc(u) = 1/cn(u)
2863 ;;
2865 coef
2869 ;; nc(4*m*K+u) = nc(u)
2870 ;; nc(0) = 1
2872 1
2875 ;; nc(4*m*K+K+u) = nc(K+u) = -ds(u)/sqrt(1-m)
2876 ;; nc(K) = infinity
2882 ;; nc(4*m*K+2*K+u) = nc(2*K+u) = -nc(u)
2883 ;; nc(2*K) = -1
2885 -1
2888 ;; nc(4*m*K+3*K+u) = nc(3*K+u) = nc(2*K+K+u) =
2889 ;; -nc(K+u) = ds(u)/sqrt(1-m)
2890 ;;
2891 ;; nc(3*K) = infinity
2902 ;; Nothing to do
2905 ;; jacobi_nd(u,m) = 1/jacobi_dn(u,m)
2908 ;; wrt u
2911 ;; wrt m
2924 m)))))))
2925 grad)
2935 args
2944 args
2951 ;; A&S 16.6.6
2954 ;; A&S 16.6.6
2957 ;; nd is even
2962 %inverse_jacobi_ns
2963 %inverse_jacobi_cn
2964 %inverse_jacobi_nc
2965 %inverse_jacobi_dn
2966 %inverse_jacobi_nd
2967 %inverse_jacobi_sc
2968 %inverse_jacobi_cs
2969 %inverse_jacobi_sd
2970 %inverse_jacobi_ds
2971 %inverse_jacobi_cd
2972 %inverse_jacobi_dc))
2977 ;; Express in terms of dn:
2978 ;; nd(x) = 1/dn(x)
2980 ;; A&S 16.20 (Jacobi's Imaginary transformation)
2982 ;; nd(i*u) = 1/dn(i*u) = 1/dc(u,1-m) = cd(u,1-m)
2985 ;; A&S 16.8.6
2986 ;;
2988 coef
2990 ;; nd has period 2K
2993 ;; nd(2*m*K+u) = nd(u)
2994 ;; nd(0) = 1
2996 1
2999 ;; nd(2*m*K+K+u) = nd(K+u) = dn(u)/sqrt(1-m)
3000 ;; nd(K) = 1/sqrt(1-m)
3008 ;; Nothing to do
3011 ;; jacobi_sc(u,m) = jacobi_sn/jacobi_cn
3014 ;; wrt u
3017 ;; wrt m
3031 m))))))
3045 m))))))))
3046 grad)
3057 args
3068 args
3076 ;; A&S 16.6.9
3079 ;; A&S 16.6.9
3082 ;; sc is odd
3087 %inverse_jacobi_ns
3088 %inverse_jacobi_cn
3089 %inverse_jacobi_nc
3090 %inverse_jacobi_dn
3091 %inverse_jacobi_nd
3092 %inverse_jacobi_sc
3093 %inverse_jacobi_cs
3094 %inverse_jacobi_sd
3095 %inverse_jacobi_ds
3096 %inverse_jacobi_cd
3097 %inverse_jacobi_dc))
3102 ;; Express in terms of sn and cn
3103 ;; sc(x) = sn(x)/cn(x)
3106 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3108 ;; sc(i*u) = sn(i*u)/cn(i*u) = i*sc(u,m1)/nc(u,m1) = i*sn(u,m1)
3112 ;; A&S 16.8.9
3113 ;; sc(2*m*K+u) = sc(u)
3115 coef
3119 ;; sc(2*m*K+ u) = sc(u)
3120 ;; sc(0) = 0
3122 1
3125 ;; sc(2*m*K + K + u) = sc(K+u)= - cs(u)/sqrt(1-m)
3126 ;; sc(K) = infinity
3134 ;; From A&S 16.3.3 and 16.5.2:
3135 ;; sc(1/2*K) = 1/(1-m)^(1/4)
3140 ;; Nothing to do
3143 ;; jacobi_sd(u,m) = jacobi_sn/jacobi_dn
3146 ;; wrt u
3149 ;; wrt m
3163 m))))))
3177 m))))))))
3178 grad)
3189 args
3200 args
3208 ;; A&S 16.6.5
3211 ;; A&S 16.6.5
3214 ;; sd is odd
3219 %inverse_jacobi_ns
3220 %inverse_jacobi_cn
3221 %inverse_jacobi_nc
3222 %inverse_jacobi_dn
3223 %inverse_jacobi_nd
3224 %inverse_jacobi_sc
3225 %inverse_jacobi_cs
3226 %inverse_jacobi_sd
3227 %inverse_jacobi_ds
3228 %inverse_jacobi_cd
3229 %inverse_jacobi_dc))
3234 ;; Express in terms of sn and dn
3237 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3239 ;; sd(i*u) = sn(i*u)/dn(i*u) = i*sc(u,m1)/dc(u,m1) = i*sd(u,m1)
3243 ;; A&S 16.8.5
3244 ;; sd(4*m*K+u) = sd(u)
3246 coef
3250 ;; sd(4*m*K+u) = sd(u)
3251 ;; sd(0) = 0
3253 0
3256 ;; sd(4*m*K+K+u) = sd(K+u) = cn(u)/sqrt(1-m)
3257 ;; sd(K) = 1/sqrt(m1)
3263 ;; sd(4*m*K+2*K+u) = sd(2*K+u) = -sd(u)
3264 ;; sd(2*K) = 0
3266 0
3269 ;; sd(4*m*K+3*K+u) = sd(3*K+u) = sd(2*K+K+u) =
3270 ;; -sd(K+u) = -cn(u)/sqrt(1-m)
3271 ;; sd(3*K) = -1/sqrt(m1)
3278 ;; jacobi_sn/jacobi_dn
3282 m)
3286 m)))
3290 ;; Nothing to do
3293 ;; jacobi_cs(u,m) = jacobi_cn/jacobi_sn
3296 ;; wrt u
3299 ;; wrt m
3314 m))))))
3327 m))))))))
3328 grad)
3339 args
3350 args
3356 ;; A&S 16.6.12
3359 ;; A&S 16.6.12
3364 ;; cs is odd
3369 %inverse_jacobi_ns
3370 %inverse_jacobi_cn
3371 %inverse_jacobi_nc
3372 %inverse_jacobi_dn
3373 %inverse_jacobi_nd
3374 %inverse_jacobi_sc
3375 %inverse_jacobi_cs
3376 %inverse_jacobi_sd
3377 %inverse_jacobi_ds
3378 %inverse_jacobi_cd
3379 %inverse_jacobi_dc))
3384 ;; Express in terms of cn an sn
3387 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3389 ;; cs(i*u) = cn(i*u)/sn(i*u) = -i*nc(u,m1)/sc(u,m1) = -i*ns(u,m1)
3393 ;; A&S 16.8.12
3394 ;;
3395 ;; cs(2*m*K + u) = cs(u)
3397 coef
3401 ;; cs(2*m*K + u) = cs(u)
3402 ;; cs(0) = infinity
3407 ;; cs(K+u) = -sqrt(1-m)*sc(u)
3408 ;; cs(K) = 0
3410 0
3415 ;; 1/jacobi_sc
3419 m)))
3423 ;; Nothing to do
3426 ;; jacobi_cd(u,m) = jacobi_cn/jacobi_dn
3429 ;; wrt u
3433 ;; wrt m
3448 m))))))
3461 m))))))))
3462 grad)
3473 args
3483 args
3490 ;; A&S 16.6.4
3493 ;; A&S 16.6.4
3496 ;; cd is even
3501 %inverse_jacobi_ns
3502 %inverse_jacobi_cn
3503 %inverse_jacobi_nc
3504 %inverse_jacobi_dn
3505 %inverse_jacobi_nd
3506 %inverse_jacobi_sc
3507 %inverse_jacobi_cs
3508 %inverse_jacobi_sd
3509 %inverse_jacobi_ds
3510 %inverse_jacobi_cd
3511 %inverse_jacobi_dc))
3516 ;; Express in terms of cn and dn
3519 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3521 ;; cd(i*u) = cn(i*u)/dn(i*u) = nc(u,m1)/dc(u,m1) = nd(u,m1)
3524 ;; A&S 16.8.4
3525 ;;
3527 coef
3531 ;; cd(4*m*K + u) = cd(u)
3532 ;; cd(0) = 1
3534 1
3537 ;; cd(4*m*K + K + u) = cd(K+u) = -sn(u)
3538 ;; cd(K) = 0
3540 0
3543 ;; cd(4*m*K + 2*K + u) = cd(2*K+u) = -cd(u)
3544 ;; cd(2*K) = -1
3546 -1
3549 ;; cd(4*m*K + 3*K + u) = cd(2*K + K + u) =
3550 ;; -cd(K+u) = sn(u)
3551 ;; cd(3*K) = 0
3553 0
3557 ;; jacobi_cn/jacobi_dn
3561 m)
3565 m)))
3567 ;; Nothing to do
3570 ;; Nothing to do
3573 ;; jacobi_ds(u,m) = jacobi_dn/jacobi_sn
3576 ;; wrt u
3579 ;; wrt m
3594 m))))))
3607 m))))))))
3608 grad)
3619 args
3630 args
3636 ;; A&S 16.6.11
3639 ;; A&S 16.6.11
3648 %inverse_jacobi_ns
3649 %inverse_jacobi_cn
3650 %inverse_jacobi_nc
3651 %inverse_jacobi_dn
3652 %inverse_jacobi_nd
3653 %inverse_jacobi_sc
3654 %inverse_jacobi_cs
3655 %inverse_jacobi_sd
3656 %inverse_jacobi_ds
3657 %inverse_jacobi_cd
3658 %inverse_jacobi_dc))
3663 ;; Express in terms of dn and sn
3666 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3668 ;; ds(i*u) = dn(i*u)/sn(i*u) = -i*dc(u,m1)/sc(u,m1) = -i*ds(u,m1)
3672 ;; A&S 16.8.11
3674 coef
3678 ;; ds(4*m*K + u) = ds(u)
3679 ;; ds(0) = infinity
3684 ;; ds(4*m*K + K + u) = ds(K+u) = sqrt(1-m)*nc(u)
3685 ;; ds(K) = sqrt(1-m)
3691 ;; ds(4*m*K + 2*K + u) = ds(2*K+u) = -ds(u)
3692 ;; ds(0) = pole
3697 ;; ds(4*m*K + 3*K + u) = ds(2*K + K + u) =
3698 ;; -ds(K+u) = -sqrt(1-m)*nc(u)
3699 ;; ds(3*K) = -sqrt(1-m)
3706 ;; jacobi_dn/jacobi_sn
3710 m)
3713 m)))
3715 ;; Nothing to do
3718 ;; Nothing to do
3721 ;; jacobi_dc(u,m) = jacobi_dn/jacobi_cn
3724 ;; wrt u
3728 ;; wrt m
3742 m))))))
3756 m))))))))
3757 grad)
3768 args
3779 args
3787 ;; A&S 16.6.7
3790 ;; A&S 16.6.7
3797 %inverse_jacobi_ns
3798 %inverse_jacobi_cn
3799 %inverse_jacobi_nc
3800 %inverse_jacobi_dn
3801 %inverse_jacobi_nd
3802 %inverse_jacobi_sc
3803 %inverse_jacobi_cs
3804 %inverse_jacobi_sd
3805 %inverse_jacobi_ds
3806 %inverse_jacobi_cd
3807 %inverse_jacobi_dc))
3812 ;; Express in terms of dn and cn
3815 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3817 ;; dc(i*u) = dn(i*u)/cn(i*u) = dc(u,m1)/nc(u,m1) = dn(u,m1)
3820 ;; See A&S 16.8.7
3822 coef
3826 ;; dc(4*m*K + u) = dc(u)
3827 ;; dc(0) = 1
3829 1
3832 ;; dc(4*m*K + K + u) = dc(K+u) = -ns(u)
3833 ;; dc(K) = pole
3838 ;; dc(4*m*K + 2*K + u) = dc(2*K+u) = -dc(u)
3839 ;; dc(2K) = -1
3841 -1
3844 ;; dc(4*m*K + 3*K + u) = dc(2*K + K + u) =
3845 ;; -dc(K+u) = ns(u)
3846 ;; dc(3*K) = ns(0) = inf
3852 ;; jacobi_dn/jacobi_cn
3856 m)
3859 m)))
3861 ;; Nothing to do
3864 ;; Nothing to do
3867 ;;; Other inverse Jacobian functions
3869 ;; inverse_jacobi_ns(x)
3870 ;;
3871 ;; Let u = inverse_jacobi_ns(x). Then jacobi_ns(u) = x or
3872 ;; 1/jacobi_sn(u) = x or
3873 ;;
3874 ;; jacobi_sn(u) = 1/x
3875 ;;
3876 ;; so u = inverse_jacobi_sn(1/x)
3879 ;; Whittaker and Watson, example in 22.122
3880 ;; inverse_jacobi_ns(u,m) = integrate(1/sqrt(t^2-1)/sqrt(t^2-m), t, u, inf)
3881 ;; -> -1/sqrt(x^2-1)/sqrt(x^2-m)
3887 ;; wrt m
3888 ; ((%derivative) ((%inverse_jacobi_ns) x m) m 1)
3889 nil)
3890 grad)
3896 ;; Numerically evaluate asn
3897 ;;
3898 ;; ans(x,m) = asn(1/x,m) = F(asin(1/x),m)
3908 args
3912 ;; ans(x,0) = F(asin(1/x),0) = asin(1/x)
3915 ;; ans(x,1) = F(asin(1/x),1) = log(tan(pi/2+asin(1/x)/2))
3925 ;; inverse_jacobi_ns(ns(u)) = u
3928 ;; Nothing to do
3931 ;; inverse_jacobi_nc(x)
3932 ;;
3933 ;; Let u = inverse_jacobi_nc(x). Then jacobi_nc(u) = x or
3934 ;; 1/jacobi_cn(u) = x or
3935 ;;
3936 ;; jacobi_cn(u) = 1/x
3937 ;;
3938 ;; so u = inverse_jacobi_cn(1/x)
3941 ;; Whittaker and Watson, example in 22.122
3942 ;; inverse_jacobi_nc(u,m) = integrate(1/sqrt(t^2-1)/sqrt((1-m)*t^2+m), t, 1, u)
3943 ;; -> 1/sqrt(x^2-1)/sqrt((1-m)*x^2+m)
3950 ;; wrt m
3951 ; ((%derivative) ((%inverse_jacobi_nc) x m) m 1)
3952 nil)
3953 grad)
3960 ;;
3970 ;; inverse_jacobi_nc(nc(u)) = u
3973 ;; Nothing to do
3976 ;; inverse_jacobi_nd(x)
3977 ;;
3978 ;; Let u = inverse_jacobi_nd(x). Then jacobi_nd(u) = x or
3979 ;; 1/jacobi_dn(u) = x or
3980 ;;
3981 ;; jacobi_dn(u) = 1/x
3982 ;;
3983 ;; so u = inverse_jacobi_dn(1/x)
3986 ;; Whittaker and Watson, example in 22.122
3987 ;; inverse_jacobi_nd(u,m) = integrate(1/sqrt(t^2-1)/sqrt(1-(1-m)*t^2), t, 1, u)
3988 ;; -> 1/sqrt(u^2-1)/sqrt(1-(1-m)*t^2)
3996 ;; wrt m
3997 ; ((%derivative) ((%inverse_jacobi_nd) x m) m 1)
3998 nil)
3999 grad)
4010 ;; jacobi_nd(1/sqrt(1-m),m) = K(m). This follows from
4011 ;; jacobi_dn(sqrt(1-m),m) = K(m).
4017 ;; inverse_jacobi_nd(nd(u)) = u
4020 ;; Nothing to do
4023 ;; inverse_jacobi_sc(x)
4024 ;;
4025 ;; Let u = inverse_jacobi_sc(x). Then jacobi_sc(u) = x or
4026 ;; x = jacobi_sn(u)/jacobi_cn(u)
4027 ;;
4028 ;; x^2 = sn^2/cn^2
4029 ;; = sn^2/(1-sn^2)
4030 ;;
4031 ;; so
4032 ;;
4033 ;; sn^2 = x^2/(1+x^2)
4034 ;;
4035 ;; sn(u) = x/sqrt(1+x^2)
4036 ;;
4037 ;; u = inverse_sn(x/sqrt(1+x^2))
4038 ;;
4041 ;; Whittaker and Watson, example in 22.122
4042 ;; inverse_jacobi_sc(u,m) = integrate(1/sqrt(1+t^2)/sqrt(1+(1-m)*t^2), t, 0, u)
4043 ;; -> 1/sqrt(1+x^2)/sqrt(1+(1-m)*x^2)
4051 ;; wrt m
4052 ; ((%derivative) ((%inverse_jacobi_sc) x m) m 1)
4053 nil)
4054 grad)
4063 m))
4065 ;; jacobi_sc(0,m) = 0
4071 ;; inverse_jacobi_sc(sc(u)) = u
4074 ;; Nothing to do
4077 ;; inverse_jacobi_sd(x)
4078 ;;
4079 ;; Let u = inverse_jacobi_sd(x). Then jacobi_sd(u) = x or
4080 ;; x = jacobi_sn(u)/jacobi_dn(u)
4081 ;;
4082 ;; x^2 = sn^2/dn^2
4083 ;; = sn^2/(1-m*sn^2)
4084 ;;
4085 ;; so
4086 ;;
4087 ;; sn^2 = x^2/(1+m*x^2)
4088 ;;
4089 ;; sn(u) = x/sqrt(1+m*x^2)
4090 ;;
4091 ;; u = inverse_sn(x/sqrt(1+m*x^2))
4092 ;;
4095 ;; Whittaker and Watson, example in 22.122
4096 ;; inverse_jacobi_sd(u,m) = integrate(1/sqrt(1-(1-m)*t^2)/sqrt(1+m*t^2), t, 0, u)
4097 ;; -> 1/sqrt(1-(1-m)*x^2)/sqrt(1+m*x^2)
4104 ;; wrt m
4105 ; ((%derivative) ((%inverse_jacobi_sd) x m) m 1)
4106 nil)
4107 grad)
4116 m))
4120 ;; inverse_jacobi_sd(1/sqrt(1-m), m) = elliptic_kc(m)
4121 ;;
4122 ;; We can see this from inverse_jacobi_sd(x,m) =
4123 ;; inverse_jacobi_sn(x/sqrt(1+m*x^2), m). So
4124 ;; inverse_jacobi_sd(1/sqrt(1-m),m) = inverse_jacobi_sn(1,m)
4130 ;; inverse_jacobi_sd(sd(u)) = u
4133 ;; Nothing to do
4136 ;; inverse_jacobi_cs(x)
4137 ;;
4138 ;; Let u = inverse_jacobi_cs(x). Then jacobi_cs(u) = x or
4139 ;; 1/x = 1/jacobi_cs(u) = jacobi_sc(u)
4140 ;;
4141 ;; u = inverse_sc(1/x)
4142 ;;
4145 ;; Whittaker and Watson, example in 22.122
4146 ;; inverse_jacobi_cs(u,m) = integrate(1/sqrt(t^2+1)/sqrt(t^2+(1-m)), t, u, inf)
4147 ;; -> -1/sqrt(x^2+1)/sqrt(x^2+(1-m))
4155 ;; wrt m
4156 ; ((%derivative) ((%inverse_jacobi_cs) x m) m 1)
4157 nil)
4158 grad)
4169 ;; Nothing to do
4172 ;; inverse_jacobi_cd(x)
4173 ;;
4174 ;; Let u = inverse_jacobi_cd(x). Then jacobi_cd(u) = x or
4175 ;; x = jacobi_cn(u)/jacobi_dn(u)
4176 ;;
4177 ;; x^2 = cn^2/dn^2
4178 ;; = (1-sn^2)/(1-m*sn^2)
4179 ;;
4180 ;; or
4181 ;;
4182 ;; sn^2 = (1-x^2)/(1-m*x^2)
4183 ;;
4184 ;; sn(u) = sqrt(1-x^2)/sqrt(1-m*x^2)
4185 ;;
4186 ;; u = inverse_sn(sqrt(1-x^2)/sqrt(1-m*x^2))
4187 ;;
4190 ;; Whittaker and Watson, example in 22.122
4191 ;; inverse_jacobi_cd(u,m) = integrate(1/sqrt(1-t^2)/sqrt(1-m*t^2), t, u, 1)
4192 ;; -> -1/sqrt(1-x^2)/sqrt(1-m*x^2)
4200 ;; wrt m
4201 ; ((%derivative) ((%inverse_jacobi_cd) x m) m 1)
4202 nil)
4203 grad)
4212 m)))
4221 ;; inverse_jacobi_cd(cd(u)) = u
4224 ;; Nothing to do
4227 ;; inverse_jacobi_ds(x)
4228 ;;
4229 ;; Let u = inverse_jacobi_ds(x). Then jacobi_ds(u) = x or
4230 ;; 1/x = 1/jacobi_ds(u) = jacobi_sd(u)
4231 ;;
4232 ;; u = inverse_sd(1/x)
4233 ;;
4236 ;; Whittaker and Watson, example in 22.122
4237 ;; inverse_jacobi_ds(u,m) = integrate(1/sqrt(t^2-(1-m))/sqrt(t^2+m), t, u, inf)
4238 ;; -> -1/sqrt(x^2-(1-m))/sqrt(x^2+m)
4246 ;; wrt m
4247 ; ((%derivative) ((%inverse_jacobi_ds) x m) m 1)
4248 nil)
4249 grad)
4260 ;; inverse_jacobi_ds(sqrt(1-m),m) = elliptic_kc(m)
4261 ;;
4262 ;; Since inverse_jacobi_ds(sqrt(1-m), m) =
4263 ;; inverse_jacobi_sd(1/sqrt(1-m),m). And we know from
4264 ;; above that this is elliptic_kc(m)
4270 ;; inverse_jacobi_ds(ds(u)) = u
4273 ;; Nothing to do
4277 ;; inverse_jacobi_dc(x)
4278 ;;
4279 ;; Let u = inverse_jacobi_dc(x). Then jacobi_dc(u) = x or
4280 ;; 1/x = 1/jacobi_dc(u) = jacobi_cd(u)
4281 ;;
4282 ;; u = inverse_cd(1/x)
4283 ;;
4286 ;; Note: Whittaker and Watson, example in 22.122 says
4287 ;; inverse_jacobi_dc(u,m) = integrate(1/sqrt(t^2-1)/sqrt(t^2-m),
4288 ;; t, u, 1) but that seems wrong. A&S 17.4.47 says
4289 ;; integrate(1/sqrt(t^2-1)/sqrt(t^2-m), t, a, u) =
4290 ;; inverse_jacobi_cd(x,m). Lawden 3.2.8 says the same.
4291 ;; functions.wolfram.com says the derivative is
4292 ;; 1/sqrt(t^2-1)/sqrt(t^2-m).
4302 ;; wrt m
4303 ; ((%derivative) ((%inverse_jacobi_dc) x m) m 1)
4304 nil)
4305 grad)
4317 ;; inverse_jacobi_dc(dc(u)) = u
4320 ;; Nothing to do
4323 ;; Convert an inverse Jacobian function into the equivalent elliptic
4324 ;; integral F.
4325 ;;
4326 ;; See A&S 17.4.41-17.4.52.
4329 e)
4331 %inverse_jacobi_nd %inverse_jacobi_dn
4332 %inverse_jacobi_sn %inverse_jacobi_cd
4333 %inverse_jacobi_dc %inverse_jacobi_ns
4334 %inverse_jacobi_nc %inverse_jacobi_ds
4335 %inverse_jacobi_sd %inverse_jacobi_cn))
4336 ;; We have some inverse Jacobi function. Convert it to the F form.
4338 e
4341 (%inverse_jacobi_sc
4342 ;; A&S 17.4.41
4344 (%inverse_jacobi_cs
4345 ;; A&S 17.4.42
4347 (%inverse_jacobi_nd
4348 ;; A&S 17.4.43
4355 m))
4356 (%inverse_jacobi_dn
4357 ;; A&S 17.4.44
4362 m))
4363 (%inverse_jacobi_sn
4364 ;; A&S 17.4.45
4366 (%inverse_jacobi_cd
4367 ;; A&S 17.4.46
4373 m))
4374 (%inverse_jacobi_dc
4375 ;; A&S 17.4.47
4381 m))
4382 (%inverse_jacobi_ns
4383 ;; A&S 17.4.48
4385 (%inverse_jacobi_nc
4386 ;; A&S 17.4.49
4388 (%inverse_jacobi_ds
4389 ;; A&S 17.4.50
4394 m))
4395 (%inverse_jacobi_sd
4396 ;; A&S 17.4.51
4402 m))
4403 (%inverse_jacobi_cn
4404 ;; A&S 17.4.52
4411 e
4425 e
4429 ;; Eu(u,m) = integrate(jacobi_dn(v,m)^2,v,0,u)
4430 ;; = integrate(sqrt((1-m*t^2)/(1-t^2)),t,0,jacobi_sn(u,m))
4431 ;;
4432 ;; Eu(u,m) = E(am(u),m)
4433 ;;
4434 ;; where E(u,m) is elliptic-e above.
4435 ;;
4436 ;; Checks.
4437 ;; Lawden gives the following relationships
4438 ;;
4439 ;; E(u+v) = E(u) + E(v) - m*sn(u)*sn(v)*sn(u+v)
4440 ;; E(u,0) = u, E(u,1) = tanh u
4441 ;;
4442 ;; E(i*u,k) = i*sc(u,k')*dn(u,k') - i*E(u,k') + i*u
4443 ;;
4444 ;; E(2*i*K') = 2*i*(K'-E')
4445 ;;
4446 ;; E(u + 2*i*K') = E(u) + 2*i*(K' - E')
4447 ;;
4448 ;; E(u+K) = E(u) + E - k^2*sn(u)*cd(u)
4451 ;; E(u + 2*n*K) = E(u) + 2*n*E
4456 ;; 0 <= u-rem < 2*K
4459 ;; 0 <= u-rem < K so
4460 ;; E(u + K) = E(u) + E - m*sn(u)*cd(u)
4463 ell-e)
4469 ;; From Lawden:
4470 ;;
4471 ;; E(u+i*v, m) = E(u,m) -i*E(v,m') + i*v + i*sc(v,m')*dn(v,m')
4472 ;; -i*m*sn(u,m)*sc(v,m')*sn(u+i*v,m)
4473 ;;
4489 ;; wrt m
4490 )
4491 grad)
4495 ;; as it stands, ELLIPTIC-EU can't handle bigfloats or complex bigfloats,
4496 ;; so handle only floats and complex floats here.
4507 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4508 ;; Integrals. At present with respect to first argument only.
4509 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4511 ;; A&S 16.24.1: integrate(jacobi_sn(u,m),u)
4512 ;; = log(jacobi_dn(u,m)-sqrt(m)*jacobi_cn(u,m))/sqrt(m)
4521 nil)
4522 integral)
4524 ;; A&S 16.24.2: integrate(jacobi_cn(u,m),u) = acos(jacobi_dn(u,m))/sqrt(m)
4529 nil)
4530 integral)
4532 ;; A&S 16.24.3: integrate(jacobi_dn(u,m),u) = asin(jacobi_sn(u,m))
4536 nil)
4537 integral)
4539 ;; A&S 16.24.4: integrate(jacobi_cd(u,m),u)
4540 ;; = log(jacobi_nd(u,m)+sqrt(m)*jacobi_sd(u,m))/sqrt(m)
4548 nil)
4549 integral)
4551 ;; integrate(jacobi_sd(u,m),u)
4552 ;;
4553 ;; A&S 16.24.5 gives
4554 ;; asin(-sqrt(m)*jacobi_cd(u,m))/sqrt(m*m_1), where m + m_1 = 1
4555 ;; but this does not pass some simple tests.
4556 ;;
4557 ;; functions.wolfram.com 09.35.21.001.01 gives
4558 ;; -asin(sqrt(m)*jacobi_cd(u,m))*sqrt(1-m*jacobi_cd(u,m)^2)*jacobi_dn(u,m)/((1-m)*sqrt(m))
4559 ;; and this does pass.
4573 nil)
4574 integral)
4576 ;; integrate(jacobi_nd(u,m),u)
4577 ;;
4578 ;; A&S 16.24.6 gives
4579 ;; acos(jacobi_cd(u,m))/sqrt(m_1), where m + m_1 = 1
4580 ;; but this does not pass some simple tests.
4581 ;;
4582 ;; functions.wolfram.com 09.32.21.0001.01 gives
4583 ;; sqrt(1-jacobi_cd(u,m)^2)*acos(jacobi_cd(u,m))/((1-m)*jacobi_sd(u,m))
4584 ;; and this does pass.
4594 nil)
4595 integral)
4597 ;; A&S 16.24.7: integrate(jacobi_dc(u,m),u) = log(jacobi_nc(u,m)+jacobi_sc(u,m))
4601 nil)
4602 integral)
4604 ;; A&S 16.24.8: integrate(jacobi_nc(u,m),u)
4605 ;; = log(jacobi_dc(u,m)+sqrt(m_1)*jacobi_sc(u,m))/sqrt(m_1), where m + m_1 = 1
4617 nil)
4618 integral)
4620 ;; A&S 16.24.9: integrate(jacobi_sc(u,m),u)
4621 ;; = log(jacobi_dc(u,m)+sqrt(m_1)*jacobi_nc(u,m))/sqrt(m_1), where m + m_1 = 1
4633 nil)
4634 integral)
4636 ;; A&S 16.24.10: integrate(jacobi_ns(u,m),u)
4637 ;; = log(jacobi_ds(u,m)-jacobi_cs(u,m))
4643 nil)
4644 integral)
4646 ;; integrate(jacobi_ds(u,m),u)
4647 ;;
4648 ;; A&S 16.24.11 gives
4649 ;; log(jacobi_ds(u,m)-jacobi_cs(u,m))
4650 ;; but this does not pass some simple tests.
4651 ;;
4652 ;; functions.wolfram.com 09.30.21.0001.01 gives
4653 ;; log((1-jacobi_cn(u,m))/jacobi_sn(u,m))
4654 ;;
4661 nil)
4662 integral)
4664 ;; A&S 16.24.12: integrate(jacobi_cs(u,m),u) = log(jacobi_ns(u,m)-jacobi_ds(u,m))
4670 nil)
4671 integral)
4673 ;; functions.wolfram.com 09.48.21.0001.01
4674 ;; integrate(inverse_jacobi_sn(u,m),u) =
4675 ;; inverse_jacobi_sn(u,m)*u
4676 ;; - log( jacobi_dn(inverse_jacobi_sn(u,m),m)
4677 ;; -sqrt(m)*jacobi_cn(inverse_jacobi_sn(u,m),m)) / sqrt(m)
4687 nil)
4688 integral)
4690 ;; functions.wolfram.com 09.38.21.0001.01
4691 ;; integrate(inverse_jacobi_cn(u,m),u) =
4692 ;; u*inverse_jacobi_cn(u,m)
4693 ;; -%i*log(%i*jacobi_dn(inverse_jacobi_cn(u,m),m)/sqrt(m)
4694 ;; -jacobi_sn(inverse_jacobi_cn(u,m),m))
4695 ;; /sqrt(m)
4706 nil)
4707 integral)
4709 ;; functions.wolfram.com 09.41.21.0001.01
4710 ;; integrate(inverse_jacobi_dn(u,m),u) =
4711 ;; u*inverse_jacobi_dn(u,m)
4712 ;; - %i*log(%i*jacobi_cn(inverse_jacobi_dn(u,m),m)
4713 ;; +jacobi_sn(inverse_jacobi_dn(u,m),m))
4723 nil)
4724 integral)
4726 \f
4727 ;; Real and imaginary part for Jacobi elliptic functions.
4735 ;; We only split on the argument, not the order
4739 ;; Pure real
4753 ;; Let s = jacobi_sn(x,m)
4754 ;; c = jacobi_cn(x,m)
4755 ;; d = jacobi_dn(x,m)
4756 ;; s1 = jacobi_sn(y,1-m)
4757 ;; c1 = jacobi_cn(y,1-m)
4758 ;; d1 = jacobi_dn(y,1-m)
4760 (%jacobi_sn
4761 ;; A&S 16.21.1
4762 ;; jacobi_sn(x+%i*y,m) =
4763 ;;
4764 ;; s*d1 + %i*c*d*s1*c1
4765 ;; -------------------
4766 ;; c1^2+m*s^2*s1^2
4767 ;;
4770 den)))
4771 (%jacobi_cn
4772 ;; A&S 16.21.2
4773 ;;
4774 ;; cn(x+%i_y, m) =
4775 ;;
4776 ;; c*c1 - %i*s*d*s1*d1
4777 ;; -------------------
4778 ;; c1^2+m*s^2*s1^2
4782 den)))
4783 (%jacobi_dn
4784 ;; A&S 16.21.3
4785 ;;
4786 ;; dn(x+%i_y, m) =
4787 ;;
4788 ;; d*c1*d1 - %i*m*s*c*s1
4789 ;; ---------------------
4790 ;; c1^2+m*s^2*s1^2
4792 den)
4794 den))))))))))
4796 \f
4797 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4798 ;; Jacobi amplitude function.
4799 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4803 ;; Arithmetic-Geometric Mean algorithm for real or complex numbers.
4804 ;; See https://dlmf.nist.gov/22.20.ii.
4805 ;;
4806 ;; Do not use this for computing jacobi sn. It loses some 7 digits of
4807 ;; accuracy for sn(1+%i,0.7).
4811 ;; Instead of allocating these array anew each time, we'll reuse
4812 ;; them and allow them to grow as needed.
4814 "Arithmetic-Geometric Mean algorithm for real or complex a0, b0, c0.
4815 Algorithm continues until |c[n]| <= tol."
4817 ;; DLMF (https://dlmf.nist.gov/22.20.ii) says for any real or
4818 ;; complex a0 and b0, b0/a0 must not be real and negative. Let's
4819 ;; check that.
4825 ;; DLMF (https://dlmf.nist.gov/22.20.ii) says that |c[n]| <=
4826 ;; C*2^(-2^n), for some constant C. Solve C*2^(-2^n) = tol to
4827 ;; get n = log(log(C/tol)/log(2))/log(2). Arbitrarily assume C
4828 ;; is one to get n = log(-(log(tol)/log(2)))/log(2). Thus, the
4829 ;; approximate number of term needed is n =
4830 ;; 1.44*log(-(1.44*log(tol))). Round to 2*log(-log2(tol)).
4845 ;; DLMF (https://dlmf.nist.gov/22.20.ii) has conditions on how
4846 ;; to choose the square root depending on the phase of a[n-1]
4847 ;; and b[n-1]. We don't check for that here.
4852 "Evaluate the jacobi_am function from real u and m with |m| <= 1. This
4853 uses the AGM method until a tolerance of TOL is reached for the
4854 error."
4858 ;; See DLMF (https://dlmf.nist.gov/22.20.ii) for the algorithm.
4861 do
4866 phi)))
4868 ;; Compute Jacobi am for real or complex values of U and M. The args
4869 ;; must be floats or bigfloat::bigfloats. TOL is the tolerance used
4870 ;; by the AGM algorithm. It is ignored if the AGM algorithm is not
4871 ;; used.
4874 ;; The case of real u and m with |m| <= 1. We can use AGM to
4875 ;; compute the result.
4878 tol))
4880 ;; Otherwise, use the formula am(u,m) = asin(jacobi_sn(u,m)).
4881 ;; (See DLMF https://dlmf.nist.gov/22.16.E1). This appears
4882 ;; to be what functions.wolfram.com is using in this case.
4890 ;; am(0,m) = 0
4893 ;; See https://dlmf.nist.gov/22.16.E4
4894 ;;
4895 ;; am(u,0) = u
4896 u)
4898 ;; See https://dlmf.nist.gov/22.16.E5. This is equivalent to
4899 ;; the Gudermannian function.
4900 ;;
4901 ;; am(u,1) = 2*atan(exp(u))-%pi/2
4907 double-float-epsilon)))
4910 args
4913 double-float-epsilon))))
4920 args
4925 ;; Nothing to do
4928 ;; Derivative of jacobi_am wrt z and m.
4931 ;; WRT z. From http://functions.wolfram.com/09.24.20.0001.01
4932 ;; jacobi_dn(z,m)
4934 ;; WRT m. From http://functions.wolfram.com/09.24.20.0003.01.
4935 ;; There are 5 different formulas listed; we chose the first,
4936 ;; arbitrarily.
4937 ;;
4938 ;; (((m-1)*z+elliptic_e(jacobi_am(z,m),m))*jacobi_dn(z,m)
4939 ;; - m*jacobi_cn(z,m)*jacobi_sn(z,m))/(2*m*(m-1))
4947 nil)
4948 grad)