?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Add more comments and clean up some typos.
[maxima.git] / src / ellipt.lisp
blob067ab17f667ba796039765ca48a794caba4707c9
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
13
14 (in-package :maxima)
15 ;;(macsyma-module ellipt)
16
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 ;;;
32
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 ;;
48
49 (in-package #:bigfloat)
50
51 (declaim (inline descending-transform ascending-transform))
52
53 (defun ascending-transform (u m)
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).
61 (let* ((root-m (sqrt m))
62 (mu (/ (* 4 root-m)
63 (expt (1+ root-m) 2)))
64 (root-mu1 (/ (- 1 root-m) (+ 1 root-m)))
65 (v (/ u (1+ root-mu1))))
66 (values v mu root-mu1)))
67
68 (defun descending-transform (u m)
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.
73 (let* ((root-m1 (sqrt (- 1 m)))
74 (root-mu (/ (- 1 root-m1) (+ 1 root-m1)))
75 (mu (* root-mu root-mu))
76 (v (/ u (1+ root-mu))))
77 (values v mu root-mu)))
78
79 ;;
80 ;; This appears to work quite well for both real and complex values
81 ;; of u.
82 (defun elliptic-sn-descending (u m)
83 (cond ((= m 1)
84 ;; A&S 16.6.1
85 (tanh u))
86 ((< (abs m) (epsilon u))
87 ;; A&S 16.6.1
88 (sin u))
89 (t
90 (multiple-value-bind (v mu root-mu)
91 (descending-transform u m)
92 (let* ((new-sn (elliptic-sn-descending v mu)))
93 (/ (* (1+ root-mu) new-sn)
94 (1+ (* root-mu new-sn new-sn))))))))
95
96 ;; AGM scale. See A&S 17.6
97 ;;
98 ;; The AGM scale is
99 ;;
100 ;; a[n] = (a[n-1]+b[n-1])/2, b[n] = sqrt(a[n-1]*b[n-1]), c[n] = (a[n-1]-b[n-1])/2.
101 ;;
102 ;; We stop when abs(c[n]) <= 10*eps
103 ;;
104 ;; A list of (n a[n] b[n] c[n]) is returned.
105 (defun agm-scale (a b c)
106 (loop for n from 0
107 while (> (abs c) (* 10 (epsilon c)))
108 collect (list n a b c)
109 do (psetf a (/ (+ a b) 2)
110 b (sqrt (* a b))
111 c (/ (- a b) 2))))
112
113 ;; WARNING: This seems to have accuracy problems when u is complex. I
114 ;; (rtoy) do not know why. For example (jacobi-agm #c(1e0 1e0) .7e0)
115 ;; returns
116 ;;
117 ;; #C(1.134045970915582 0.3522523454566013)
118 ;; #C(0.57149659007575 -0.6989899153338323)
119 ;; #C(0.6229715431044184 -0.4488635962149656)
120 ;;
121 ;; But the actual value of sn(1+%i, .7) is .3522523469224946 %i +
122 ;; 1.134045971912365. We've lost about 7 digits of accuracy!
123 (defun jacobi-agm (u m)
124 ;; A&S 16.4.
125 ;;
126 ;; Compute the AGM scale with a = 1, b = sqrt(1-m), c = sqrt(m).
127 ;;
128 ;; Then phi[N] = 2^N*a[N]*u and compute phi[n] from
129 ;;
130 ;; sin(2*phi[n-1] - phi[n]) = c[n]/a[n]*sin(phi[n])
131 ;;
132 ;; Finally,
133 ;;
134 ;; sn(u|m) = sin(phi[0]), cn(u|m) = cos(phi[0])
135 ;; dn(u|m) = cos(phi[0])/cos(phi[1]-phi[0])
136 ;;
137 ;; Returns the three values sn, cn, dn.
138 (let* ((agm-data (nreverse (rest (agm-scale 1 (sqrt (- 1 m)) (sqrt m)))))
139 (phi (destructuring-bind (n a b c)
140 (first agm-data)
141 (declare (ignore b c))
142 (* a u (ash 1 n))))
143 (phi1 0e0))
144 (dolist (agm agm-data)
145 (destructuring-bind (n a b c)
146 agm
147 (declare (ignore n b))
148 (setf phi1 phi
149 phi (/ (+ phi (asin (* (/ c a) (sin phi)))) 2))))
150 (values (sin phi) (cos phi) (/ (cos phi) (cos (- phi1 phi))))))
151
152 (defun sn (u m)
153 (cond ((zerop m)
154 ;; jacobi_sn(u,0) = sin(u). Should we use A&S 16.13.1 if m
155 ;; is small enough?
156 ;;
157 ;; sn(u,m) = sin(u) - 1/4*m(u-sin(u)*cos(u))*cos(u)
158 (sin u))
159 ((= m 1)
160 ;; jacobi_sn(u,1) = tanh(u). Should we use A&S 16.15.1 if m
161 ;; is close enough to 1?
162 ;;
163 ;; sn(u,m) = tanh(u) + 1/4*(1-m)*(sinh(u)*cosh(u)-u)*sech(u)^2
164 (tanh u))
165 (t
166 ;; Use the ascending Landen transformation to compute sn.
167 (let ((s (elliptic-sn-descending u m)))
168 (if (and (realp u) (realp m))
169 (realpart s)
170 s)))))
171
172 (defun dn (u m)
173 (cond ((zerop m)
174 ;; jacobi_dn(u,0) = 1. Should we use A&S 16.13.3 for small m?
175 ;;
176 ;; dn(u,m) = 1 - 1/2*m*sin(u)^2
177 1)
178 ((= m 1)
179 ;; jacobi_dn(u,1) = sech(u). Should we use A&S 16.15.3 if m
180 ;; is close enough to 1?
181 ;;
182 ;; dn(u,m) = sech(u) + 1/4*(1-m)*(sinh(u)*cosh(u)+u)*tanh(u)*sech(u)
183 (/ (cosh u)))
184 (t
185 ;; Use the Gauss transformation from
186 ;; http://functions.wolfram.com/09.29.16.0013.01:
187 ;;
188 ;;
189 ;; dn((1+sqrt(m))*z, 4*sqrt(m)/(1+sqrt(m))^2)
190 ;; = (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2)
191 ;;
192 ;; So
193 ;;
194 ;; dn(y, mu) = (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2)
195 ;;
196 ;; where z = y/(1+sqrt(m)) and mu=4*sqrt(m)/(1+sqrt(m))^2.
197 ;;
198 ;; Solve for m, and we get
199 ;;
200 ;; sqrt(m) = -(mu+2*sqrt(1-mu)-2)/mu or (-mu+2*sqrt(1-mu)+2)/mu.
201 ;;
202 ;; I don't think it matters which sqrt we use, so I (rtoy)
203 ;; arbitrarily choose the first one above.
204 ;;
205 ;; Note that (1-sqrt(1-mu))/(1+sqrt(1-mu)) is the same as
206 ;; -(mu+2*sqrt(1-mu)-2)/mu. Also, the former is more
207 ;; accurate for small mu.
208 (let* ((root (let ((root-1-m (sqrt (- 1 m))))
209 (/ (- 1 root-1-m)
210 (+ 1 root-1-m))))
211 (z (/ u (+ 1 root)))
212 (s (elliptic-sn-descending z (* root root)))
213 (p (* root s s )))
214 (/ (- 1 p)
215 (+ 1 p))))))
216
217 (defun cn (u m)
218 (cond ((zerop m)
219 ;; jacobi_cn(u,0) = cos(u). Should we use A&S 16.13.2 for
220 ;; small m?
221 ;;
222 ;; cn(u,m) = cos(u) + 1/4*m*(u-sin(u)*cos(u))*sin(u)
223 (cos u))
224 ((= m 1)
225 ;; jacobi_cn(u,1) = sech(u). Should we use A&S 16.15.3 if m
226 ;; is close enough to 1?
227 ;;
228 ;; cn(u,m) = sech(u) - 1/4*(1-m)*(sinh(u)*cosh(u)-u)*tanh(u)*sech(u)
229 (/ (cosh u)))
230 (t
231 ;; Use the ascending Landen transformation, A&S 16.14.3.
232 (multiple-value-bind (v mu root-mu1)
233 (ascending-transform u m)
234 (let ((d (dn v mu)))
235 (* (/ (+ 1 root-mu1) mu)
236 (/ (- (* d d) root-mu1)
237 d)))))))
238
239 (in-package :maxima)
240
241 ;; Tell maxima what the derivatives are.
242 ;;
243 ;; Lawden says the derivative wrt to k but that's not what we want.
244 ;;
245 ;; Here's the derivation we used, based on how Lawden get's his results.
246 ;;
247 ;; Let
248 ;;
249 ;; diff(sn(u,m),m) = s
250 ;; diff(cn(u,m),m) = p
251 ;; diff(dn(u,m),m) = q
252 ;;
253 ;; From the derivatives of sn, cn, dn wrt to u, we have
254 ;;
255 ;; diff(sn(u,m),u) = cn(u)*dn(u)
256 ;; diff(cn(u,m),u) = -cn(u)*dn(u)
257 ;; diff(dn(u,m),u) = -m*sn(u)*cn(u)
258 ;;
259
260 ;; Differentiate these wrt to m:
261 ;;
262 ;; diff(s,u) = p*dn + cn*q
263 ;; diff(p,u) = -p*dn - q*dn
264 ;; diff(q,u) = -sn*cn - m*s*cn - m*sn*q
265 ;;
266 ;; Also recall that
267 ;;
268 ;; sn(u)^2 + cn(u)^2 = 1
269 ;; dn(u)^2 + m*sn(u)^2 = 1
270 ;;
271 ;; Differentiate these wrt to m:
272 ;;
273 ;; sn*s + cn*p = 0
274 ;; 2*dn*q + sn^2 + 2*m*sn*s = 0
275 ;;
276 ;; Thus,
277 ;;
278 ;; p = -s*sn/cn
279 ;; q = -m*s*sn/dn - sn^2/dn/2
280 ;;
281 ;; So
282 ;; diff(s,u) = -s*sn*dn/cn - m*s*sn*cn/dn - sn^2*cn/dn/2
283 ;;
284 ;; or
285 ;;
286 ;; diff(s,u) + s*(sn*dn/cn + m*sn*cn/dn) = -1/2*sn^2*cn/dn
287 ;;
288 ;; diff(s,u) + s*sn/cn/dn*(dn^2 + m*cn^2) = -1/2*sn^2*cn/dn
289 ;;
290 ;; Multiply through by the integrating factor 1/cn/dn:
291 ;;
292 ;; diff(s/cn/dn, u) = -1/2*sn^2/dn^2 = -1/2*sd^2.
293 ;;
294 ;; Integrate this to get
295 ;;
296 ;; s/cn/dn = C + -1/2*int sd^2
297 ;;
298 ;; It can be shown that C is zero.
299 ;;
300 ;; We know that (by differentiating this expression)
301 ;;
302 ;; int dn^2 = (1-m)*u+m*sn*cd + m*(1-m)*int sd^2
303 ;;
304 ;; or
305 ;;
306 ;; int sd^2 = 1/m/(1-m)*int dn^2 - u/m - sn*cd/(1-m)
307 ;;
308 ;; Thus, we get
309 ;;
310 ;; s/cn/dn = u/(2*m) + sn*cd/(2*(1-m)) - 1/2/m/(1-m)*int dn^2
311 ;;
312 ;; or
313 ;;
314 ;; s = 1/(2*m)*u*cn*dn + 1/(2*(1-m))*sn*cn^2 - 1/2/(m*(1-m))*cn*dn*E(u)
315 ;;
316 ;; where E(u) = int dn^2 = elliptic_e(am(u)) = elliptic_e(asin(sn(u)))
317 ;;
318 ;; This is our desired result:
319 ;;
320 ;; s = 1/(2*m)*cn*dn*[u - elliptic_e(asin(sn(u)),m)/(1-m)] + sn*cn^2/2/(1-m)
321 ;;
322 ;;
323 ;; Since diff(cn(u,m),m) = p = -s*sn/cn, we have
324 ;;
325 ;; p = -1/(2*m)*sn*dn[u - elliptic_e(asin(sn(u)),m)/(1-m)] - sn^2*cn/2/(1-m)
326 ;;
327 ;; diff(dn(u,m),m) = q = -m*s*sn/dn - sn^2/dn/2
328 ;;
329 ;; q = -1/2*sn*cn*[u-elliptic_e(asin(sn),m)/(1-m)] - m*sn^2*cn^2/dn/2/(1-m)
330 ;;
331 ;; - sn^2/dn/2
332 ;;
333 ;; = -1/2*sn*cn*[u-elliptic_e(asin(sn),m)/(1-m)] + dn*sn^2/2/(m-1)
334 ;;
335 (defprop %jacobi_sn
336 ((u m)
337 ((mtimes) ((%jacobi_cn) u m) ((%jacobi_dn) u m))
338 ((mplus simp)
339 ((mtimes simp) ((rat simp) 1 2)
340 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
341 ((mexpt simp) ((%jacobi_cn simp) u m) 2) ((%jacobi_sn simp) u m))
342 ((mtimes simp) ((rat simp) 1 2) ((mexpt simp) m -1)
343 ((%jacobi_cn simp) u m) ((%jacobi_dn simp) u m)
344 ((mplus simp) u
345 ((mtimes simp) -1 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
346 ((%elliptic_e simp) ((%asin simp) ((%jacobi_sn simp) u m)) m))))))
347 grad)
348
349 (defprop %jacobi_cn
350 ((u m)
351 ((mtimes simp) -1 ((%jacobi_sn simp) u m) ((%jacobi_dn simp) u m))
352 ((mplus simp)
353 ((mtimes simp) ((rat simp) -1 2)
354 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
355 ((%jacobi_cn simp) u m) ((mexpt simp) ((%jacobi_sn simp) u m) 2))
356 ((mtimes simp) ((rat simp) -1 2) ((mexpt simp) m -1)
357 ((%jacobi_dn simp) u m) ((%jacobi_sn simp) u m)
358 ((mplus simp) u
359 ((mtimes simp) -1 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
360 ((%elliptic_e simp) ((%asin simp) ((%jacobi_sn simp) u m)) m))))))
361 grad)
362
363 (defprop %jacobi_dn
364 ((u m)
365 ((mtimes) -1 m ((%jacobi_sn) u m) ((%jacobi_cn) u m))
366 ((mplus simp)
367 ((mtimes simp) ((rat simp) -1 2)
368 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
369 ((%jacobi_dn simp) u m) ((mexpt simp) ((%jacobi_sn simp) u m) 2))
370 ((mtimes simp) ((rat simp) -1 2) ((%jacobi_cn simp) u m)
371 ((%jacobi_sn simp) u m)
372 ((mplus simp) u
373 ((mtimes simp) -1
374 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
375 ((%elliptic_e simp) ((%asin simp) ((%jacobi_sn simp) u m)) m))))))
376 grad)
377
378 ;; The inverse elliptic functions.
379 ;;
380 ;; F(phi|m) = asn(sin(phi),m)
381 ;;
382 ;; so asn(u,m) = F(asin(u)|m)
383 (defprop %inverse_jacobi_sn
384 ((x m)
385 ;; Lawden 3.1.2:
386 ;; inverse_jacobi_sn(x) = integrate(1/sqrt(1-t^2)/sqrt(1-m*t^2),t,0,x)
387 ;; -> 1/sqrt(1-x^2)/sqrt(1-m*x^2)
388 ((mtimes simp)
389 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
390 ((rat simp) -1 2))
391 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) x 2)))
392 ((rat simp) -1 2)))
393 ;; diff(F(asin(u)|m),m)
394 ((mtimes simp) ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
395 ((mplus simp)
396 ((mtimes simp) -1 x
397 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
398 ((rat simp) 1 2))
399 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) x 2)))
400 ((rat simp) -1 2)))
401 ((mtimes simp) ((mexpt simp) m -1)
402 ((mplus simp) ((%elliptic_e simp) ((%asin simp) x) m)
403 ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
404 ((%elliptic_f simp) ((%asin simp) x) m)))))))
405 grad)
406
407 ;; Let u = inverse_jacobi_cn(x). Then jacobi_cn(u) = x or
408 ;; sqrt(1-jacobi_sn(u)^2) = x. Or
409 ;;
410 ;; jacobi_sn(u) = sqrt(1-x^2)
411 ;;
412 ;; So u = inverse_jacobi_sn(sqrt(1-x^2),m) = inverse_jacob_cn(x,m)
413 ;;
414 (defprop %inverse_jacobi_cn
415 ((x m)
416 ;; Whittaker and Watson, 22.121
417 ;; inverse_jacobi_cn(u,m) = integrate(1/sqrt(1-t^2)/sqrt(1-m+m*t^2), t, u, 1)
418 ;; -> -1/sqrt(1-x^2)/sqrt(1-m+m*x^2)
419 ((mtimes simp) -1
420 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
421 ((rat simp) -1 2))
422 ((mexpt simp)
423 ((mplus simp) 1 ((mtimes simp) -1 m)
424 ((mtimes simp) m ((mexpt simp) x 2)))
425 ((rat simp) -1 2)))
426 ((mtimes simp) ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
427 ((mplus simp)
428 ((mtimes simp) -1
429 ((mexpt simp)
430 ((mplus simp) 1
431 ((mtimes simp) -1 m ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))))
432 ((rat simp) -1 2))
433 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2))) ((rat simp) 1 2))
434 ((mabs simp) x))
435 ((mtimes simp) ((mexpt simp) m -1)
436 ((mplus simp)
437 ((%elliptic_e simp)
438 ((%asin simp)
439 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2))) ((rat simp) 1 2)))
440 m)
441 ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
442 ((%elliptic_f simp)
443 ((%asin simp)
444 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2))) ((rat simp) 1 2)))
445 m)))))))
446 grad)
447
448 ;; Let u = inverse_jacobi_dn(x). Then
449 ;;
450 ;; jacobi_dn(u) = x or
451 ;;
452 ;; x^2 = jacobi_dn(u)^2 = 1 - m*jacobi_sn(u)^2
453 ;;
454 ;; so jacobi_sn(u) = sqrt(1-x^2)/sqrt(m)
455 ;;
456 ;; or u = inverse_jacobi_sn(sqrt(1-x^2)/sqrt(m))
457 (defprop %inverse_jacobi_dn
458 ((x m)
459 ;; Whittaker and Watson, 22.121
460 ;; inverse_jacobi_dn(u,m) = integrate(1/sqrt(1-t^2)/sqrt(t^2-(1-m)), t, u, 1)
461 ;; -> -1/sqrt(1-x^2)/sqrt(x^2+m-1)
462 ((mtimes simp)
463 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
464 ((rat simp) -1 2))
465 ((mexpt simp) ((mplus simp) -1 m ((mexpt simp) x 2)) ((rat simp) -1 2)))
466 ((mplus simp)
467 ((mtimes simp) ((rat simp) -1 2) ((mexpt simp) m ((rat simp) -3 2))
468 ((mexpt simp)
469 ((mplus simp) 1
470 ((mtimes simp) -1 ((mexpt simp) m -1)
471 ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))))
472 ((rat simp) -1 2))
473 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
474 ((rat simp) 1 2))
475 ((mexpt simp) ((mabs simp) x) -1))
476 ((mtimes simp) ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
477 ((mplus simp)
478 ((mtimes simp) -1 ((mexpt simp) m ((rat simp) -1 2))
479 ((mexpt simp)
480 ((mplus simp) 1
481 ((mtimes simp) -1 ((mexpt simp) m -1)
482 ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))))
483 ((rat simp) 1 2))
484 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 ((mexpt simp) x 2)))
485 ((rat simp) 1 2))
486 ((mexpt simp) ((mabs simp) x) -1))
487 ((mtimes simp) ((mexpt simp) m -1)
488 ((mplus simp)
489 ((%elliptic_e simp)
490 ((%asin simp)
491 ((mtimes simp) ((mexpt simp) m ((rat simp) -1 2))
492 ((mexpt simp) ((mplus simp) 1
493 ((mtimes simp) -1 ((mexpt simp) x 2)))
494 ((rat simp) 1 2))))
495 m)
496 ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
497 ((%elliptic_f simp)
498 ((%asin simp)
499 ((mtimes simp) ((mexpt simp) m ((rat simp) -1 2))
500 ((mexpt simp) ((mplus simp) 1
501 ((mtimes simp) -1 ((mexpt simp) x 2)))
502 ((rat simp) 1 2))))
503 m))))))))
504 grad)
505
506
507 ;; Possible forms of a complex number:
508 ;;
509 ;; 2.3
510 ;; $%i
511 ;; ((mplus simp) 2.3 ((mtimes simp) 2.3 $%i))
512 ;; ((mplus simp) 2.3 $%i))
513 ;; ((mtimes simp) 2.3 $%i)
514 ;;
515
516
517 ;; Is argument u a complex number with real and imagpart satisfying predicate ntypep?
518 (defun complex-number-p (u &optional (ntypep 'numberp))
519 (let ((R 0) (I 0))
520 (labels ((a1 (x) (cadr x))
521 (a2 (x) (caddr x))
522 (a3+ (x) (cdddr x))
523 (N (x) (funcall ntypep x)) ; N
524 (i (x) (and (eq x '$%i) (N 1))) ; %i
525 (N+i (x) (and (null (a3+ x)) ; mplus test is precondition
526 (N (setq R (a1 x)))
527 (or (and (i (a2 x)) (setq I 1) t)
528 (and (mtimesp (a2 x)) (N*i (a2 x))))))
529 (N*i (x) (and (null (a3+ x)) ; mtimes test is precondition
530 (N (setq I (a1 x)))
531 (eq (a2 x) '$%i))))
532 (declare (inline a1 a2 a3+ N i N+i N*i))
533 (cond ((N u) (values t u 0)) ;2.3
534 ((atom u) (if (i u) (values t 0 1))) ;%i
535 ((mplusp u) (if (N+i u) (values t R I))) ;N+%i, N+N*%i
536 ((mtimesp u) (if (N*i u) (values t R I))) ;N*%i
537 (t nil)))))
538
539 (defun complexify (x)
540 ;; Convert a Lisp number to a maxima number
541 (cond ((realp x) x)
542 ((complexp x) (add (realpart x) (mul '$%i (imagpart x))))
543 (t (merror (intl:gettext "COMPLEXIFY: argument must be a Lisp real or complex number.~%COMPLEXIFY: found: ~:M") x))))
544
545 (defun kc-arg (exp m)
546 ;; Replace elliptic_kc(m) in the expression with sym. Check to see
547 ;; if the resulting expression is linear in sym and the constant
548 ;; term is zero. If so, return the coefficient of sym, i.e, the
549 ;; coefficient of elliptic_kc(m).
550 (let* ((sym (gensym))
551 (arg (maxima-substitute sym `((%elliptic_kc) ,m) exp)))
552 (if (and (not (equalp arg exp))
553 (linearp arg sym)
554 (zerop1 (coefficient arg sym 0)))
555 (coefficient arg sym 1)
556 nil)))
557
558 (defun kc-arg2 (exp m)
559 ;; Replace elliptic_kc(m) in the expression with sym. Check to see
560 ;; if the resulting expression is linear in sym and the constant
561 ;; term is zero. If so, return the coefficient of sym, i.e, the
562 ;; coefficient of elliptic_kc(m), and the constant term. Otherwise,
563 ;; return NIL.
564 (let* ((sym (gensym))
565 (arg (maxima-substitute sym `((%elliptic_kc) ,m) exp)))
566 (if (and (not (equalp arg exp))
567 (linearp arg sym))
568 (list (coefficient arg sym 1)
569 (coefficient arg sym 0))
570 nil)))
571
572 ;; Tell maxima how to simplify the functions
573
574 (def-simplifier jacobi_sn (u m)
575 (let (coef args)
576 (cond
577 ((float-numerical-eval-p u m)
578 (to (bigfloat::sn (bigfloat:to ($float u)) (bigfloat:to ($float m)))))
579 ((setf args (complex-float-numerical-eval-p u m))
580 (destructuring-bind (u m)
581 args
582 (to (bigfloat::sn (bigfloat:to ($float u)) (bigfloat:to ($float m))))))
583 ((bigfloat-numerical-eval-p u m)
584 (to (bigfloat::sn (bigfloat:to ($bfloat u)) (bigfloat:to ($bfloat m)))))
585 ((setf args (complex-bigfloat-numerical-eval-p u m))
586 (destructuring-bind (u m)
587 args
588 (to (bigfloat::sn (bigfloat:to ($bfloat u)) (bigfloat:to ($bfloat m))))))
589 ((zerop1 u)
590 ;; A&S 16.5.1
591 0)
592 ((zerop1 m)
593 ;; A&S 16.6.1
594 (ftake '%sin u))
595 ((onep1 m)
596 ;; A&S 16.6.1
597 (ftake '%tanh u))
598 ((and $trigsign (mminusp* u))
599 (neg (ftake* '%jacobi_sn (neg u) m)))
600 ((and $triginverses
601 (listp u)
602 (member (caar u) '(%inverse_jacobi_sn
603 %inverse_jacobi_ns
604 %inverse_jacobi_cn
605 %inverse_jacobi_nc
606 %inverse_jacobi_dn
607 %inverse_jacobi_nd
608 %inverse_jacobi_sc
609 %inverse_jacobi_cs
610 %inverse_jacobi_sd
611 %inverse_jacobi_ds
612 %inverse_jacobi_cd
613 %inverse_jacobi_dc))
614 (alike1 (third u) m))
615 (let ((inv-arg (second u)))
616 (ecase (caar u)
617 (%inverse_jacobi_sn
618 ;; jacobi_sn(inverse_jacobi_sn(u,m), m) = u
619 inv-arg)
620 (%inverse_jacobi_ns
621 ;; inverse_jacobi_ns(u,m) = inverse_jacobi_sn(1/u,m)
622 (div 1 inv-arg))
623 (%inverse_jacobi_cn
624 ;; sn(x)^2 + cn(x)^2 = 1 so sn(x) = sqrt(1-cn(x)^2)
625 (power (sub 1 (mul inv-arg inv-arg)) 1//2))
626 (%inverse_jacobi_nc
627 ;; inverse_jacobi_nc(u) = inverse_jacobi_cn(1/u)
628 (ftake '%jacobi_sn (ftake '%inverse_jacobi_cn (div 1 inv-arg) m)
629 m))
630 (%inverse_jacobi_dn
631 ;; dn(x)^2 + m*sn(x)^2 = 1 so
632 ;; sn(x) = 1/sqrt(m)*sqrt(1-dn(x)^2)
633 (mul (div 1 (power m 1//2))
634 (power (sub 1 (mul inv-arg inv-arg)) 1//2)))
635 (%inverse_jacobi_nd
636 ;; inverse_jacobi_nd(u) = inverse_jacobi_dn(1/u)
637 (ftake '%jacobi_sn (ftake '%inverse_jacobi_dn (div 1 inv-arg) m)
638 m))
639 (%inverse_jacobi_sc
640 ;; See below for inverse_jacobi_sc.
641 (div inv-arg (power (add 1 (mul inv-arg inv-arg)) 1//2)))
642 (%inverse_jacobi_cs
643 ;; inverse_jacobi_cs(u) = inverse_jacobi_sc(1/u)
644 (ftake '%jacobi_sn (ftake '%inverse_jacobi_sc (div 1 inv-arg) m)
645 m))
646 (%inverse_jacobi_sd
647 ;; See below for inverse_jacobi_sd
648 (div inv-arg (power (add 1 (mul m (mul inv-arg inv-arg))) 1//2)))
649 (%inverse_jacobi_ds
650 ;; inverse_jacobi_ds(u) = inverse_jacobi_sd(1/u)
651 (ftake '%jacobi_sn (ftake '%inverse_jacobi_sd (div 1 inv-arg) m)
652 m))
653 (%inverse_jacobi_cd
654 ;; See below
655 (div (power (sub 1 (mul inv-arg inv-arg)) 1//2)
656 (power (sub 1 (mul m (mul inv-arg inv-arg))) 1//2)))
657 (%inverse_jacobi_dc
658 (ftake '%jacobi_sn (ftake '%inverse_jacobi_cd (div 1 inv-arg) m) m)))))
659 ;; A&S 16.20.1 (Jacobi's Imaginary transformation)
660 ((and $%iargs (multiplep u '$%i))
661 (mul '$%i
662 (ftake* '%jacobi_sc (coeff u '$%i 1) (add 1 (neg m)))))
663 ((setq coef (kc-arg2 u m))
664 ;; sn(m*K+u)
665 ;;
666 ;; A&S 16.8.1
667 (destructuring-bind (lin const)
668 coef
669 (cond ((integerp lin)
670 (ecase (mod lin 4)
671 (0
672 ;; sn(4*m*K + u) = sn(u), sn(0) = 0
673 (if (zerop1 const)
674 0
675 (ftake '%jacobi_sn const m)))
676 (1
677 ;; sn(4*m*K + K + u) = sn(K+u) = cd(u)
678 ;; sn(K) = 1
679 (if (zerop1 const)
680 1
681 (ftake '%jacobi_cd const m)))
682 (2
683 ;; sn(4*m*K+2*K + u) = sn(2*K+u) = -sn(u)
684 ;; sn(2*K) = 0
685 (if (zerop1 const)
686 0
687 (neg (ftake '%jacobi_sn const m))))
688 (3
689 ;; sn(4*m*K+3*K+u) = sn(2*K + K + u) = -sn(K+u) = -cd(u)
690 ;; sn(3*K) = -1
691 (if (zerop1 const)
692 -1
693 (neg (ftake '%jacobi_cd const m))))))
694 ((and (alike1 lin 1//2)
695 (zerop1 const))
696 ;; A&S 16.5.2
697 ;;
698 ;; sn(1/2*K) = 1/sqrt(1+sqrt(1-m))
699 (div 1
700 (power (add 1 (power (sub 1 m) 1//2))
701 1//2)))
702 ((and (alike1 lin 3//2)
703 (zerop1 const))
704 ;; A&S 16.5.2
705 ;;
706 ;; sn(1/2*K + K) = cd(1/2*K,m)
707 (ftake '%jacobi_cd (mul 1//2
708 (ftake '%elliptic_kc m))
709 m))
710 (t
711 (give-up)))))
712 (t
713 ;; Nothing to do
714 (give-up)))))
715
716 (def-simplifier jacobi_cn (u m)
717 (let (coef args)
718 (cond
719 ((float-numerical-eval-p u m)
720 (to (bigfloat::cn (bigfloat:to ($float u)) (bigfloat:to ($float m)))))
721 ((setf args (complex-float-numerical-eval-p u m))
722 (destructuring-bind (u m)
723 args
724 (to (bigfloat::cn (bigfloat:to ($float u)) (bigfloat:to ($float m))))))
725 ((bigfloat-numerical-eval-p u m)
726 (to (bigfloat::cn (bigfloat:to ($bfloat u)) (bigfloat:to ($bfloat m)))))
727 ((setf args (complex-bigfloat-numerical-eval-p u m))
728 (destructuring-bind (u m)
729 args
730 (to (bigfloat::cn (bigfloat:to ($bfloat u)) (bigfloat:to ($bfloat m))))))
731 ((zerop1 u)
732 ;; A&S 16.5.1
733 1)
734 ((zerop1 m)
735 ;; A&S 16.6.2
736 (ftake '%cos u))
737 ((onep1 m)
738 ;; A&S 16.6.2
739 (ftake '%sech u))
740 ((and $trigsign (mminusp* u))
741 (ftake* '%jacobi_cn (neg u) m))
742 ((and $triginverses
743 (listp u)
744 (member (caar u) '(%inverse_jacobi_sn
745 %inverse_jacobi_ns
746 %inverse_jacobi_cn
747 %inverse_jacobi_nc
748 %inverse_jacobi_dn
749 %inverse_jacobi_nd
750 %inverse_jacobi_sc
751 %inverse_jacobi_cs
752 %inverse_jacobi_sd
753 %inverse_jacobi_ds
754 %inverse_jacobi_cd
755 %inverse_jacobi_dc))
756 (alike1 (third u) m))
757 (cond ((eq (caar u) '%inverse_jacobi_cn)
758 (second u))
759 (t
760 ;; I'm lazy. Use cn(x) = sqrt(1-sn(x)^2). Hope
761 ;; this is right.
762 (power (sub 1 (power (ftake '%jacobi_sn u (third u)) 2))
763 1//2))))
764 ;; A&S 16.20.2 (Jacobi's Imaginary transformation)
765 ((and $%iargs (multiplep u '$%i))
766 (ftake* '%jacobi_nc (coeff u '$%i 1) (add 1 (neg m))))
767 ((setq coef (kc-arg2 u m))
768 ;; cn(m*K+u)
769 ;;
770 ;; A&S 16.8.2
771 (destructuring-bind (lin const)
772 coef
773 (cond ((integerp lin)
774 (ecase (mod lin 4)
775 (0
776 ;; cn(4*m*K + u) = cn(u),
777 ;; cn(0) = 1
778 (if (zerop1 const)
779 1
780 (ftake '%jacobi_cn const m)))
781 (1
782 ;; cn(4*m*K + K + u) = cn(K+u) = -sqrt(m1)*sd(u)
783 ;; cn(K) = 0
784 (if (zerop1 const)
785 0
786 (neg (mul (power (sub 1 m) 1//2)
787 (ftake '%jacobi_sd const m)))))
788 (2
789 ;; cn(4*m*K + 2*K + u) = cn(2*K+u) = -cn(u)
790 ;; cn(2*K) = -1
791 (if (zerop1 const)
792 -1
793 (neg (ftake '%jacobi_cn const m))))
794 (3
795 ;; cn(4*m*K + 3*K + u) = cn(2*K + K + u) =
796 ;; -cn(K+u) = sqrt(m1)*sd(u)
797 ;;
798 ;; cn(3*K) = 0
799 (if (zerop1 const)
800 0
801 (mul (power (sub 1 m) 1//2)
802 (ftake '%jacobi_sd const m))))))
803 ((and (alike1 lin 1//2)
804 (zerop1 const))
805 ;; A&S 16.5.2
806 ;; cn(1/2*K) = (1-m)^(1/4)/sqrt(1+sqrt(1-m))
807 (mul (power (sub 1 m) (div 1 4))
808 (power (add 1
809 (power (sub 1 m)
810 1//2))
811 1//2)))
812 (t
813 (give-up)))))
814 (t
815 (give-up)))))
816
817 (def-simplifier jacobi_dn (u m)
818 (let (coef args)
819 (cond
820 ((float-numerical-eval-p u m)
821 (to (bigfloat::dn (bigfloat:to ($float u)) (bigfloat:to ($float m)))))
822 ((setf args (complex-float-numerical-eval-p u m))
823 (destructuring-bind (u m)
824 args
825 (to (bigfloat::dn (bigfloat:to ($float u)) (bigfloat:to ($float m))))))
826 ((bigfloat-numerical-eval-p u m)
827 (to (bigfloat::dn (bigfloat:to ($bfloat u)) (bigfloat:to ($bfloat m)))))
828 ((setf args (complex-bigfloat-numerical-eval-p u m))
829 (destructuring-bind (u m)
830 args
831 (to (bigfloat::dn (bigfloat:to ($bfloat u)) (bigfloat:to ($bfloat m))))))
832 ((zerop1 u)
833 ;; A&S 16.5.1
834 1)
835 ((zerop1 m)
836 ;; A&S 16.6.3
837 1)
838 ((onep1 m)
839 ;; A&S 16.6.3
840 (ftake '%sech u))
841 ((and $trigsign (mminusp* u))
842 (ftake* '%jacobi_dn (neg u) m))
843 ((and $triginverses
844 (listp u)
845 (member (caar u) '(%inverse_jacobi_sn
846 %inverse_jacobi_ns
847 %inverse_jacobi_cn
848 %inverse_jacobi_nc
849 %inverse_jacobi_dn
850 %inverse_jacobi_nd
851 %inverse_jacobi_sc
852 %inverse_jacobi_cs
853 %inverse_jacobi_sd
854 %inverse_jacobi_ds
855 %inverse_jacobi_cd
856 %inverse_jacobi_dc))
857 (alike1 (third u) m))
858 (cond ((eq (caar u) '%inverse_jacobi_dn)
859 ;; jacobi_dn(inverse_jacobi_dn(u,m), m) = u
860 (second u))
861 (t
862 ;; Express in terms of sn:
863 ;; dn(x) = sqrt(1-m*sn(x)^2)
864 (power (sub 1 (mul m
865 (power (ftake '%jacobi_sn u m) 2)))
866 1//2))))
867 ((zerop1 ($ratsimp (sub u (power (sub 1 m) 1//2))))
868 ;; A&S 16.5.3
869 ;; dn(sqrt(1-m),m) = K(m)
870 (ftake '%elliptic_kc m))
871 ;; A&S 16.20.2 (Jacobi's Imaginary transformation)
872 ((and $%iargs (multiplep u '$%i))
873 (ftake* '%jacobi_dc (coeff u '$%i 1)
874 (add 1 (neg m))))
875 ((setq coef (kc-arg2 u m))
876 ;; A&S 16.8.3
877 ;;
878 ;; dn(m*K+u) has period 2K
879 ;;
880 (destructuring-bind (lin const)
881 coef
882 (cond ((integerp lin)
883 (ecase (mod lin 2)
884 (0
885 ;; dn(2*m*K + u) = dn(u)
886 ;; dn(0) = 1
887 (if (zerop1 const)
888 1
889 ;; dn(4*m*K+2*K + u) = dn(2*K+u) = dn(u)
890 (ftake '%jacobi_dn const m)))
891 (1
892 ;; dn(2*m*K + K + u) = dn(K + u) = sqrt(1-m)*nd(u)
893 ;; dn(K) = sqrt(1-m)
894 (if (zerop1 const)
895 (power (sub 1 m) 1//2)
896 (mul (power (sub 1 m) 1//2)
897 (ftake '%jacobi_nd const m))))))
898 ((and (alike1 lin 1//2)
899 (zerop1 const))
900 ;; A&S 16.5.2
901 ;; dn(1/2*K) = (1-m)^(1/4)
902 (power (sub 1 m)
903 (div 1 4)))
904 (t
905 (give-up)))))
906 (t (give-up)))))
907
908 ;; Should we simplify the inverse elliptic functions into the
909 ;; appropriate incomplete elliptic integral? I think we should leave
910 ;; it, but perhaps allow some way to do that transformation if
911 ;; desired.
912
913 (def-simplifier inverse_jacobi_sn (u m)
914 (let (args)
915 ;; To numerically evaluate inverse_jacobi_sn (asn), use
916 ;;
917 ;; asn(x,m) = F(asin(x),m)
918 ;;
919 ;; But F(phi,m) = sin(phi)*rf(cos(phi)^2, 1-m*sin(phi)^2,1). Thus
920 ;;
921 ;; asn(x,m) = F(asin(x),m)
922 ;; = x*rf(1-x^2,1-m*x^2,1)
923 ;;
924 ;; I (rtoy) am not 100% about the first identity above for all
925 ;; complex values of x and m, but tests seem to indicate that it
926 ;; produces the correct value as verified by verifying
927 ;; jacobi_sn(inverse_jacobi_sn(x,m),m) = x.
928 (cond ((float-numerical-eval-p u m)
929 (let ((uu (bigfloat:to ($float u)))
930 (mm (bigfloat:to ($float m))))
931 (complexify
932 (* uu
933 (bigfloat::bf-rf (bigfloat:to (- 1 (* uu uu)))
934 (bigfloat:to (- 1 (* mm uu uu)))
935 1)))))
936 ((setf args (complex-float-numerical-eval-p u m))
937 (destructuring-bind (u m)
938 args
939 (let ((uu (bigfloat:to ($float u)))
940 (mm (bigfloat:to ($float m))))
941 (complexify (* uu (bigfloat::bf-rf (- 1 (* uu uu))
942 (- 1 (* mm uu uu))
943 1))))))
944 ((bigfloat-numerical-eval-p u m)
945 (let ((uu (bigfloat:to u))
946 (mm (bigfloat:to m)))
947 (to (bigfloat:* uu
948 (bigfloat::bf-rf (bigfloat:- 1 (bigfloat:* uu uu))
949 (bigfloat:- 1 (bigfloat:* mm uu uu))
950 1)))))
951 ((setf args (complex-bigfloat-numerical-eval-p u m))
952 (destructuring-bind (u m)
953 args
954 (let ((uu (bigfloat:to u))
955 (mm (bigfloat:to m)))
956 (to (bigfloat:* uu
957 (bigfloat::bf-rf (bigfloat:- 1 (bigfloat:* uu uu))
958 (bigfloat:- 1 (bigfloat:* mm uu uu))
959 1))))))
960 ((zerop1 u)
961 ;; asn(0,m) = 0
962 0)
963 ((onep1 u)
964 ;; asn(1,m) = elliptic_kc(m)
965 (ftake '%elliptic_kc m))
966 ((and (numberp u) (onep1 (- u)))
967 ;; asn(-1,m) = -elliptic_kc(m)
968 (mul -1 (ftake '%elliptic_kc m)))
969 ((zerop1 m)
970 ;; asn(x,0) = F(asin(x),0) = asin(x)
971 (ftake '%asin u))
972 ((onep1 m)
973 ;; asn(x,1) = F(asin(x),1) = log(tan(pi/4+asin(x)/2))
974 (ftake '%elliptic_f (ftake '%asin u) 1))
975 ((and (eq $triginverses '$all)
976 (listp u)
977 (eq (caar u) '%jacobi_sn)
978 (alike1 (third u) m))
979 ;; inverse_jacobi_sn(sn(u)) = u
980 (second u))
981 (t
982 ;; Nothing to do
983 (give-up)))))
984
985 (def-simplifier inverse_jacobi_cn (u m)
986 (let (args)
987 (cond ((float-numerical-eval-p u m)
988 ;; Numerically evaluate acn
989 ;;
990 ;; acn(x,m) = F(acos(x),m)
991 (to (elliptic-f (cl:acos ($float u)) ($float m))))
992 ((setf args (complex-float-numerical-eval-p u m))
993 (destructuring-bind (u m)
994 args
995 (to (elliptic-f (cl:acos (bigfloat:to ($float u)))
996 (bigfloat:to ($float m))))))
997 ((bigfloat-numerical-eval-p u m)
998 (to (bigfloat::bf-elliptic-f (bigfloat:acos (bigfloat:to u))
999 (bigfloat:to m))))
1000 ((setf args (complex-bigfloat-numerical-eval-p u m))
1001 (destructuring-bind (u m)
1002 args
1003 (to (bigfloat::bf-elliptic-f (bigfloat:acos (bigfloat:to u))
1004 (bigfloat:to m)))))
1005 ((zerop1 m)
1006 ;; asn(x,0) = F(acos(x),0) = acos(x)
1007 (ftake '%elliptic_f (ftake '%acos u) 0))
1008 ((onep1 m)
1009 ;; asn(x,1) = F(asin(x),1) = log(tan(pi/2+asin(x)/2))
1010 (ftake '%elliptic_f (ftake '%acos u) 1))
1011 ((zerop1 u)
1012 (ftake '%elliptic_kc m))
1013 ((onep1 u)
1014 0)
1015 ((and (eq $triginverses '$all)
1016 (listp u)
1017 (eq (caar u) '%jacobi_cn)
1018 (alike1 (third u) m))
1019 ;; inverse_jacobi_cn(cn(u)) = u
1020 (second u))
1021 (t
1022 ;; Nothing to do
1023 (give-up)))))
1024
1025 (def-simplifier inverse_jacobi_dn (u m)
1026 (let (args)
1027 (cond ((float-numerical-eval-p u m)
1028 (to (bigfloat::bf-inverse-jacobi-dn (bigfloat:to (float u))
1029 (bigfloat:to (float m)))))
1030 ((setf args (complex-float-numerical-eval-p u m))
1031 (destructuring-bind (u m)
1032 args
1033 (let ((uu (bigfloat:to ($float u)))
1034 (mm (bigfloat:to ($float m))))
1035 (to (bigfloat::bf-inverse-jacobi-dn uu mm)))))
1036 ((bigfloat-numerical-eval-p u m)
1037 (let ((uu (bigfloat:to u))
1038 (mm (bigfloat:to m)))
1039 (to (bigfloat::bf-inverse-jacobi-dn uu mm))))
1040 ((setf args (complex-bigfloat-numerical-eval-p u m))
1041 (destructuring-bind (u m)
1042 args
1043 (to (bigfloat::bf-inverse-jacobi-dn (bigfloat:to u) (bigfloat:to m)))))
1044 ((onep1 m)
1045 ;; x = dn(u,1) = sech(u). so u = asech(x)
1046 (ftake '%asech u))
1047 ((onep1 u)
1048 ;; jacobi_dn(0,m) = 1
1049 0)
1050 ((zerop1 ($ratsimp (sub u (power (sub 1 m) 1//2))))
1051 ;; jacobi_dn(K(m),m) = sqrt(1-m) so
1052 ;; inverse_jacobi_dn(sqrt(1-m),m) = K(m)
1053 (ftake '%elliptic_kc m))
1054 ((and (eq $triginverses '$all)
1055 (listp u)
1056 (eq (caar u) '%jacobi_dn)
1057 (alike1 (third u) m))
1058 ;; inverse_jacobi_dn(dn(u)) = u
1059 (second u))
1060 (t
1061 ;; Nothing to do
1062 (give-up)))))
1063
1064 ;;;; Elliptic integrals
1065
1066 (let ((errtol (expt (* 4 +flonum-epsilon+) 1/6))
1067 (c1 (float 1/24))
1068 (c2 (float 3/44))
1069 (c3 (float 1/14)))
1070 (declare (type flonum errtol c1 c2 c3))
1071 (defun crf (x y z)
1072 "Compute Carlson's incomplete or complete elliptic integral of the
1073 first kind:
1074
1075 INF
1076 /
1077 [ 1
1078 RF(x, y, z) = I ----------------------------------- dt
1079 ] SQRT(x + t) SQRT(y + t) SQRT(z + t)
1080 /
1081 0
1082
1083 x, y, and z may be complex.
1084 "
1085 (declare (number x y z))
1086 (let ((x (coerce x '(complex flonum)))
1087 (y (coerce y '(complex flonum)))
1088 (z (coerce z '(complex flonum))))
1089 (declare (type (complex flonum) x y z))
1090 (loop
1091 (let* ((mu (/ (+ x y z) 3))
1092 (x-dev (- 2 (/ (+ mu x) mu)))
1093 (y-dev (- 2 (/ (+ mu y) mu)))
1094 (z-dev (- 2 (/ (+ mu z) mu))))
1095 (when (< (max (abs x-dev) (abs y-dev) (abs z-dev)) errtol)
1096 (let ((e2 (- (* x-dev y-dev) (* z-dev z-dev)))
1097 (e3 (* x-dev y-dev z-dev)))
1098 (return (/ (+ 1
1099 (* e2 (- (* c1 e2)
1100 1/10
1101 (* c2 e3)))
1102 (* c3 e3))
1103 (sqrt mu)))))
1104 (let* ((x-root (sqrt x))
1105 (y-root (sqrt y))
1106 (z-root (sqrt z))
1107 (lam (+ (* x-root (+ y-root z-root)) (* y-root z-root))))
1108 (setf x (* (+ x lam) 1/4))
1109 (setf y (* (+ y lam) 1/4))
1110 (setf z (* (+ z lam) 1/4))))))))
1111
1112 ;; Elliptic integral of the first kind (Legendre's form):
1113 ;;
1114 ;;
1115 ;; phi
1116 ;; /
1117 ;; [ 1
1118 ;; I ------------------- ds
1119 ;; ] 2
1120 ;; / SQRT(1 - m SIN (s))
1121 ;; 0
1122
1123 (defun elliptic-f (phi-arg m-arg)
1124 (flet ((base (phi-arg m-arg)
1125 (cond ((and (realp m-arg) (realp phi-arg))
1126 (let ((phi (float phi-arg))
1127 (m (float m-arg)))
1128 (cond ((> m 1)
1129 ;; A&S 17.4.15
1130 ;;
1131 ;; F(phi|m) = 1/sqrt(m)*F(theta|1/m)
1132 ;;
1133 ;; with sin(theta) = sqrt(m)*sin(phi)
1134 (/ (elliptic-f (cl:asin (* (sqrt m) (sin phi))) (/ m))
1135 (sqrt m)))
1136 ((< m 0)
1137 ;; A&S 17.4.17
1138 (let* ((m (- m))
1139 (m+1 (+ 1 m))
1140 (root (sqrt m+1))
1141 (m/m+1 (/ m m+1)))
1142 (- (/ (elliptic-f (float (/ pi 2)) m/m+1)
1143 root)
1144 (/ (elliptic-f (- (float (/ pi 2)) phi) m/m+1)
1145 root))))
1146 ((= m 0)
1147 ;; A&S 17.4.19
1148 phi)
1149 ((= m 1)
1150 ;; A&S 17.4.21
1151 ;;
1152 1 ;; F(phi,1) = log(sec(phi)+tan(phi))
1153 ;; = log(tan(pi/4+pi/2))
1154 (log (cl:tan (+ (/ phi 2) (float (/ pi 4))))))
1155 ((minusp phi)
1156 (- (elliptic-f (- phi) m)))
1157 ((> phi pi)
1158 ;; A&S 17.4.3
1159 (multiple-value-bind (s phi-rem)
1160 (truncate phi (float pi))
1161 (+ (* 2 s (elliptic-k m))
1162 (elliptic-f phi-rem m))))
1163 ((<= phi (/ pi 2))
1164 (let ((sin-phi (sin phi))
1165 (cos-phi (cos phi))
1166 (k (sqrt m)))
1167 (* sin-phi
1168 (bigfloat::bf-rf (* cos-phi cos-phi)
1169 (* (- 1 (* k sin-phi))
1170 (+ 1 (* k sin-phi)))
1171 1.0))))
1172 ((< phi pi)
1173 (+ (* 2 (elliptic-k m))
1174 (elliptic-f (- phi (float pi)) m)))
1175 (t
1176 (error "Shouldn't happen! Unhandled case in elliptic-f: ~S ~S~%"
1177 phi-arg m-arg)))))
1178 (t
1179 (let ((phi (coerce phi-arg '(complex flonum)))
1180 (m (coerce m-arg '(complex flonum))))
1181 (let ((sin-phi (sin phi))
1182 (cos-phi (cos phi))
1183 (k (sqrt m)))
1184 (* sin-phi
1185 (crf (* cos-phi cos-phi)
1186 (* (- 1 (* k sin-phi))
1187 (+ 1 (* k sin-phi)))
1188 1.0))))))))
1189 ;; Elliptic F is quasi-periodic wrt to z:
1190 ;;
1191 ;; F(z|m) = F(z - pi*round(Re(z)/pi)|m) + 2*round(Re(z)/pi)*K(m)
1192 (let ((period (round (realpart phi-arg) pi)))
1193 (+ (base (- phi-arg (* pi period)) m-arg)
1194 (if (zerop period)
1195 0
1196 (* 2 period
1197 (bigfloat:to (elliptic-k m-arg))))))))
1198
1199 ;; Complete elliptic integral of the first kind
1200 (defun elliptic-k (m)
1201 (cond ((realp m)
1202 (cond ((< m 0)
1203 ;; A&S 17.4.17
1204 (let* ((m (- m))
1205 (m+1 (+ 1 m))
1206 (root (sqrt m+1))
1207 (m/m+1 (/ m m+1)))
1208 (- (/ (elliptic-k m/m+1)
1209 root)
1210 (/ (elliptic-f 0.0 m/m+1)
1211 root))))
1212 ((= m 0)
1213 ;; A&S 17.4.19
1214 (float (/ pi 2)))
1215 ((= m 1)
1216 (maxima::merror
1217 (intl:gettext "elliptic_kc: elliptic_kc(1) is undefined.")))
1218 (t
1219 (bigfloat::bf-rf 0.0 (- 1 m)
1220 1.0))))
1221 (t
1222 (bigfloat::bf-rf 0.0 (- 1 m)
1223 1.0))))
1224
1225 ;; Elliptic integral of the second kind (Legendre's form):
1226 ;;
1227 ;;
1228 ;; phi
1229 ;; /
1230 ;; [ 2
1231 ;; I SQRT(1 - m SIN (s)) ds
1232 ;; ]
1233 ;; /
1234 ;; 0
1235
1236 (defun elliptic-e (phi m)
1237 (declare (type flonum phi m))
1238 (flet ((base (phi m)
1239 (cond ((= m 0)
1240 ;; A&S 17.4.23
1241 phi)
1242 ((= m 1)
1243 ;; A&S 17.4.25
1244 (sin phi))
1245 (t
1246 (let* ((sin-phi (sin phi))
1247 (cos-phi (cos phi))
1248 (k (sqrt m))
1249 (y (* (- 1 (* k sin-phi))
1250 (+ 1 (* k sin-phi)))))
1251 (to (- (* sin-phi
1252 (bigfloat::bf-rf (* cos-phi cos-phi) y 1.0))
1253 (* (/ m 3)
1254 (expt sin-phi 3)
1255 (bigfloat::bf-rd (* cos-phi cos-phi) y 1.0)))))))))
1256 ;; Elliptic E is quasi-periodic wrt to phi:
1257 ;;
1258 ;; E(z|m) = E(z - %pi*round(Re(z)/%pi)|m) + 2*round(Re(z)/%pi)*E(m)
1259 (let ((period (round (realpart phi) pi)))
1260 (+ (base (- phi (* pi period)) m)
1261 (* 2 period (elliptic-ec m))))))
1262
1263 ;; Complete version
1264 (defun elliptic-ec (m)
1265 (declare (type flonum m))
1266 (cond ((= m 0)
1267 ;; A&S 17.4.23
1268 (float (/ pi 2)))
1269 ((= m 1)
1270 ;; A&S 17.4.25
1271 1.0)
1272 (t
1273 (let* ((y (- 1 m)))
1274 (to (- (bigfloat::bf-rf 0.0 y 1.0)
1275 (* (/ m 3)
1276 (bigfloat::bf-rd 0.0 y 1.0))))))))
1277
1278
1279 ;; Define the elliptic integrals for maxima
1280 ;;
1281 ;; We use the definitions given in A&S 17.2.6 and 17.2.8. In particular:
1282 ;;
1283 ;; phi
1284 ;; /
1285 ;; [ 1
1286 ;; F(phi|m) = I ------------------- ds
1287 ;; ] 2
1288 ;; / SQRT(1 - m SIN (s))
1289 ;; 0
1290 ;;
1291 ;; and
1292 ;;
1293 ;; phi
1294 ;; /
1295 ;; [ 2
1296 ;; E(phi|m) = I SQRT(1 - m SIN (s)) ds
1297 ;; ]
1298 ;; /
1299 ;; 0
1300 ;;
1301 ;; That is, we do not use the modular angle, alpha, as the second arg;
1302 ;; the parameter m = sin(alpha)^2 is used.
1303 ;;
1304
1305
1306 ;; The derivative of F(phi|m) wrt to phi is easy. The derivative wrt
1307 ;; to m is harder. Here is a derivation. Hope I got it right.
1308 ;;
1309 ;; diff(integrate(1/sqrt(1-m*sin(x)^2),x,0,phi), m);
1310 ;;
1311 ;; PHI
1312 ;; / 2
1313 ;; [ SIN (x)
1314 ;; I ------------------ dx
1315 ;; ] 2 3/2
1316 ;; / (1 - m SIN (x))
1317 ;; 0
1318 ;; --------------------------
1319 ;; 2
1320 ;;
1321 ;;
1322 ;; Now use the following relationship that is easily verified:
1323 ;;
1324 ;; 2 2
1325 ;; (1 - m) SIN (x) COS (x) COS(x) SIN(x)
1326 ;; ------------------- = ------------------- - DIFF(-------------------, x)
1327 ;; 2 2 2
1328 ;; SQRT(1 - m SIN (x)) SQRT(1 - m SIN (x)) SQRT(1 - m SIN (x))
1329 ;;
1330 ;;
1331 ;; Now integrate this to get:
1332 ;;
1333 ;;
1334 ;; PHI
1335 ;; / 2
1336 ;; [ SIN (x)
1337 ;; (1 - m) I ------------------- dx =
1338 ;; ] 2
1339 ;; / SQRT(1 - m SIN (x))
1340 ;; 0
1341
1342 ;;
1343 ;; PHI
1344 ;; / 2
1345 ;; [ COS (x)
1346 ;; + I ------------------- dx
1347 ;; ] 2
1348 ;; / SQRT(1 - m SIN (x))
1349 ;; 0
1350 ;; COS(PHI) SIN(PHI)
1351 ;; - ---------------------
1352 ;; 2
1353 ;; SQRT(1 - m SIN (PHI))
1354 ;;
1355 ;; Use the fact that cos(x)^2 = 1 - sin(x)^2 to show that this
1356 ;; integral on the RHS is:
1357 ;;
1358 ;;
1359 ;; (1 - m) elliptic_F(PHI, m) + elliptic_E(PHI, m)
1360 ;; -------------------------------------------
1361 ;; m
1362 ;; So, finally, we have
1363 ;;
1364 ;;
1365 ;;
1366 ;; d
1367 ;; 2 -- (elliptic_F(PHI, m)) =
1368 ;; dm
1369 ;;
1370 ;; elliptic_E(PHI, m) - (1 - m) elliptic_F(PHI, m) COS(PHI) SIN(PHI)
1371 ;; ---------------------------------------------- - ---------------------
1372 ;; m 2
1373 ;; SQRT(1 - m SIN (PHI))
1374 ;; ----------------------------------------------------------------------
1375 ;; 1 - m
1376
1377 (defprop %elliptic_f
1378 ((phi m)
1379 ;; diff wrt phi
1380 ;; 1/sqrt(1-m*sin(phi)^2)
1381 ((mexpt simp)
1382 ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) ((%sin simp) phi) 2)))
1383 ((rat simp) -1 2))
1384 ;; diff wrt m
1385 ((mtimes simp) ((rat simp) 1 2)
1386 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
1387 ((mplus simp)
1388 ((mtimes simp) ((mexpt simp) m -1)
1389 ((mplus simp) ((%elliptic_e simp) phi m)
1390 ((mtimes simp) -1 ((mplus simp) 1 ((mtimes simp) -1 m))
1391 ((%elliptic_f simp) phi m))))
1392 ((mtimes simp) -1 ((%cos simp) phi) ((%sin simp) phi)
1393 ((mexpt simp)
1394 ((mplus simp) 1
1395 ((mtimes simp) -1 m ((mexpt simp) ((%sin simp) phi) 2)))
1396 ((rat simp) -1 2))))))
1397 grad)
1398
1399 ;;
1400 ;; The derivative of E(phi|m) wrt to m is much simpler to derive than F(phi|m).
1401 ;;
1402 ;; Take the derivative of the definition to get
1403 ;;
1404 ;; PHI
1405 ;; / 2
1406 ;; [ SIN (x)
1407 ;; I ------------------- dx
1408 ;; ] 2
1409 ;; / SQRT(1 - m SIN (x))
1410 ;; 0
1411 ;; - ---------------------------
1412 ;; 2
1413 ;;
1414 ;; It is easy to see that
1415 ;;
1416 ;; PHI
1417 ;; / 2
1418 ;; [ SIN (x)
1419 ;; elliptic_F(PHI, m) - m I ------------------- dx = elliptic_E(PHI, m)
1420 ;; ] 2
1421 ;; / SQRT(1 - m SIN (x))
1422 ;; 0
1423 ;;
1424 ;; So we finally have
1425 ;;
1426 ;; d elliptic_E(PHI, m) - elliptic_F(PHI, m)
1427 ;; -- (elliptic_E(PHI, m)) = ---------------------------------------
1428 ;; dm 2 m
1429
1430 (defprop %elliptic_e
1431 ((phi m)
1432 ;; sqrt(1-m*sin(phi)^2)
1433 ((mexpt simp)
1434 ((mplus simp) 1 ((mtimes simp) -1 m ((mexpt simp) ((%sin simp) phi) 2)))
1435 ((rat simp) 1 2))
1436 ;; diff wrt m
1437 ((mtimes simp) ((rat simp) 1 2) ((mexpt simp) m -1)
1438 ((mplus simp) ((%elliptic_e simp) phi m)
1439 ((mtimes simp) -1 ((%elliptic_f simp) phi m)))))
1440 grad)
1441
1442 (def-simplifier elliptic_f (phi m)
1443 (let (args)
1444 (cond ((float-numerical-eval-p phi m)
1445 ;; Numerically evaluate it
1446 (to (elliptic-f ($float phi) ($float m))))
1447 ((setf args (complex-float-numerical-eval-p phi m))
1448 (destructuring-bind (phi m)
1449 args
1450 (to (elliptic-f (bigfloat:to ($float phi))
1451 (bigfloat:to ($float m))))))
1452 ((bigfloat-numerical-eval-p phi m)
1453 (to (bigfloat::bf-elliptic-f (bigfloat:to ($bfloat phi))
1454 (bigfloat:to ($bfloat m)))))
1455 ((setf args (complex-bigfloat-numerical-eval-p phi m))
1456 (destructuring-bind (phi m)
1457 args
1458 (to (bigfloat::bf-elliptic-f (bigfloat:to ($bfloat phi))
1459 (bigfloat:to ($bfloat m))))))
1460 ((zerop1 phi)
1461 0)
1462 ((zerop1 m)
1463 ;; A&S 17.4.19
1464 phi)
1465 ((onep1 m)
1466 ;; A&S 17.4.21. Let's pick the log tan form. But this
1467 ;; isn't right if we know that abs(phi) > %pi/2, where
1468 ;; elliptic_f is undefined (or infinity).
1469 (cond ((not (eq '$pos (csign (sub ($abs phi) (div '$%pi 2)))))
1470 (ftake '%log
1471 (ftake '%tan
1472 (add (mul '$%pi (div 1 4))
1473 (mul 1//2 phi)))))
1474 (t
1475 (merror (intl:gettext "elliptic_f: elliptic_f(~:M, ~:M) is undefined.")
1476 phi m))))
1477 ((alike1 phi (div '$%pi 2))
1478 ;; Complete elliptic integral
1479 (ftake '%elliptic_kc m))
1480 (t
1481 ;; Nothing to do
1482 (give-up)))))
1483
1484 (def-simplifier elliptic_e (phi m)
1485 (let (args)
1486 (cond ((float-numerical-eval-p phi m)
1487 ;; Numerically evaluate it
1488 (elliptic-e ($float phi) ($float m)))
1489 ((complex-float-numerical-eval-p phi m)
1490 (complexify (bigfloat::bf-elliptic-e (complex ($float ($realpart phi)) ($float ($imagpart phi)))
1491 (complex ($float ($realpart m)) ($float ($imagpart m))))))
1492 ((bigfloat-numerical-eval-p phi m)
1493 (to (bigfloat::bf-elliptic-e (bigfloat:to ($bfloat phi))
1494 (bigfloat:to ($bfloat m)))))
1495 ((setf args (complex-bigfloat-numerical-eval-p phi m))
1496 (destructuring-bind (phi m)
1497 args
1498 (to (bigfloat::bf-elliptic-e (bigfloat:to ($bfloat phi))
1499 (bigfloat:to ($bfloat m))))))
1500 ((zerop1 phi)
1501 0)
1502 ((zerop1 m)
1503 ;; A&S 17.4.23
1504 phi)
1505 ((onep1 m)
1506 ;; A&S 17.4.25, but handle periodicity:
1507 ;; elliptic_e(x,m) = elliptic_e(x-%pi*round(x/%pi), m)
1508 ;; + 2*round(x/%pi)*elliptic_ec(m)
1509 ;;
1510 ;; Or
1511 ;;
1512 ;; elliptic_e(x,1) = sin(x-%pi*round(x/%pi)) + 2*round(x/%pi)*elliptic_ec(m)
1513 ;;
1514 (let ((mult-pi (ftake '%round (div phi '$%pi))))
1515 (add (ftake '%sin (sub phi
1516 (mul '$%pi
1517 mult-pi)))
1518 (mul 2
1519 (mul mult-pi
1520 (ftake '%elliptic_ec m))))))
1521 ((alike1 phi (div '$%pi 2))
1522 ;; Complete elliptic integral
1523 (ftake '%elliptic_ec m))
1524 ((and ($numberp phi)
1525 (let ((r ($round (div phi '$%pi))))
1526 (and ($numberp r)
1527 (not (zerop1 r)))))
1528 ;; Handle the case where phi is a number where we can apply
1529 ;; the periodicity property without blowing up the
1530 ;; expression.
1531 (add (ftake '%elliptic_e
1532 (add phi
1533 (mul (mul -1 '$%pi)
1534 (ftake '%round (div phi '$%pi))))
1535 m)
1536 (mul 2
1537 (mul (ftake '%round (div phi '$%pi))
1538 (ftake '%elliptic_ec m)))))
1539 (t
1540 ;; Nothing to do
1541 (give-up)))))
1542
1543 ;; Complete elliptic integrals
1544 ;;
1545 ;; elliptic_kc(m) = elliptic_f(%pi/2, m)
1546 ;;
1547 ;; elliptic_ec(m) = elliptic_e(%pi/2, m)
1548 ;;
1549
1550 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1551
1552 ;;; We support a simplim%function. The function is looked up in simplimit and
1553 ;;; handles specific values of the function.
1554
1555 (defprop %elliptic_kc simplim%elliptic_kc simplim%function)
1556
1557 (defun simplim%elliptic_kc (expr var val)
1558 ;; Look for the limit of the argument
1559 (let ((m (limit (cadr expr) var val 'think)))
1560 (cond ((onep1 m)
1561 ;; For an argument 1 return $infinity.
1562 '$infinity)
1563 (t
1564 ;; All other cases are handled by the simplifier of the function.
1565 (simplify (list '(%elliptic_kc) m))))))
1566
1567 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1568
1569 (def-simplifier elliptic_kc (m)
1570 (let (args)
1571 (cond ((onep1 m)
1572 ;; elliptic_kc(1) is complex infinity. Maxima can not handle
1573 ;; infinities correctly, throw a Maxima error.
1574 (merror
1575 (intl:gettext "elliptic_kc: elliptic_kc(~:M) is undefined.")
1576 m))
1577 ((float-numerical-eval-p m)
1578 ;; Numerically evaluate it
1579 (to (elliptic-k ($float m))))
1580 ((complex-float-numerical-eval-p m)
1581 (complexify (bigfloat::bf-elliptic-k (complex ($float ($realpart m)) ($float ($imagpart m))))))
1582 ((setf args (complex-bigfloat-numerical-eval-p m))
1583 (destructuring-bind (m)
1584 args
1585 (to (bigfloat::bf-elliptic-k (bigfloat:to ($bfloat m))))))
1586 ((zerop1 m)
1587 (mul 1//2 '$%pi))
1588 ((alike1 m 1//2)
1589 ;; http://functions.wolfram.com/EllipticIntegrals/EllipticK/03/01/
1590 ;;
1591 ;; elliptic_kc(1/2) = 8*%pi^(3/2)/gamma(-1/4)^2
1592 (div (mul 8 (power '$%pi (div 3 2)))
1593 (power (gm (div -1 4)) 2)))
1594 ((eql -1 m)
1595 ;; elliptic_kc(-1) = gamma(1/4)^2/(4*sqrt(2*%pi))
1596 (div (power (gm (div 1 4)) 2)
1597 (mul 4 (power (mul 2 '$%pi) 1//2))))
1598 ((alike1 m (add 17 (mul -12 (power 2 1//2))))
1599 ;; elliptic_kc(17-12*sqrt(2)) = 2*(2+sqrt(2))*%pi^(3/2)/gamma(-1/4)^2
1600 (div (mul 2 (mul (add 2 (power 2 1//2))
1601 (power '$%pi (div 3 2))))
1602 (power (gm (div -1 4)) 2)))
1603 ((or (alike1 m (div (add 2 (power 3 1//2))
1604 4))
1605 (alike1 m (add (div (power 3 1//2)
1606 4)
1607 1//2)))
1608 ;; elliptic_kc((sqrt(3)+2)/4) = sqrt(%pi)*gamma(1/3)/gamma(5/6).
1609 ;;
1610 ;; First evaluate this integral, where y = sqrt(1+t^3).
1611 ;;
1612 ;; integrate(1/y,t,-1,inf) = integrate(1/y,t,-1,0) + integrate(1/y,t,0,inf).
1613 ;;
1614 ;; The second integral, maxima gives beta(1/6,1/3)/3.
1615 ;;
1616 ;; For the first, we can use the change of variable x=-u^(1/3) to get
1617 ;;
1618 ;; integrate(1/sqrt(1-u)/u^(2/3),u,0,1)
1619 ;;
1620 ;; which is a beta integral that maxima can evaluate to
1621 ;; beta(1/3,1/2)/3. Then we see the value of the initial
1622 ;; integral is
1623 ;;
1624 ;; beta(1/6,1/3)/3 + beta(1/3,1/2)/3
1625 ;;
1626 ;; (Thanks to Guilherme Namen for this derivation on the mailing list, 2023-03-09.)
1627 ;;
1628 ;; We can simplify this expression by converting to gamma functions:
1629 ;;
1630 ;; beta(1/6,1/3)/3 + beta(1/3,1/2)/3 =
1631 ;; (gamma(1/3)*(gamma(1/6)*gamma(5/6)+%pi))/(3*sqrt(%pi)*gamma(5/6));
1632 ;;
1633 ;; Using the reflection formula gamma(1-z)*gamma(z) =
1634 ;; %pi/sin(%pi*z), we can write gamma(1/6)*gamma(5/6) =
1635 ;; %pi/sin(%pi*1/6) = 2*%pi. Finally, we have
1636 ;;
1637 ;; sqrt(%pi)*gamma(1/3)/gamma(5/6);
1638 ;;
1639 ;; All that remains is to show that integrate(1/y,t) can be
1640 ;; written as an inverse_jacobi_cn function with modulus
1641 ;; (sqrt(3)+2)/4.
1642 ;;
1643 ;; First apply the substitution
1644 ;;
1645 ;; s = (t+sqrt(3)+1)/(t-sqrt(3)+1). We then have the integral
1646 ;;
1647 ;; C*integrate(1/sqrt(s^2-1)/sqrt(s^2+4*sqrt(3)+7),s)
1648 ;;
1649 ;; where C is some constant. From A&S 14.4.49, we can see
1650 ;; this integral is the inverse_jacobi_nc function with
1651 ;; modulus of (4*sqrt(3)+7)/(4*sqrt(3)+7+1) =
1652 ;; (sqrt(3)+2)/4.
1653 (div (mul (power '$%pi 1//2)
1654 (ftake '%gamma (div 1 3)))
1655 (ftake '%gamma (div 5 6))))
1656 ($hypergeometric_representation
1657 ;; See http://functions.wolfram.com/08.02.26.0001.01
1658 ;;
1659 ;; elliptic_kc(z) = %pi/2*%f[2,1]([1/2,1/2],[1], z)
1660 ;;
1661 (mul (div '$%pi 2)
1662 (ftake '%hypergeometric
1663 (make-mlist 1//2 1//2)
1664 (make-mlist 1)
1665 m)))
1666 (t
1667 ;; Nothing to do
1668 (give-up)))))
1669
1670 (defprop %elliptic_kc
1671 ((m)
1672 ;; diff wrt m
1673 ((mtimes)
1674 ((rat) 1 2)
1675 ((mplus) ((%elliptic_ec) m)
1676 ((mtimes) -1
1677 ((%elliptic_kc) m)
1678 ((mplus) 1 ((mtimes) -1 m))))
1679 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
1680 ((mexpt) m -1)))
1681 grad)
1682
1683 (def-simplifier elliptic_ec (m)
1684 (let (args)
1685 (cond ((float-numerical-eval-p m)
1686 ;; Numerically evaluate it
1687 (elliptic-ec ($float m)))
1688 ((setf args (complex-float-numerical-eval-p m))
1689 (destructuring-bind (m)
1690 args
1691 (complexify (bigfloat::bf-elliptic-ec (bigfloat:to ($float m))))))
1692 ((setf args (complex-bigfloat-numerical-eval-p m))
1693 (destructuring-bind (m)
1694 args
1695 (to (bigfloat::bf-elliptic-ec (bigfloat:to ($bfloat m))))))
1696 ;; Some special cases we know about.
1697 ((zerop1 m)
1698 (div '$%pi 2))
1699 ((onep1 m)
1700 1)
1701 ((alike1 m 1//2)
1702 ;; elliptic_ec(1/2). Use the identity
1703 ;;
1704 ;; elliptic_ec(z)*elliptic_kc(1-z) - elliptic_kc(z)*elliptic_kc(1-z)
1705 ;; + elliptic_ec(1-z)*elliptic_kc(z) = %pi/2;
1706 ;;
1707 ;; Let z = 1/2 to get
1708 ;;
1709 ;; %pi^(3/2)*'elliptic_ec(1/2)/gamma(3/4)^2-%pi^3/(4*gamma(3/4)^4) = %pi/2
1710 ;;
1711 ;; since we know that elliptic_kc(1/2) = %pi^(3/2)/(2*gamma(3/4)^2). Hence
1712 ;;
1713 ;; elliptic_ec(1/2)
1714 ;; = (2*%pi*gamma(3/4)^4+%pi^3)/(4*%pi^(3/2)*gamma(3/4)^2)
1715 ;; = gamma(3/4)^2/(2*sqrt(%pi))+%pi^(3/2)/(4*gamma(3/4)^2)
1716 ;;
1717 (add (div (power (ftake '%gamma (div 3 4)) 2)
1718 (mul 2 (power '$%pi 1//2)))
1719 (div (power '$%pi (div 3 2))
1720 (mul 4 (power (ftake '%gamma (div 3 4)) 2)))))
1721 ((zerop1 (add 1 m))
1722 ;; elliptic_ec(-1). Use the identity
1723 ;; http://functions.wolfram.com/08.01.17.0002.01
1724 ;;
1725 ;;
1726 ;; elliptic_ec(z) = sqrt(1 - z)*elliptic_ec(z/(z-1))
1727 ;;
1728 ;; Let z = -1 to get
1729 ;;
1730 ;; elliptic_ec(-1) = sqrt(2)*elliptic_ec(1/2)
1731 ;;
1732 ;; Should we expand out elliptic_ec(1/2) using the above result?
1733 (mul (power 2 1//2)
1734 (ftake '%elliptic_ec 1//2)))
1735 ($hypergeometric_representation
1736 ;; See http://functions.wolfram.com/08.01.26.0001.01
1737 ;;
1738 ;; elliptic_ec(z) = %pi/2*%f[2,1]([-1/2,1/2],[1], z)
1739 ;;
1740 (mul (div '$%pi 2)
1741 (ftake '%hypergeometric
1742 (make-mlist -1//2 1//2)
1743 (make-mlist 1)
1744 m)))
1745 (t
1746 ;; Nothing to do
1747 (give-up)))))
1748
1749 (defprop %elliptic_ec
1750 ((m)
1751 ((mtimes) ((rat) 1 2)
1752 ((mplus) ((%elliptic_ec) m)
1753 ((mtimes) -1 ((%elliptic_kc)
1754 m)))
1755 ((mexpt) m -1)))
1756 grad)
1757
1758 ;;
1759 ;; Elliptic integral of the third kind:
1760 ;;
1761 ;; (A&S 17.2.14)
1762 ;;
1763 ;; phi
1764 ;; /
1765 ;; [ 1
1766 ;; PI(n;phi|m) = I ----------------------------------- ds
1767 ;; ] 2 2
1768 ;; / SQRT(1 - m SIN (s)) (1 - n SIN (s))
1769 ;; 0
1770 ;;
1771 ;; As with E and F, we do not use the modular angle alpha but the
1772 ;; parameter m = sin(alpha)^2.
1773 ;;
1774 (def-simplifier elliptic_pi (n phi m)
1775 (let (args)
1776 (cond
1777 ((float-numerical-eval-p n phi m)
1778 ;; Numerically evaluate it
1779 (elliptic-pi ($float n) ($float phi) ($float m)))
1780 ((setf args (complex-float-numerical-eval-p n phi m))
1781 (destructuring-bind (n phi m)
1782 args
1783 (elliptic-pi (bigfloat:to ($float n))
1784 (bigfloat:to ($float phi))
1785 (bigfloat:to ($float m)))))
1786 ((bigfloat-numerical-eval-p n phi m)
1787 (to (bigfloat::bf-elliptic-pi (bigfloat:to n)
1788 (bigfloat:to phi)
1789 (bigfloat:to m))))
1790 ((setq args (complex-bigfloat-numerical-eval-p n phi m))
1791 (destructuring-bind (n phi m)
1792 args
1793 (to (bigfloat::bf-elliptic-pi (bigfloat:to n)
1794 (bigfloat:to phi)
1795 (bigfloat:to m)))))
1796 ((zerop1 n)
1797 (ftake '%elliptic_f phi m))
1798 ((zerop1 m)
1799 ;; 3 cases depending on n < 1, n > 1, or n = 1.
1800 (let ((s (asksign (add -1 n))))
1801 (case s
1802 ($positive
1803 (div (ftake '%atanh (mul (power (add n -1) 1//2)
1804 (ftake '%tan phi)))
1805 (power (add n -1) 1//2)))
1806 ($negative
1807 (div (ftake '%atan (mul (power (sub 1 n) 1//2)
1808 (ftake '%tan phi)))
1809 (power (sub 1 n) 1//2)))
1810 ($zero
1811 (ftake '%tan phi)))))
1812 (t
1813 ;; Nothing to do
1814 (give-up)))))
1815
1816 ;; Complete elliptic-pi. That is phi = %pi/2. Then
1817 ;; elliptic_pi(n,m)
1818 ;; = Rf(0, 1-m,1) + Rj(0,1-m,1-n)*n/3;
1819 (defun elliptic-pi-complete (n m)
1820 (to (bigfloat:+ (bigfloat::bf-rf 0 (- 1 m) 1)
1821 (bigfloat:* 1/3 n (bigfloat::bf-rj 0 (- 1 m) 1 (- 1 n))))))
1822
1823 ;; To compute elliptic_pi for all z, we use the property
1824 ;; (http://functions.wolfram.com/08.06.16.0002.01)
1825 ;;
1826 ;; elliptic_pi(n, z + %pi*k, m)
1827 ;; = 2*k*elliptic_pi(n, %pi/2, m) + elliptic_pi(n, z, m)
1828 ;;
1829 ;; So we are left with computing the integral for 0 <= z < %pi. Using
1830 ;; Carlson's formulation produces the wrong values for %pi/2 < z <
1831 ;; %pi. How to do that?
1832 ;;
1833 ;; Let
1834 ;;
1835 ;; I(a,b) = integrate(1/(1-n*sin(x)^2)/sqrt(1 - m*sin(x)^2), x, a, b)
1836 ;;
1837 ;; That is, I(a,b) is the integral for the elliptic_pi function but
1838 ;; with a lower limit of a and an upper limit of b.
1839 ;;
1840 ;; Then, we want to compute I(0, z), with %pi <= z < %pi. Let w = z +
1841 ;; %pi/2, 0 <= w < %pi/2. Then
1842 ;;
1843 ;; I(0, w+%pi/2) = I(0, %pi/2) + I(%pi/2, w+%pi/2)
1844 ;;
1845 ;; To evaluate I(%pi/2, w+%pi/2), use a change of variables:
1846 ;;
1847 ;; changevar('integrate(1/(1-n*sin(x)^2)/sqrt(1 - m*sin(x)^2), x, %pi/2, w + %pi/2),
1848 ;; x-%pi+u,u,x)
1849 ;;
1850 ;; = integrate(-1/(sqrt(1-m*sin(u)^2)*(1-n*sin(u)^2)),u,%pi/2-w,%pi/2)
1851 ;; = I(%pi/2-w,%pi/2)
1852 ;; = I(0,%pi/2) - I(0,%pi/2-w)
1853 ;;
1854 ;; Thus,
1855 ;;
1856 ;; I(0,%pi/2+w) = 2*I(0,%pi/2) - I(0,%pi/2-w)
1857 ;;
1858 ;; This allows us to compute the general result with 0 <= z < %pi
1859 ;;
1860 ;; I(0, k*%pi + z) = 2*k*I(0,%pi/2) + I(0,z);
1861 ;;
1862 ;; If 0 <= z < %pi/2, then the we are done. If %pi/2 <= z < %pi, let
1863 ;; z = w+%pi/2. Then
1864 ;;
1865 ;; I(0,z) = 2*I(0,%pi/2) - I(0,%pi/2-w)
1866 ;;
1867 ;; Or, since w = z-%pi/2:
1868 ;;
1869 ;; I(0,z) = 2*I(0,%pi/2) - I(0,%pi-z)
1870
1871 (defun elliptic-pi (n phi m)
1872 ;; elliptic_pi(n, -phi, m) = -elliptic_pi(n, phi, m). That is, it
1873 ;; is an odd function of phi.
1874 (when (minusp (realpart phi))
1875 (return-from elliptic-pi (- (elliptic-pi n (- phi) m))))
1876
1877 ;; Note: Carlson's DRJ has n defined as the negative of the n given
1878 ;; in A&S.
1879 (flet ((base (n phi m)
1880 ;; elliptic_pi(n,phi,m) =
1881 ;; sin(phi)*Rf(cos(phi)^2, 1-m*sin(phi)^2, 1)
1882 ;; - (-n / 3) * sin(phi)^3
1883 ;; * Rj(cos(phi)^2, 1-m*sin(phi)^2, 1, 1-n*sin(phi)^2)
1884 (let* ((nn (- n))
1885 (sin-phi (sin phi))
1886 (cos-phi (cos phi))
1887 (k (sqrt m))
1888 (k2sin (* (- 1 (* k sin-phi))
1889 (+ 1 (* k sin-phi)))))
1890 (- (* sin-phi (bigfloat::bf-rf (expt cos-phi 2) k2sin 1.0))
1891 (* (/ nn 3) (expt sin-phi 3)
1892 (bigfloat::bf-rj (expt cos-phi 2) k2sin 1.0
1893 (- 1 (* n (expt sin-phi 2)))))))))
1894 ;; FIXME: Reducing the arg by pi has significant round-off.
1895 ;; Consider doing something better.
1896 (let* ((cycles (round (realpart phi) pi))
1897 (rem (- phi (* cycles pi))))
1898 (let ((complete (elliptic-pi-complete n m)))
1899 (to (+ (* 2 cycles complete)
1900 (base n rem m)))))))
1901
1902 ;;; Deriviatives from functions.wolfram.com
1903 ;;; http://functions.wolfram.com/EllipticIntegrals/EllipticPi3/20/
1904 (defprop %elliptic_pi
1905 ((n z m)
1906 ;Derivative wrt first argument
1907 ((mtimes) ((rat) 1 2)
1908 ((mexpt) ((mplus) m ((mtimes) -1 n)) -1)
1909 ((mexpt) ((mplus) -1 n) -1)
1910 ((mplus)
1911 ((mtimes) ((mexpt) n -1)
1912 ((mplus) ((mtimes) -1 m) ((mexpt) n 2))
1913 ((%elliptic_pi) n z m))
1914 ((%elliptic_e) z m)
1915 ((mtimes) ((mplus) m ((mtimes) -1 n)) ((mexpt) n -1)
1916 ((%elliptic_f) z m))
1917 ((mtimes) ((rat) -1 2) n
1918 ((mexpt)
1919 ((mplus) 1 ((mtimes) -1 m ((mexpt) ((%sin) z) 2)))
1920 ((rat) 1 2))
1921 ((mexpt)
1922 ((mplus) 1 ((mtimes) -1 n ((mexpt) ((%sin) z) 2)))
1923 -1)
1924 ((%sin) ((mtimes) 2 z)))))
1925 ;derivative wrt second argument
1926 ((mtimes)
1927 ((mexpt)
1928 ((mplus) 1 ((mtimes) -1 m ((mexpt) ((%sin) z) 2)))
1929 ((rat) -1 2))
1930 ((mexpt)
1931 ((mplus) 1 ((mtimes) -1 n ((mexpt) ((%sin) z) 2))) -1))
1932 ;Derivative wrt third argument
1933 ((mtimes) ((rat) 1 2)
1934 ((mexpt) ((mplus) ((mtimes) -1 m) n) -1)
1935 ((mplus) ((%elliptic_pi) n z m)
1936 ((mtimes) ((mexpt) ((mplus) -1 m) -1)
1937 ((%elliptic_e) z m))
1938 ((mtimes) ((rat) -1 2) ((mexpt) ((mplus) -1 m) -1) m
1939 ((mexpt)
1940 ((mplus) 1 ((mtimes) -1 m ((mexpt) ((%sin) z) 2)))
1941 ((rat) -1 2))
1942 ((%sin) ((mtimes) 2 z))))))
1943 grad)
1944
1945 (in-package #:bigfloat)
1946 ;; Translation of Jim FitzSimons' bigfloat implementation of elliptic
1947 ;; integrals from http://www.getnet.com/~cherry/elliptbf3.mac.
1948 ;;
1949 ;; The algorithms are based on B.C. Carlson's "Numerical Computation
1950 ;; of Real or Complex Elliptic Integrals". These are updated to the
1951 ;; algorithms in Journal of Computational and Applied Mathematics 118
1952 ;; (2000) 71-85 "Reduction Theorems for Elliptic Integrands with the
1953 ;; Square Root of two quadritic factors"
1954 ;;
1955
1956 ;; NOTE: Despite the names indicating these are for bigfloat numbers,
1957 ;; the algorithms and routines are generic and will work with floats
1958 ;; and bigfloats.
1959
1960 (defun bferrtol (&rest args)
1961 ;; Compute error tolerance as sqrt(2^(-fpprec)). Not sure this is
1962 ;; quite right, but it makes the routines more accurate as fpprec
1963 ;; increases.
1964 (sqrt (reduce #'min (mapcar #'(lambda (x)
1965 (if (rationalp (realpart x))
1966 maxima::+flonum-epsilon+
1967 (epsilon x)))
1968 args))))
1969
1970 ;; rc(x,y) = integrate(1/2*(t+x)^(-1/2)/(t+y), t, 0, inf)
1971 ;;
1972 ;; log(x) = (x-1)*rc(((1+x)/2)^2, x), x > 0
1973 ;; asin(x) = x * rc(1-x^2, 1), |x|<= 1
1974 ;; acos(x) = sqrt(1-x^2)*rc(x^2,1), 0 <= x <=1
1975 ;; atan(x) = x * rc(1,1+x^2)
1976 ;; asinh(x) = x * rc(1+x^2,1)
1977 ;; acosh(x) = sqrt(x^2-1) * rc(x^2,1), x >= 1
1978 ;; atanh(x) = x * rc(1,1-x^2), |x|<=1
1979
1980 (defun bf-rc (x y)
1981 (let ((yn y)
1982 xn z w a an pwr4 n epslon lambda sn s)
1983 (cond ((and (zerop (imagpart yn))
1984 (minusp (realpart yn)))
1985 (setf xn (- x y))
1986 (setf yn (- yn))
1987 (setf z yn)
1988 (setf w (sqrt (/ x xn))))
1989 (t
1990 (setf xn x)
1991 (setf z yn)
1992 (setf w 1)))
1993 (setf a (/ (+ xn yn yn) 3))
1994 (setf epslon (/ (abs (- a xn)) (bferrtol x y)))
1995 (setf an a)
1996 (setf pwr4 1)
1997 (setf n 0)
1998 (loop while (> (* epslon pwr4) (abs an))
1999 do
2000 (setf pwr4 (/ pwr4 4))
2001 (setf lambda (+ (* 2 (sqrt xn) (sqrt yn)) yn))
2002 (setf an (/ (+ an lambda) 4))
2003 (setf xn (/ (+ xn lambda) 4))
2004 (setf yn (/ (+ yn lambda) 4))
2005 (incf n))
2006 ;; c2=3/10,c3=1/7,c4=3/8,c5=9/22,c6=159/208,c7=9/8
2007 (setf sn (/ (* pwr4 (- z a)) an))
2008 (setf s (* sn sn (+ 3/10
2009 (* sn (+ 1/7
2010 (* sn (+ 3/8
2011 (* sn (+ 9/22
2012 (* sn (+ 159/208
2013 (* sn 9/8))))))))))))
2014 (/ (* w (+ 1 s))
2015 (sqrt an))))
2016
2017
2018
2019 ;; See https://dlmf.nist.gov/19.16.E5:
2020 ;;
2021 ;; rd(x,y,z) = integrate(3/2/sqrt(t+x)/sqrt(t+y)/sqrt(t+z)/(t+z), t, 0, inf)
2022 ;;
2023 ;; rd(1,1,1) = 1
2024 ;; E(K) = rf(0, 1-K^2, 1) - (K^2/3)*rd(0,1-K^2,1)
2025 ;;
2026 ;; B = integrate(s^2/sqrt(1-s^4), s, 0 ,1)
2027 ;; = beta(3/4,1/2)/4
2028 ;; = sqrt(%pi)*gamma(3/4)/gamma(1/4)
2029 ;; = 1/3*rd(0,2,1)
2030
2031 (defun bf-rd (x y z)
2032 (let* ((xn x)
2033 (yn y)
2034 (zn z)
2035 (a (/ (+ xn yn (* 3 zn)) 5))
2036 (epslon (/ (max (abs (- a xn))
2037 (abs (- a yn))
2038 (abs (- a zn)))
2039 (bferrtol x y z)))
2040 (an a)
2041 (sigma 0)
2042 (power4 1)
2043 (n 0)
2044 xnroot ynroot znroot lam)
2045 (loop while (> (* power4 epslon) (abs an))
2046 do
2047 (setf xnroot (sqrt xn))
2048 (setf ynroot (sqrt yn))
2049 (setf znroot (sqrt zn))
2050 (setf lam (+ (* xnroot ynroot)
2051 (* xnroot znroot)
2052 (* ynroot znroot)))
2053 (setf sigma (+ sigma (/ power4
2054 (* znroot (+ zn lam)))))
2055 (setf power4 (* power4 1/4))
2056 (setf xn (* (+ xn lam) 1/4))
2057 (setf yn (* (+ yn lam) 1/4))
2058 (setf zn (* (+ zn lam) 1/4))
2059 (setf an (* (+ an lam) 1/4))
2060 (incf n))
2061 ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26
2062 (let* ((xndev (/ (* (- a x) power4) an))
2063 (yndev (/ (* (- a y) power4) an))
2064 (zndev (- (* (+ xndev yndev) 1/3)))
2065 (ee2 (- (* xndev yndev) (* 6 zndev zndev)))
2066 (ee3 (* (- (* 3 xndev yndev)
2067 (* 8 zndev zndev))
2068 zndev))
2069 (ee4 (* 3 (- (* xndev yndev) (* zndev zndev)) zndev zndev))
2070 (ee5 (* xndev yndev zndev zndev zndev))
2071 (s (+ 1
2072 (* -3/14 ee2)
2073 (* 1/6 ee3)
2074 (* 9/88 ee2 ee2)
2075 (* -3/22 ee4)
2076 (* -9/52 ee2 ee3)
2077 (* 3/26 ee5)
2078 (* -1/16 ee2 ee2 ee2)
2079 (* 3/10 ee3 ee3)
2080 (* 3/20 ee2 ee4)
2081 (* 45/272 ee2 ee2 ee3)
2082 (* -9/68 (+ (* ee2 ee5) (* ee3 ee4))))))
2083 (+ (* 3 sigma)
2084 (/ (* power4 s)
2085 (expt an 3/2))))))
2086
2087 ;; See https://dlmf.nist.gov/19.16.E1
2088 ;;
2089 ;; rf(x,y,z) = 1/2*integrate(1/(sqrt(t+x)*sqrt(t+y)*sqrt(t+z)), t, 0, inf);
2090 ;;
2091 ;; rf(1,1,1) = 1
2092 (defun bf-rf (x y z)
2093 (let* ((xn x)
2094 (yn y)
2095 (zn z)
2096 (a (/ (+ xn yn zn) 3))
2097 (epslon (/ (max (abs (- a xn))
2098 (abs (- a yn))
2099 (abs (- a zn)))
2100 (bferrtol x y z)))
2101 (an a)
2102 (power4 1)
2103 (n 0)
2104 xnroot ynroot znroot lam)
2105 (loop while (> (* power4 epslon) (abs an))
2106 do
2107 (setf xnroot (sqrt xn))
2108 (setf ynroot (sqrt yn))
2109 (setf znroot (sqrt zn))
2110 (setf lam (+ (* xnroot ynroot)
2111 (* xnroot znroot)
2112 (* ynroot znroot)))
2113 (setf power4 (* power4 1/4))
2114 (setf xn (* (+ xn lam) 1/4))
2115 (setf yn (* (+ yn lam) 1/4))
2116 (setf zn (* (+ zn lam) 1/4))
2117 (setf an (* (+ an lam) 1/4))
2118 (incf n))
2119 ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26
2120 (let* ((xndev (/ (* (- a x) power4) an))
2121 (yndev (/ (* (- a y) power4) an))
2122 (zndev (- (+ xndev yndev)))
2123 (ee2 (- (* xndev yndev) (* 6 zndev zndev)))
2124 (ee3 (* xndev yndev zndev))
2125 (s (+ 1
2126 (* -1/10 ee2)
2127 (* 1/14 ee3)
2128 (* 1/24 ee2 ee2)
2129 (* -3/44 ee2 ee3))))
2130 (/ s (sqrt an)))))
2131
2132 (defun bf-rj1 (x y z p)
2133 (let* ((xn x)
2134 (yn y)
2135 (zn z)
2136 (pn p)
2137 (en (* (- pn xn)
2138 (- pn yn)
2139 (- pn zn)))
2140 (sigma 0)
2141 (power4 1)
2142 (k 0)
2143 (a (/ (+ xn yn zn pn pn) 5))
2144 (epslon (/ (max (abs (- a xn))
2145 (abs (- a yn))
2146 (abs (- a zn))
2147 (abs (- a pn)))
2148 (bferrtol x y z p)))
2149 (an a)
2150 xnroot ynroot znroot pnroot lam dn)
2151 (loop while (> (* power4 epslon) (abs an))
2152 do
2153 (setf xnroot (sqrt xn))
2154 (setf ynroot (sqrt yn))
2155 (setf znroot (sqrt zn))
2156 (setf pnroot (sqrt pn))
2157 (setf lam (+ (* xnroot ynroot)
2158 (* xnroot znroot)
2159 (* ynroot znroot)))
2160 (setf dn (* (+ pnroot xnroot)
2161 (+ pnroot ynroot)
2162 (+ pnroot znroot)))
2163 (setf sigma (+ sigma
2164 (/ (* power4
2165 (bf-rc 1 (+ 1 (/ en (* dn dn)))))
2166 dn)))
2167 (setf power4 (* power4 1/4))
2168 (setf en (/ en 64))
2169 (setf xn (* (+ xn lam) 1/4))
2170 (setf yn (* (+ yn lam) 1/4))
2171 (setf zn (* (+ zn lam) 1/4))
2172 (setf pn (* (+ pn lam) 1/4))
2173 (setf an (* (+ an lam) 1/4))
2174 (incf k))
2175 (let* ((xndev (/ (* (- a x) power4) an))
2176 (yndev (/ (* (- a y) power4) an))
2177 (zndev (/ (* (- a z) power4) an))
2178 (pndev (* -0.5 (+ xndev yndev zndev)))
2179 (ee2 (+ (* xndev yndev)
2180 (* xndev zndev)
2181 (* yndev zndev)
2182 (* -3 pndev pndev)))
2183 (ee3 (+ (* xndev yndev zndev)
2184 (* 2 ee2 pndev)
2185 (* 4 pndev pndev pndev)))
2186 (ee4 (* (+ (* 2 xndev yndev zndev)
2187 (* ee2 pndev)
2188 (* 3 pndev pndev pndev))
2189 pndev))
2190 (ee5 (* xndev yndev zndev pndev pndev))
2191 (s (+ 1
2192 (* -3/14 ee2)
2193 (* 1/6 ee3)
2194 (* 9/88 ee2 ee2)
2195 (* -3/22 ee4)
2196 (* -9/52 ee2 ee3)
2197 (* 3/26 ee5)
2198 (* -1/16 ee2 ee2 ee2)
2199 (* 3/10 ee3 ee3)
2200 (* 3/20 ee2 ee4)
2201 (* 45/272 ee2 ee2 ee3)
2202 (* -9/68 (+ (* ee2 ee5) (* ee3 ee4))))))
2203 (+ (* 6 sigma)
2204 (/ (* power4 s)
2205 (sqrt (* an an an)))))))
2206
2207 (defun bf-rj (x y z p)
2208 (let* ((xn x)
2209 (yn y)
2210 (zn z)
2211 (qn (- p)))
2212 (cond ((and (and (zerop (imagpart xn)) (>= (realpart xn) 0))
2213 (and (zerop (imagpart yn)) (>= (realpart yn) 0))
2214 (and (zerop (imagpart zn)) (>= (realpart zn) 0))
2215 (and (zerop (imagpart qn)) (> (realpart qn) 0)))
2216 (destructuring-bind (xn yn zn)
2217 (sort (list xn yn zn) #'<)
2218 (let* ((pn (+ yn (* (- zn yn) (/ (- yn xn) (+ yn qn)))))
2219 (s (- (* (- pn yn) (bf-rj1 xn yn zn pn))
2220 (* 3 (bf-rf xn yn zn)))))
2221 (setf s (+ s (* 3 (sqrt (/ (* xn yn zn)
2222 (+ (* xn zn) (* pn qn))))
2223 (bf-rc (+ (* xn zn) (* pn qn)) (* pn qn)))))
2224 (/ s (+ yn qn)))))
2225 (t
2226 (bf-rj1 x y z p)))))
2227
2228 (defun bf-rg (x y z)
2229 (* 0.5
2230 (+ (* z (bf-rf x y z))
2231 (* (- z x)
2232 (- y z)
2233 (bf-rd x y z)
2234 1/3)
2235 (sqrt (/ (* x y) z)))))
2236
2237 ;; elliptic_f(phi,m) = sin(phi)*rf(cos(phi)^2, 1-m*sin(phi)^2,1)
2238 (defun bf-elliptic-f (phi m)
2239 (flet ((base (phi m)
2240 (cond ((= m 1)
2241 ;; F(z|1) = log(tan(z/2+%pi/4))
2242 (log (tan (+ (/ phi 2) (/ (%pi phi) 4)))))
2243 (t
2244 (let ((s (sin phi))
2245 (c (cos phi)))
2246 (* s (bf-rf (* c c) (- 1 (* m s s)) 1)))))))
2247 ;; Handle periodicity (see elliptic-f)
2248 (let* ((bfpi (%pi phi))
2249 (period (round (realpart phi) bfpi)))
2250 (+ (base (- phi (* bfpi period)) m)
2251 (if (zerop period)
2252 0
2253 (* 2 period (bf-elliptic-k m)))))))
2254
2255 ;; elliptic_kc(k) = rf(0, 1-k^2,1)
2256 ;;
2257 ;; or
2258 ;; elliptic_kc(m) = rf(0, 1-m,1)
2259
2260 (defun bf-elliptic-k (m)
2261 (cond ((= m 0)
2262 (if (maxima::$bfloatp m)
2263 (maxima::$bfloat (maxima::div 'maxima::$%pi 2))
2264 (float (/ pi 2) 1e0)))
2265 ((= m 1)
2266 (maxima::merror
2267 (intl:gettext "elliptic_kc: elliptic_kc(1) is undefined.")))
2268 (t
2269 (bf-rf 0 (- 1 m) 1))))
2270
2271 ;; elliptic_e(phi, k) = sin(phi)*rf(cos(phi)^2,1-k^2*sin(phi)^2,1)
2272 ;; - (k^2/3)*sin(phi)^3*rd(cos(phi)^2, 1-k^2*sin(phi)^2,1)
2273 ;;
2274 ;;
2275 ;; or
2276 ;; elliptic_e(phi, m) = sin(phi)*rf(cos(phi)^2,1-m*sin(phi)^2,1)
2277 ;; - (m/3)*sin(phi)^3*rd(cos(phi)^2, 1-m*sin(phi)^2,1)
2278 ;;
2279 (defun bf-elliptic-e (phi m)
2280 (flet ((base (phi m)
2281 (let* ((s (sin phi))
2282 (c (cos phi))
2283 (c2 (* c c))
2284 (s2 (- 1 (* m s s))))
2285 (- (* s (bf-rf c2 s2 1))
2286 (* (/ m 3) (* s s s) (bf-rd c2 s2 1))))))
2287 ;; Elliptic E is quasi-periodic wrt to phi:
2288 ;;
2289 ;; E(z|m) = E(z - %pi*round(Re(z)/%pi)|m) + 2*round(Re(z)/%pi)*E(m)
2290 (let* ((bfpi (%pi phi))
2291 (period (round (realpart phi) bfpi)))
2292 (+ (base (- phi (* bfpi period)) m)
2293 (* 2 period (bf-elliptic-ec m))))))
2294
2295
2296 ;; elliptic_ec(k) = rf(0,1-k^2,1) - (k^2/3)*rd(0,1-k^2,1);
2297 ;;
2298 ;; or
2299 ;; elliptic_ec(m) = rf(0,1-m,1) - (m/3)*rd(0,1-m,1);
2300
2301 (defun bf-elliptic-ec (m)
2302 (cond ((= m 0)
2303 (if (typep m 'bigfloat)
2304 (bigfloat (maxima::$bfloat (maxima::div 'maxima::$%pi 2)))
2305 (float (/ pi 2) 1e0)))
2306 ((= m 1)
2307 (if (typep m 'bigfloat)
2308 (bigfloat 1)
2309 1e0))
2310 (t
2311 (let ((m1 (- 1 m)))
2312 (- (bf-rf 0 m1 1)
2313 (* m 1/3 (bf-rd 0 m1 1)))))))
2314
2315 (defun bf-elliptic-pi-complete (n m)
2316 (+ (bf-rf 0 (- 1 m) 1)
2317 (* 1/3 n (bf-rj 0 (- 1 m) 1 (- 1 n)))))
2318
2319 (defun bf-elliptic-pi (n phi m)
2320 ;; Note: Carlson's DRJ has n defined as the negative of the n given
2321 ;; in A&S.
2322 (flet ((base (n phi m)
2323 (let* ((nn (- n))
2324 (sin-phi (sin phi))
2325 (cos-phi (cos phi))
2326 (k (sqrt m))
2327 (k2sin (* (- 1 (* k sin-phi))
2328 (+ 1 (* k sin-phi)))))
2329 (- (* sin-phi (bf-rf (expt cos-phi 2) k2sin 1.0))
2330 (* (/ nn 3) (expt sin-phi 3)
2331 (bf-rj (expt cos-phi 2) k2sin 1.0
2332 (- 1 (* n (expt sin-phi 2)))))))))
2333 ;; FIXME: Reducing the arg by pi has significant round-off.
2334 ;; Consider doing something better.
2335 (let* ((bf-pi (%pi (realpart phi)))
2336 (cycles (round (realpart phi) bf-pi))
2337 (rem (- phi (* cycles bf-pi))))
2338 (let ((complete (bf-elliptic-pi-complete n m)))
2339 (+ (* 2 cycles complete)
2340 (base n rem m))))))
2341
2342 ;; Compute inverse_jacobi_sn, for float or bigfloat args.
2343 (defun bf-inverse-jacobi-sn (u m)
2344 (* u (bf-rf (- 1 (* u u))
2345 (- 1 (* m u u))
2346 1)))
2347
2348 ;; Compute inverse_jacobi_dn. We use the following identity
2349 ;; from Gradshteyn & Ryzhik, 8.153.6
2350 ;;
2351 ;; w = dn(z|m) = cn(sqrt(m)*z, 1/m)
2352 ;;
2353 ;; Solve for z to get
2354 ;;
2355 ;; z = inverse_jacobi_dn(w,m)
2356 ;; = 1/sqrt(m) * inverse_jacobi_cn(w, 1/m)
2357 (defun bf-inverse-jacobi-dn (w m)
2358 (cond ((= w 1)
2359 (float 0 w))
2360 ((= m 1)
2361 ;; jacobi_dn(x,1) = sech(x) so the inverse is asech(x)
2362 (maxima::take '(maxima::%asech) (maxima::to w)))
2363 (t
2364 ;; We should do something better to make sure that things
2365 ;; that should be real are real.
2366 (/ (to (maxima::take '(maxima::%inverse_jacobi_cn)
2367 (maxima::to w)
2368 (maxima::to (/ m))))
2369 (sqrt m)))))
2370
2371 (in-package :maxima)
2372
2373 ;; Define Carlson's elliptic integrals.
2374
2375 (def-simplifier carlson_rc (x y)
2376 (let (args)
2377 (flet ((calc (x y)
2378 (flet ((floatify (z)
2379 ;; If z is a complex rational, convert to a
2380 ;; complex double-float. Otherwise, leave it as
2381 ;; is. If we don't do this, %i is handled as
2382 ;; #c(0 1), which makes bf-rc use single-float
2383 ;; arithmetic instead of the desired
2384 ;; double-float.
2385 (if (and (complexp z) (rationalp (realpart z)))
2386 (complex (float (realpart z))
2387 (float (imagpart z)))
2388 z)))
2389 (to (bigfloat::bf-rc (floatify (bigfloat:to x))
2390 (floatify (bigfloat:to y)))))))
2391 ;; See comments from bf-rc
2392 (cond ((float-numerical-eval-p x y)
2393 (calc ($float x) ($float y)))
2394 ((bigfloat-numerical-eval-p x y)
2395 (calc ($bfloat x) ($bfloat y)))
2396 ((setf args (complex-float-numerical-eval-p x y))
2397 (destructuring-bind (x y)
2398 args
2399 (calc ($float x) ($float y))))
2400 ((setf args (complex-bigfloat-numerical-eval-p x y))
2401 (destructuring-bind (x y)
2402 args
2403 (calc ($bfloat x) ($bfloat y))))
2404 ((and (zerop1 x)
2405 (onep1 y))
2406 ;; rc(0, 1) = %pi/2
2407 (div '$%pi 2))
2408 ((and (zerop1 x)
2409 (alike1 y (div 1 4)))
2410 ;; rc(0,1/4) = %pi
2411 '$%pi)
2412 ((and (eql x 2)
2413 (onep1 y))
2414 ;; rc(2,1) = 1/2*integrate(1/sqrt(t+2)/(t+1), t, 0, inf)
2415 ;; = (log(sqrt(2)+1)-log(sqrt(2)-1))/2
2416 ;; ratsimp(logcontract(%)),algebraic:
2417 ;; = -log(3-2^(3/2))/2
2418 ;; = -log(sqrt(3-2^(3/2)))
2419 ;; = -log(sqrt(2)-1)
2420 ;; = log(1/(sqrt(2)-1))
2421 ;; ratsimp(%),algebraic;
2422 ;; = log(sqrt(2)+1)
2423 (ftake '%log (add 1 (power 2 1//2))))
2424 ((and (alike x '$%i)
2425 (alike y (add 1 '$%i)))
2426 ;; rc(%i, %i+1) = 1/2*integrate(1/sqrt(t+%i)/(t+%i+1), t, 0, inf)
2427 ;; = %pi/2-atan((-1)^(1/4))
2428 ;; ratsimp(logcontract(ratsimp(rectform(%o42)))),algebraic;
2429 ;; = (%i*log(3-2^(3/2))+%pi)/4
2430 ;; = (%i*log(3-2^(3/2)))/4+%pi/4
2431 ;; = %i*log(sqrt(3-2^(3/2)))/2+%pi/4
2432 ;; sqrtdenest(%);
2433 ;; = %pi/4 + %i*log(sqrt(2)-1)/2
2434 (add (div '$%pi 4)
2435 (mul '$%i
2436 1//2
2437 (ftake '%log (sub (power 2 1//2) 1)))))
2438 ((and (zerop1 x)
2439 (alike1 y '$%i))
2440 ;; rc(0,%i) = 1/2*integrate(1/(sqrt(t)*(t+%i)), t, 0, inf)
2441 ;; = -((sqrt(2)*%i-sqrt(2))*%pi)/4
2442 ;; = ((1-%i)*%pi)/2^(3/2)
2443 (div (mul (sub 1 '$%i)
2444 '$%pi)
2445 (power 2 3//2)))
2446 ((and (alike1 x y)
2447 (eq ($sign ($realpart x)) '$pos))
2448 ;; carlson_rc(x,x) = 1/2*integrate(1/sqrt(t+x)/(t+x), t, 0, inf)
2449 ;; = 1/sqrt(x)
2450 (power x -1//2))
2451 ((and (alike1 x (power (div (add 1 y) 2) 2))
2452 (eq ($sign ($realpart y)) '$pos))
2453 ;; Rc(((1+x)/2)^2,x) = log(x)/(x-1) for x > 0.
2454 ;;
2455 ;; This is done by looking at Rc(x,y) and seeing if
2456 ;; ((1+y)/2)^2 is the same as x.
2457 (div (ftake '%log y)
2458 (sub y 1)))
2459 (t
2460 (give-up))))))
2461
2462 (def-simplifier carlson_rd (x y z)
2463 (let (args)
2464 (flet ((calc (x y z)
2465 (to (bigfloat::bf-rd (bigfloat:to x)
2466 (bigfloat:to y)
2467 (bigfloat:to z)))))
2468 ;; See https://dlmf.nist.gov/19.20.E18
2469 (cond ((and (eql x 1)
2470 (eql x y)
2471 (eql y z))
2472 ;; Rd(1,1,1) = 1
2473 1)
2474 ((and (alike1 x y)
2475 (alike1 y z))
2476 ;; Rd(x,x,x) = x^(-3/2)
2477 (power x (div -3 2)))
2478 ((and (zerop1 x)
2479 (alike1 y z))
2480 ;; Rd(0,y,y) = 3/4*%pi*y^(-3/2)
2481 (mul (div 3 4)
2482 '$%pi
2483 (power y (div -3 2))))
2484 ((alike1 y z)
2485 ;; Rd(x,y,y) = 3/(2*(y-x))*(Rc(x, y) - sqrt(x)/y)
2486 (mul (div 3 (mul 2 (sub y x)))
2487 (sub (ftake '%carlson_rc x y)
2488 (div (power x 1//2)
2489 y))))
2490 ((alike1 x y)
2491 ;; Rd(x,x,z) = 3/(z-x)*(Rc(z,x) - 1/sqrt(z))
2492 (mul (div 3 (sub z x))
2493 (sub (ftake '%carlson_rc z x)
2494 (div 1 (power z 1//2)))))
2495 ((and (eql z 1)
2496 (or (and (eql x 0)
2497 (eql y 2))
2498 (and (eql x 2)
2499 (eql y 0))))
2500 ;; Rd(0,2,1) = 3*(gamma(3/4)^2)/sqrt(2*%pi)
2501 ;; See https://dlmf.nist.gov/19.20.E22.
2502 ;;
2503 ;; But that's the same as
2504 ;; 3*sqrt(%pi)*gamma(3/4)/gamma(1/4). We can see this by
2505 ;; taking the ratio to get
2506 ;; gamma(1/4)*gamma(3/4)/sqrt(2)*%pi. But
2507 ;; gamma(1/4)*gamma(3/4) = beta(1/4,3/4) = sqrt(2)*%pi.
2508 ;; Hence, the ratio is 1.
2509 ;;
2510 ;; Note also that Rd(x,y,z) = Rd(y,x,z)
2511 (mul 3
2512 (power '$%pi 1//2)
2513 (div (ftake '%gamma (div 3 4))
2514 (ftake '%gamma (div 1 4)))))
2515 ((and (or (eql x 0) (eql y 0))
2516 (eql z 1))
2517 ;; 1/3*m*Rd(0,1-m,1) = K(m) - E(m).
2518 ;; See https://dlmf.nist.gov/19.25.E1
2519 ;;
2520 ;; Thus, Rd(0,y,1) = 3/(1-y)*(K(1-y) - E(1-y))
2521 ;;
2522 ;; Note that Rd(x,y,z) = Rd(y,x,z).
2523 (let ((m (sub 1 y)))
2524 (mul (div 3 m)
2525 (sub (ftake '%elliptic_kc m)
2526 (ftake '%elliptic_ec m)))))
2527 ((or (and (eql x 0)
2528 (eql y 1))
2529 (and (eql x 1)
2530 (eql y 0)))
2531 ;; 1/3*m*(1-m)*Rd(0,1,1-m) = E(m) - (1-m)*K(m)
2532 ;; See https://dlmf.nist.gov/19.25.E1
2533 ;;
2534 ;; Thus
2535 ;; Rd(0,1,z) = 3/(z*(1-z))*(E(1-z) - z*K(1-z))
2536 ;; Recall that Rd(x,y,z) = Rd(y,x,z).
2537 (mul (div 3 (mul z (sub 1 z)))
2538 (sub (ftake '%elliptic_ec (sub 1 z))
2539 (mul z
2540 (ftake '%elliptic_kc (sub 1 z))))))
2541 ((float-numerical-eval-p x y z)
2542 (calc ($float x) ($float y) ($float z)))
2543 ((bigfloat-numerical-eval-p x y z)
2544 (calc ($bfloat x) ($bfloat y) ($bfloat z)))
2545 ((setf args (complex-float-numerical-eval-p x y z))
2546 (destructuring-bind (x y z)
2547 args
2548 (calc ($float x) ($float y) ($float z))))
2549 ((setf args (complex-bigfloat-numerical-eval-p x y z))
2550 (destructuring-bind (x y z)
2551 args
2552 (calc ($bfloat x) ($bfloat y) ($bfloat z))))
2553 (t
2554 (give-up))))))
2555
2556 (def-simplifier carlson_rf (x y z)
2557 (let (args)
2558 (flet ((calc (x y z)
2559 (to (bigfloat::bf-rf (bigfloat:to x)
2560 (bigfloat:to y)
2561 (bigfloat:to z)))))
2562 ;; See https://dlmf.nist.gov/19.20.i
2563 (cond ((and (alike1 x y)
2564 (alike1 y z))
2565 ;; Rf(x,x,x) = x^(-1/2)
2566 (power x -1//2))
2567 ((and (zerop1 x)
2568 (alike1 y z))
2569 ;; Rf(0,y,y) = 1/2*%pi*y^(-1/2)
2570 (mul 1//2 '$%pi
2571 (power y -1//2)))
2572 ((alike1 y z)
2573 (ftake '%carlson_rc x y))
2574 ((some #'(lambda (args)
2575 (destructuring-bind (x y z)
2576 args
2577 (and (zerop1 x)
2578 (eql y 1)
2579 (eql z 2))))
2580 (list (list x y z)
2581 (list x z y)
2582 (list y x z)
2583 (list y z x)
2584 (list z x y)
2585 (list z y x)))
2586 ;; Rf(0,1,2) = (gamma(1/4))^2/(4*sqrt(2*%pi))
2587 ;;
2588 ;; And Rf is symmetric in all the args, so check every
2589 ;; permutation too. This could probably be simplified
2590 ;; without consing all the lists, but I'm lazy.
2591 (div (power (ftake '%gamma (div 1 4)) 2)
2592 (mul 4 (power (mul 2 '$%pi) 1//2))))
2593 ((some #'(lambda (args)
2594 (destructuring-bind (x y z)
2595 args
2596 (and (alike1 x '$%i)
2597 (alike1 y (mul -1 '$%i))
2598 (eql z 0))))
2599 (list (list x y z)
2600 (list x z y)
2601 (list y x z)
2602 (list y z x)
2603 (list z x y)
2604 (list z y x)))
2605 ;; rf(%i, -%i, 0)
2606 ;; = 1/2*integrate(1/sqrt(t^2+1)/sqrt(t),t,0,inf)
2607 ;; = beta(1/4,1/4)/4;
2608 ;; makegamma(%)
2609 ;; = gamma(1/4)^2/(4*sqrt(%pi))
2610 ;;
2611 ;; Rf is symmetric, so check all the permutations too.
2612 (div (power (ftake '%gamma (div 1 4)) 2)
2613 (mul 4 (power '$%pi 1//2))))
2614 ((setf args
2615 (some #'(lambda (args)
2616 (destructuring-bind (x y z)
2617 args
2618 ;; Check that x = 0 and z = 1, and
2619 ;; return y.
2620 (and (zerop1 x)
2621 (eql z 1)
2622 y)))
2623 (list (list x y z)
2624 (list x z y)
2625 (list y x z)
2626 (list y z x)
2627 (list z x y)
2628 (list z y x))))
2629 ;; Rf(0,1-m,1) = elliptic_kc(m).
2630 ;; See https://dlmf.nist.gov/19.25.E1
2631 (ftake '%elliptic_kc (sub 1 args)))
2632 ((some #'(lambda (args)
2633 (destructuring-bind (x y z)
2634 args
2635 (and (alike1 x '$%i)
2636 (alike1 y (mul -1 '$%i))
2637 (eql z 0))))
2638 (list (list x y z)
2639 (list x z y)
2640 (list y x z)
2641 (list y z x)
2642 (list z x y)
2643 (list z y x)))
2644 ;; rf(%i, -%i, 0)
2645 ;; = 1/2*integrate(1/sqrt(t^2+1)/sqrt(t),t,0,inf)
2646 ;; = beta(1/4,1/4)/4;
2647 ;; makegamma(%)
2648 ;; = gamma(1/4)^2/(4*sqrt(%pi))
2649 ;;
2650 ;; Rf is symmetric, so check all the permutations too.
2651 (div (power (ftake '%gamma (div 1 4)) 2)
2652 (mul 4 (power '$%pi 1//2))))
2653 ((float-numerical-eval-p x y z)
2654 (calc ($float x) ($float y) ($float z)))
2655 ((bigfloat-numerical-eval-p x y z)
2656 (calc ($bfloat x) ($bfloat y) ($bfloat z)))
2657 ((setf args (complex-float-numerical-eval-p x y z))
2658 (destructuring-bind (x y z)
2659 args
2660 (calc ($float x) ($float y) ($float z))))
2661 ((setf args (complex-bigfloat-numerical-eval-p x y z))
2662 (destructuring-bind (x y z)
2663 args
2664 (calc ($bfloat x) ($bfloat y) ($bfloat z))))
2665 (t
2666 (give-up))))))
2667
2668 (def-simplifier carlson_rj (x y z p)
2669 (let (args)
2670 (flet ((calc (x y z p)
2671 (to (bigfloat::bf-rj (bigfloat:to x)
2672 (bigfloat:to y)
2673 (bigfloat:to z)
2674 (bigfloat:to p)))))
2675 ;; See https://dlmf.nist.gov/19.20.iii
2676 (cond ((and (alike1 x y)
2677 (alike1 y z)
2678 (alike1 z p))
2679 ;; Rj(x,x,x,x) = x^(-3/2)
2680 (power x (div -3 2)))
2681 ((alike1 z p)
2682 ;; Rj(x,y,z,z) = Rd(x,y,z)
2683 (ftake '%carlson_rd x y z))
2684 ((and (zerop1 x)
2685 (alike1 y z))
2686 ;; Rj(0,y,y,p) = 3*%pi/(2*(y*sqrt(p)+p*sqrt(y)))
2687 (div (mul 3 '$%pi)
2688 (mul 2
2689 (add (mul y (power p 1//2))
2690 (mul p (power y 1//2))))))
2691 ((alike1 y z)
2692 ;; Rj(x,y,y,p) = 3/(p-y)*(Rc(x,y) - Rc(x,p))
2693 (mul (div 3 (sub p y))
2694 (sub (ftake '%carlson_rc x y)
2695 (ftake '%carlson_rc x p))))
2696 ((and (alike1 y z)
2697 (alike1 y p))
2698 ;; Rj(x,y,y,y) = Rd(x,y,y)
2699 (ftake '%carlson_rd x y y))
2700 ((float-numerical-eval-p x y z p)
2701 (calc ($float x) ($float y) ($float z) ($float p)))
2702 ((bigfloat-numerical-eval-p x y z p)
2703 (calc ($bfloat x) ($bfloat y) ($bfloat z) ($bfloat p)))
2704 ((setf args (complex-float-numerical-eval-p x y z p))
2705 (destructuring-bind (x y z p)
2706 args
2707 (calc ($float x) ($float y) ($float z) ($float p))))
2708 ((setf args (complex-bigfloat-numerical-eval-p x y z p))
2709 (destructuring-bind (x y z p)
2710 args
2711 (calc ($bfloat x) ($bfloat y) ($bfloat z) ($bfloat p))))
2712 (t
2713 (give-up))))))
2714
2715 ;;; Other Jacobian elliptic functions
2716
2717 ;; jacobi_ns(u,m) = 1/jacobi_sn(u,m)
2718 (defprop %jacobi_ns
2719 ((u m)
2720 ;; diff wrt u
2721 ((mtimes) -1 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
2722 ((mexpt) ((%jacobi_sn) u m) -2))
2723 ;; diff wrt m
2724 ((mtimes) -1 ((mexpt) ((%jacobi_sn) u m) -2)
2725 ((mplus)
2726 ((mtimes) ((rat) 1 2)
2727 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
2728 ((mexpt) ((%jacobi_cn) u m) 2)
2729 ((%jacobi_sn) u m))
2730 ((mtimes) ((rat) 1 2) ((mexpt) m -1)
2731 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
2732 ((mplus) u
2733 ((mtimes) -1
2734 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
2735 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
2736 m)))))))
2737 grad)
2738
2739 (def-simplifier jacobi_ns (u m)
2740 (let (coef args)
2741 (cond
2742 ((float-numerical-eval-p u m)
2743 (to (bigfloat:/ (bigfloat::sn (bigfloat:to ($float u))
2744 (bigfloat:to ($float m))))))
2745 ((setf args (complex-float-numerical-eval-p u m))
2746 (destructuring-bind (u m)
2747 args
2748 (to (bigfloat:/ (bigfloat::sn (bigfloat:to ($float u))
2749 (bigfloat:to ($float m)))))))
2750 ((bigfloat-numerical-eval-p u m)
2751 (let ((uu (bigfloat:to ($bfloat u)))
2752 (mm (bigfloat:to ($bfloat m))))
2753 (to (bigfloat:/ (bigfloat::sn uu mm)))))
2754 ((setf args (complex-bigfloat-numerical-eval-p u m))
2755 (destructuring-bind (u m)
2756 args
2757 (let ((uu (bigfloat:to ($bfloat u)))
2758 (mm (bigfloat:to ($bfloat m))))
2759 (to (bigfloat:/ (bigfloat::sn uu mm))))))
2760 ((zerop1 m)
2761 ;; A&S 16.6.10
2762 (ftake '%csc u))
2763 ((onep1 m)
2764 ;; A&S 16.6.10
2765 (ftake '%coth u))
2766 ((zerop1 u)
2767 (dbz-err1 'jacobi_ns))
2768 ((and $trigsign (mminusp* u))
2769 ;; ns is odd
2770 (neg (ftake* '%jacobi_ns (neg u) m)))
2771 ((and $triginverses
2772 (listp u)
2773 (member (caar u) '(%inverse_jacobi_sn
2774 %inverse_jacobi_ns
2775 %inverse_jacobi_cn
2776 %inverse_jacobi_nc
2777 %inverse_jacobi_dn
2778 %inverse_jacobi_nd
2779 %inverse_jacobi_sc
2780 %inverse_jacobi_cs
2781 %inverse_jacobi_sd
2782 %inverse_jacobi_ds
2783 %inverse_jacobi_cd
2784 %inverse_jacobi_dc))
2785 (alike1 (third u) m))
2786 (cond ((eq (caar u) '%inverse_jacobi_ns)
2787 (second u))
2788 (t
2789 ;; Express in terms of sn:
2790 ;; ns(x) = 1/sn(x)
2791 (div 1 (ftake '%jacobi_sn u m)))))
2792 ;; A&S 16.20 (Jacobi's Imaginary transformation)
2793 ((and $%iargs (multiplep u '$%i))
2794 ;; ns(i*u) = 1/sn(i*u) = -i/sc(u,m1) = -i*cs(u,m1)
2795 (neg (mul '$%i
2796 (ftake* '%jacobi_cs (coeff u '$%i 1) (add 1 (neg m))))))
2797 ((setq coef (kc-arg2 u m))
2798 ;; A&S 16.8.10
2799 ;;
2800 ;; ns(m*K+u) = 1/sn(m*K+u)
2801 ;;
2802 (destructuring-bind (lin const)
2803 coef
2804 (cond ((integerp lin)
2805 (ecase (mod lin 4)
2806 (0
2807 ;; ns(4*m*K+u) = ns(u)
2808 ;; ns(0) = infinity
2809 (if (zerop1 const)
2810 (dbz-err1 'jacobi_ns)
2811 (ftake '%jacobi_ns const m)))
2812 (1
2813 ;; ns(4*m*K + K + u) = ns(K+u) = dc(u)
2814 ;; ns(K) = 1
2815 (if (zerop1 const)
2816 1
2817 (ftake '%jacobi_dc const m)))
2818 (2
2819 ;; ns(4*m*K+2*K + u) = ns(2*K+u) = -ns(u)
2820 ;; ns(2*K) = infinity
2821 (if (zerop1 const)
2822 (dbz-err1 'jacobi_ns)
2823 (neg (ftake '%jacobi_ns const m))))
2824 (3
2825 ;; ns(4*m*K+3*K+u) = ns(2*K + K + u) = -ns(K+u) = -dc(u)
2826 ;; ns(3*K) = -1
2827 (if (zerop1 const)
2828 -1
2829 (neg (ftake '%jacobi_dc const m))))))
2830 ((and (alike1 lin 1//2)
2831 (zerop1 const))
2832 (div 1 (ftake '%jacobi_sn u m)))
2833 (t
2834 (give-up)))))
2835 (t
2836 ;; Nothing to do
2837 (give-up)))))
2838
2839 ;; jacobi_nc(u,m) = 1/jacobi_cn(u,m)
2840 (defprop %jacobi_nc
2841 ((u m)
2842 ;; wrt u
2843 ((mtimes) ((mexpt) ((%jacobi_cn) u m) -2)
2844 ((%jacobi_dn) u m) ((%jacobi_sn) u m))
2845 ;; wrt m
2846 ((mtimes) -1 ((mexpt) ((%jacobi_cn) u m) -2)
2847 ((mplus)
2848 ((mtimes) ((rat) -1 2)
2849 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
2850 ((%jacobi_cn) u m) ((mexpt) ((%jacobi_sn) u m) 2))
2851 ((mtimes) ((rat) -1 2) ((mexpt) m -1)
2852 ((%jacobi_dn) u m) ((%jacobi_sn) u m)
2853 ((mplus) u
2854 ((mtimes) -1 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
2855 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m)) m)))))))
2856 grad)
2857
2858 (def-simplifier jacobi_nc (u m)
2859 (let (coef args)
2860 (cond
2861 ((float-numerical-eval-p u m)
2862 (to (bigfloat:/ (bigfloat::cn (bigfloat:to ($float u))
2863 (bigfloat:to ($float m))))))
2864 ((setf args (complex-float-numerical-eval-p u m))
2865 (destructuring-bind (u m)
2866 args
2867 (to (bigfloat:/ (bigfloat::cn (bigfloat:to ($float u))
2868 (bigfloat:to ($float m)))))))
2869 ((bigfloat-numerical-eval-p u m)
2870 (let ((uu (bigfloat:to ($bfloat u)))
2871 (mm (bigfloat:to ($bfloat m))))
2872 (to (bigfloat:/ (bigfloat::cn uu mm)))))
2873 ((setf args (complex-bigfloat-numerical-eval-p u m))
2874 (destructuring-bind (u m)
2875 args
2876 (let ((uu (bigfloat:to ($bfloat u)))
2877 (mm (bigfloat:to ($bfloat m))))
2878 (to (bigfloat:/ (bigfloat::cn uu mm))))))
2879 ((zerop1 u)
2880 1)
2881 ((zerop1 m)
2882 ;; A&S 16.6.8
2883 (ftake '%sec u))
2884 ((onep1 m)
2885 ;; A&S 16.6.8
2886 (ftake '%cosh u))
2887 ((and $trigsign (mminusp* u))
2888 ;; nc is even
2889 (ftake* '%jacobi_nc (neg u) m))
2890 ((and $triginverses
2891 (listp u)
2892 (member (caar u) '(%inverse_jacobi_sn
2893 %inverse_jacobi_ns
2894 %inverse_jacobi_cn
2895 %inverse_jacobi_nc
2896 %inverse_jacobi_dn
2897 %inverse_jacobi_nd
2898 %inverse_jacobi_sc
2899 %inverse_jacobi_cs
2900 %inverse_jacobi_sd
2901 %inverse_jacobi_ds
2902 %inverse_jacobi_cd
2903 %inverse_jacobi_dc))
2904 (alike1 (third u) m))
2905 (cond ((eq (caar u) '%inverse_jacobi_nc)
2906 (second u))
2907 (t
2908 ;; Express in terms of cn:
2909 ;; nc(x) = 1/cn(x)
2910 (div 1 (ftake '%jacobi_cn u m)))))
2911 ;; A&S 16.20 (Jacobi's Imaginary transformation)
2912 ((and $%iargs (multiplep u '$%i))
2913 ;; nc(i*u) = 1/cn(i*u) = 1/nc(u,1-m) = cn(u,1-m)
2914 (ftake* '%jacobi_cn (coeff u '$%i 1) (add 1 (neg m))))
2915 ((setq coef (kc-arg2 u m))
2916 ;; A&S 16.8.8
2917 ;;
2918 ;; nc(u) = 1/cn(u)
2919 ;;
2920 (destructuring-bind (lin const)
2921 coef
2922 (cond ((integerp lin)
2923 (ecase (mod lin 4)
2924 (0
2925 ;; nc(4*m*K+u) = nc(u)
2926 ;; nc(0) = 1
2927 (if (zerop1 const)
2928 1
2929 (ftake '%jacobi_nc const m)))
2930 (1
2931 ;; nc(4*m*K+K+u) = nc(K+u) = -ds(u)/sqrt(1-m)
2932 ;; nc(K) = infinity
2933 (if (zerop1 const)
2934 (dbz-err1 'jacobi_nc)
2935 (neg (div (ftake '%jacobi_ds const m)
2936 (power (sub 1 m) 1//2)))))
2937 (2
2938 ;; nc(4*m*K+2*K+u) = nc(2*K+u) = -nc(u)
2939 ;; nc(2*K) = -1
2940 (if (zerop1 const)
2941 -1
2942 (neg (ftake '%jacobi_nc const m))))
2943 (3
2944 ;; nc(4*m*K+3*K+u) = nc(3*K+u) = nc(2*K+K+u) =
2945 ;; -nc(K+u) = ds(u)/sqrt(1-m)
2946 ;;
2947 ;; nc(3*K) = infinity
2948 (if (zerop1 const)
2949 (dbz-err1 'jacobi_nc)
2950 (div (ftake '%jacobi_ds const m)
2951 (power (sub 1 m) 1//2))))))
2952 ((and (alike1 1//2 lin)
2953 (zerop1 const))
2954 (div 1 (ftake '%jacobi_cn u m)))
2955 (t
2956 (give-up)))))
2957 (t
2958 ;; Nothing to do
2959 (give-up)))))
2960
2961 ;; jacobi_nd(u,m) = 1/jacobi_dn(u,m)
2962 (defprop %jacobi_nd
2963 ((u m)
2964 ;; wrt u
2965 ((mtimes) m ((%jacobi_cn) u m)
2966 ((mexpt) ((%jacobi_dn) u m) -2) ((%jacobi_sn) u m))
2967 ;; wrt m
2968 ((mtimes) -1 ((mexpt) ((%jacobi_dn) u m) -2)
2969 ((mplus)
2970 ((mtimes) ((rat) -1 2)
2971 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
2972 ((%jacobi_dn) u m)
2973 ((mexpt) ((%jacobi_sn) u m) 2))
2974 ((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
2975 ((%jacobi_sn) u m)
2976 ((mplus) u
2977 ((mtimes) -1
2978 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
2979 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
2980 m)))))))
2981 grad)
2982
2983 (def-simplifier jacobi_nd (u m)
2984 (let (coef args)
2985 (cond
2986 ((float-numerical-eval-p u m)
2987 (to (bigfloat:/ (bigfloat::dn (bigfloat:to ($float u))
2988 (bigfloat:to ($float m))))))
2989 ((setf args (complex-float-numerical-eval-p u m))
2990 (destructuring-bind (u m)
2991 args
2992 (to (bigfloat:/ (bigfloat::dn (bigfloat:to ($float u))
2993 (bigfloat:to ($float m)))))))
2994 ((bigfloat-numerical-eval-p u m)
2995 (let ((uu (bigfloat:to ($bfloat u)))
2996 (mm (bigfloat:to ($bfloat m))))
2997 (to (bigfloat:/ (bigfloat::dn uu mm)))))
2998 ((setf args (complex-bigfloat-numerical-eval-p u m))
2999 (destructuring-bind (u m)
3000 args
3001 (let ((uu (bigfloat:to ($bfloat u)))
3002 (mm (bigfloat:to ($bfloat m))))
3003 (to (bigfloat:/ (bigfloat::dn uu mm))))))
3004 ((zerop1 u)
3005 1)
3006 ((zerop1 m)
3007 ;; A&S 16.6.6
3008 1)
3009 ((onep1 m)
3010 ;; A&S 16.6.6
3011 (ftake '%cosh u))
3012 ((and $trigsign (mminusp* u))
3013 ;; nd is even
3014 (ftake* '%jacobi_nd (neg u) m))
3015 ((and $triginverses
3016 (listp u)
3017 (member (caar u) '(%inverse_jacobi_sn
3018 %inverse_jacobi_ns
3019 %inverse_jacobi_cn
3020 %inverse_jacobi_nc
3021 %inverse_jacobi_dn
3022 %inverse_jacobi_nd
3023 %inverse_jacobi_sc
3024 %inverse_jacobi_cs
3025 %inverse_jacobi_sd
3026 %inverse_jacobi_ds
3027 %inverse_jacobi_cd
3028 %inverse_jacobi_dc))
3029 (alike1 (third u) m))
3030 (cond ((eq (caar u) '%inverse_jacobi_nd)
3031 (second u))
3032 (t
3033 ;; Express in terms of dn:
3034 ;; nd(x) = 1/dn(x)
3035 (div 1 (ftake '%jacobi_dn u m)))))
3036 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3037 ((and $%iargs (multiplep u '$%i))
3038 ;; nd(i*u) = 1/dn(i*u) = 1/dc(u,1-m) = cd(u,1-m)
3039 (ftake* '%jacobi_cd (coeff u '$%i 1) (add 1 (neg m))))
3040 ((setq coef (kc-arg2 u m))
3041 ;; A&S 16.8.6
3042 ;;
3043 (destructuring-bind (lin const)
3044 coef
3045 (cond ((integerp lin)
3046 ;; nd has period 2K
3047 (ecase (mod lin 2)
3048 (0
3049 ;; nd(2*m*K+u) = nd(u)
3050 ;; nd(0) = 1
3051 (if (zerop1 const)
3052 1
3053 (ftake '%jacobi_nd const m)))
3054 (1
3055 ;; nd(2*m*K+K+u) = nd(K+u) = dn(u)/sqrt(1-m)
3056 ;; nd(K) = 1/sqrt(1-m)
3057 (if (zerop1 const)
3058 (power (sub 1 m) -1//2)
3059 (div (ftake '%jacobi_nd const m)
3060 (power (sub 1 m) 1//2))))))
3061 (t
3062 (give-up)))))
3063 (t
3064 ;; Nothing to do
3065 (give-up)))))
3066
3067 ;; jacobi_sc(u,m) = jacobi_sn/jacobi_cn
3068 (defprop %jacobi_sc
3069 ((u m)
3070 ;; wrt u
3071 ((mtimes) ((mexpt) ((%jacobi_cn) u m) -2)
3072 ((%jacobi_dn) u m))
3073 ;; wrt m
3074 ((mplus)
3075 ((mtimes) ((mexpt) ((%jacobi_cn) u m) -1)
3076 ((mplus)
3077 ((mtimes) ((rat) 1 2)
3078 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3079 ((mexpt) ((%jacobi_cn) u m) 2)
3080 ((%jacobi_sn) u m))
3081 ((mtimes) ((rat) 1 2) ((mexpt) m -1)
3082 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
3083 ((mplus) u
3084 ((mtimes) -1
3085 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3086 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3087 m))))))
3088 ((mtimes) -1 ((mexpt) ((%jacobi_cn) u m) -2)
3089 ((%jacobi_sn) u m)
3090 ((mplus)
3091 ((mtimes) ((rat) -1 2)
3092 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3093 ((%jacobi_cn) u m)
3094 ((mexpt) ((%jacobi_sn) u m) 2))
3095 ((mtimes) ((rat) -1 2) ((mexpt) m -1)
3096 ((%jacobi_dn) u m) ((%jacobi_sn) u m)
3097 ((mplus) u
3098 ((mtimes) -1
3099 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3100 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3101 m))))))))
3102 grad)
3103
3104 (def-simplifier jacobi_sc (u m)
3105 (let (coef args)
3106 (cond
3107 ((float-numerical-eval-p u m)
3108 (let ((fu (bigfloat:to ($float u)))
3109 (fm (bigfloat:to ($float m))))
3110 (to (bigfloat:/ (bigfloat::sn fu fm) (bigfloat::cn fu fm)))))
3111 ((setf args (complex-float-numerical-eval-p u m))
3112 (destructuring-bind (u m)
3113 args
3114 (let ((fu (bigfloat:to ($float u)))
3115 (fm (bigfloat:to ($float m))))
3116 (to (bigfloat:/ (bigfloat::sn fu fm) (bigfloat::cn fu fm))))))
3117 ((bigfloat-numerical-eval-p u m)
3118 (let ((uu (bigfloat:to ($bfloat u)))
3119 (mm (bigfloat:to ($bfloat m))))
3120 (to (bigfloat:/ (bigfloat::sn uu mm)
3121 (bigfloat::cn uu mm)))))
3122 ((setf args (complex-bigfloat-numerical-eval-p u m))
3123 (destructuring-bind (u m)
3124 args
3125 (let ((uu (bigfloat:to ($bfloat u)))
3126 (mm (bigfloat:to ($bfloat m))))
3127 (to (bigfloat:/ (bigfloat::sn uu mm)
3128 (bigfloat::cn uu mm))))))
3129 ((zerop1 u)
3130 0)
3131 ((zerop1 m)
3132 ;; A&S 16.6.9
3133 (ftake '%tan u))
3134 ((onep1 m)
3135 ;; A&S 16.6.9
3136 (ftake '%sinh u))
3137 ((and $trigsign (mminusp* u))
3138 ;; sc is odd
3139 (neg (ftake* '%jacobi_sc (neg u) m)))
3140 ((and $triginverses
3141 (listp u)
3142 (member (caar u) '(%inverse_jacobi_sn
3143 %inverse_jacobi_ns
3144 %inverse_jacobi_cn
3145 %inverse_jacobi_nc
3146 %inverse_jacobi_dn
3147 %inverse_jacobi_nd
3148 %inverse_jacobi_sc
3149 %inverse_jacobi_cs
3150 %inverse_jacobi_sd
3151 %inverse_jacobi_ds
3152 %inverse_jacobi_cd
3153 %inverse_jacobi_dc))
3154 (alike1 (third u) m))
3155 (cond ((eq (caar u) '%inverse_jacobi_sc)
3156 (second u))
3157 (t
3158 ;; Express in terms of sn and cn
3159 ;; sc(x) = sn(x)/cn(x)
3160 (div (ftake '%jacobi_sn u m)
3161 (ftake '%jacobi_cn u m)))))
3162 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3163 ((and $%iargs (multiplep u '$%i))
3164 ;; sc(i*u) = sn(i*u)/cn(i*u) = i*sc(u,m1)/nc(u,m1) = i*sn(u,m1)
3165 (mul '$%i
3166 (ftake* '%jacobi_sn (coeff u '$%i 1) (add 1 (neg m)))))
3167 ((setq coef (kc-arg2 u m))
3168 ;; A&S 16.8.9
3169 ;; sc(2*m*K+u) = sc(u)
3170 (destructuring-bind (lin const)
3171 coef
3172 (cond ((integerp lin)
3173 (ecase (mod lin 2)
3174 (0
3175 ;; sc(2*m*K+ u) = sc(u)
3176 ;; sc(0) = 0
3177 (if (zerop1 const)
3178 1
3179 (ftake '%jacobi_sc const m)))
3180 (1
3181 ;; sc(2*m*K + K + u) = sc(K+u)= - cs(u)/sqrt(1-m)
3182 ;; sc(K) = infinity
3183 (if (zerop1 const)
3184 (dbz-err1 'jacobi_sc)
3185 (mul -1
3186 (div (ftake* '%jacobi_cs const m)
3187 (power (sub 1 m) 1//2)))))))
3188 ((and (alike1 lin 1//2)
3189 (zerop1 const))
3190 ;; From A&S 16.3.3 and 16.5.2:
3191 ;; sc(1/2*K) = 1/(1-m)^(1/4)
3192 (power (sub 1 m) (div -1 4)))
3193 (t
3194 (give-up)))))
3195 (t
3196 ;; Nothing to do
3197 (give-up)))))
3198
3199 ;; jacobi_sd(u,m) = jacobi_sn/jacobi_dn
3200 (defprop %jacobi_sd
3201 ((u m)
3202 ;; wrt u
3203 ((mtimes) ((%jacobi_cn) u m)
3204 ((mexpt) ((%jacobi_dn) u m) -2))
3205 ;; wrt m
3206 ((mplus)
3207 ((mtimes) ((mexpt) ((%jacobi_dn) u m) -1)
3208 ((mplus)
3209 ((mtimes) ((rat) 1 2)
3210 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3211 ((mexpt) ((%jacobi_cn) u m) 2)
3212 ((%jacobi_sn) u m))
3213 ((mtimes) ((rat) 1 2) ((mexpt) m -1)
3214 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
3215 ((mplus) u
3216 ((mtimes) -1
3217 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3218 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3219 m))))))
3220 ((mtimes) -1 ((mexpt) ((%jacobi_dn) u m) -2)
3221 ((%jacobi_sn) u m)
3222 ((mplus)
3223 ((mtimes) ((rat) -1 2)
3224 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3225 ((%jacobi_dn) u m)
3226 ((mexpt) ((%jacobi_sn) u m) 2))
3227 ((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
3228 ((%jacobi_sn) u m)
3229 ((mplus) u
3230 ((mtimes) -1
3231 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3232 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3233 m))))))))
3234 grad)
3235
3236 (def-simplifier jacobi_sd (u m)
3237 (let (coef args)
3238 (cond
3239 ((float-numerical-eval-p u m)
3240 (let ((fu (bigfloat:to ($float u)))
3241 (fm (bigfloat:to ($float m))))
3242 (to (bigfloat:/ (bigfloat::sn fu fm) (bigfloat::dn fu fm)))))
3243 ((setf args (complex-float-numerical-eval-p u m))
3244 (destructuring-bind (u m)
3245 args
3246 (let ((fu (bigfloat:to ($float u)))
3247 (fm (bigfloat:to ($float m))))
3248 (to (bigfloat:/ (bigfloat::sn fu fm) (bigfloat::dn fu fm))))))
3249 ((bigfloat-numerical-eval-p u m)
3250 (let ((uu (bigfloat:to ($bfloat u)))
3251 (mm (bigfloat:to ($bfloat m))))
3252 (to (bigfloat:/ (bigfloat::sn uu mm)
3253 (bigfloat::dn uu mm)))))
3254 ((setf args (complex-bigfloat-numerical-eval-p u m))
3255 (destructuring-bind (u m)
3256 args
3257 (let ((uu (bigfloat:to ($bfloat u)))
3258 (mm (bigfloat:to ($bfloat m))))
3259 (to (bigfloat:/ (bigfloat::sn uu mm)
3260 (bigfloat::dn uu mm))))))
3261 ((zerop1 u)
3262 0)
3263 ((zerop1 m)
3264 ;; A&S 16.6.5
3265 (ftake '%sin u))
3266 ((onep1 m)
3267 ;; A&S 16.6.5
3268 (ftake '%sinh u))
3269 ((and $trigsign (mminusp* u))
3270 ;; sd is odd
3271 (neg (ftake* '%jacobi_sd (neg u) m)))
3272 ((and $triginverses
3273 (listp u)
3274 (member (caar u) '(%inverse_jacobi_sn
3275 %inverse_jacobi_ns
3276 %inverse_jacobi_cn
3277 %inverse_jacobi_nc
3278 %inverse_jacobi_dn
3279 %inverse_jacobi_nd
3280 %inverse_jacobi_sc
3281 %inverse_jacobi_cs
3282 %inverse_jacobi_sd
3283 %inverse_jacobi_ds
3284 %inverse_jacobi_cd
3285 %inverse_jacobi_dc))
3286 (alike1 (third u) m))
3287 (cond ((eq (caar u) '%inverse_jacobi_sd)
3288 (second u))
3289 (t
3290 ;; Express in terms of sn and dn
3291 (div (ftake '%jacobi_sn u m)
3292 (ftake '%jacobi_dn u m)))))
3293 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3294 ((and $%iargs (multiplep u '$%i))
3295 ;; sd(i*u) = sn(i*u)/dn(i*u) = i*sc(u,m1)/dc(u,m1) = i*sd(u,m1)
3296 (mul '$%i
3297 (ftake* '%jacobi_sd (coeff u '$%i 1) (add 1 (neg m)))))
3298 ((setq coef (kc-arg2 u m))
3299 ;; A&S 16.8.5
3300 ;; sd(4*m*K+u) = sd(u)
3301 (destructuring-bind (lin const)
3302 coef
3303 (cond ((integerp lin)
3304 (ecase (mod lin 4)
3305 (0
3306 ;; sd(4*m*K+u) = sd(u)
3307 ;; sd(0) = 0
3308 (if (zerop1 const)
3309 0
3310 (ftake '%jacobi_sd const m)))
3311 (1
3312 ;; sd(4*m*K+K+u) = sd(K+u) = cn(u)/sqrt(1-m)
3313 ;; sd(K) = 1/sqrt(m1)
3314 (if (zerop1 const)
3315 (power (sub 1 m) 1//2)
3316 (div (ftake '%jacobi_cn const m)
3317 (power (sub 1 m) 1//2))))
3318 (2
3319 ;; sd(4*m*K+2*K+u) = sd(2*K+u) = -sd(u)
3320 ;; sd(2*K) = 0
3321 (if (zerop1 const)
3322 0
3323 (neg (ftake '%jacobi_sd const m))))
3324 (3
3325 ;; sd(4*m*K+3*K+u) = sd(3*K+u) = sd(2*K+K+u) =
3326 ;; -sd(K+u) = -cn(u)/sqrt(1-m)
3327 ;; sd(3*K) = -1/sqrt(m1)
3328 (if (zerop1 const)
3329 (neg (power (sub 1 m) -1//2))
3330 (neg (div (ftake '%jacobi_cn const m)
3331 (power (sub 1 m) 1//2)))))))
3332 ((and (alike1 lin 1//2)
3333 (zerop1 const))
3334 ;; jacobi_sn/jacobi_dn
3335 (div (ftake '%jacobi_sn
3336 (mul 1//2
3337 (ftake '%elliptic_kc m))
3338 m)
3339 (ftake '%jacobi_dn
3340 (mul 1//2
3341 (ftake '%elliptic_kc m))
3342 m)))
3343 (t
3344 (give-up)))))
3345 (t
3346 ;; Nothing to do
3347 (give-up)))))
3348
3349 ;; jacobi_cs(u,m) = jacobi_cn/jacobi_sn
3350 (defprop %jacobi_cs
3351 ((u m)
3352 ;; wrt u
3353 ((mtimes) -1 ((%jacobi_dn) u m)
3354 ((mexpt) ((%jacobi_sn) u m) -2))
3355 ;; wrt m
3356 ((mplus)
3357 ((mtimes) -1 ((%jacobi_cn) u m)
3358 ((mexpt) ((%jacobi_sn) u m) -2)
3359 ((mplus)
3360 ((mtimes) ((rat) 1 2)
3361 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3362 ((mexpt) ((%jacobi_cn) u m) 2)
3363 ((%jacobi_sn) u m))
3364 ((mtimes) ((rat) 1 2) ((mexpt) m -1)
3365 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
3366 ((mplus) u
3367 ((mtimes) -1
3368 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3369 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3370 m))))))
3371 ((mtimes) ((mexpt) ((%jacobi_sn) u m) -1)
3372 ((mplus)
3373 ((mtimes) ((rat) -1 2)
3374 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3375 ((%jacobi_cn) u m)
3376 ((mexpt) ((%jacobi_sn) u m) 2))
3377 ((mtimes) ((rat) -1 2) ((mexpt) m -1)
3378 ((%jacobi_dn) u m) ((%jacobi_sn) u m)
3379 ((mplus) u
3380 ((mtimes) -1
3381 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3382 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3383 m))))))))
3384 grad)
3385
3386 (def-simplifier jacobi_cs (u m)
3387 (let (coef args)
3388 (cond
3389 ((float-numerical-eval-p u m)
3390 (let ((fu (bigfloat:to ($float u)))
3391 (fm (bigfloat:to ($float m))))
3392 (to (bigfloat:/ (bigfloat::cn fu fm) (bigfloat::sn fu fm)))))
3393 ((setf args (complex-float-numerical-eval-p u m))
3394 (destructuring-bind (u m)
3395 args
3396 (let ((fu (bigfloat:to ($float u)))
3397 (fm (bigfloat:to ($float m))))
3398 (to (bigfloat:/ (bigfloat::cn fu fm) (bigfloat::sn fu fm))))))
3399 ((bigfloat-numerical-eval-p u m)
3400 (let ((uu (bigfloat:to ($bfloat u)))
3401 (mm (bigfloat:to ($bfloat m))))
3402 (to (bigfloat:/ (bigfloat::cn uu mm)
3403 (bigfloat::sn uu mm)))))
3404 ((setf args (complex-bigfloat-numerical-eval-p u m))
3405 (destructuring-bind (u m)
3406 args
3407 (let ((uu (bigfloat:to ($bfloat u)))
3408 (mm (bigfloat:to ($bfloat m))))
3409 (to (bigfloat:/ (bigfloat::cn uu mm)
3410 (bigfloat::sn uu mm))))))
3411 ((zerop1 m)
3412 ;; A&S 16.6.12
3413 (ftake '%cot u))
3414 ((onep1 m)
3415 ;; A&S 16.6.12
3416 (ftake '%csch u))
3417 ((zerop1 u)
3418 (dbz-err1 'jacobi_cs))
3419 ((and $trigsign (mminusp* u))
3420 ;; cs is odd
3421 (neg (ftake* '%jacobi_cs (neg u) m)))
3422 ((and $triginverses
3423 (listp u)
3424 (member (caar u) '(%inverse_jacobi_sn
3425 %inverse_jacobi_ns
3426 %inverse_jacobi_cn
3427 %inverse_jacobi_nc
3428 %inverse_jacobi_dn
3429 %inverse_jacobi_nd
3430 %inverse_jacobi_sc
3431 %inverse_jacobi_cs
3432 %inverse_jacobi_sd
3433 %inverse_jacobi_ds
3434 %inverse_jacobi_cd
3435 %inverse_jacobi_dc))
3436 (alike1 (third u) m))
3437 (cond ((eq (caar u) '%inverse_jacobi_cs)
3438 (second u))
3439 (t
3440 ;; Express in terms of cn an sn
3441 (div (ftake '%jacobi_cn u m)
3442 (ftake '%jacobi_sn u m)))))
3443 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3444 ((and $%iargs (multiplep u '$%i))
3445 ;; cs(i*u) = cn(i*u)/sn(i*u) = -i*nc(u,m1)/sc(u,m1) = -i*ns(u,m1)
3446 (neg (mul '$%i
3447 (ftake* '%jacobi_ns (coeff u '$%i 1) (add 1 (neg m))))))
3448 ((setq coef (kc-arg2 u m))
3449 ;; A&S 16.8.12
3450 ;;
3451 ;; cs(2*m*K + u) = cs(u)
3452 (destructuring-bind (lin const)
3453 coef
3454 (cond ((integerp lin)
3455 (ecase (mod lin 2)
3456 (0
3457 ;; cs(2*m*K + u) = cs(u)
3458 ;; cs(0) = infinity
3459 (if (zerop1 const)
3460 (dbz-err1 'jacobi_cs)
3461 (ftake '%jacobi_cs const m)))
3462 (1
3463 ;; cs(K+u) = -sqrt(1-m)*sc(u)
3464 ;; cs(K) = 0
3465 (if (zerop1 const)
3466 0
3467 (neg (mul (power (sub 1 m) 1//2)
3468 (ftake '%jacobi_sc const m)))))))
3469 ((and (alike1 lin 1//2)
3470 (zerop1 const))
3471 ;; 1/jacobi_sc
3472 (div 1
3473 (ftake '%jacobi_sc (mul 1//2
3474 (ftake '%elliptic_kc m))
3475 m)))
3476 (t
3477 (give-up)))))
3478 (t
3479 ;; Nothing to do
3480 (give-up)))))
3481
3482 ;; jacobi_cd(u,m) = jacobi_cn/jacobi_dn
3483 (defprop %jacobi_cd
3484 ((u m)
3485 ;; wrt u
3486 ((mtimes) ((mplus) -1 m)
3487 ((mexpt) ((%jacobi_dn) u m) -2)
3488 ((%jacobi_sn) u m))
3489 ;; wrt m
3490 ((mplus)
3491 ((mtimes) -1 ((%jacobi_cn) u m)
3492 ((mexpt) ((%jacobi_dn) u m) -2)
3493 ((mplus)
3494 ((mtimes) ((rat) -1 2)
3495 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3496 ((%jacobi_dn) u m)
3497 ((mexpt) ((%jacobi_sn) u m) 2))
3498 ((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
3499 ((%jacobi_sn) u m)
3500 ((mplus) u
3501 ((mtimes) -1
3502 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3503 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3504 m))))))
3505 ((mtimes) ((mexpt) ((%jacobi_dn) u m) -1)
3506 ((mplus)
3507 ((mtimes) ((rat) -1 2)
3508 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3509 ((%jacobi_cn) u m)
3510 ((mexpt) ((%jacobi_sn) u m) 2))
3511 ((mtimes) ((rat) -1 2) ((mexpt) m -1)
3512 ((%jacobi_dn) u m) ((%jacobi_sn) u m)
3513 ((mplus) u
3514 ((mtimes) -1
3515 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3516 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3517 m))))))))
3518 grad)
3519
3520 (def-simplifier jacobi_cd (u m)
3521 (let (coef args)
3522 (cond
3523 ((float-numerical-eval-p u m)
3524 (let ((fu (bigfloat:to ($float u)))
3525 (fm (bigfloat:to ($float m))))
3526 (to (bigfloat:/ (bigfloat::cn fu fm) (bigfloat::dn fu fm)))))
3527 ((setf args (complex-float-numerical-eval-p u m))
3528 (destructuring-bind (u m)
3529 args
3530 (let ((fu (bigfloat:to ($float u)))
3531 (fm (bigfloat:to ($float m))))
3532 (to (bigfloat:/ (bigfloat::cn fu fm) (bigfloat::dn fu fm))))))
3533 ((bigfloat-numerical-eval-p u m)
3534 (let ((uu (bigfloat:to ($bfloat u)))
3535 (mm (bigfloat:to ($bfloat m))))
3536 (to (bigfloat:/ (bigfloat::cn uu mm) (bigfloat::dn uu mm)))))
3537 ((setf args (complex-bigfloat-numerical-eval-p u m))
3538 (destructuring-bind (u m)
3539 args
3540 (let ((uu (bigfloat:to ($bfloat u)))
3541 (mm (bigfloat:to ($bfloat m))))
3542 (to (bigfloat:/ (bigfloat::cn uu mm) (bigfloat::dn uu mm))))))
3543 ((zerop1 u)
3544 1)
3545 ((zerop1 m)
3546 ;; A&S 16.6.4
3547 (ftake '%cos u))
3548 ((onep1 m)
3549 ;; A&S 16.6.4
3550 1)
3551 ((and $trigsign (mminusp* u))
3552 ;; cd is even
3553 (ftake* '%jacobi_cd (neg u) m))
3554 ((and $triginverses
3555 (listp u)
3556 (member (caar u) '(%inverse_jacobi_sn
3557 %inverse_jacobi_ns
3558 %inverse_jacobi_cn
3559 %inverse_jacobi_nc
3560 %inverse_jacobi_dn
3561 %inverse_jacobi_nd
3562 %inverse_jacobi_sc
3563 %inverse_jacobi_cs
3564 %inverse_jacobi_sd
3565 %inverse_jacobi_ds
3566 %inverse_jacobi_cd
3567 %inverse_jacobi_dc))
3568 (alike1 (third u) m))
3569 (cond ((eq (caar u) '%inverse_jacobi_cd)
3570 (second u))
3571 (t
3572 ;; Express in terms of cn and dn
3573 (div (ftake '%jacobi_cn u m)
3574 (ftake '%jacobi_dn u m)))))
3575 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3576 ((and $%iargs (multiplep u '$%i))
3577 ;; cd(i*u) = cn(i*u)/dn(i*u) = nc(u,m1)/dc(u,m1) = nd(u,m1)
3578 (ftake* '%jacobi_nd (coeff u '$%i 1) (add 1 (neg m))))
3579 ((setf coef (kc-arg2 u m))
3580 ;; A&S 16.8.4
3581 ;;
3582 (destructuring-bind (lin const)
3583 coef
3584 (cond ((integerp lin)
3585 (ecase (mod lin 4)
3586 (0
3587 ;; cd(4*m*K + u) = cd(u)
3588 ;; cd(0) = 1
3589 (if (zerop1 const)
3590 1
3591 (ftake '%jacobi_cd const m)))
3592 (1
3593 ;; cd(4*m*K + K + u) = cd(K+u) = -sn(u)
3594 ;; cd(K) = 0
3595 (if (zerop1 const)
3596 0
3597 (neg (ftake '%jacobi_sn const m))))
3598 (2
3599 ;; cd(4*m*K + 2*K + u) = cd(2*K+u) = -cd(u)
3600 ;; cd(2*K) = -1
3601 (if (zerop1 const)
3602 -1
3603 (neg (ftake '%jacobi_cd const m))))
3604 (3
3605 ;; cd(4*m*K + 3*K + u) = cd(2*K + K + u) =
3606 ;; -cd(K+u) = sn(u)
3607 ;; cd(3*K) = 0
3608 (if (zerop1 const)
3609 0
3610 (ftake '%jacobi_sn const m)))))
3611 ((and (alike1 lin 1//2)
3612 (zerop1 const))
3613 ;; jacobi_cn/jacobi_dn
3614 (div (ftake '%jacobi_cn
3615 (mul 1//2
3616 (ftake '%elliptic_kc m))
3617 m)
3618 (ftake '%jacobi_dn
3619 (mul 1//2
3620 (ftake '%elliptic_kc m))
3621 m)))
3622 (t
3623 ;; Nothing to do
3624 (give-up)))))
3625 (t
3626 ;; Nothing to do
3627 (give-up)))))
3628
3629 ;; jacobi_ds(u,m) = jacobi_dn/jacobi_sn
3630 (defprop %jacobi_ds
3631 ((u m)
3632 ;; wrt u
3633 ((mtimes) -1 ((%jacobi_cn) u m)
3634 ((mexpt) ((%jacobi_sn) u m) -2))
3635 ;; wrt m
3636 ((mplus)
3637 ((mtimes) -1 ((%jacobi_dn) u m)
3638 ((mexpt) ((%jacobi_sn) u m) -2)
3639 ((mplus)
3640 ((mtimes) ((rat) 1 2)
3641 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3642 ((mexpt) ((%jacobi_cn) u m) 2)
3643 ((%jacobi_sn) u m))
3644 ((mtimes) ((rat) 1 2) ((mexpt) m -1)
3645 ((%jacobi_cn) u m) ((%jacobi_dn) u m)
3646 ((mplus) u
3647 ((mtimes) -1
3648 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3649 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3650 m))))))
3651 ((mtimes) ((mexpt) ((%jacobi_sn) u m) -1)
3652 ((mplus)
3653 ((mtimes) ((rat) -1 2)
3654 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3655 ((%jacobi_dn) u m)
3656 ((mexpt) ((%jacobi_sn) u m) 2))
3657 ((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
3658 ((%jacobi_sn) u m)
3659 ((mplus) u
3660 ((mtimes) -1
3661 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3662 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3663 m))))))))
3664 grad)
3665
3666 (def-simplifier jacobi_ds (u m)
3667 (let (coef args)
3668 (cond
3669 ((float-numerical-eval-p u m)
3670 (let ((fu (bigfloat:to ($float u)))
3671 (fm (bigfloat:to ($float m))))
3672 (to (bigfloat:/ (bigfloat::dn fu fm) (bigfloat::sn fu fm)))))
3673 ((setf args (complex-float-numerical-eval-p u m))
3674 (destructuring-bind (u m)
3675 args
3676 (let ((fu (bigfloat:to ($float u)))
3677 (fm (bigfloat:to ($float m))))
3678 (to (bigfloat:/ (bigfloat::dn fu fm) (bigfloat::sn fu fm))))))
3679 ((bigfloat-numerical-eval-p u m)
3680 (let ((uu (bigfloat:to ($bfloat u)))
3681 (mm (bigfloat:to ($bfloat m))))
3682 (to (bigfloat:/ (bigfloat::dn uu mm)
3683 (bigfloat::sn uu mm)))))
3684 ((setf args (complex-bigfloat-numerical-eval-p u m))
3685 (destructuring-bind (u m)
3686 args
3687 (let ((uu (bigfloat:to ($bfloat u)))
3688 (mm (bigfloat:to ($bfloat m))))
3689 (to (bigfloat:/ (bigfloat::dn uu mm)
3690 (bigfloat::sn uu mm))))))
3691 ((zerop1 m)
3692 ;; A&S 16.6.11
3693 (ftake '%csc u))
3694 ((onep1 m)
3695 ;; A&S 16.6.11
3696 (ftake '%csch u))
3697 ((zerop1 u)
3698 (dbz-err1 'jacobi_ds))
3699 ((and $trigsign (mminusp* u))
3700 (neg (ftake* '%jacobi_ds (neg u) m)))
3701 ((and $triginverses
3702 (listp u)
3703 (member (caar u) '(%inverse_jacobi_sn
3704 %inverse_jacobi_ns
3705 %inverse_jacobi_cn
3706 %inverse_jacobi_nc
3707 %inverse_jacobi_dn
3708 %inverse_jacobi_nd
3709 %inverse_jacobi_sc
3710 %inverse_jacobi_cs
3711 %inverse_jacobi_sd
3712 %inverse_jacobi_ds
3713 %inverse_jacobi_cd
3714 %inverse_jacobi_dc))
3715 (alike1 (third u) m))
3716 (cond ((eq (caar u) '%inverse_jacobi_ds)
3717 (second u))
3718 (t
3719 ;; Express in terms of dn and sn
3720 (div (ftake '%jacobi_dn u m)
3721 (ftake '%jacobi_sn u m)))))
3722 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3723 ((and $%iargs (multiplep u '$%i))
3724 ;; ds(i*u) = dn(i*u)/sn(i*u) = -i*dc(u,m1)/sc(u,m1) = -i*ds(u,m1)
3725 (neg (mul '$%i
3726 (ftake* '%jacobi_ds (coeff u '$%i 1) (add 1 (neg m))))))
3727 ((setf coef (kc-arg2 u m))
3728 ;; A&S 16.8.11
3729 (destructuring-bind (lin const)
3730 coef
3731 (cond ((integerp lin)
3732 (ecase (mod lin 4)
3733 (0
3734 ;; ds(4*m*K + u) = ds(u)
3735 ;; ds(0) = infinity
3736 (if (zerop1 const)
3737 (dbz-err1 'jacobi_ds)
3738 (ftake '%jacobi_ds const m)))
3739 (1
3740 ;; ds(4*m*K + K + u) = ds(K+u) = sqrt(1-m)*nc(u)
3741 ;; ds(K) = sqrt(1-m)
3742 (if (zerop1 const)
3743 (power (sub 1 m) 1//2)
3744 (mul (power (sub 1 m) 1//2)
3745 (ftake '%jacobi_nc const m))))
3746 (2
3747 ;; ds(4*m*K + 2*K + u) = ds(2*K+u) = -ds(u)
3748 ;; ds(0) = pole
3749 (if (zerop1 const)
3750 (dbz-err1 'jacobi_ds)
3751 (neg (ftake '%jacobi_ds const m))))
3752 (3
3753 ;; ds(4*m*K + 3*K + u) = ds(2*K + K + u) =
3754 ;; -ds(K+u) = -sqrt(1-m)*nc(u)
3755 ;; ds(3*K) = -sqrt(1-m)
3756 (if (zerop1 const)
3757 (neg (power (sub 1 m) 1//2))
3758 (neg (mul (power (sub 1 m) 1//2)
3759 (ftake '%jacobi_nc u m)))))))
3760 ((and (alike1 lin 1//2)
3761 (zerop1 const))
3762 ;; jacobi_dn/jacobi_sn
3763 (div
3764 (ftake '%jacobi_dn
3765 (mul 1//2 (ftake '%elliptic_kc m))
3766 m)
3767 (ftake '%jacobi_sn
3768 (mul 1//2 (ftake '%elliptic_kc m))
3769 m)))
3770 (t
3771 ;; Nothing to do
3772 (give-up)))))
3773 (t
3774 ;; Nothing to do
3775 (give-up)))))
3776
3777 ;; jacobi_dc(u,m) = jacobi_dn/jacobi_cn
3778 (defprop %jacobi_dc
3779 ((u m)
3780 ;; wrt u
3781 ((mtimes) ((mplus) 1 ((mtimes) -1 m))
3782 ((mexpt) ((%jacobi_cn) u m) -2)
3783 ((%jacobi_sn) u m))
3784 ;; wrt m
3785 ((mplus)
3786 ((mtimes) ((mexpt) ((%jacobi_cn) u m) -1)
3787 ((mplus)
3788 ((mtimes) ((rat) -1 2)
3789 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3790 ((%jacobi_dn) u m)
3791 ((mexpt) ((%jacobi_sn) u m) 2))
3792 ((mtimes) ((rat) -1 2) ((%jacobi_cn) u m)
3793 ((%jacobi_sn) u m)
3794 ((mplus) u
3795 ((mtimes) -1
3796 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3797 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3798 m))))))
3799 ((mtimes) -1 ((mexpt) ((%jacobi_cn) u m) -2)
3800 ((%jacobi_dn) u m)
3801 ((mplus)
3802 ((mtimes) ((rat) -1 2)
3803 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3804 ((%jacobi_cn) u m)
3805 ((mexpt) ((%jacobi_sn) u m) 2))
3806 ((mtimes) ((rat) -1 2) ((mexpt) m -1)
3807 ((%jacobi_dn) u m) ((%jacobi_sn) u m)
3808 ((mplus) u
3809 ((mtimes) -1
3810 ((mexpt) ((mplus) 1 ((mtimes) -1 m)) -1)
3811 ((%elliptic_e) ((%asin) ((%jacobi_sn) u m))
3812 m))))))))
3813 grad)
3814
3815 (def-simplifier jacobi_dc (u m)
3816 (let (coef args)
3817 (cond
3818 ((float-numerical-eval-p u m)
3819 (let ((fu (bigfloat:to ($float u)))
3820 (fm (bigfloat:to ($float m))))
3821 (to (bigfloat:/ (bigfloat::dn fu fm) (bigfloat::cn fu fm)))))
3822 ((setf args (complex-float-numerical-eval-p u m))
3823 (destructuring-bind (u m)
3824 args
3825 (let ((fu (bigfloat:to ($float u)))
3826 (fm (bigfloat:to ($float m))))
3827 (to (bigfloat:/ (bigfloat::dn fu fm) (bigfloat::cn fu fm))))))
3828 ((bigfloat-numerical-eval-p u m)
3829 (let ((uu (bigfloat:to ($bfloat u)))
3830 (mm (bigfloat:to ($bfloat m))))
3831 (to (bigfloat:/ (bigfloat::dn uu mm)
3832 (bigfloat::cn uu mm)))))
3833 ((setf args (complex-bigfloat-numerical-eval-p u m))
3834 (destructuring-bind (u m)
3835 args
3836 (let ((uu (bigfloat:to ($bfloat u)))
3837 (mm (bigfloat:to ($bfloat m))))
3838 (to (bigfloat:/ (bigfloat::dn uu mm)
3839 (bigfloat::cn uu mm))))))
3840 ((zerop1 u)
3841 1)
3842 ((zerop1 m)
3843 ;; A&S 16.6.7
3844 (ftake '%sec u))
3845 ((onep1 m)
3846 ;; A&S 16.6.7
3847 1)
3848 ((and $trigsign (mminusp* u))
3849 (ftake* '%jacobi_dc (neg u) m))
3850 ((and $triginverses
3851 (listp u)
3852 (member (caar u) '(%inverse_jacobi_sn
3853 %inverse_jacobi_ns
3854 %inverse_jacobi_cn
3855 %inverse_jacobi_nc
3856 %inverse_jacobi_dn
3857 %inverse_jacobi_nd
3858 %inverse_jacobi_sc
3859 %inverse_jacobi_cs
3860 %inverse_jacobi_sd
3861 %inverse_jacobi_ds
3862 %inverse_jacobi_cd
3863 %inverse_jacobi_dc))
3864 (alike1 (third u) m))
3865 (cond ((eq (caar u) '%inverse_jacobi_dc)
3866 (second u))
3867 (t
3868 ;; Express in terms of dn and cn
3869 (div (ftake '%jacobi_dn u m)
3870 (ftake '%jacobi_cn u m)))))
3871 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3872 ((and $%iargs (multiplep u '$%i))
3873 ;; dc(i*u) = dn(i*u)/cn(i*u) = dc(u,m1)/nc(u,m1) = dn(u,m1)
3874 (ftake* '%jacobi_dn (coeff u '$%i 1) (add 1 (neg m))))
3875 ((setf coef (kc-arg2 u m))
3876 ;; See A&S 16.8.7
3877 (destructuring-bind (lin const)
3878 coef
3879 (cond ((integerp lin)
3880 (ecase (mod lin 4)
3881 (0
3882 ;; dc(4*m*K + u) = dc(u)
3883 ;; dc(0) = 1
3884 (if (zerop1 const)
3885 1
3886 (ftake '%jacobi_dc const m)))
3887 (1
3888 ;; dc(4*m*K + K + u) = dc(K+u) = -ns(u)
3889 ;; dc(K) = pole
3890 (if (zerop1 const)
3891 (dbz-err1 'jacobi_dc)
3892 (neg (ftake '%jacobi_ns const m))))
3893 (2
3894 ;; dc(4*m*K + 2*K + u) = dc(2*K+u) = -dc(u)
3895 ;; dc(2K) = -1
3896 (if (zerop1 const)
3897 -1
3898 (neg (ftake '%jacobi_dc const m))))
3899 (3
3900 ;; dc(4*m*K + 3*K + u) = dc(2*K + K + u) =
3901 ;; -dc(K+u) = ns(u)
3902 ;; dc(3*K) = ns(0) = inf
3903 (if (zerop1 const)
3904 (dbz-err1 'jacobi_dc)
3905 (ftake '%jacobi_dc const m)))))
3906 ((and (alike1 lin 1//2)
3907 (zerop1 const))
3908 ;; jacobi_dn/jacobi_cn
3909 (div
3910 (ftake '%jacobi_dn
3911 (mul 1//2 (ftake '%elliptic_kc m))
3912 m)
3913 (ftake '%jacobi_cn
3914 (mul 1//2 (ftake '%elliptic_kc m))
3915 m)))
3916 (t
3917 ;; Nothing to do
3918 (give-up)))))
3919 (t
3920 ;; Nothing to do
3921 (give-up)))))
3922
3923 ;;; Other inverse Jacobian functions
3924
3925 ;; inverse_jacobi_ns(x)
3926 ;;
3927 ;; Let u = inverse_jacobi_ns(x). Then jacobi_ns(u) = x or
3928 ;; 1/jacobi_sn(u) = x or
3929 ;;
3930 ;; jacobi_sn(u) = 1/x
3931 ;;
3932 ;; so u = inverse_jacobi_sn(1/x)
3933 (defprop %inverse_jacobi_ns
3934 ((x m)
3935 ;; Whittaker and Watson, example in 22.122
3936 ;; inverse_jacobi_ns(u,m) = integrate(1/sqrt(t^2-1)/sqrt(t^2-m), t, u, inf)
3937 ;; -> -1/sqrt(x^2-1)/sqrt(x^2-m)
3938 ((mtimes) -1
3939 ((mexpt) ((mplus) -1 ((mexpt) x 2)) ((rat) -1 2))
3940 ((mexpt)
3941 ((mplus) ((mtimes simp ratsimp) -1 m) ((mexpt) x 2))
3942 ((rat) -1 2)))
3943 ;; wrt m
3944 ; ((%derivative) ((%inverse_jacobi_ns) x m) m 1)
3945 nil)
3946 grad)
3947
3948 (def-simplifier inverse_jacobi_ns (u m)
3949 (let (args)
3950 (cond
3951 ((float-numerical-eval-p u m)
3952 ;; Numerically evaluate asn
3953 ;;
3954 ;; ans(x,m) = asn(1/x,m) = F(asin(1/x),m)
3955 (to (elliptic-f (cl:asin (/ ($float u))) ($float m))))
3956 ((complex-float-numerical-eval-p u m)
3957 (to (elliptic-f (cl:asin (/ (complex ($realpart ($float u)) ($imagpart ($float u)))))
3958 (complex ($realpart ($float m)) ($imagpart ($float m))))))
3959 ((bigfloat-numerical-eval-p u m)
3960 (to (bigfloat::bf-elliptic-f (bigfloat:asin (bigfloat:/ (bigfloat:to ($bfloat u))))
3961 (bigfloat:to ($bfloat m)))))
3962 ((setf args (complex-bigfloat-numerical-eval-p u m))
3963 (destructuring-bind (u m)
3964 args
3965 (to (bigfloat::bf-elliptic-f (bigfloat:asin (bigfloat:/ (bigfloat:to ($bfloat u))))
3966 (bigfloat:to ($bfloat m))))))
3967 ((zerop1 m)
3968 ;; ans(x,0) = F(asin(1/x),0) = asin(1/x)
3969 (ftake '%elliptic_f (ftake '%asin (div 1 u)) 0))
3970 ((onep1 m)
3971 ;; ans(x,1) = F(asin(1/x),1) = log(tan(pi/2+asin(1/x)/2))
3972 (ftake '%elliptic_f (ftake '%asin (div 1 u)) 1))
3973 ((onep1 u)
3974 (ftake '%elliptic_kc m))
3975 ((alike1 u -1)
3976 (neg (ftake '%elliptic_kc m)))
3977 ((and (eq $triginverses '$all)
3978 (listp u)
3979 (eq (caar u) '%jacobi_ns)
3980 (alike1 (third u) m))
3981 ;; inverse_jacobi_ns(ns(u)) = u
3982 (second u))
3983 (t
3984 ;; Nothing to do
3985 (give-up)))))
3986
3987 ;; inverse_jacobi_nc(x)
3988 ;;
3989 ;; Let u = inverse_jacobi_nc(x). Then jacobi_nc(u) = x or
3990 ;; 1/jacobi_cn(u) = x or
3991 ;;
3992 ;; jacobi_cn(u) = 1/x
3993 ;;
3994 ;; so u = inverse_jacobi_cn(1/x)
3995 (defprop %inverse_jacobi_nc
3996 ((x m)
3997 ;; Whittaker and Watson, example in 22.122
3998 ;; inverse_jacobi_nc(u,m) = integrate(1/sqrt(t^2-1)/sqrt((1-m)*t^2+m), t, 1, u)
3999 ;; -> 1/sqrt(x^2-1)/sqrt((1-m)*x^2+m)
4000 ((mtimes)
4001 ((mexpt) ((mplus) -1 ((mexpt) x 2)) ((rat) -1 2))
4002 ((mexpt)
4003 ((mplus) m
4004 ((mtimes) -1 ((mplus) -1 m) ((mexpt) x 2)))
4005 ((rat) -1 2)))
4006 ;; wrt m
4007 ; ((%derivative) ((%inverse_jacobi_nc) x m) m 1)
4008 nil)
4009 grad)
4010
4011 (def-simplifier inverse_jacobi_nc (u m)
4012 (cond ((or (float-numerical-eval-p u m)
4013 (complex-float-numerical-eval-p u m)
4014 (bigfloat-numerical-eval-p u m)
4015 (complex-bigfloat-numerical-eval-p u m))
4016 ;;
4017 (ftake '%inverse_jacobi_cn ($rectform (div 1 u)) m))
4018 ((onep1 u)
4019 0)
4020 ((alike1 u -1)
4021 (mul 2 (ftake '%elliptic_kc m)))
4022 ((and (eq $triginverses '$all)
4023 (listp u)
4024 (eq (caar u) '%jacobi_nc)
4025 (alike1 (third u) m))
4026 ;; inverse_jacobi_nc(nc(u)) = u
4027 (second u))
4028 (t
4029 ;; Nothing to do
4030 (give-up))))
4031
4032 ;; inverse_jacobi_nd(x)
4033 ;;
4034 ;; Let u = inverse_jacobi_nd(x). Then jacobi_nd(u) = x or
4035 ;; 1/jacobi_dn(u) = x or
4036 ;;
4037 ;; jacobi_dn(u) = 1/x
4038 ;;
4039 ;; so u = inverse_jacobi_dn(1/x)
4040 (defprop %inverse_jacobi_nd
4041 ((x m)
4042 ;; Whittaker and Watson, example in 22.122
4043 ;; inverse_jacobi_nd(u,m) = integrate(1/sqrt(t^2-1)/sqrt(1-(1-m)*t^2), t, 1, u)
4044 ;; -> 1/sqrt(u^2-1)/sqrt(1-(1-m)*t^2)
4045 ((mtimes)
4046 ((mexpt) ((mplus) -1 ((mexpt simp ratsimp) x 2))
4047 ((rat) -1 2))
4048 ((mexpt)
4049 ((mplus) 1
4050 ((mtimes) ((mplus) -1 m) ((mexpt simp ratsimp) x 2)))
4051 ((rat) -1 2)))
4052 ;; wrt m
4053 ; ((%derivative) ((%inverse_jacobi_nd) x m) m 1)
4054 nil)
4055 grad)
4056
4057 (def-simplifier inverse_jacobi_nd (u m)
4058 (cond ((or (float-numerical-eval-p u m)
4059 (complex-float-numerical-eval-p u m)
4060 (bigfloat-numerical-eval-p u m)
4061 (complex-bigfloat-numerical-eval-p u m))
4062 (ftake '%inverse_jacobi_dn ($rectform (div 1 u)) m))
4063 ((onep1 u)
4064 0)
4065 ((onep1 ($ratsimp (mul (power (sub 1 m) 1//2) u)))
4066 ;; jacobi_nd(1/sqrt(1-m),m) = K(m). This follows from
4067 ;; jacobi_dn(sqrt(1-m),m) = K(m).
4068 (ftake '%elliptic_kc m))
4069 ((and (eq $triginverses '$all)
4070 (listp u)
4071 (eq (caar u) '%jacobi_nd)
4072 (alike1 (third u) m))
4073 ;; inverse_jacobi_nd(nd(u)) = u
4074 (second u))
4075 (t
4076 ;; Nothing to do
4077 (give-up))))
4078
4079 ;; inverse_jacobi_sc(x)
4080 ;;
4081 ;; Let u = inverse_jacobi_sc(x). Then jacobi_sc(u) = x or
4082 ;; x = jacobi_sn(u)/jacobi_cn(u)
4083 ;;
4084 ;; x^2 = sn^2/cn^2
4085 ;; = sn^2/(1-sn^2)
4086 ;;
4087 ;; so
4088 ;;
4089 ;; sn^2 = x^2/(1+x^2)
4090 ;;
4091 ;; sn(u) = x/sqrt(1+x^2)
4092 ;;
4093 ;; u = inverse_sn(x/sqrt(1+x^2))
4094 ;;
4095 (defprop %inverse_jacobi_sc
4096 ((x m)
4097 ;; Whittaker and Watson, example in 22.122
4098 ;; inverse_jacobi_sc(u,m) = integrate(1/sqrt(1+t^2)/sqrt(1+(1-m)*t^2), t, 0, u)
4099 ;; -> 1/sqrt(1+x^2)/sqrt(1+(1-m)*x^2)
4100 ((mtimes)
4101 ((mexpt) ((mplus) 1 ((mexpt) x 2))
4102 ((rat) -1 2))
4103 ((mexpt)
4104 ((mplus) 1
4105 ((mtimes) -1 ((mplus) -1 m) ((mexpt) x 2)))
4106 ((rat) -1 2)))
4107 ;; wrt m
4108 ; ((%derivative) ((%inverse_jacobi_sc) x m) m 1)
4109 nil)
4110 grad)
4111
4112 (def-simplifier inverse_jacobi_sc (u m)
4113 (cond ((or (float-numerical-eval-p u m)
4114 (complex-float-numerical-eval-p u m)
4115 (bigfloat-numerical-eval-p u m)
4116 (complex-bigfloat-numerical-eval-p u m))
4117 (ftake '%inverse_jacobi_sn
4118 ($rectform (div u (power (add 1 (mul u u)) 1//2)))
4119 m))
4120 ((zerop1 u)
4121 ;; jacobi_sc(0,m) = 0
4122 0)
4123 ((and (eq $triginverses '$all)
4124 (listp u)
4125 (eq (caar u) '%jacobi_sc)
4126 (alike1 (third u) m))
4127 ;; inverse_jacobi_sc(sc(u)) = u
4128 (second u))
4129 (t
4130 ;; Nothing to do
4131 (give-up))))
4132
4133 ;; inverse_jacobi_sd(x)
4134 ;;
4135 ;; Let u = inverse_jacobi_sd(x). Then jacobi_sd(u) = x or
4136 ;; x = jacobi_sn(u)/jacobi_dn(u)
4137 ;;
4138 ;; x^2 = sn^2/dn^2
4139 ;; = sn^2/(1-m*sn^2)
4140 ;;
4141 ;; so
4142 ;;
4143 ;; sn^2 = x^2/(1+m*x^2)
4144 ;;
4145 ;; sn(u) = x/sqrt(1+m*x^2)
4146 ;;
4147 ;; u = inverse_sn(x/sqrt(1+m*x^2))
4148 ;;
4149 (defprop %inverse_jacobi_sd
4150 ((x m)
4151 ;; Whittaker and Watson, example in 22.122
4152 ;; inverse_jacobi_sd(u,m) = integrate(1/sqrt(1-(1-m)*t^2)/sqrt(1+m*t^2), t, 0, u)
4153 ;; -> 1/sqrt(1-(1-m)*x^2)/sqrt(1+m*x^2)
4154 ((mtimes)
4155 ((mexpt)
4156 ((mplus) 1 ((mtimes) ((mplus) -1 m) ((mexpt) x 2)))
4157 ((rat) -1 2))
4158 ((mexpt) ((mplus) 1 ((mtimes) m ((mexpt) x 2)))
4159 ((rat) -1 2)))
4160 ;; wrt m
4161 ; ((%derivative) ((%inverse_jacobi_sd) x m) m 1)
4162 nil)
4163 grad)
4164
4165 (def-simplifier inverse_jacobi_sd (u m)
4166 (cond ((or (float-numerical-eval-p u m)
4167 (complex-float-numerical-eval-p u m)
4168 (bigfloat-numerical-eval-p u m)
4169 (complex-bigfloat-numerical-eval-p u m))
4170 (ftake '%inverse_jacobi_sn
4171 ($rectform (div u (power (add 1 (mul m (mul u u))) 1//2)))
4172 m))
4173 ((zerop1 u)
4174 0)
4175 ((eql 0 ($ratsimp (sub u (div 1 (power (sub 1 m) 1//2)))))
4176 ;; inverse_jacobi_sd(1/sqrt(1-m), m) = elliptic_kc(m)
4177 ;;
4178 ;; We can see this from inverse_jacobi_sd(x,m) =
4179 ;; inverse_jacobi_sn(x/sqrt(1+m*x^2), m). So
4180 ;; inverse_jacobi_sd(1/sqrt(1-m),m) = inverse_jacobi_sn(1,m)
4181 (ftake '%elliptic_kc m))
4182 ((and (eq $triginverses '$all)
4183 (listp u)
4184 (eq (caar u) '%jacobi_sd)
4185 (alike1 (third u) m))
4186 ;; inverse_jacobi_sd(sd(u)) = u
4187 (second u))
4188 (t
4189 ;; Nothing to do
4190 (give-up))))
4191
4192 ;; inverse_jacobi_cs(x)
4193 ;;
4194 ;; Let u = inverse_jacobi_cs(x). Then jacobi_cs(u) = x or
4195 ;; 1/x = 1/jacobi_cs(u) = jacobi_sc(u)
4196 ;;
4197 ;; u = inverse_sc(1/x)
4198 ;;
4199 (defprop %inverse_jacobi_cs
4200 ((x m)
4201 ;; Whittaker and Watson, example in 22.122
4202 ;; inverse_jacobi_cs(u,m) = integrate(1/sqrt(t^2+1)/sqrt(t^2+(1-m)), t, u, inf)
4203 ;; -> -1/sqrt(x^2+1)/sqrt(x^2+(1-m))
4204 ((mtimes) -1
4205 ((mexpt) ((mplus) 1 ((mexpt simp ratsimp) x 2))
4206 ((rat) -1 2))
4207 ((mexpt) ((mplus) 1
4208 ((mtimes simp ratsimp) -1 m)
4209 ((mexpt simp ratsimp) x 2))
4210 ((rat) -1 2)))
4211 ;; wrt m
4212 ; ((%derivative) ((%inverse_jacobi_cs) x m) m 1)
4213 nil)
4214 grad)
4215
4216 (def-simplifier inverse_jacobi_cs (u m)
4217 (cond ((or (float-numerical-eval-p u m)
4218 (complex-float-numerical-eval-p u m)
4219 (bigfloat-numerical-eval-p u m)
4220 (complex-bigfloat-numerical-eval-p u m))
4221 (ftake '%inverse_jacobi_sc ($rectform (div 1 u)) m))
4222 ((zerop1 u)
4223 (ftake '%elliptic_kc m))
4224 (t
4225 ;; Nothing to do
4226 (give-up))))
4227
4228 ;; inverse_jacobi_cd(x)
4229 ;;
4230 ;; Let u = inverse_jacobi_cd(x). Then jacobi_cd(u) = x or
4231 ;; x = jacobi_cn(u)/jacobi_dn(u)
4232 ;;
4233 ;; x^2 = cn^2/dn^2
4234 ;; = (1-sn^2)/(1-m*sn^2)
4235 ;;
4236 ;; or
4237 ;;
4238 ;; sn^2 = (1-x^2)/(1-m*x^2)
4239 ;;
4240 ;; sn(u) = sqrt(1-x^2)/sqrt(1-m*x^2)
4241 ;;
4242 ;; u = inverse_sn(sqrt(1-x^2)/sqrt(1-m*x^2))
4243 ;;
4244 (defprop %inverse_jacobi_cd
4245 ((x m)
4246 ;; Whittaker and Watson, example in 22.122
4247 ;; inverse_jacobi_cd(u,m) = integrate(1/sqrt(1-t^2)/sqrt(1-m*t^2), t, u, 1)
4248 ;; -> -1/sqrt(1-x^2)/sqrt(1-m*x^2)
4249 ((mtimes) -1
4250 ((mexpt)
4251 ((mplus) 1 ((mtimes) -1 ((mexpt) x 2)))
4252 ((rat) -1 2))
4253 ((mexpt)
4254 ((mplus) 1 ((mtimes) -1 m ((mexpt) x 2)))
4255 ((rat) -1 2)))
4256 ;; wrt m
4257 ; ((%derivative) ((%inverse_jacobi_cd) x m) m 1)
4258 nil)
4259 grad)
4260
4261 (def-simplifier inverse_jacobi_cd (u m)
4262 (cond ((or (complex-float-numerical-eval-p u m)
4263 (complex-bigfloat-numerical-eval-p u m))
4264 (let (($numer t))
4265 (ftake '%inverse_jacobi_sn
4266 ($rectform (div (power (mul (sub 1 u) (add 1 u)) 1//2)
4267 (power (sub 1 (mul m (mul u u))) 1//2)))
4268 m)))
4269 ((onep1 u)
4270 0)
4271 ((zerop1 u)
4272 (ftake '%elliptic_kc m))
4273 ((and (eq $triginverses '$all)
4274 (listp u)
4275 (eq (caar u) '%jacobi_cd)
4276 (alike1 (third u) m))
4277 ;; inverse_jacobi_cd(cd(u)) = u
4278 (second u))
4279 (t
4280 ;; Nothing to do
4281 (give-up))))
4282
4283 ;; inverse_jacobi_ds(x)
4284 ;;
4285 ;; Let u = inverse_jacobi_ds(x). Then jacobi_ds(u) = x or
4286 ;; 1/x = 1/jacobi_ds(u) = jacobi_sd(u)
4287 ;;
4288 ;; u = inverse_sd(1/x)
4289 ;;
4290 (defprop %inverse_jacobi_ds
4291 ((x m)
4292 ;; Whittaker and Watson, example in 22.122
4293 ;; inverse_jacobi_ds(u,m) = integrate(1/sqrt(t^2-(1-m))/sqrt(t^2+m), t, u, inf)
4294 ;; -> -1/sqrt(x^2-(1-m))/sqrt(x^2+m)
4295 ((mtimes) -1
4296 ((mexpt)
4297 ((mplus) -1 m ((mexpt simp ratsimp) x 2))
4298 ((rat) -1 2))
4299 ((mexpt)
4300 ((mplus) m ((mexpt simp ratsimp) x 2))
4301 ((rat) -1 2)))
4302 ;; wrt m
4303 ; ((%derivative) ((%inverse_jacobi_ds) x m) m 1)
4304 nil)
4305 grad)
4306
4307 (def-simplifier inverse_jacobi_ds (u m)
4308 (cond ((or (float-numerical-eval-p u m)
4309 (complex-float-numerical-eval-p u m)
4310 (bigfloat-numerical-eval-p u m)
4311 (complex-bigfloat-numerical-eval-p u m))
4312 (ftake '%inverse_jacobi_sd ($rectform (div 1 u)) m))
4313 ((and $trigsign (mminusp* u))
4314 (neg (ftake* '%inverse_jacobi_ds (neg u) m)))
4315 ((eql 0 ($ratsimp (sub u (power (sub 1 m) 1//2))))
4316 ;; inverse_jacobi_ds(sqrt(1-m),m) = elliptic_kc(m)
4317 ;;
4318 ;; Since inverse_jacobi_ds(sqrt(1-m), m) =
4319 ;; inverse_jacobi_sd(1/sqrt(1-m),m). And we know from
4320 ;; above that this is elliptic_kc(m)
4321 (ftake '%elliptic_kc m))
4322 ((and (eq $triginverses '$all)
4323 (listp u)
4324 (eq (caar u) '%jacobi_ds)
4325 (alike1 (third u) m))
4326 ;; inverse_jacobi_ds(ds(u)) = u
4327 (second u))
4328 (t
4329 ;; Nothing to do
4330 (give-up))))
4331
4332
4333 ;; inverse_jacobi_dc(x)
4334 ;;
4335 ;; Let u = inverse_jacobi_dc(x). Then jacobi_dc(u) = x or
4336 ;; 1/x = 1/jacobi_dc(u) = jacobi_cd(u)
4337 ;;
4338 ;; u = inverse_cd(1/x)
4339 ;;
4340 (defprop %inverse_jacobi_dc
4341 ((x m)
4342 ;; Note: Whittaker and Watson, example in 22.122 says
4343 ;; inverse_jacobi_dc(u,m) = integrate(1/sqrt(t^2-1)/sqrt(t^2-m),
4344 ;; t, u, 1) but that seems wrong. A&S 17.4.47 says
4345 ;; integrate(1/sqrt(t^2-1)/sqrt(t^2-m), t, a, u) =
4346 ;; inverse_jacobi_cd(x,m). Lawden 3.2.8 says the same.
4347 ;; functions.wolfram.com says the derivative is
4348 ;; 1/sqrt(t^2-1)/sqrt(t^2-m).
4349 ((mtimes)
4350 ((mexpt)
4351 ((mplus) -1 ((mexpt simp ratsimp) x 2))
4352 ((rat) -1 2))
4353 ((mexpt)
4354 ((mplus)
4355 ((mtimes simp ratsimp) -1 m)
4356 ((mexpt simp ratsimp) x 2))
4357 ((rat) -1 2)))
4358 ;; wrt m
4359 ; ((%derivative) ((%inverse_jacobi_dc) x m) m 1)
4360 nil)
4361 grad)
4362
4363 (def-simplifier inverse_jacobi_dc (u m)
4364 (cond ((or (complex-float-numerical-eval-p u m)
4365 (complex-bigfloat-numerical-eval-p u m))
4366 (ftake '%inverse_jacobi_cd ($rectform (div 1 u)) m))
4367 ((onep1 u)
4368 0)
4369 ((and (eq $triginverses '$all)
4370 (listp u)
4371 (eq (caar u) '%jacobi_dc)
4372 (alike1 (third u) m))
4373 ;; inverse_jacobi_dc(dc(u)) = u
4374 (second u))
4375 (t
4376 ;; Nothing to do
4377 (give-up))))
4378
4379 ;; Convert an inverse Jacobian function into the equivalent elliptic
4380 ;; integral F.
4381 ;;
4382 ;; See A&S 17.4.41-17.4.52.
4383 (defun make-elliptic-f (e)
4384 (cond ((atom e)
4385 e)
4386 ((member (caar e) '(%inverse_jacobi_sc %inverse_jacobi_cs
4387 %inverse_jacobi_nd %inverse_jacobi_dn
4388 %inverse_jacobi_sn %inverse_jacobi_cd
4389 %inverse_jacobi_dc %inverse_jacobi_ns
4390 %inverse_jacobi_nc %inverse_jacobi_ds
4391 %inverse_jacobi_sd %inverse_jacobi_cn))
4392 ;; We have some inverse Jacobi function. Convert it to the F form.
4393 (destructuring-bind ((fn &rest ops) u m)
4394 e
4395 (declare (ignore ops))
4396 (ecase fn
4397 (%inverse_jacobi_sc
4398 ;; A&S 17.4.41
4399 (ftake '%elliptic_f (ftake '%atan u) m))
4400 (%inverse_jacobi_cs
4401 ;; A&S 17.4.42
4402 (ftake '%elliptic_f (ftake '%atan (div 1 u)) m))
4403 (%inverse_jacobi_nd
4404 ;; A&S 17.4.43
4405 (ftake '%elliptic_f
4406 (ftake '%asin
4407 (mul (power m -1//2)
4408 (div 1 u)
4409 (power (add -1 (mul u u))
4410 1//2)))
4411 m))
4412 (%inverse_jacobi_dn
4413 ;; A&S 17.4.44
4414 (ftake '%elliptic_f
4415 (ftake '%asin (mul
4416 (power m -1//2)
4417 (power (sub 1 (power u 2)) 1//2)))
4418 m))
4419 (%inverse_jacobi_sn
4420 ;; A&S 17.4.45
4421 (ftake '%elliptic_f (ftake '%asin u) m))
4422 (%inverse_jacobi_cd
4423 ;; A&S 17.4.46
4424 (ftake '%elliptic_f
4425 (ftake '%asin
4426 (power (mul (sub 1 (mul u u))
4427 (sub 1 (mul m u u)))
4428 1//2))
4429 m))
4430 (%inverse_jacobi_dc
4431 ;; A&S 17.4.47
4432 (ftake '%elliptic_f
4433 (ftake '%asin
4434 (power (mul (sub (mul u u) 1)
4435 (sub (mul u u) m))
4436 1//2))
4437 m))
4438 (%inverse_jacobi_ns
4439 ;; A&S 17.4.48
4440 (ftake '%elliptic_f (ftake '%asin (div 1 u)) m))
4441 (%inverse_jacobi_nc
4442 ;; A&S 17.4.49
4443 (ftake '%elliptic_f (ftake '%acos (div 1 u)) m))
4444 (%inverse_jacobi_ds
4445 ;; A&S 17.4.50
4446 (ftake '%elliptic_f
4447 (ftake '%asin
4448 (power (add m (mul u u))
4449 1//2))
4450 m))
4451 (%inverse_jacobi_sd
4452 ;; A&S 17.4.51
4453 (ftake '%elliptic_f
4454 (ftake '%asin
4455 (div u
4456 (power (add 1 (mul m u u))
4457 1//2)))
4458 m))
4459 (%inverse_jacobi_cn
4460 ;; A&S 17.4.52
4461 (ftake '%elliptic_f (ftake '%acos u) m)))))
4462 (t
4463 (recur-apply #'make-elliptic-f e))))
4464
4465 (defmfun $make_elliptic_f (e)
4466 (if (atom e)
4467 e
4468 (simplify (make-elliptic-f e))))
4469
4470 (defun make-elliptic-e (e)
4471 (cond ((atom e) e)
4472 ((eq (caar e) '$elliptic_eu)
4473 (destructuring-bind ((ffun &rest ops) u m) e
4474 (declare (ignore ffun ops))
4475 (ftake '%elliptic_e (ftake '%asin (ftake '%jacobi_sn u m)) m)))
4476 (t
4477 (recur-apply #'make-elliptic-e e))))
4478
4479 (defmfun $make_elliptic_e (e)
4480 (if (atom e)
4481 e
4482 (simplify (make-elliptic-e e))))
4483
4484
4485 ;; Eu(u,m) = integrate(jacobi_dn(v,m)^2,v,0,u)
4486 ;; = integrate(sqrt((1-m*t^2)/(1-t^2)),t,0,jacobi_sn(u,m))
4487 ;;
4488 ;; Eu(u,m) = E(am(u),m)
4489 ;;
4490 ;; where E(u,m) is elliptic-e above.
4491 ;;
4492 ;; Checks.
4493 ;; Lawden gives the following relationships
4494 ;;
4495 ;; E(u+v) = E(u) + E(v) - m*sn(u)*sn(v)*sn(u+v)
4496 ;; E(u,0) = u, E(u,1) = tanh u
4497 ;;
4498 ;; E(i*u,k) = i*sc(u,k')*dn(u,k') - i*E(u,k') + i*u
4499 ;;
4500 ;; E(2*i*K') = 2*i*(K'-E')
4501 ;;
4502 ;; E(u + 2*i*K') = E(u) + 2*i*(K' - E')
4503 ;;
4504 ;; E(u+K) = E(u) + E - k^2*sn(u)*cd(u)
4505 (defun elliptic-eu (u m)
4506 (cond ((realp u)
4507 ;; E(u + 2*n*K) = E(u) + 2*n*E
4508 (let ((ell-k (to (elliptic-k m)))
4509 (ell-e (elliptic-ec m)))
4510 (multiple-value-bind (n u-rem)
4511 (floor u (* 2 ell-k))
4512 ;; 0 <= u-rem < 2*K
4513 (+ (* 2 n ell-e)
4514 (cond ((>= u-rem ell-k)
4515 ;; 0 <= u-rem < K so
4516 ;; E(u + K) = E(u) + E - m*sn(u)*cd(u)
4517 (let ((u-k (- u ell-k)))
4518 (- (+ (elliptic-e (cl:asin (bigfloat::sn u-k m)) m)
4519 ell-e)
4520 (/ (* m (bigfloat::sn u-k m) (bigfloat::cn u-k m))
4521 (bigfloat::dn u-k m)))))
4522 (t
4523 (elliptic-e (cl:asin (bigfloat::sn u m)) m)))))))
4524 ((complexp u)
4525 ;; From Lawden:
4526 ;;
4527 ;; E(u+i*v, m) = E(u,m) -i*E(v,m') + i*v + i*sc(v,m')*dn(v,m')
4528 ;; -i*m*sn(u,m)*sc(v,m')*sn(u+i*v,m)
4529 ;;
4530 (let ((u-r (realpart u))
4531 (u-i (imagpart u))
4532 (m1 (- 1 m)))
4533 (+ (elliptic-eu u-r m)
4534 (* #c(0 1)
4535 (- (+ u-i
4536 (/ (* (bigfloat::sn u-i m1) (bigfloat::dn u-i m1))
4537 (bigfloat::cn u-i m1)))
4538 (+ (elliptic-eu u-i m1)
4539 (/ (* m (bigfloat::sn u-r m) (bigfloat::sn u-i m1) (bigfloat::sn u m))
4540 (bigfloat::cn u-i m1))))))))))
4541
4542 (defprop $elliptic_eu
4543 ((u m)
4544 ((mexpt) ((%jacobi_dn) u m) 2)
4545 ;; wrt m
4546 )
4547 grad)
4548
4549 (def-simplifier elliptic_eu (u m)
4550 (cond
4551 ;; as it stands, ELLIPTIC-EU can't handle bigfloats or complex bigfloats,
4552 ;; so handle only floats and complex floats here.
4553 ((float-numerical-eval-p u m)
4554 (elliptic-eu ($float u) ($float m)))
4555 ((complex-float-numerical-eval-p u m)
4556 (let ((u-r ($realpart u))
4557 (u-i ($imagpart u))
4558 (m ($float m)))
4559 (complexify (elliptic-eu (complex u-r u-i) m))))
4560 (t
4561 (give-up))))
4562
4563 #+nil
4564 (def-simplifier jacobi_am (u m)
4565 (cond
4566 ;; as it stands, BIGFLOAT::SN can't handle bigfloats or complex bigfloats,
4567 ;; so handle only floats and complex floats here.
4568 ((float-numerical-eval-p u m)
4569 (cl:asin (bigfloat::sn ($float u) ($float m))))
4570 ((complex-float-numerical-eval-p u m)
4571 (let ((u-r ($realpart ($float u)))
4572 (u-i ($imagpart ($float u)))
4573 (m ($float m)))
4574 (complexify (cl:asin (bigfloat::sn (complex u-r u-i) m)))))
4575 (t
4576 ;; Nothing to do
4577 (give-up))))
4578
4579 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4580 ;; Integrals. At present with respect to first argument only.
4581 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4582
4583 ;; A&S 16.24.1: integrate(jacobi_sn(u,m),u)
4584 ;; = log(jacobi_dn(u,m)-sqrt(m)*jacobi_cn(u,m))/sqrt(m)
4585 (defprop %jacobi_sn
4586 ((u m)
4587 ((mtimes simp) ((mexpt simp) m ((rat simp) -1 2))
4588 ((%log simp)
4589 ((mplus simp)
4590 ((mtimes simp) -1 ((mexpt simp) m ((rat simp) 1 2))
4591 ((%jacobi_cn simp) u m))
4592 ((%jacobi_dn simp) u m))))
4593 nil)
4594 integral)
4595
4596 ;; A&S 16.24.2: integrate(jacobi_cn(u,m),u) = acos(jacobi_dn(u,m))/sqrt(m)
4597 (defprop %jacobi_cn
4598 ((u m)
4599 ((mtimes simp) ((mexpt simp) m ((rat simp) -1 2))
4600 ((%acos simp) ((%jacobi_dn simp) u m)))
4601 nil)
4602 integral)
4603
4604 ;; A&S 16.24.3: integrate(jacobi_dn(u,m),u) = asin(jacobi_sn(u,m))
4605 (defprop %jacobi_dn
4606 ((u m)
4607 ((%asin simp) ((%jacobi_sn simp) u m))
4608 nil)
4609 integral)
4610
4611 ;; A&S 16.24.4: integrate(jacobi_cd(u,m),u)
4612 ;; = log(jacobi_nd(u,m)+sqrt(m)*jacobi_sd(u,m))/sqrt(m)
4613 (defprop %jacobi_cd
4614 ((u m)
4615 ((mtimes simp) ((mexpt simp) m ((rat simp) -1 2))
4616 ((%log simp)
4617 ((mplus simp) ((%jacobi_nd simp) u m)
4618 ((mtimes simp) ((mexpt simp) m ((rat simp) 1 2))
4619 ((%jacobi_sd simp) u m)))))
4620 nil)
4621 integral)
4622
4623 ;; integrate(jacobi_sd(u,m),u)
4624 ;;
4625 ;; A&S 16.24.5 gives
4626 ;; asin(-sqrt(m)*jacobi_cd(u,m))/sqrt(m*m_1), where m + m_1 = 1
4627 ;; but this does not pass some simple tests.
4628 ;;
4629 ;; functions.wolfram.com 09.35.21.001.01 gives
4630 ;; -asin(sqrt(m)*jacobi_cd(u,m))*sqrt(1-m*jacobi_cd(u,m)^2)*jacobi_dn(u,m)/((1-m)*sqrt(m))
4631 ;; and this does pass.
4632 (defprop %jacobi_sd
4633 ((u m)
4634 ((mtimes simp) -1
4635 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
4636 ((mexpt simp) m ((rat simp) -1 2))
4637 ((mexpt simp)
4638 ((mplus simp) 1
4639 ((mtimes simp) -1 $m ((mexpt simp) ((%jacobi_cd simp) u m) 2)))
4640 ((rat simp) 1 2))
4641 ((%jacobi_dn simp) u m)
4642 ((%asin simp)
4643 ((mtimes simp) ((mexpt simp) m ((rat simp) 1 2))
4644 ((%jacobi_cd simp) u m))))
4645 nil)
4646 integral)
4647
4648 ;; integrate(jacobi_nd(u,m),u)
4649 ;;
4650 ;; A&S 16.24.6 gives
4651 ;; acos(jacobi_cd(u,m))/sqrt(m_1), where m + m_1 = 1
4652 ;; but this does not pass some simple tests.
4653 ;;
4654 ;; functions.wolfram.com 09.32.21.0001.01 gives
4655 ;; sqrt(1-jacobi_cd(u,m)^2)*acos(jacobi_cd(u,m))/((1-m)*jacobi_sd(u,m))
4656 ;; and this does pass.
4657 (defprop %jacobi_nd
4658 ((u m)
4659 ((mtimes simp) ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m)) -1)
4660 ((mexpt simp)
4661 ((mplus simp) 1
4662 ((mtimes simp) -1 ((mexpt simp) ((%jacobi_cd simp) u m) 2)))
4663 ((rat simp) 1 2))
4664 ((mexpt simp) ((%jacobi_sd simp) u m) -1)
4665 ((%acos simp) ((%jacobi_cd simp) u m)))
4666 nil)
4667 integral)
4668
4669 ;; A&S 16.24.7: integrate(jacobi_dc(u,m),u) = log(jacobi_nc(u,m)+jacobi_sc(u,m))
4670 (defprop %jacobi_dc
4671 ((u m)
4672 ((%log simp) ((mplus simp) ((%jacobi_nc simp) u m) ((%jacobi_sc simp) u m)))
4673 nil)
4674 integral)
4675
4676 ;; A&S 16.24.8: integrate(jacobi_nc(u,m),u)
4677 ;; = log(jacobi_dc(u,m)+sqrt(m_1)*jacobi_sc(u,m))/sqrt(m_1), where m + m_1 = 1
4678 (defprop %jacobi_nc
4679 ((u m)
4680 ((mtimes simp)
4681 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m))
4682 ((rat simp) -1 2))
4683 ((%log simp)
4684 ((mplus simp) ((%jacobi_dc simp) u m)
4685 ((mtimes simp)
4686 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m))
4687 ((rat simp) 1 2))
4688 ((%jacobi_sc simp) u m)))))
4689 nil)
4690 integral)
4691
4692 ;; A&S 16.24.9: integrate(jacobi_sc(u,m),u)
4693 ;; = log(jacobi_dc(u,m)+sqrt(m_1)*jacobi_nc(u,m))/sqrt(m_1), where m + m_1 = 1
4694 (defprop %jacobi_sc
4695 ((u m)
4696 ((mtimes simp)
4697 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m))
4698 ((rat simp) -1 2))
4699 ((%log simp)
4700 ((mplus simp) ((%jacobi_dc simp) u m)
4701 ((mtimes simp)
4702 ((mexpt simp) ((mplus simp) 1 ((mtimes simp) -1 m))
4703 ((rat simp) 1 2))
4704 ((%jacobi_nc simp) u m)))))
4705 nil)
4706 integral)
4707
4708 ;; A&S 16.24.10: integrate(jacobi_ns(u,m),u)
4709 ;; = log(jacobi_ds(u,m)-jacobi_cs(u,m))
4710 (defprop %jacobi_ns
4711 ((u m)
4712 ((%log simp)
4713 ((mplus simp) ((mtimes simp) -1 ((%jacobi_cs simp) u m))
4714 ((%jacobi_ds simp) u m)))
4715 nil)
4716 integral)
4717
4718 ;; integrate(jacobi_ds(u,m),u)
4719 ;;
4720 ;; A&S 16.24.11 gives
4721 ;; log(jacobi_ds(u,m)-jacobi_cs(u,m))
4722 ;; but this does not pass some simple tests.
4723 ;;
4724 ;; functions.wolfram.com 09.30.21.0001.01 gives
4725 ;; log((1-jacobi_cn(u,m))/jacobi_sn(u,m))
4726 ;;
4727 (defprop %jacobi_ds
4728 ((u m)
4729 ((%log simp)
4730 ((mtimes simp)
4731 ((mplus simp) 1 ((mtimes simp) -1 ((%jacobi_cn simp) u m)))
4732 ((mexpt simp) ((%jacobi_sn simp) u m) -1)))
4733 nil)
4734 integral)
4735
4736 ;; A&S 16.24.12: integrate(jacobi_cs(u,m),u) = log(jacobi_ns(u,m)-jacobi_ds(u,m))
4737 (defprop %jacobi_cs
4738 ((u m)
4739 ((%log simp)
4740 ((mplus simp) ((mtimes simp) -1 ((%jacobi_ds simp) u m))
4741 ((%jacobi_ns simp) u m)))
4742 nil)
4743 integral)
4744
4745 ;; functions.wolfram.com 09.48.21.0001.01
4746 ;; integrate(inverse_jacobi_sn(u,m),u) =
4747 ;; inverse_jacobi_sn(u,m)*u
4748 ;; - log( jacobi_dn(inverse_jacobi_sn(u,m),m)
4749 ;; -sqrt(m)*jacobi_cn(inverse_jacobi_sn(u,m),m)) / sqrt(m)
4750 (defprop %inverse_jacobi_sn
4751 ((u m)
4752 ((mplus simp) ((mtimes simp) u ((%inverse_jacobi_sn simp) u m))
4753 ((mtimes simp) -1 ((mexpt simp) m ((rat simp) -1 2))
4754 ((%log simp)
4755 ((mplus simp)
4756 ((mtimes simp) -1 ((mexpt simp) m ((rat simp) 1 2))
4757 ((%jacobi_cn simp) ((%inverse_jacobi_sn simp) u m) m))
4758 ((%jacobi_dn simp) ((%inverse_jacobi_sn simp) u m) m)))))
4759 nil)
4760 integral)
4761
4762 ;; functions.wolfram.com 09.38.21.0001.01
4763 ;; integrate(inverse_jacobi_cn(u,m),u) =
4764 ;; u*inverse_jacobi_cn(u,m)
4765 ;; -%i*log(%i*jacobi_dn(inverse_jacobi_cn(u,m),m)/sqrt(m)
4766 ;; -jacobi_sn(inverse_jacobi_cn(u,m),m))
4767 ;; /sqrt(m)
4768 (defprop %inverse_jacobi_cn
4769 ((u m)
4770 ((mplus simp) ((mtimes simp) u ((%inverse_jacobi_cn simp) u m))
4771 ((mtimes simp) -1 $%i ((mexpt simp) m ((rat simp) -1 2))
4772 ((%log simp)
4773 ((mplus simp)
4774 ((mtimes simp) $%i ((mexpt simp) m ((rat simp) -1 2))
4775 ((%jacobi_dn simp) ((%inverse_jacobi_cn simp) u m) m))
4776 ((mtimes simp) -1
4777 ((%jacobi_sn simp) ((%inverse_jacobi_cn simp) u m) m))))))
4778 nil)
4779 integral)
4780
4781 ;; functions.wolfram.com 09.41.21.0001.01
4782 ;; integrate(inverse_jacobi_dn(u,m),u) =
4783 ;; u*inverse_jacobi_dn(u,m)
4784 ;; - %i*log(%i*jacobi_cn(inverse_jacobi_dn(u,m),m)
4785 ;; +jacobi_sn(inverse_jacobi_dn(u,m),m))
4786 (defprop %inverse_jacobi_dn
4787 ((u m)
4788 ((mplus simp) ((mtimes simp) u ((%inverse_jacobi_dn simp) u m))
4789 ((mtimes simp) -1 $%i
4790 ((%log simp)
4791 ((mplus simp)
4792 ((mtimes simp) $%i
4793 ((%jacobi_cn simp) ((%inverse_jacobi_dn simp) u m) m))
4794 ((%jacobi_sn simp) ((%inverse_jacobi_dn simp) u m) m)))))
4795 nil)
4796 integral)
4797
4798 \f
4799 ;; Real and imaginary part for Jacobi elliptic functions.
4800 (defprop %jacobi_sn risplit-sn-cn-dn risplit-function)
4801 (defprop %jacobi_cn risplit-sn-cn-dn risplit-function)
4802 (defprop %jacobi_dn risplit-sn-cn-dn risplit-function)
4803
4804 (defun risplit-sn-cn-dn (expr)
4805 (let* ((arg (second expr))
4806 (param (third expr)))
4807 ;; We only split on the argument, not the order
4808 (destructuring-bind (arg-r . arg-i)
4809 (risplit arg)
4810 (cond ((=0 arg-i)
4811 ;; Pure real
4812 (cons (take (first expr) arg-r param)
4813 0))
4814 (t
4815 (let* ((s (ftake '%jacobi_sn arg-r param))
4816 (c (ftake '%jacobi_cn arg-r param))
4817 (d (ftake '%jacobi_dn arg-r param))
4818 (s1 (ftake '%jacobi_sn arg-i (sub 1 param)))
4819 (c1 (ftake '%jacobi_cn arg-i (sub 1 param)))
4820 (d1 (ftake '%jacobi_dn arg-i (sub 1 param)))
4821 (den (add (mul c1 c1)
4822 (mul param
4823 (mul (mul s s)
4824 (mul s1 s1))))))
4825 ;; Let s = jacobi_sn(x,m)
4826 ;; c = jacobi_cn(x,m)
4827 ;; d = jacobi_dn(x,m)
4828 ;; s1 = jacobi_sn(y,1-m)
4829 ;; c1 = jacobi_cn(y,1-m)
4830 ;; d1 = jacobi_dn(y,1-m)
4831 (case (caar expr)
4832 (%jacobi_sn
4833 ;; A&S 16.21.1
4834 ;; jacobi_sn(x+%i*y,m) =
4835 ;;
4836 ;; s*d1 + %i*c*d*s1*c1
4837 ;; -------------------
4838 ;; c1^2+m*s^2*s1^2
4839 ;;
4840 (cons (div (mul s d1) den)
4841 (div (mul c (mul d (mul s1 c1)))
4842 den)))
4843 (%jacobi_cn
4844 ;; A&S 16.21.2
4845 ;;
4846 ;; cn(x+%i_y, m) =
4847 ;;
4848 ;; c*c1 - %i*s*d*s1*d1
4849 ;; -------------------
4850 ;; c1^2+m*s^2*s1^2
4851 (cons (div (mul c c1) den)
4852 (div (mul -1
4853 (mul s (mul d (mul s1 d1))))
4854 den)))
4855 (%jacobi_dn
4856 ;; A&S 16.21.3
4857 ;;
4858 ;; dn(x+%i_y, m) =
4859 ;;
4860 ;; d*c1*d1 - %i*m*s*c*s1
4861 ;; ---------------------
4862 ;; c1^2+m*s^2*s1^2
4863 (cons (div (mul d (mul c1 d1))
4864 den)
4865 (div (mul -1 (mul param (mul s (mul c s1))))
4866 den))))))))))
4867
4868 \f
4869 ;; Jacobi amplitude function.
4870
4871 (in-package #:bigfloat)
4872
4873 ;; Arithmetic-Geometric Mean algorithm for real or complex numbers.
4874 ;; See https://dlmf.nist.gov/22.20.ii.
4875 (let ((an (make-array 100 :fill-pointer 0))
4876 (bn (make-array 100 :fill-pointer 0))
4877 (cn (make-array 100 :fill-pointer 0)))
4878 ;; Instead of allocating these array anew each time, we'll reuse
4879 ;; them and allow them to grow as needed.
4880 (defun agm (a0 b0 c0 tol)
4881 "Arithmetic-Geometric Mean algorithm for real or complex a0, b0, c0.
4882 Algorithm continues until |c[n]| <= tol."
4883
4884 ;; DLMF (https://dlmf.nist.gov/22.20.ii) says for any real or
4885 ;; complex a0 and b0, b0/a0 must not be real and negative. Let's
4886 ;; check that.
4887 (let ((q (/ b0 a0)))
4888 (when (and (= (imagpart q) 0)
4889 (minusp (realpart q)))
4890 (error "Invalid arguments for AGM: ~A ~A~%" a0 b0)))
4891 (let ((nd (max (* 2 (ceiling (log (- (log tol 2))))) 8)))
4892 ;; DLMF (https://dlmf.nist.gov/22.20.ii) says that |c[n]| <=
4893 ;; C*2^(-2^n), for some constant C. Solve C*2^(-2^n) = tol to
4894 ;; get n = log(log(C/tol)/log(2))/log(2). Arbitrarily assume C
4895 ;; is one to get n = log(-(log(tol)/log(2)))/log(2). Thus, the
4896 ;; approximate number of term needed is n =
4897 ;; 1.44*log(-(1.44*log(tol))). Round to 2*log(-log2(tol)).
4898 (setf (fill-pointer an) 0
4899 (fill-pointer bn) 0
4900 (fill-pointer cn) 0)
4901 (vector-push-extend a0 an)
4902 (vector-push-extend b0 bn)
4903 (vector-push-extend c0 cn)
4904
4905 (do ((k 0 (1+ k)))
4906 ((or (<= (abs (aref cn k)) tol)
4907 (>= k nd))
4908 (if (>= k nd)
4909 (error "Failed to converge")
4910 (values k an bn cn)))
4911 (vector-push-extend (/ (+ (aref an k) (aref bn k)) 2) an)
4912 ;; DLMF (https://dlmf.nist.gov/22.20.ii) has conditions on how
4913 ;; to choose the square root depending on the phase of a[n-1]
4914 ;; and b[n-1]. We don't check for that here.
4915 (vector-push-extend (sqrt (* (aref an k) (aref bn k))) bn)
4916 (vector-push-extend (/ (- (aref an k) (aref bn k)) 2) cn)))))
4917
4918 (defun jacobi-am-agm (u m tol)
4919 "Evaluate the jacobi_am function from real u and m with |m| <= 1. This
4920 uses the AGM method until a tolerance of TOL is reached for the
4921 error."
4922 (multiple-value-bind (n an bn cn)
4923 (agm 1 (sqrt (- 1 m)) (sqrt m) tol)
4924 (declare (ignore bn))
4925 ;; See DLMF (https://dlmf.nist.gov/22.20.ii) for the algorithm.
4926 (let ((phi (* u (aref an n) (expt 2 n))))
4927 (loop for k from n downto 1
4928 do
4929 (setf phi (/ (+ phi (asin (* (/ (aref cn k)
4930 (aref an k))
4931 (sin phi))))
4932 2)))
4933 phi)))
4934
4935 ;; For real z and real m with |m| < 1, this appears to be less
4936 ;; accurate than using AGM. For am(3,0.7), we get 2.1263556865337185,
4937 ;; but AGM gives 2.126355671062825, which matches what Wolfram gives.
4938 #+nil
4939 (defun am-q-series (z m limit)
4940 (let* ((K (bf-elliptic-k m))
4941 (K-prime (bf-elliptic-k (- 1 m)))
4942 (2-arg (/ (* (%pi z) z)
4943 K))
4944 (q (exp (- (* (%pi z)
4945 (/ K-prime K))))))
4946 (format t "K = ~A~%" K)
4947 (format t "K-prime = ~A~%" K-prime)
4948 (format t "2-arg = ~A~%" 2-arg)
4949 (format t "q = ~A~%" q)
4950
4951 (do* ((n 1 (1+ n))
4952 (term (* (sin (* n 2-arg))
4953 (/ (* (expt q n))
4954 (* n (1+ (expt q (* 2 n))))))
4955 (* (sin (* n 2-arg))
4956 (/ (* (expt q n))
4957 (* n (1+ (expt q (* 2 n)))))))
4958 (sum term
4959 (+ sum term)))
4960 ((<= (abs term) limit)
4961 (+ (* 2 sum)
4962 (/ (* (%pi z) z)
4963 (* 2 K))))
4964 #+nil
4965 (format t "~4d: term = ~A~%" n term))))
4966
4967 ;; Compute Jacobi am for real or complex values of U and M. The args
4968 ;; must be floats or bigfloat::bigfloats. TOL is the tolerance used
4969 ;; by the AGM algorithm. It is ignored if the AGM algorithm is not
4970 ;; used.
4971 (defun bf-jacobi-am (u m tol)
4972 (cond ((and (realp u) (realp m) (<= (abs m) 1))
4973 ;; The case of real u and m with |m| <= 1. We can use AGM to
4974 ;; compute the result.
4975 (jacobi-am-agm (to u)
4976 (to m)
4977 tol))
4978 (t
4979 ;; Otherwise, use the formula am(u,m) = asin(jacobi_sn(u,m)).
4980 ;; (See DLMF https://dlmf.nist.gov/22.16.E1). This appears
4981 ;; to be what functions.wolfram.com is using in this case.
4982 (asin (sn (to u) (to m))))))
4983
4984 (in-package :maxima)
4985 (def-simplifier jacobi_am (u m)
4986 (let (args)
4987 (cond
4988 ((zerop1 u)
4989 ;; am(0,m) = 0
4990 0)
4991 ((zerop1 m)
4992 ;; am(u,0) = u
4993 u)
4994 ((eql m 1)
4995 ;; am(u,1) = 2*atan(exp(u))-%pi/2
4996 (sub (mul 2 (ftake '%atan (ftake '%exp u)))
4997 (div '$%pi 2)))
4998 ((float-numerical-eval-p u m)
4999 ;; For |m| <= 1, we want to use AGM. But for |m| > 1, we want
5000 ;; to use am(z,m) = asin(jacobi_sn(z,m)).
5001 (to (bigfloat::bf-jacobi-am ($float u)
5002 ($float m)
5003 double-float-epsilon)))
5004 ((setf args (complex-float-numerical-eval-p u m))
5005 (destructuring-bind (u m)
5006 args
5007 (to (bigfloat::bf-jacobi-am ($float u)
5008 ($float m)
5009 double-float-epsilon))))
5010 ((bigfloat-numerical-eval-p u m)
5011 (to (bigfloat::bf-jacobi-am (bigfloat:to ($bfloat u))
5012 (bigfloat:to ($bfloat m))
5013 (expt 2 (- fpprec)))))
5014 ((setf args (complex-bigfloat-numerical-eval-p u m))
5015 (destructuring-bind (u m)
5016 args
5017 (to (bigfloat::bf-jacobi-am (bigfloat:to ($bfloat u))
5018 (bigfloat:to ($bfloat m))
5019 (expt 2 (- fpprec))))))
5020 (t
5021 ;; Nothing to do
5022 (give-up)))))
Maxima, a Computer Algebra System