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