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