1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10-*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancments. ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
8 ;;; Copyright (c) 2001 by Raymond Toy. Replaced everything and added ;;;;;
9 ;;; support for symbolic manipulation of all 12 Jacobian elliptic ;;;;;
10 ;;; functions and the complete and incomplete elliptic integral ;;;;;
11 ;;; of the first, second and third kinds. ;;;;;
12 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
15 ;;(macsyma-module ellipt)
17 (defvar 3//2 '((rat simp
) 3 2))
18 (defvar 1//2 '((rat simp
) 1 2))
19 (defvar -
1//2 '((rat simp
) -
1 2))
22 ;;; Jacobian elliptic functions and elliptic integrals.
26 ;;; [1] Abramowitz and Stegun
27 ;;; [2] Lawden, Elliptic Functions and Applications, Springer-Verlag, 1989
28 ;;; [3] Whittaker & Watson, A Course of Modern Analysis
30 ;;; We use the definitions from Abramowitz and Stegun where our
31 ;;; sn(u,m) is sn(u|m). That is, the second arg is the parameter,
32 ;;; instead of the modulus k or modular angle alpha.
34 ;;; Note that m = k^2 and k = sin(alpha).
38 ;; Routines for computing the basic elliptic functions sn, cn, and dn.
41 ;; A&S gives several methods for computing elliptic functions
42 ;; including the AGM method (16.4) and ascending and descending Landen
43 ;; transformations (16.12 and 16.14). The latter are actually quite
44 ;; fast, only requiring simple arithmetic and square roots for the
45 ;; transformation until the last step. The AGM requires evaluation of
46 ;; several trigonometric functions at each stage.
48 ;; However, the Landen transformations appear to have some round-off
49 ;; issues. For example, using the ascending transform to compute cn,
50 ;; cn(100,.7) > 1e10. This is clearly not right since |cn| <= 1.
53 (in-package #-gcl
#:bigfloat
#+gcl
"BIGFLOAT")
55 (declaim (inline descending-transform ascending-transform
))
57 (defun ascending-transform (u m
)
60 ;; Take care in computing this transform. For the case where
61 ;; m is complex, we should compute sqrt(mu1) first as
62 ;; (1-sqrt(m))/(1+sqrt(m)), and then square this to get mu1.
63 ;; If not, we may choose the wrong branch when computing
65 (let* ((root-m (sqrt m
))
67 (expt (1+ root-m
) 2)))
68 (root-mu1 (/ (- 1 root-m
) (+ 1 root-m
)))
69 (v (/ u
(1+ root-mu1
))))
70 (values v mu root-mu1
)))
72 (defun descending-transform (u m
)
73 ;; Note: Don't calculate mu first, as given in 16.12.1. We
74 ;; should calculate sqrt(mu) = (1-sqrt(m1)/(1+sqrt(m1)), and
75 ;; then compute mu = sqrt(mu)^2. If we calculate mu first,
76 ;; sqrt(mu) loses information when m or m1 is complex.
77 (let* ((root-m1 (sqrt (- 1 m
)))
78 (root-mu (/ (- 1 root-m1
) (+ 1 root-m1
)))
79 (mu (* root-mu root-mu
))
80 (v (/ u
(1+ root-mu
))))
81 (values v mu root-mu
)))
84 ;; This appears to work quite well for both real and complex values
86 (defun elliptic-sn-descending (u m
)
90 ((< (abs m
) (epsilon u
))
94 (multiple-value-bind (v mu root-mu
)
95 (descending-transform u m
)
96 (let* ((new-sn (elliptic-sn-descending v mu
)))
97 (/ (* (1+ root-mu
) new-sn
)
98 (1+ (* root-mu new-sn new-sn
))))))))
100 ;; AGM scale. See A&S 17.6
104 ;; a[n] = (a[n-1]+b[n-1])/2, b[n] = sqrt(a[n-1]*b[n-1]), c[n] = (a[n-1]-b[n-1])/2.
106 ;; We stop when abs(c[n]) <= 10*eps
108 ;; A list of (n a[n] b[n] c[n]) is returned.
109 (defun agm-scale (a b c
)
111 while
(> (abs c
) (* 10 (epsilon c
)))
112 collect
(list n a b c
)
113 do
(psetf a
(/ (+ a b
) 2)
117 ;; WARNING: This seems to have accuracy problems when u is complex. I
118 ;; (rtoy) do not know why. For example (jacobi-agm #c(1e0 1e0) .7e0)
121 ;; #C(1.134045970915582 0.3522523454566013)
122 ;; #C(0.57149659007575 -0.6989899153338323)
123 ;; #C(0.6229715431044184 -0.4488635962149656)
125 ;; But the actual value of sn(1+%i, .7) is .3522523469224946 %i +
126 ;; 1.134045971912365. We've lost about 7 digits of accuracy!
127 (defun jacobi-agm (u m
)
130 ;; Compute the AGM scale with a = 1, b = sqrt(1-m), c = sqrt(m).
132 ;; Then phi[N] = 2^N*a[N]*u and compute phi[n] from
134 ;; sin(2*phi[n-1] - phi[n]) = c[n]/a[n]*sin(phi[n])
138 ;; sn(u|m) = sin(phi[0]), cn(u|m) = cos(phi[0])
139 ;; dn(u|m) = cos(phi[0])/cos(phi[1]-phi[0])
141 ;; Returns the three values sn, cn, dn.
142 (let* ((agm-data (nreverse (rest (agm-scale 1 (sqrt (- 1 m
)) (sqrt m
)))))
143 (phi (destructuring-bind (n a b c
)
145 (declare (ignore b c
))
148 (dolist (agm agm-data
)
149 (destructuring-bind (n a b c
)
151 (declare (ignore n b
))
153 phi
(/ (+ phi
(asin (* (/ c a
) (sin phi
)))) 2))))
154 (values (sin phi
) (cos phi
) (/ (cos phi
) (cos (- phi1 phi
))))))
158 ;; jacobi_sn(u,0) = sin(u). Should we use A&S 16.13.1 if m
161 ;; sn(u,m) = sin(u) - 1/4*m(u-sin(u)*cos(u))*cos(u)
164 ;; jacobi_sn(u,1) = tanh(u). Should we use A&S 16.15.1 if m
165 ;; is close enough to 1?
167 ;; sn(u,m) = tanh(u) + 1/4*(1-m)*(sinh(u)*cosh(u)-u)*sech(u)^2
170 ;; Use the ascending Landen transformation to compute sn.
171 (let ((s (elliptic-sn-descending u m
)))
172 (if (and (realp u
) (realp m
))
178 ;; jacobi_dn(u,0) = 1. Should we use A&S 16.13.3 for small m?
180 ;; dn(u,m) = 1 - 1/2*m*sin(u)^2
183 ;; jacobi_dn(u,1) = sech(u). Should we use A&S 16.15.3 if m
184 ;; is close enough to 1?
186 ;; dn(u,m) = sech(u) + 1/4*(1-m)*(sinh(u)*cosh(u)+u)*tanh(u)*sech(u)
189 ;; Use the Gauss transformation from
190 ;; http://functions.wolfram.com/09.29.16.0013.01:
193 ;; dn((1+sqrt(m))*z, 4*sqrt(m)/(1+sqrt(m))^2)
194 ;; = (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2)
198 ;; dn(y, mu) = (1-sqrt(m)*sn(z, m)^2)/(1+sqrt(m)*sn(z,m)^2)
200 ;; where z = y/(1+sqrt(m)) and mu=4*sqrt(m)/(1+sqrt(m))^2.
202 ;; Solve for m, and we get
204 ;; sqrt(m) = -(mu+2*sqrt(1-mu)-2)/mu or (-mu+2*sqrt(1-mu)+2)/mu.
206 ;; I don't think it matters which sqrt we use, so I (rtoy)
207 ;; arbitrarily choose the first one above.
209 ;; Note that (1-sqrt(1-mu))/(1+sqrt(1-mu)) is the same as
210 ;; -(mu+2*sqrt(1-mu)-2)/mu. Also, the former is more
211 ;; accurate for small mu.
212 (let* ((root (let ((root-1-m (sqrt (- 1 m
))))
216 (s (elliptic-sn-descending z
(* root root
)))
223 ;; jacobi_cn(u,0) = cos(u). Should we use A&S 16.13.2 for
226 ;; cn(u,m) = cos(u) + 1/4*m*(u-sin(u)*cos(u))*sin(u)
229 ;; jacobi_cn(u,1) = sech(u). Should we use A&S 16.15.3 if m
230 ;; is close enough to 1?
232 ;; cn(u,m) = sech(u) - 1/4*(1-m)*(sinh(u)*cosh(u)-u)*tanh(u)*sech(u)
235 ;; Use the ascending Landen transformation, A&S 16.14.3.
236 (multiple-value-bind (v mu root-mu1
)
237 (ascending-transform u m
)
239 (* (/ (+ 1 root-mu1
) mu
)
240 (/ (- (* d d
) root-mu1
)
245 ;; Tell maxima what the derivatives are.
247 ;; Lawden says the derivative wrt to k but that's not what we want.
249 ;; Here's the derivation we used, based on how Lawden get's his results.
253 ;; diff(sn(u,m),m) = s
254 ;; diff(cn(u,m),m) = p
255 ;; diff(dn(u,m),m) = q
257 ;; From the derivatives of sn, cn, dn wrt to u, we have
259 ;; diff(sn(u,m),u) = cn(u)*dn(u)
260 ;; diff(cn(u,m),u) = -cn(u)*dn(u)
261 ;; diff(dn(u,m),u) = -m*sn(u)*cn(u)
264 ;; Differentiate these wrt to m:
266 ;; diff(s,u) = p*dn + cn*q
267 ;; diff(p,u) = -p*dn - q*dn
268 ;; diff(q,u) = -sn*cn - m*s*cn - m*sn*q
272 ;; sn(u)^2 + cn(u)^2 = 1
273 ;; dn(u)^2 + m*sn(u)^2 = 1
275 ;; Differentiate these wrt to m:
278 ;; 2*dn*q + sn^2 + 2*m*sn*s = 0
283 ;; q = -m*s*sn/dn - sn^2/dn/2
286 ;; diff(s,u) = -s*sn*dn/cn - m*s*sn*cn/dn - sn^2*cn/dn/2
290 ;; diff(s,u) + s*(sn*dn/cn + m*sn*cn/dn) = -1/2*sn^2*cn/dn
292 ;; diff(s,u) + s*sn/cn/dn*(dn^2 + m*cn^2) = -1/2*sn^2*cn/dn
294 ;; Multiply through by the integrating factor 1/cn/dn:
296 ;; diff(s/cn/dn, u) = -1/2*sn^2/dn^2 = -1/2*sd^2.
298 ;; Interate this to get
300 ;; s/cn/dn = C + -1/2*int sd^2
302 ;; It can be shown that C is zero.
304 ;; We know that (by differentiating this expression)
306 ;; int dn^2 = (1-m)*u+m*sn*cd + m*(1-m)*int sd^2
310 ;; int sd^2 = 1/m/(1-m)*int dn^2 - u/m - sn*cd/(1-m)
314 ;; s/cn/dn = u/(2*m) + sn*cd/(2*(1-m)) - 1/2/m/(1-m)*int dn^2
318 ;; s = 1/(2*m)*u*cn*dn + 1/(2*(1-m))*sn*cn^2 - 1/2/(m*(1-m))*cn*dn*E(u)
320 ;; where E(u) = int dn^2 = elliptic_e(am(u)) = elliptic_e(asin(sn(u)))
322 ;; This is our desired result:
324 ;; s = 1/(2*m)*cn*dn*[u - elliptic_e(asin(sn(u)),m)/(1-m)] + sn*cn^2/2/(1-m)
327 ;; Since diff(cn(u,m),m) = p = -s*sn/cn, we have
329 ;; p = -1/(2*m)*sn*dn[u - elliptic_e(asin(sn(u)),m)/(1-m)] - sn^2*cn/2/(1-m)
331 ;; diff(dn(u,m),m) = q = -m*s*sn/dn - sn^2/dn/2
333 ;; q = -1/2*sn*cn*[u-elliptic_e(asin(sn),m)/(1-m)] - m*sn^2*cn^2/dn/2/(1-m)
337 ;; = -1/2*sn*cn*[u-elliptic_e(asin(sn),m)/(1-m)] + dn*sn^2/2/(m-1)
341 ((mtimes) ((%jacobi_cn
) u m
) ((%jacobi_dn
) u m
))
343 ((mtimes simp
) ((rat simp
) 1 2)
344 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
345 ((mexpt simp
) ((%jacobi_cn simp
) u m
) 2) ((%jacobi_sn simp
) u m
))
346 ((mtimes simp
) ((rat simp
) 1 2) ((mexpt simp
) m -
1)
347 ((%jacobi_cn simp
) u m
) ((%jacobi_dn simp
) u m
)
349 ((mtimes simp
) -
1 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
350 ((%elliptic_e simp
) ((%asin simp
) ((%jacobi_sn simp
) u m
)) m
))))))
355 ((mtimes simp
) -
1 ((%jacobi_sn simp
) u m
) ((%jacobi_dn simp
) u m
))
357 ((mtimes simp
) ((rat simp
) -
1 2)
358 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
359 ((%jacobi_cn simp
) u m
) ((mexpt simp
) ((%jacobi_sn simp
) u m
) 2))
360 ((mtimes simp
) ((rat simp
) -
1 2) ((mexpt simp
) m -
1)
361 ((%jacobi_dn simp
) u m
) ((%jacobi_sn simp
) u m
)
363 ((mtimes simp
) -
1 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
364 ((%elliptic_e simp
) ((%asin simp
) ((%jacobi_sn simp
) u m
)) m
))))))
369 ((mtimes) -
1 m
((%jacobi_sn
) u m
) ((%jacobi_cn
) u m
))
371 ((mtimes simp
) ((rat simp
) -
1 2)
372 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
373 ((%jacobi_dn simp
) u m
) ((mexpt simp
) ((%jacobi_sn simp
) u m
) 2))
374 ((mtimes simp
) ((rat simp
) -
1 2) ((%jacobi_cn simp
) u m
)
375 ((%jacobi_sn simp
) u m
)
378 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
379 ((%elliptic_e simp
) ((%asin simp
) ((%jacobi_sn simp
) u m
)) m
))))))
382 ;; The inverse elliptic functions.
384 ;; F(phi|m) = asn(sin(phi),m)
386 ;; so asn(u,m) = F(asin(u)|m)
387 (defprop %inverse_jacobi_sn
390 ;; inverse_jacobi_sn(x) = integrate(1/sqrt(1-t^2)/sqrt(1-m*t^2),t,0,x)
391 ;; -> 1/sqrt(1-x^2)/sqrt(1-m*x^2)
393 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2)))
395 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
((mexpt simp
) x
2)))
397 ;; diff(F(asin(u)|m),m)
398 ((mtimes simp
) ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
401 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2)))
403 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
((mexpt simp
) x
2)))
405 ((mtimes simp
) ((mexpt simp
) m -
1)
406 ((mplus simp
) ((%elliptic_e simp
) ((%asin simp
) x
) m
)
407 ((mtimes simp
) -
1 ((mplus simp
) 1 ((mtimes simp
) -
1 m
))
408 ((%elliptic_f simp
) ((%asin simp
) x
) m
)))))))
411 ;; Let u = inverse_jacobi_cn(x). Then jacobi_cn(u) = x or
412 ;; sqrt(1-jacobi_sn(u)^2) = x. Or
414 ;; jacobi_sn(u) = sqrt(1-x^2)
416 ;; So u = inverse_jacobi_sn(sqrt(1-x^2),m) = inverse_jacob_cn(x,m)
418 (defprop %inverse_jacobi_cn
420 ;; Whittaker and Watson, 22.121
421 ;; inverse_jacobi_cn(u,m) = integrate(1/sqrt(1-t^2)/sqrt(1-m+m*t^2), t, u, 1)
422 ;; -> -1/sqrt(1-x^2)/sqrt(1-m+m*x^2)
424 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2)))
427 ((mplus simp
) 1 ((mtimes simp
) -
1 m
)
428 ((mtimes simp
) m
((mexpt simp
) x
2)))
430 ((mtimes simp
) ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
435 ((mtimes simp
) -
1 m
((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2)))))
437 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2))) ((rat simp
) 1 2))
439 ((mtimes simp
) ((mexpt simp
) m -
1)
443 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2))) ((rat simp
) 1 2)))
445 ((mtimes simp
) -
1 ((mplus simp
) 1 ((mtimes simp
) -
1 m
))
448 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2))) ((rat simp
) 1 2)))
452 ;; Let u = inverse_jacobi_dn(x). Then
454 ;; jacobi_dn(u) = x or
456 ;; x^2 = jacobi_dn(u)^2 = 1 - m*jacobi_sn(u)^2
458 ;; so jacobi_sn(u) = sqrt(1-x^2)/sqrt(m)
460 ;; or u = inverse_jacobi_sn(sqrt(1-x^2)/sqrt(m))
461 (defprop %inverse_jacobi_dn
463 ;; Whittaker and Watson, 22.121
464 ;; inverse_jacobi_dn(u,m) = integrate(1/sqrt(1-t^2)/sqrt(t^2-(1-m)), t, u, 1)
465 ;; -> -1/sqrt(1-x^2)/sqrt(x^2+m-1)
467 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2)))
469 ((mexpt simp
) ((mplus simp
) -
1 m
((mexpt simp
) x
2)) ((rat simp
) -
1 2)))
471 ((mtimes simp
) ((rat simp
) -
1 2) ((mexpt simp
) m
((rat simp
) -
3 2))
474 ((mtimes simp
) -
1 ((mexpt simp
) m -
1)
475 ((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2)))))
477 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2)))
479 ((mexpt simp
) ((mabs simp
) x
) -
1))
480 ((mtimes simp
) ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
482 ((mtimes simp
) -
1 ((mexpt simp
) m
((rat simp
) -
1 2))
485 ((mtimes simp
) -
1 ((mexpt simp
) m -
1)
486 ((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2)))))
488 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 ((mexpt simp
) x
2)))
490 ((mexpt simp
) ((mabs simp
) x
) -
1))
491 ((mtimes simp
) ((mexpt simp
) m -
1)
495 ((mtimes simp
) ((mexpt simp
) m
((rat simp
) -
1 2))
496 ((mexpt simp
) ((mplus simp
) 1
497 ((mtimes simp
) -
1 ((mexpt simp
) x
2)))
500 ((mtimes simp
) -
1 ((mplus simp
) 1 ((mtimes simp
) -
1 m
))
503 ((mtimes simp
) ((mexpt simp
) m
((rat simp
) -
1 2))
504 ((mexpt simp
) ((mplus simp
) 1
505 ((mtimes simp
) -
1 ((mexpt simp
) x
2)))
511 ;; Possible forms of a complex number:
515 ;; ((mplus simp) 2.3 ((mtimes simp) 2.3 $%i))
516 ;; ((mplus simp) 2.3 $%i))
517 ;; ((mtimes simp) 2.3 $%i)
521 ;; Is argument u a complex number with real and imagpart satisfying predicate ntypep?
522 (defun complex-number-p (u &optional
(ntypep 'numberp
))
524 (labels ((a1 (x) (cadr x
))
527 (N (x) (funcall ntypep x
)) ; N
528 (i (x) (and (eq x
'$%i
) (N 1))) ; %i
529 (N+i
(x) (and (null (a3+ x
)) ; mplus test is precondition
531 (or (and (i (a2 x
)) (setq I
1) t
)
532 (and (mtimesp (a2 x
)) (N*i
(a2 x
))))))
533 (N*i
(x) (and (null (a3+ x
)) ; mtimes test is precondition
536 (declare (inline a1 a2 a3
+ N i N
+i N
*i
))
537 (cond ((N u
) (values t u
0)) ;2.3
538 ((atom u
) (if (i u
) (values t
0 1))) ;%i
539 ((mplusp u
) (if (N+i u
) (values t R I
))) ;N+%i, N+N*%i
540 ((mtimesp u
) (if (N*i u
) (values t R I
))) ;N*%i
543 (defun complexify (x)
544 ;; Convert a Lisp number to a maxima number
546 ((complexp x
) (add (realpart x
) (mul '$%i
(imagpart x
))))
547 (t (merror (intl:gettext
"COMPLEXIFY: argument must be a Lisp real or complex number.~%COMPLEXIFY: found: ~:M") x
))))
549 (defun kc-arg (exp m
)
550 ;; Replace elliptic_kc(m) in the expression with sym. Check to see
551 ;; if the resulting expression is linear in sym and the constant
552 ;; term is zero. If so, return the coefficient of sym, i.e, the
553 ;; coefficient of elliptic_kc(m).
554 (let* ((sym (gensym))
555 (arg (maxima-substitute sym
`((%elliptic_kc
) ,m
) exp
)))
556 (if (and (not (equalp arg exp
))
558 (zerop1 (coefficient arg sym
0)))
559 (coefficient arg sym
1)
562 (defun kc-arg2 (exp m
)
563 ;; Replace elliptic_kc(m) in the expression with sym. Check to see
564 ;; if the resulting expression is linear in sym and the constant
565 ;; term is zero. If so, return the coefficient of sym, i.e, the
566 ;; coefficient of elliptic_kc(m), and the constant term. Otherwise,
568 (let* ((sym (gensym))
569 (arg (maxima-substitute sym
`((%elliptic_kc
) ,m
) exp
)))
570 (if (and (not (equalp arg exp
))
572 (list (coefficient arg sym
1)
573 (coefficient arg sym
0))
576 ;; Tell maxima how to simplify the functions
578 (def-simplifier jacobi_sn
(u m
)
581 ((float-numerical-eval-p u m
)
582 (to (bigfloat::sn
(bigfloat:to
($float u
)) (bigfloat:to
($float m
)))))
583 ((setf args
(complex-float-numerical-eval-p u m
))
584 (destructuring-bind (u m
)
586 (to (bigfloat::sn
(bigfloat:to
($float u
)) (bigfloat:to
($float m
))))))
587 ((bigfloat-numerical-eval-p u m
)
588 (to (bigfloat::sn
(bigfloat:to
($bfloat u
)) (bigfloat:to
($bfloat m
)))))
589 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
590 (destructuring-bind (u m
)
592 (to (bigfloat::sn
(bigfloat:to
($bfloat u
)) (bigfloat:to
($bfloat m
))))))
602 ((and $trigsign
(mminusp* u
))
603 (neg (ftake* '%jacobi_sn
(neg u
) m
)))
606 (member (caar u
) '(%inverse_jacobi_sn
618 (alike1 (third u
) m
))
619 (let ((inv-arg (second u
)))
622 ;; jacobi_sn(inverse_jacobi_sn(u,m), m) = u
625 ;; inverse_jacobi_ns(u,m) = inverse_jacobi_sn(1/u,m)
628 ;; sn(x)^2 + cn(x)^2 = 1 so sn(x) = sqrt(1-cn(x)^2)
629 (power (sub 1 (mul inv-arg inv-arg
)) 1//2))
631 ;; inverse_jacobi_nc(u) = inverse_jacobi_cn(1/u)
632 (ftake '%jacobi_sn
(ftake '%inverse_jacobi_cn
(div 1 inv-arg
) m
)
635 ;; dn(x)^2 + m*sn(x)^2 = 1 so
636 ;; sn(x) = 1/sqrt(m)*sqrt(1-dn(x)^2)
637 (mul (div 1 (power m
1//2))
638 (power (sub 1 (mul inv-arg inv-arg
)) 1//2)))
640 ;; inverse_jacobi_nd(u) = inverse_jacobi_dn(1/u)
641 (ftake '%jacobi_sn
(ftake '%inverse_jacobi_dn
(div 1 inv-arg
) m
)
644 ;; See below for inverse_jacobi_sc.
645 (div inv-arg
(power (add 1 (mul inv-arg inv-arg
)) 1//2)))
647 ;; inverse_jacobi_cs(u) = inverse_jacobi_sc(1/u)
648 (ftake '%jacobi_sn
(ftake '%inverse_jacobi_sc
(div 1 inv-arg
) m
)
651 ;; See below for inverse_jacobi_sd
652 (div inv-arg
(power (add 1 (mul m
(mul inv-arg inv-arg
))) 1//2)))
654 ;; inverse_jacobi_ds(u) = inverse_jacobi_sd(1/u)
655 (ftake '%jacobi_sn
(ftake '%inverse_jacobi_sd
(div 1 inv-arg
) m
)
659 (div (power (sub 1 (mul inv-arg inv-arg
)) 1//2)
660 (power (sub 1 (mul m
(mul inv-arg inv-arg
))) 1//2)))
662 (ftake '%jacobi_sn
(ftake '%inverse_jacobi_cd
(div 1 inv-arg
) m
) m
)))))
663 ;; A&S 16.20.1 (Jacobi's Imaginary transformation)
664 ((and $%iargs
(multiplep u
'$%i
))
666 (ftake* '%jacobi_sc
(coeff u
'$%i
1) (add 1 (neg m
)))))
667 ((setq coef
(kc-arg2 u m
))
671 (destructuring-bind (lin const
)
673 (cond ((integerp lin
)
676 ;; sn(4*m*K + u) = sn(u), sn(0) = 0
679 (ftake '%jacobi_sn const m
)))
681 ;; sn(4*m*K + K + u) = sn(K+u) = cd(u)
685 (ftake '%jacobi_cd const m
)))
687 ;; sn(4*m*K+2*K + u) = sn(2*K+u) = -sn(u)
691 (neg (ftake '%jacobi_sn const m
))))
693 ;; sn(4*m*K+3*K+u) = sn(2*K + K + u) = -sn(K+u) = -cd(u)
697 (neg (ftake '%jacobi_cd const m
))))))
698 ((and (alike1 lin
1//2)
702 ;; sn(1/2*K) = 1/sqrt(1+sqrt(1-m))
704 (power (add 1 (power (sub 1 m
) 1//2))
706 ((and (alike1 lin
3//2)
710 ;; sn(1/2*K + K) = cd(1/2*K,m)
711 (ftake '%jacobi_cd
(mul 1//2
712 (ftake '%elliptic_kc m
))
720 (def-simplifier jacobi_cn
(u m
)
723 ((float-numerical-eval-p u m
)
724 (to (bigfloat::cn
(bigfloat:to
($float u
)) (bigfloat:to
($float m
)))))
725 ((setf args
(complex-float-numerical-eval-p u m
))
726 (destructuring-bind (u m
)
728 (to (bigfloat::cn
(bigfloat:to
($float u
)) (bigfloat:to
($float m
))))))
729 ((bigfloat-numerical-eval-p u m
)
730 (to (bigfloat::cn
(bigfloat:to
($bfloat u
)) (bigfloat:to
($bfloat m
)))))
731 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
732 (destructuring-bind (u m
)
734 (to (bigfloat::cn
(bigfloat:to
($bfloat u
)) (bigfloat:to
($bfloat m
))))))
744 ((and $trigsign
(mminusp* u
))
745 (ftake* '%jacobi_cn
(neg u
) m
))
748 (member (caar u
) '(%inverse_jacobi_sn
760 (alike1 (third u
) m
))
761 (cond ((eq (caar u
) '%inverse_jacobi_cn
)
764 ;; I'm lazy. Use cn(x) = sqrt(1-sn(x)^2). Hope
766 (power (sub 1 (power (ftake '%jacobi_sn u
(third u
)) 2))
768 ;; A&S 16.20.2 (Jacobi's Imaginary transformation)
769 ((and $%iargs
(multiplep u
'$%i
))
770 (ftake* '%jacobi_nc
(coeff u
'$%i
1) (add 1 (neg m
))))
771 ((setq coef
(kc-arg2 u m
))
775 (destructuring-bind (lin const
)
777 (cond ((integerp lin
)
780 ;; cn(4*m*K + u) = cn(u),
784 (ftake '%jacobi_cn const m
)))
786 ;; cn(4*m*K + K + u) = cn(K+u) = -sqrt(m1)*sd(u)
790 (neg (mul (power (sub 1 m
) 1//2)
791 (ftake '%jacobi_sd const m
)))))
793 ;; cn(4*m*K + 2*K + u) = cn(2*K+u) = -cn(u)
797 (neg (ftake '%jacobi_cn const m
))))
799 ;; cn(4*m*K + 3*K + u) = cn(2*K + K + u) =
800 ;; -cn(K+u) = sqrt(m1)*sd(u)
805 (mul (power (sub 1 m
) 1//2)
806 (ftake '%jacobi_sd const m
))))))
807 ((and (alike1 lin
1//2)
810 ;; cn(1/2*K) = (1-m)^(1/4)/sqrt(1+sqrt(1-m))
811 (mul (power (sub 1 m
) (div 1 4))
821 (def-simplifier jacobi_dn
(u m
)
824 ((float-numerical-eval-p u m
)
825 (to (bigfloat::dn
(bigfloat:to
($float u
)) (bigfloat:to
($float m
)))))
826 ((setf args
(complex-float-numerical-eval-p u m
))
827 (destructuring-bind (u m
)
829 (to (bigfloat::dn
(bigfloat:to
($float u
)) (bigfloat:to
($float m
))))))
830 ((bigfloat-numerical-eval-p u m
)
831 (to (bigfloat::dn
(bigfloat:to
($bfloat u
)) (bigfloat:to
($bfloat m
)))))
832 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
833 (destructuring-bind (u m
)
835 (to (bigfloat::dn
(bigfloat:to
($bfloat u
)) (bigfloat:to
($bfloat m
))))))
845 ((and $trigsign
(mminusp* u
))
846 (ftake* '%jacobi_dn
(neg u
) m
))
849 (member (caar u
) '(%inverse_jacobi_sn
861 (alike1 (third u
) m
))
862 (cond ((eq (caar u
) '%inverse_jacobi_dn
)
863 ;; jacobi_dn(inverse_jacobi_dn(u,m), m) = u
866 ;; Express in terms of sn:
867 ;; dn(x) = sqrt(1-m*sn(x)^2)
869 (power (ftake '%jacobi_sn u m
) 2)))
871 ((zerop1 ($ratsimp
(sub u
(power (sub 1 m
) 1//2))))
873 ;; dn(sqrt(1-m),m) = K(m)
874 (ftake '%elliptic_kc m
))
875 ;; A&S 16.20.2 (Jacobi's Imaginary transformation)
876 ((and $%iargs
(multiplep u
'$%i
))
877 (ftake* '%jacobi_dc
(coeff u
'$%i
1)
879 ((setq coef
(kc-arg2 u m
))
882 ;; dn(m*K+u) has period 2K
884 (destructuring-bind (lin const
)
886 (cond ((integerp lin
)
889 ;; dn(2*m*K + u) = dn(u)
893 ;; dn(4*m*K+2*K + u) = dn(2*K+u) = dn(u)
894 (ftake '%jacobi_dn const m
)))
896 ;; dn(2*m*K + K + u) = dn(K + u) = sqrt(1-m)*nd(u)
899 (power (sub 1 m
) 1//2)
900 (mul (power (sub 1 m
) 1//2)
901 (ftake '%jacobi_nd const m
))))))
902 ((and (alike1 lin
1//2)
905 ;; dn(1/2*K) = (1-m)^(1/4)
912 ;; Should we simplify the inverse elliptic functions into the
913 ;; appropriate incomplete elliptic integral? I think we should leave
914 ;; it, but perhaps allow some way to do that transformation if
917 (def-simplifier inverse_jacobi_sn
(u m
)
919 ;; To numerically evaluate inverse_jacobi_sn (asn), use
921 ;; asn(x,m) = F(asin(x),m)
923 ;; But F(phi,m) = sin(phi)*rf(cos(phi)^2, 1-m*sin(phi)^2,1). Thus
925 ;; asn(x,m) = F(asin(x),m)
926 ;; = x*rf(1-x^2,1-m*x^2,1)
928 ;; I (rtoy) am not 100% about the first identity above for all
929 ;; complex values of x and m, but tests seem to indicate that it
930 ;; produces the correct value as verified by verifying
931 ;; jacobi_sn(inverse_jacobi_sn(x,m),m) = x.
932 (cond ((float-numerical-eval-p u m
)
933 (let ((uu (bigfloat:to
($float u
)))
934 (mm (bigfloat:to
($float m
))))
937 (bigfloat::bf-rf
(bigfloat:to
(- 1 (* uu uu
)))
938 (bigfloat:to
(- 1 (* mm uu uu
)))
940 ((setf args
(complex-float-numerical-eval-p u m
))
941 (destructuring-bind (u m
)
943 (let ((uu (bigfloat:to
($float u
)))
944 (mm (bigfloat:to
($float m
))))
945 (complexify (* uu
(bigfloat::bf-rf
(- 1 (* uu uu
))
948 ((bigfloat-numerical-eval-p u m
)
949 (let ((uu (bigfloat:to u
))
950 (mm (bigfloat:to m
)))
952 (bigfloat::bf-rf
(bigfloat:-
1 (bigfloat:* uu uu
))
953 (bigfloat:-
1 (bigfloat:* mm uu uu
))
955 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
956 (destructuring-bind (u m
)
958 (let ((uu (bigfloat:to u
))
959 (mm (bigfloat:to m
)))
961 (bigfloat::bf-rf
(bigfloat:-
1 (bigfloat:* uu uu
))
962 (bigfloat:-
1 (bigfloat:* mm uu uu
))
968 ;; asn(1,m) = elliptic_kc(m)
969 (ftake '%elliptic_kc m
))
970 ((and (numberp u
) (onep1 (- u
)))
971 ;; asn(-1,m) = -elliptic_kc(m)
972 (mul -
1 (ftake '%elliptic_kc m
)))
974 ;; asn(x,0) = F(asin(x),0) = asin(x)
977 ;; asn(x,1) = F(asin(x),1) = log(tan(pi/4+asin(x)/2))
978 (ftake '%elliptic_f
(ftake '%asin u
) 1))
979 ((and (eq $triginverses
'$all
)
981 (eq (caar u
) '%jacobi_sn
)
982 (alike1 (third u
) m
))
983 ;; inverse_jacobi_sn(sn(u)) = u
989 (def-simplifier inverse_jacobi_cn
(u m
)
991 (cond ((float-numerical-eval-p u m
)
992 ;; Numerically evaluate acn
994 ;; acn(x,m) = F(acos(x),m)
995 (to (elliptic-f (cl:acos
($float u
)) ($float m
))))
996 ((setf args
(complex-float-numerical-eval-p u m
))
997 (destructuring-bind (u m
)
999 (to (elliptic-f (cl:acos
(bigfloat:to
($float u
)))
1000 (bigfloat:to
($float m
))))))
1001 ((bigfloat-numerical-eval-p u m
)
1002 (to (bigfloat::bf-elliptic-f
(bigfloat:acos
(bigfloat:to u
))
1004 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
1005 (destructuring-bind (u m
)
1007 (to (bigfloat::bf-elliptic-f
(bigfloat:acos
(bigfloat:to u
))
1010 ;; asn(x,0) = F(acos(x),0) = acos(x)
1011 (ftake '%elliptic_f
(ftake '%acos u
) 0))
1013 ;; asn(x,1) = F(asin(x),1) = log(tan(pi/2+asin(x)/2))
1014 (ftake '%elliptic_f
(ftake '%acos u
) 1))
1016 (ftake '%elliptic_kc m
))
1019 ((and (eq $triginverses
'$all
)
1021 (eq (caar u
) '%jacobi_cn
)
1022 (alike1 (third u
) m
))
1023 ;; inverse_jacobi_cn(cn(u)) = u
1029 (def-simplifier inverse_jacobi_dn
(u m
)
1031 (cond ((float-numerical-eval-p u m
)
1032 (to (bigfloat::bf-inverse-jacobi-dn
(bigfloat:to
(float u
))
1033 (bigfloat:to
(float m
)))))
1034 ((setf args
(complex-float-numerical-eval-p u m
))
1035 (destructuring-bind (u m
)
1037 (let ((uu (bigfloat:to
($float u
)))
1038 (mm (bigfloat:to
($float m
))))
1039 (to (bigfloat::bf-inverse-jacobi-dn uu mm
)))))
1040 ((bigfloat-numerical-eval-p u m
)
1041 (let ((uu (bigfloat:to u
))
1042 (mm (bigfloat:to m
)))
1043 (to (bigfloat::bf-inverse-jacobi-dn uu mm
))))
1044 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
1045 (destructuring-bind (u m
)
1047 (to (bigfloat::bf-inverse-jacobi-dn
(bigfloat:to u
) (bigfloat:to m
)))))
1049 ;; x = dn(u,1) = sech(u). so u = asech(x)
1052 ;; jacobi_dn(0,m) = 1
1054 ((zerop1 ($ratsimp
(sub u
(power (sub 1 m
) 1//2))))
1055 ;; jacobi_dn(K(m),m) = sqrt(1-m) so
1056 ;; inverse_jacobi_dn(sqrt(1-m),m) = K(m)
1057 (ftake '%elliptic_kc m
))
1058 ((and (eq $triginverses
'$all
)
1060 (eq (caar u
) '%jacobi_dn
)
1061 (alike1 (third u
) m
))
1062 ;; inverse_jacobi_dn(dn(u)) = u
1068 ;;;; Elliptic integrals
1070 (let ((errtol (expt (* 4 flonum-epsilon
) 1/6))
1074 (declare (type flonum errtol c1 c2 c3
))
1076 "Compute Carlson's incomplete or complete elliptic integral of the
1082 RF(x, y, z) = I ----------------------------------- dt
1083 ] SQRT(x + t) SQRT(y + t) SQRT(z + t)
1087 x, y, and z may be complex.
1089 (declare (number x y z
))
1090 (let ((x (coerce x
'(complex flonum
)))
1091 (y (coerce y
'(complex flonum
)))
1092 (z (coerce z
'(complex flonum
))))
1093 (declare (type (complex flonum
) x y z
))
1095 (let* ((mu (/ (+ x y z
) 3))
1096 (x-dev (- 2 (/ (+ mu x
) mu
)))
1097 (y-dev (- 2 (/ (+ mu y
) mu
)))
1098 (z-dev (- 2 (/ (+ mu z
) mu
))))
1099 (when (< (max (abs x-dev
) (abs y-dev
) (abs z-dev
)) errtol
)
1100 (let ((e2 (- (* x-dev y-dev
) (* z-dev z-dev
)))
1101 (e3 (* x-dev y-dev z-dev
)))
1108 (let* ((x-root (sqrt x
))
1111 (lam (+ (* x-root
(+ y-root z-root
)) (* y-root z-root
))))
1112 (setf x
(* (+ x lam
) 1/4))
1113 (setf y
(* (+ y lam
) 1/4))
1114 (setf z
(* (+ z lam
) 1/4))))))))
1116 ;; Elliptic integral of the first kind (Legendre's form):
1122 ;; I ------------------- ds
1124 ;; / SQRT(1 - m SIN (s))
1127 (defun elliptic-f (phi-arg m-arg
)
1128 (flet ((base (phi-arg m-arg
)
1129 (cond ((and (realp m-arg
) (realp phi-arg
))
1130 (let ((phi (float phi-arg
))
1135 ;; F(phi|m) = 1/sqrt(m)*F(theta|1/m)
1137 ;; with sin(theta) = sqrt(m)*sin(phi)
1138 (/ (elliptic-f (cl:asin
(* (sqrt m
) (sin phi
))) (/ m
))
1146 (- (/ (elliptic-f (float (/ pi
2)) m
/m
+1)
1148 (/ (elliptic-f (- (float (/ pi
2)) phi
) m
/m
+1)
1156 1 ;; F(phi,1) = log(sec(phi)+tan(phi))
1157 ;; = log(tan(pi/4+pi/2))
1158 (log (cl:tan
(+ (/ phi
2) (float (/ pi
4))))))
1160 (- (elliptic-f (- phi
) m
)))
1163 (multiple-value-bind (s phi-rem
)
1164 (truncate phi
(float pi
))
1165 (+ (* 2 s
(elliptic-k m
))
1166 (elliptic-f phi-rem m
))))
1168 (let ((sin-phi (sin phi
))
1172 (bigfloat::bf-rf
(* cos-phi cos-phi
)
1173 (* (- 1 (* k sin-phi
))
1174 (+ 1 (* k sin-phi
)))
1177 (+ (* 2 (elliptic-k m
))
1178 (elliptic-f (- phi
(float pi
)) m
)))
1180 (error "Shouldn't happen! Unhandled case in elliptic-f: ~S ~S~%"
1183 (let ((phi (coerce phi-arg
'(complex flonum
)))
1184 (m (coerce m-arg
'(complex flonum
))))
1185 (let ((sin-phi (sin phi
))
1189 (crf (* cos-phi cos-phi
)
1190 (* (- 1 (* k sin-phi
))
1191 (+ 1 (* k sin-phi
)))
1193 ;; Elliptic F is quasi-periodic wrt to z:
1195 ;; F(z|m) = F(z - pi*round(Re(z)/pi)|m) + 2*round(Re(z)/pi)*K(m)
1196 (let ((period (round (realpart phi-arg
) pi
)))
1197 (+ (base (- phi-arg
(* pi period
)) m-arg
)
1201 (bigfloat:to
(elliptic-k m-arg
))))))))
1203 ;; Complete elliptic integral of the first kind
1204 (defun elliptic-k (m)
1212 (- (/ (elliptic-k m
/m
+1)
1214 (/ (elliptic-f 0.0 m
/m
+1)
1221 (intl:gettext
"elliptic_kc: elliptic_kc(1) is undefined.")))
1223 (bigfloat::bf-rf
0.0 (- 1 m
)
1226 (bigfloat::bf-rf
0.0 (- 1 m
)
1229 ;; Elliptic integral of the second kind (Legendre's form):
1235 ;; I SQRT(1 - m SIN (s)) ds
1240 (defun elliptic-e (phi m
)
1241 (declare (type flonum phi m
))
1242 (flet ((base (phi m
)
1250 (let* ((sin-phi (sin phi
))
1253 (y (* (- 1 (* k sin-phi
))
1254 (+ 1 (* k sin-phi
)))))
1256 (bigfloat::bf-rf
(* cos-phi cos-phi
) y
1.0))
1259 (bigfloat::bf-rd
(* cos-phi cos-phi
) y
1.0)))))))))
1260 ;; Elliptic E is quasi-periodic wrt to phi:
1262 ;; E(z|m) = E(z - %pi*round(Re(z)/%pi)|m) + 2*round(Re(z)/%pi)*E(m)
1263 (let ((period (round (realpart phi
) pi
)))
1264 (+ (base (- phi
(* pi period
)) m
)
1265 (* 2 period
(elliptic-ec m
))))))
1268 (defun elliptic-ec (m)
1269 (declare (type flonum m
))
1278 (to (- (bigfloat::bf-rf
0.0 y
1.0)
1280 (bigfloat::bf-rd
0.0 y
1.0))))))))
1283 ;; Define the elliptic integrals for maxima
1285 ;; We use the definitions given in A&S 17.2.6 and 17.2.8. In particular:
1290 ;; F(phi|m) = I ------------------- ds
1292 ;; / SQRT(1 - m SIN (s))
1300 ;; E(phi|m) = I SQRT(1 - m SIN (s)) ds
1305 ;; That is, we do not use the modular angle, alpha, as the second arg;
1306 ;; the parameter m = sin(alpha)^2 is used.
1310 ;; The derivative of F(phi|m) wrt to phi is easy. The derivative wrt
1311 ;; to m is harder. Here is a derivation. Hope I got it right.
1313 ;; diff(integrate(1/sqrt(1-m*sin(x)^2),x,0,phi), m);
1318 ;; I ------------------ dx
1320 ;; / (1 - m SIN (x))
1322 ;; --------------------------
1326 ;; Now use the following relationship that is easily verified:
1329 ;; (1 - m) SIN (x) COS (x) COS(x) SIN(x)
1330 ;; ------------------- = ------------------- - DIFF(-------------------, x)
1332 ;; SQRT(1 - m SIN (x)) SQRT(1 - m SIN (x)) SQRT(1 - m SIN (x))
1335 ;; Now integrate this to get:
1341 ;; (1 - m) I ------------------- dx =
1343 ;; / SQRT(1 - m SIN (x))
1350 ;; + I ------------------- dx
1352 ;; / SQRT(1 - m SIN (x))
1354 ;; COS(PHI) SIN(PHI)
1355 ;; - ---------------------
1357 ;; SQRT(1 - m SIN (PHI))
1359 ;; Use the fact that cos(x)^2 = 1 - sin(x)^2 to show that this
1360 ;; integral on the RHS is:
1363 ;; (1 - m) elliptic_F(PHI, m) + elliptic_E(PHI, m)
1364 ;; -------------------------------------------
1366 ;; So, finally, we have
1371 ;; 2 -- (elliptic_F(PHI, m)) =
1374 ;; elliptic_E(PHI, m) - (1 - m) elliptic_F(PHI, m) COS(PHI) SIN(PHI)
1375 ;; ---------------------------------------------- - ---------------------
1377 ;; SQRT(1 - m SIN (PHI))
1378 ;; ----------------------------------------------------------------------
1381 (defprop %elliptic_f
1384 ;; 1/sqrt(1-m*sin(phi)^2)
1386 ((mplus simp
) 1 ((mtimes simp
) -
1 m
((mexpt simp
) ((%sin simp
) phi
) 2)))
1389 ((mtimes simp
) ((rat simp
) 1 2)
1390 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
1392 ((mtimes simp
) ((mexpt simp
) m -
1)
1393 ((mplus simp
) ((%elliptic_e simp
) phi m
)
1394 ((mtimes simp
) -
1 ((mplus simp
) 1 ((mtimes simp
) -
1 m
))
1395 ((%elliptic_f simp
) phi m
))))
1396 ((mtimes simp
) -
1 ((%cos simp
) phi
) ((%sin simp
) phi
)
1399 ((mtimes simp
) -
1 m
((mexpt simp
) ((%sin simp
) phi
) 2)))
1400 ((rat simp
) -
1 2))))))
1404 ;; The derivative of E(phi|m) wrt to m is much simpler to derive than F(phi|m).
1406 ;; Take the derivative of the definition to get
1411 ;; I ------------------- dx
1413 ;; / SQRT(1 - m SIN (x))
1415 ;; - ---------------------------
1418 ;; It is easy to see that
1423 ;; elliptic_F(PHI, m) - m I ------------------- dx = elliptic_E(PHI, m)
1425 ;; / SQRT(1 - m SIN (x))
1428 ;; So we finally have
1430 ;; d elliptic_E(PHI, m) - elliptic_F(PHI, m)
1431 ;; -- (elliptic_E(PHI, m)) = ---------------------------------------
1434 (defprop %elliptic_e
1436 ;; sqrt(1-m*sin(phi)^2)
1438 ((mplus simp
) 1 ((mtimes simp
) -
1 m
((mexpt simp
) ((%sin simp
) phi
) 2)))
1441 ((mtimes simp
) ((rat simp
) 1 2) ((mexpt simp
) m -
1)
1442 ((mplus simp
) ((%elliptic_e simp
) phi m
)
1443 ((mtimes simp
) -
1 ((%elliptic_f simp
) phi m
)))))
1446 (def-simplifier elliptic_f
(phi m
)
1448 (cond ((float-numerical-eval-p phi m
)
1449 ;; Numerically evaluate it
1450 (to (elliptic-f ($float phi
) ($float m
))))
1451 ((setf args
(complex-float-numerical-eval-p phi m
))
1452 (destructuring-bind (phi m
)
1454 (to (elliptic-f (bigfloat:to
($float phi
))
1455 (bigfloat:to
($float m
))))))
1456 ((bigfloat-numerical-eval-p phi m
)
1457 (to (bigfloat::bf-elliptic-f
(bigfloat:to
($bfloat phi
))
1458 (bigfloat:to
($bfloat m
)))))
1459 ((setf args
(complex-bigfloat-numerical-eval-p phi m
))
1460 (destructuring-bind (phi m
)
1462 (to (bigfloat::bf-elliptic-f
(bigfloat:to
($bfloat phi
))
1463 (bigfloat:to
($bfloat m
))))))
1470 ;; A&S 17.4.21. Let's pick the log tan form. But this
1471 ;; isn't right if we know that abs(phi) > %pi/2, where
1472 ;; elliptic_f is undefined (or infinity).
1473 (cond ((not (eq '$pos
(csign (sub ($abs phi
) (div '$%pi
2)))))
1476 (add (mul '$%pi
(div 1 4))
1479 (merror (intl:gettext
"elliptic_f: elliptic_f(~:M, ~:M) is undefined.")
1481 ((alike1 phi
'((mtimes) ((rat) 1 2) $%pi
))
1482 ;; Complete elliptic integral
1483 (ftake '%elliptic_kc m
))
1488 (def-simplifier elliptic_e
(phi m
)
1490 (cond ((float-numerical-eval-p phi m
)
1491 ;; Numerically evaluate it
1492 (elliptic-e ($float phi
) ($float m
)))
1493 ((complex-float-numerical-eval-p phi m
)
1494 (complexify (bigfloat::bf-elliptic-e
(complex ($float
($realpart phi
)) ($float
($imagpart phi
)))
1495 (complex ($float
($realpart m
)) ($float
($imagpart m
))))))
1496 ((bigfloat-numerical-eval-p phi m
)
1497 (to (bigfloat::bf-elliptic-e
(bigfloat:to
($bfloat phi
))
1498 (bigfloat:to
($bfloat m
)))))
1499 ((setf args
(complex-bigfloat-numerical-eval-p phi m
))
1500 (destructuring-bind (phi m
)
1502 (to (bigfloat::bf-elliptic-e
(bigfloat:to
($bfloat phi
))
1503 (bigfloat:to
($bfloat m
))))))
1510 ;; A&S 17.4.25, but handle periodicity:
1511 ;; elliptic_e(x,m) = elliptic_e(x-%pi*round(x/%pi), m)
1512 ;; + 2*round(x/%pi)*elliptic_ec(m)
1516 ;; elliptic_e(x,1) = sin(phi) + 2*round(x/%pi)*elliptic_ec(m)
1518 (add (ftake '%sin phi
)
1520 (mul (ftake '%round
(div phi
'$%pi
))
1521 (ftake '%elliptic_ec m
)))))
1522 ((alike1 phi
'((mtimes) ((rat) 1 2) $%pi
))
1523 ;; Complete elliptic integral
1524 (ftake '%elliptic_ec m
))
1525 ((and ($numberp phi
)
1526 (let ((r ($round
(div phi
'$%pi
))))
1529 ;; Handle the case where phi is a number where we can apply
1530 ;; the periodicity property without blowing up the
1532 (add (ftake '%elliptic_e
1535 (ftake '%round
(div phi
'$%pi
))))
1538 (mul (ftake '%round
(div phi
'$%pi
))
1539 (ftake '%elliptic_ec m
)))))
1544 ;; Complete elliptic integrals
1546 ;; elliptic_kc(m) = elliptic_f(%pi/2, m)
1548 ;; elliptic_ec(m) = elliptic_e(%pi/2, m)
1551 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1553 ;;; We support a simplim%function. The function is looked up in simplimit and
1554 ;;; handles specific values of the function.
1556 (defprop %elliptic_kc simplim%elliptic_kc simplim%function
)
1558 (defun simplim%elliptic_kc
(expr var val
)
1559 ;; Look for the limit of the argument
1560 (let ((m (limit (cadr expr
) var val
'think
)))
1562 ;; For an argument 1 return $infinity.
1565 ;; All other cases are handled by the simplifier of the function.
1566 (simplify (list '(%elliptic_kc
) m
))))))
1568 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1570 (def-simplifier elliptic_kc
(m)
1573 ;; elliptic_kc(1) is complex infinity. Maxima can not handle
1574 ;; infinities correctly, throw a Maxima error.
1576 (intl:gettext
"elliptic_kc: elliptic_kc(~:M) is undefined.")
1578 ((float-numerical-eval-p m
)
1579 ;; Numerically evaluate it
1580 (to (elliptic-k ($float m
))))
1581 ((complex-float-numerical-eval-p m
)
1582 (complexify (bigfloat::bf-elliptic-k
(complex ($float
($realpart m
)) ($float
($imagpart m
))))))
1583 ((setf args
(complex-bigfloat-numerical-eval-p m
))
1584 (destructuring-bind (m)
1586 (to (bigfloat::bf-elliptic-k
(bigfloat:to
($bfloat m
))))))
1588 '((mtimes) ((rat) 1 2) $%pi
))
1590 ;; http://functions.wolfram.com/EllipticIntegrals/EllipticK/03/01/
1592 ;; elliptic_kc(1/2) = 8*%pi^(3/2)/gamma(-1/4)^2
1593 (div (mul 8 (power '$%pi
(div 3 2)))
1594 (power (gm (div -
1 4)) 2)))
1596 ;; elliptic_kc(-1) = gamma(1/4)^2/(4*sqrt(2*%pi))
1597 (div (power (gm (div 1 4)) 2)
1598 (mul 4 (power (mul 2 '$%pi
) 1//2))))
1599 ((alike1 m
(add 17 (mul -
12 (power 2 1//2))))
1600 ;; elliptic_kc(17-12*sqrt(2)) = 2*(2+sqrt(2))*%pi^(3/2)/gamma(-1/4)^2
1601 (div (mul 2 (mul (add 2 (power 2 1//2))
1602 (power '$%pi
(div 3 2))))
1603 (power (gm (div -
1 4)) 2)))
1608 (defprop %elliptic_kc
1613 ((mplus) ((%elliptic_ec
) m
)
1616 ((mplus) 1 ((mtimes) -
1 m
))))
1617 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
1621 (def-simplifier elliptic_ec
(m)
1623 (cond ((float-numerical-eval-p m
)
1624 ;; Numerically evaluate it
1625 (elliptic-ec ($float m
)))
1626 ((setf args
(complex-float-numerical-eval-p m
))
1627 (destructuring-bind (m)
1629 (complexify (bigfloat::bf-elliptic-ec
(bigfloat:to
($float m
))))))
1630 ((setf args
(complex-bigfloat-numerical-eval-p m
))
1631 (destructuring-bind (m)
1633 (to (bigfloat::bf-elliptic-ec
(bigfloat:to
($bfloat m
))))))
1634 ;; Some special cases we know about.
1636 '((mtimes) ((rat) 1 2) $%pi
))
1640 ;; elliptic_ec(1/2). Use the identity
1642 ;; elliptic_ec(z)*elliptic_kc(1-z) - elliptic_kc(z)*elliptic_kc(1-z)
1643 ;; + elliptic_ec(1-z)*elliptic_kc(z) = %pi/2;
1645 ;; Let z = 1/2 to get
1647 ;; %pi^(3/2)*'elliptic_ec(1/2)/gamma(3/4)^2-%pi^3/(4*gamma(3/4)^4) = %pi/2
1649 ;; since we know that elliptic_kc(1/2) = %pi^(3/2)/(2*gamma(3/4)^2). Hence
1652 ;; = (2*%pi*gamma(3/4)^4+%pi^3)/(4*%pi^(3/2)*gamma(3/4)^2)
1653 ;; = gamma(3/4)^2/(2*sqrt(%pi))+%pi^(3/2)/(4*gamma(3/4)^2)
1655 (add (div (power (ftake '%gamma
(div 3 4)) 2)
1656 (mul 2 (power '$%pi
1//2)))
1657 (div (power '$%pi
(div 3 2))
1658 (mul 4 (power (ftake '%gamma
(div 3 4)) 2)))))
1660 ;; elliptic_ec(-1). Use the identity
1661 ;; http://functions.wolfram.com/08.01.17.0002.01
1664 ;; elliptic_ec(z) = sqrt(1 - z)*elliptic_ec(z/(z-1))
1666 ;; Let z = -1 to get
1668 ;; elliptic_ec(-1) = sqrt(2)*elliptic_ec(1/2)
1670 ;; Should we expand out elliptic_ec(1/2) using the above result?
1672 (ftake '%elliptic_ec
1//2)))
1677 (defprop %elliptic_ec
1679 ((mtimes) ((rat) 1 2)
1680 ((mplus) ((%elliptic_ec
) m
)
1681 ((mtimes) -
1 ((%elliptic_kc
)
1687 ;; Elliptic integral of the third kind:
1694 ;; PI(n;phi|m) = I ----------------------------------- ds
1696 ;; / SQRT(1 - m SIN (s)) (1 - n SIN (s))
1699 ;; As with E and F, we do not use the modular angle alpha but the
1700 ;; parameter m = sin(alpha)^2.
1702 (def-simplifier elliptic_pi
(n phi m
)
1705 ((float-numerical-eval-p n phi m
)
1706 ;; Numerically evaluate it
1707 (elliptic-pi ($float n
) ($float phi
) ($float m
)))
1708 ((setf args
(complex-float-numerical-eval-p n phi m
))
1709 (destructuring-bind (n phi m
)
1711 (elliptic-pi (bigfloat:to
($float n
))
1712 (bigfloat:to
($float phi
))
1713 (bigfloat:to
($float m
)))))
1714 ((bigfloat-numerical-eval-p n phi m
)
1715 (to (bigfloat::bf-elliptic-pi
(bigfloat:to n
)
1718 ((setq args
(complex-bigfloat-numerical-eval-p n phi m
))
1719 (destructuring-bind (n phi m
)
1721 (to (bigfloat::bf-elliptic-pi
(bigfloat:to n
)
1725 (ftake '%elliptic_f phi m
))
1727 ;; 3 cases depending on n < 1, n > 1, or n = 1.
1728 (let ((s (asksign (add -
1 n
))))
1731 (div (ftake '%atanh
(mul (power (add n -
1) 1//2)
1733 (power (add n -
1) 1//2)))
1735 (div (ftake '%atan
(mul (power (sub 1 n
) 1//2)
1737 (power (sub 1 n
) 1//2)))
1739 (ftake '%tan phi
)))))
1744 ;; Complete elliptic-pi. That is phi = %pi/2. Then
1746 ;; = Rf(0, 1-m,1) + Rj(0,1-m,1-n)*n/3;
1747 (defun elliptic-pi-complete (n m
)
1748 (to (bigfloat:+ (bigfloat::bf-rf
0 (- 1 m
) 1)
1749 (bigfloat:* 1/3 n
(bigfloat::bf-rj
0 (- 1 m
) 1 (- 1 n
))))))
1751 ;; To compute elliptic_pi for all z, we use the property
1752 ;; (http://functions.wolfram.com/08.06.16.0002.01)
1754 ;; elliptic_pi(n, z + %pi*k, m)
1755 ;; = 2*k*elliptic_pi(n, %pi/2, m) + elliptic_pi(n, z, m)
1757 ;; So we are left with computing the integral for 0 <= z < %pi. Using
1758 ;; Carlson's formulation produces the wrong values for %pi/2 < z <
1759 ;; %pi. How to do that?
1763 ;; I(a,b) = integrate(1/(1-n*sin(x)^2)/sqrt(1 - m*sin(x)^2), x, a, b)
1765 ;; That is, I(a,b) is the integral for the elliptic_pi function but
1766 ;; with a lower limit of a and an upper limit of b.
1768 ;; Then, we want to compute I(0, z), with %pi <= z < %pi. Let w = z +
1769 ;; %pi/2, 0 <= w < %pi/2. Then
1771 ;; I(0, w+%pi/2) = I(0, %pi/2) + I(%pi/2, w+%pi/2)
1773 ;; To evaluate I(%pi/2, w+%pi/2), use a change of variables:
1775 ;; changevar('integrate(1/(1-n*sin(x)^2)/sqrt(1 - m*sin(x)^2), x, %pi/2, w + %pi/2),
1778 ;; = integrate(-1/(sqrt(1-m*sin(u)^2)*(1-n*sin(u)^2)),u,%pi/2-w,%pi/2)
1779 ;; = I(%pi/2-w,%pi/2)
1780 ;; = I(0,%pi/2) - I(0,%pi/2-w)
1784 ;; I(0,%pi/2+w) = 2*I(0,%pi/2) - I(0,%pi/2-w)
1786 ;; This allows us to compute the general result with 0 <= z < %pi
1788 ;; I(0, k*%pi + z) = 2*k*I(0,%pi/2) + I(0,z);
1790 ;; If 0 <= z < %pi/2, then the we are done. If %pi/2 <= z < %pi, let
1791 ;; z = w+%pi/2. Then
1793 ;; I(0,z) = 2*I(0,%pi/2) - I(0,%pi/2-w)
1795 ;; Or, since w = z-%pi/2:
1797 ;; I(0,z) = 2*I(0,%pi/2) - I(0,%pi-z)
1799 (defun elliptic-pi (n phi m
)
1800 ;; elliptic_pi(n, -phi, m) = -elliptic_pi(n, phi, m). That is, it
1801 ;; is an odd function of phi.
1802 (when (minusp (realpart phi
))
1803 (return-from elliptic-pi
(- (elliptic-pi n
(- phi
) m
))))
1805 ;; Note: Carlson's DRJ has n defined as the negative of the n given
1807 (flet ((base (n phi m
)
1808 ;; elliptic_pi(n,phi,m) =
1809 ;; sin(phi)*Rf(cos(phi)^2, 1-m*sin(phi)^2, 1)
1810 ;; - (-n / 3) * sin(phi)^3
1811 ;; * Rj(cos(phi)^2, 1-m*sin(phi)^2, 1, 1-n*sin(phi)^2)
1816 (k2sin (* (- 1 (* k sin-phi
))
1817 (+ 1 (* k sin-phi
)))))
1818 (- (* sin-phi
(bigfloat::bf-rf
(expt cos-phi
2) k2sin
1.0))
1819 (* (/ nn
3) (expt sin-phi
3)
1820 (bigfloat::bf-rj
(expt cos-phi
2) k2sin
1.0
1821 (- 1 (* n
(expt sin-phi
2)))))))))
1822 ;; FIXME: Reducing the arg by pi has significant round-off.
1823 ;; Consider doing something better.
1824 (let* ((cycles (round (realpart phi
) pi
))
1825 (rem (- phi
(* cycles pi
))))
1826 (let ((complete (elliptic-pi-complete n m
)))
1827 (to (+ (* 2 cycles complete
)
1828 (base n rem m
)))))))
1830 ;;; Deriviatives from functions.wolfram.com
1831 ;;; http://functions.wolfram.com/EllipticIntegrals/EllipticPi3/20/
1832 (defprop %elliptic_pi
1834 ;Derivative wrt first argument
1835 ((mtimes) ((rat) 1 2)
1836 ((mexpt) ((mplus) m
((mtimes) -
1 n
)) -
1)
1837 ((mexpt) ((mplus) -
1 n
) -
1)
1839 ((mtimes) ((mexpt) n -
1)
1840 ((mplus) ((mtimes) -
1 m
) ((mexpt) n
2))
1841 ((%elliptic_pi
) n z m
))
1843 ((mtimes) ((mplus) m
((mtimes) -
1 n
)) ((mexpt) n -
1)
1844 ((%elliptic_f
) z m
))
1845 ((mtimes) ((rat) -
1 2) n
1847 ((mplus) 1 ((mtimes) -
1 m
((mexpt) ((%sin
) z
) 2)))
1850 ((mplus) 1 ((mtimes) -
1 n
((mexpt) ((%sin
) z
) 2)))
1852 ((%sin
) ((mtimes) 2 z
)))))
1853 ;derivative wrt second argument
1856 ((mplus) 1 ((mtimes) -
1 m
((mexpt) ((%sin
) z
) 2)))
1859 ((mplus) 1 ((mtimes) -
1 n
((mexpt) ((%sin
) z
) 2))) -
1))
1860 ;Derivative wrt third argument
1861 ((mtimes) ((rat) 1 2)
1862 ((mexpt) ((mplus) ((mtimes) -
1 m
) n
) -
1)
1863 ((mplus) ((%elliptic_pi
) n z m
)
1864 ((mtimes) ((mexpt) ((mplus) -
1 m
) -
1)
1865 ((%elliptic_e
) z m
))
1866 ((mtimes) ((rat) -
1 2) ((mexpt) ((mplus) -
1 m
) -
1) m
1868 ((mplus) 1 ((mtimes) -
1 m
((mexpt) ((%sin
) z
) 2)))
1870 ((%sin
) ((mtimes) 2 z
))))))
1873 (in-package #-gcl
#:bigfloat
#+gcl
"BIGFLOAT")
1874 ;; Translation of Jim FitzSimons' bigfloat implementation of elliptic
1875 ;; integrals from http://www.getnet.com/~cherry/elliptbf3.mac.
1877 ;; The algorithms are based on B.C. Carlson's "Numerical Computation
1878 ;; of Real or Complex Elliptic Integrals". These are updated to the
1879 ;; algorithms in Journal of Computational and Applied Mathematics 118
1880 ;; (2000) 71-85 "Reduction Theorems for Elliptic Integrands with the
1881 ;; Square Root of two quadritic factors"
1884 ;; NOTE: Despite the names indicating these are for bigfloat numbers,
1885 ;; the algorithms and routines are generic and will work with floats
1888 (defun bferrtol (&rest args
)
1889 ;; Compute error tolerance as sqrt(2^(-fpprec)). Not sure this is
1890 ;; quite right, but it makes the routines more accurate as fpprec
1892 (sqrt (reduce #'min
(mapcar #'(lambda (x)
1893 (if (rationalp (realpart x
))
1894 maxima
::flonum-epsilon
1898 ;; rc(x,y) = integrate(1/2*(t+x)^(-1/2)/(t+y), t, 0, inf)
1900 ;; log(x) = (x-1)*rc(((1+x)/2)^2, x), x > 0
1901 ;; asin(x) = x * rc(1-x^2, 1), |x|<= 1
1902 ;; acos(x) = sqrt(1-x^2)*rc(x^2,1), 0 <= x <=1
1903 ;; atan(x) = x * rc(1,1+x^2)
1904 ;; asinh(x) = x * rc(1+x^2,1)
1905 ;; acosh(x) = sqrt(x^2-1) * rc(x^2,1), x >= 1
1906 ;; atanh(x) = x * rc(1,1-x^2), |x|<=1
1910 xn z w a an pwr4 n epslon lambda sn s
)
1911 (cond ((and (zerop (imagpart yn
))
1912 (minusp (realpart yn
)))
1916 (setf w
(sqrt (/ x xn
))))
1921 (setf a
(/ (+ xn yn yn
) 3))
1922 (setf epslon
(/ (abs (- a xn
)) (bferrtol x y
)))
1926 (loop while
(> (* epslon pwr4
) (abs an
))
1928 (setf pwr4
(/ pwr4
4))
1929 (setf lambda
(+ (* 2 (sqrt xn
) (sqrt yn
)) yn
))
1930 (setf an
(/ (+ an lambda
) 4))
1931 (setf xn
(/ (+ xn lambda
) 4))
1932 (setf yn
(/ (+ yn lambda
) 4))
1934 ;; c2=3/10,c3=1/7,c4=3/8,c5=9/22,c6=159/208,c7=9/8
1935 (setf sn
(/ (* pwr4
(- z a
)) an
))
1936 (setf s
(* sn sn
(+ 3/10
1941 (* sn
9/8))))))))))))
1947 ;; See https://dlmf.nist.gov/19.16.E5:
1949 ;; rd(x,y,z) = integrate(3/2/sqrt(t+x)/sqrt(t+y)/sqrt(t+z)/(t+z), t, 0, inf)
1952 ;; E(K) = rf(0, 1-K^2, 1) - (K^2/3)*rd(0,1-K^2,1)
1954 ;; B = integrate(s^2/sqrt(1-s^4), s, 0 ,1)
1955 ;; = beta(3/4,1/2)/4
1956 ;; = sqrt(%pi)*gamma(3/4)/gamma(1/4)
1959 (defun bf-rd (x y z
)
1963 (a (/ (+ xn yn
(* 3 zn
)) 5))
1964 (epslon (/ (max (abs (- a xn
))
1972 xnroot ynroot znroot lam
)
1973 (loop while
(> (* power4 epslon
) (abs an
))
1975 (setf xnroot
(sqrt xn
))
1976 (setf ynroot
(sqrt yn
))
1977 (setf znroot
(sqrt zn
))
1978 (setf lam
(+ (* xnroot ynroot
)
1981 (setf sigma
(+ sigma
(/ power4
1982 (* znroot
(+ zn lam
)))))
1983 (setf power4
(* power4
1/4))
1984 (setf xn
(* (+ xn lam
) 1/4))
1985 (setf yn
(* (+ yn lam
) 1/4))
1986 (setf zn
(* (+ zn lam
) 1/4))
1987 (setf an
(* (+ an lam
) 1/4))
1989 ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26
1990 (let* ((xndev (/ (* (- a x
) power4
) an
))
1991 (yndev (/ (* (- a y
) power4
) an
))
1992 (zndev (- (* (+ xndev yndev
) 1/3)))
1993 (ee2 (- (* xndev yndev
) (* 6 zndev zndev
)))
1994 (ee3 (* (- (* 3 xndev yndev
)
1997 (ee4 (* 3 (- (* xndev yndev
) (* zndev zndev
)) zndev zndev
))
1998 (ee5 (* xndev yndev zndev zndev zndev
))
2006 (* -
1/16 ee2 ee2 ee2
)
2009 (* 45/272 ee2 ee2 ee3
)
2010 (* -
9/68 (+ (* ee2 ee5
) (* ee3 ee4
))))))
2015 ;; See https://dlmf.nist.gov/19.16.E1
2017 ;; rf(x,y,z) = 1/2*integrate(1/(sqrt(t+x)*sqrt(t+y)*sqrt(t+z)), t, 0, inf);
2020 (defun bf-rf (x y z
)
2024 (a (/ (+ xn yn zn
) 3))
2025 (epslon (/ (max (abs (- a xn
))
2032 xnroot ynroot znroot lam
)
2033 (loop while
(> (* power4 epslon
) (abs an
))
2035 (setf xnroot
(sqrt xn
))
2036 (setf ynroot
(sqrt yn
))
2037 (setf znroot
(sqrt zn
))
2038 (setf lam
(+ (* xnroot ynroot
)
2041 (setf power4
(* power4
1/4))
2042 (setf xn
(* (+ xn lam
) 1/4))
2043 (setf yn
(* (+ yn lam
) 1/4))
2044 (setf zn
(* (+ zn lam
) 1/4))
2045 (setf an
(* (+ an lam
) 1/4))
2047 ;; c1=-3/14,c2=1/6,c3=9/88,c4=9/22,c5=-3/22,c6=-9/52,c7=3/26
2048 (let* ((xndev (/ (* (- a x
) power4
) an
))
2049 (yndev (/ (* (- a y
) power4
) an
))
2050 (zndev (- (+ xndev yndev
)))
2051 (ee2 (- (* xndev yndev
) (* 6 zndev zndev
)))
2052 (ee3 (* xndev yndev zndev
))
2057 (* -
3/44 ee2 ee3
))))
2060 (defun bf-rj1 (x y z p
)
2071 (a (/ (+ xn yn zn pn pn
) 5))
2072 (epslon (/ (max (abs (- a xn
))
2076 (bferrtol x y z p
)))
2078 xnroot ynroot znroot pnroot lam dn
)
2079 (loop while
(> (* power4 epslon
) (abs an
))
2081 (setf xnroot
(sqrt xn
))
2082 (setf ynroot
(sqrt yn
))
2083 (setf znroot
(sqrt zn
))
2084 (setf pnroot
(sqrt pn
))
2085 (setf lam
(+ (* xnroot ynroot
)
2088 (setf dn
(* (+ pnroot xnroot
)
2091 (setf sigma
(+ sigma
2093 (bf-rc 1 (+ 1 (/ en
(* dn dn
)))))
2095 (setf power4
(* power4
1/4))
2097 (setf xn
(* (+ xn lam
) 1/4))
2098 (setf yn
(* (+ yn lam
) 1/4))
2099 (setf zn
(* (+ zn lam
) 1/4))
2100 (setf pn
(* (+ pn lam
) 1/4))
2101 (setf an
(* (+ an lam
) 1/4))
2103 (let* ((xndev (/ (* (- a x
) power4
) an
))
2104 (yndev (/ (* (- a y
) power4
) an
))
2105 (zndev (/ (* (- a z
) power4
) an
))
2106 (pndev (* -
0.5 (+ xndev yndev zndev
)))
2107 (ee2 (+ (* xndev yndev
)
2110 (* -
3 pndev pndev
)))
2111 (ee3 (+ (* xndev yndev zndev
)
2113 (* 4 pndev pndev pndev
)))
2114 (ee4 (* (+ (* 2 xndev yndev zndev
)
2116 (* 3 pndev pndev pndev
))
2118 (ee5 (* xndev yndev zndev pndev pndev
))
2126 (* -
1/16 ee2 ee2 ee2
)
2129 (* 45/272 ee2 ee2 ee3
)
2130 (* -
9/68 (+ (* ee2 ee5
) (* ee3 ee4
))))))
2133 (sqrt (* an an an
)))))))
2135 (defun bf-rj (x y z p
)
2140 (cond ((and (and (zerop (imagpart xn
)) (>= (realpart xn
) 0))
2141 (and (zerop (imagpart yn
)) (>= (realpart yn
) 0))
2142 (and (zerop (imagpart zn
)) (>= (realpart zn
) 0))
2143 (and (zerop (imagpart qn
)) (> (realpart qn
) 0)))
2144 (destructuring-bind (xn yn zn
)
2145 (sort (list xn yn zn
) #'<)
2146 (let* ((pn (+ yn
(* (- zn yn
) (/ (- yn xn
) (+ yn qn
)))))
2147 (s (- (* (- pn yn
) (bf-rj1 xn yn zn pn
))
2148 (* 3 (bf-rf xn yn zn
)))))
2149 (setf s
(+ s
(* 3 (sqrt (/ (* xn yn zn
)
2150 (+ (* xn zn
) (* pn qn
))))
2151 (bf-rc (+ (* xn zn
) (* pn qn
)) (* pn qn
)))))
2154 (bf-rj1 x y z p
)))))
2156 (defun bf-rg (x y z
)
2158 (+ (* z
(bf-rf x y z
))
2163 (sqrt (/ (* x y
) z
)))))
2165 ;; elliptic_f(phi,m) = sin(phi)*rf(cos(phi)^2, 1-m*sin(phi)^2,1)
2166 (defun bf-elliptic-f (phi m
)
2167 (flet ((base (phi m
)
2169 ;; F(z|1) = log(tan(z/2+%pi/4))
2170 (log (tan (+ (/ phi
2) (/ (%pi phi
) 4)))))
2174 (* s
(bf-rf (* c c
) (- 1 (* m s s
)) 1)))))))
2175 ;; Handle periodicity (see elliptic-f)
2176 (let* ((bfpi (%pi phi
))
2177 (period (round (realpart phi
) bfpi
)))
2178 (+ (base (- phi
(* bfpi period
)) m
)
2181 (* 2 period
(bf-elliptic-k m
)))))))
2183 ;; elliptic_kc(k) = rf(0, 1-k^2,1)
2186 ;; elliptic_kc(m) = rf(0, 1-m,1)
2188 (defun bf-elliptic-k (m)
2190 (if (maxima::$bfloatp m
)
2191 (maxima::$bfloat
(maxima::div
'maxima
::$%pi
2))
2192 (float (/ pi
2) 1e0
)))
2195 (intl:gettext
"elliptic_kc: elliptic_kc(1) is undefined.")))
2197 (bf-rf 0 (- 1 m
) 1))))
2199 ;; elliptic_e(phi, k) = sin(phi)*rf(cos(phi)^2,1-k^2*sin(phi)^2,1)
2200 ;; - (k^2/3)*sin(phi)^3*rd(cos(phi)^2, 1-k^2*sin(phi)^2,1)
2204 ;; elliptic_e(phi, m) = sin(phi)*rf(cos(phi)^2,1-m*sin(phi)^2,1)
2205 ;; - (m/3)*sin(phi)^3*rd(cos(phi)^2, 1-m*sin(phi)^2,1)
2207 (defun bf-elliptic-e (phi m
)
2208 (flet ((base (phi m
)
2209 (let* ((s (sin phi
))
2212 (s2 (- 1 (* m s s
))))
2213 (- (* s
(bf-rf c2 s2
1))
2214 (* (/ m
3) (* s s s
) (bf-rd c2 s2
1))))))
2215 ;; Elliptic E is quasi-periodic wrt to phi:
2217 ;; E(z|m) = E(z - %pi*round(Re(z)/%pi)|m) + 2*round(Re(z)/%pi)*E(m)
2218 (let* ((bfpi (%pi phi
))
2219 (period (round (realpart phi
) bfpi
)))
2220 (+ (base (- phi
(* bfpi period
)) m
)
2221 (* 2 period
(bf-elliptic-ec m
))))))
2224 ;; elliptic_ec(k) = rf(0,1-k^2,1) - (k^2/3)*rd(0,1-k^2,1);
2227 ;; elliptic_ec(m) = rf(0,1-m,1) - (m/3)*rd(0,1-m,1);
2229 (defun bf-elliptic-ec (m)
2231 (if (typep m
'bigfloat
)
2232 (bigfloat (maxima::$bfloat
(maxima::div
'maxima
::$%pi
2)))
2233 (float (/ pi
2) 1e0
)))
2235 (if (typep m
'bigfloat
)
2241 (* m
1/3 (bf-rd 0 m1
1)))))))
2243 (defun bf-elliptic-pi-complete (n m
)
2244 (+ (bf-rf 0 (- 1 m
) 1)
2245 (* 1/3 n
(bf-rj 0 (- 1 m
) 1 (- 1 n
)))))
2247 (defun bf-elliptic-pi (n phi m
)
2248 ;; Note: Carlson's DRJ has n defined as the negative of the n given
2250 (flet ((base (n phi m
)
2255 (k2sin (* (- 1 (* k sin-phi
))
2256 (+ 1 (* k sin-phi
)))))
2257 (- (* sin-phi
(bf-rf (expt cos-phi
2) k2sin
1.0))
2258 (* (/ nn
3) (expt sin-phi
3)
2259 (bf-rj (expt cos-phi
2) k2sin
1.0
2260 (- 1 (* n
(expt sin-phi
2)))))))))
2261 ;; FIXME: Reducing the arg by pi has significant round-off.
2262 ;; Consider doing something better.
2263 (let* ((bf-pi (%pi
(realpart phi
)))
2264 (cycles (round (realpart phi
) bf-pi
))
2265 (rem (- phi
(* cycles bf-pi
))))
2266 (let ((complete (bf-elliptic-pi-complete n m
)))
2267 (+ (* 2 cycles complete
)
2270 ;; Compute inverse_jacobi_sn, for float or bigfloat args.
2271 (defun bf-inverse-jacobi-sn (u m
)
2272 (* u
(bf-rf (- 1 (* u u
))
2276 ;; Compute inverse_jacobi_dn. We use the following identity
2277 ;; from Gradshteyn & Ryzhik, 8.153.6
2279 ;; w = dn(z|m) = cn(sqrt(m)*z, 1/m)
2281 ;; Solve for z to get
2283 ;; z = inverse_jacobi_dn(w,m)
2284 ;; = 1/sqrt(m) * inverse_jacobi_cn(w, 1/m)
2285 (defun bf-inverse-jacobi-dn (w m
)
2289 ;; jacobi_dn(x,1) = sech(x) so the inverse is asech(x)
2290 (maxima::take
'(maxima::%asech
) (maxima::to w
)))
2292 ;; We should do something better to make sure that things
2293 ;; that should be real are real.
2294 (/ (to (maxima::take
'(maxima::%inverse_jacobi_cn
)
2296 (maxima::to
(/ m
))))
2299 (in-package :maxima
)
2301 ;; Define Carlson's elliptic integrals.
2303 (def-simplifier carlson_rc
(x y
)
2306 (flet ((floatify (z)
2307 ;; If z is a complex rational, convert to a
2308 ;; complex double-float. Otherwise, leave it as
2309 ;; is. If we don't do this, %i is handled as
2310 ;; #c(0 1), which makes bf-rc use single-float
2311 ;; arithmetic instead of the desired
2313 (if (and (complexp z
) (rationalp (realpart z
)))
2314 (complex (float (realpart z
))
2315 (float (imagpart z
)))
2317 (to (bigfloat::bf-rc
(floatify (bigfloat:to x
))
2318 (floatify (bigfloat:to y
)))))))
2319 ;; See comments from bf-rc
2320 (cond ((float-numerical-eval-p x y
)
2321 (calc ($float x
) ($float y
)))
2322 ((bigfloat-numerical-eval-p x y
)
2323 (calc ($bfloat x
) ($bfloat y
)))
2324 ((setf args
(complex-float-numerical-eval-p x y
))
2325 (destructuring-bind (x y
)
2327 (calc ($float x
) ($float y
))))
2328 ((setf args
(complex-bigfloat-numerical-eval-p x y
))
2329 (destructuring-bind (x y
)
2331 (calc ($bfloat x
) ($bfloat y
))))
2337 (alike1 y
(div 1 4)))
2342 ;; rc(2,1) = 1/2*integrate(1/sqrt(t+2)/(t+1), t, 0, inf)
2343 ;; = (log(sqrt(2)+1)-log(sqrt(2)-1))/2
2344 ;; ratsimp(logcontract(%)),algebraic:
2345 ;; = -log(3-2^(3/2))/2
2346 ;; = -log(sqrt(3-2^(3/2)))
2347 ;; = -log(sqrt(2)-1)
2348 ;; = log(1/(sqrt(2)-1))
2349 ;; ratsimp(%),algebraic;
2351 (ftake '%log
(add 1 (power 2 1//2))))
2352 ((and (alike x
'$%i
)
2353 (alike y
(add 1 '$%i
)))
2354 ;; rc(%i, %i+1) = 1/2*integrate(1/sqrt(t+%i)/(t+%i+1), t, 0, inf)
2355 ;; = %pi/2-atan((-1)^(1/4))
2356 ;; ratsimp(logcontract(ratsimp(rectform(%o42)))),algebraic;
2357 ;; = (%i*log(3-2^(3/2))+%pi)/4
2358 ;; = (%i*log(3-2^(3/2)))/4+%pi/4
2359 ;; = %i*log(sqrt(3-2^(3/2)))/2+%pi/4
2361 ;; = %pi/4 + %i*log(sqrt(2)-1)/2
2365 (ftake '%log
(sub (power 2 1//2) 1)))))
2368 ;; rc(0,%i) = 1/2*integrate(1/(sqrt(t)*(t+%i)), t, 0, inf)
2369 ;; = -((sqrt(2)*%i-sqrt(2))*%pi)/4
2370 ;; = ((1-%i)*%pi)/2^(3/2)
2371 (div (mul (sub 1 '$%i
)
2375 (eq ($sign
($realpart x
)) '$pos
))
2376 ;; carlson_rc(x,x) = 1/2*integrate(1/sqrt(t+x)/(t+x), t, 0, inf)
2379 ((and (alike1 x
(power (div (add 1 y
) 2) 2))
2380 (eq ($sign
($realpart y
)) '$pos
))
2381 ;; Rc(((1+x)/2)^2,x) = log(x)/(x-1) for x > 0.
2383 ;; This is done by looking at Rc(x,y) and seeing if
2384 ;; ((1+y)/2)^2 is the same as x.
2385 (div (ftake '%log y
)
2390 (def-simplifier carlson_rd
(x y z
)
2392 (flet ((calc (x y z
)
2393 (to (bigfloat::bf-rd
(bigfloat:to x
)
2396 ;; See https://dlmf.nist.gov/19.20.E18
2397 (cond ((and (eql x
1)
2404 ;; Rd(x,x,x) = x^(-3/2)
2405 (power x
(div -
3 2)))
2408 ;; Rd(0,y,y) = 3/4*%pi*y^(-3/2)
2411 (power y
(div -
3 2))))
2413 ;; Rd(x,y,y) = 3/(2*(y-x))*(Rc(x, y) - sqrt(x)/y)
2414 (mul (div 3 (mul 2 (sub y x
)))
2415 (sub (ftake '%carlson_rc x y
)
2419 ;; Rd(x,x,z) = 3/(z-x)*(Rc(z,x) - 1/sqrt(z))
2420 (mul (div 3 (sub z x
))
2421 (sub (ftake '%carlson_rc z x
)
2422 (div 1 (power z
1//2)))))
2428 ;; Rd(0,2,1) = 3*(gamma(3/4)^2)/sqrt(2*%pi)
2429 ;; See https://dlmf.nist.gov/19.20.E22.
2431 ;; But that's the same as
2432 ;; 3*sqrt(%pi)*gamma(3/4)/gamma(1/4). We can see this by
2433 ;; taking the ratio to get
2434 ;; gamma(1/4)*gamma(3/4)/sqrt(2)*%pi. But
2435 ;; gamma(1/4)*gamma(3/4) = beta(1/4,3/4) = sqrt(2)*%pi.
2436 ;; Hence, the ratio is 1.
2438 ;; Note also that Rd(x,y,z) = Rd(y,x,z)
2441 (div (ftake '%gamma
(div 3 4))
2442 (ftake '%gamma
(div 1 4)))))
2443 ((and (or (eql x
0) (eql y
0))
2445 ;; 1/3*m*Rd(0,1-m,1) = K(m) - E(m).
2446 ;; See https://dlmf.nist.gov/19.25.E1
2448 ;; Thus, Rd(0,y,1) = 3/(1-y)*(K(1-y) - E(1-y))
2450 ;; Note that Rd(x,y,z) = Rd(y,x,z).
2451 (let ((m (sub 1 y
)))
2453 (sub (ftake '%elliptic_kc m
)
2454 (ftake '%elliptic_ec m
)))))
2459 ;; 1/3*m*(1-m)*Rd(0,1,1-m) = E(m) - (1-m)*K(m)
2460 ;; See https://dlmf.nist.gov/19.25.E1
2463 ;; Rd(0,1,z) = 3/(z*(1-z))*(E(1-z) - z*K(1-z))
2464 ;; Recall that Rd(x,y,z) = Rd(y,x,z).
2465 (mul (div 3 (mul z
(sub 1 z
)))
2466 (sub (ftake '%elliptic_ec
(sub 1 z
))
2468 (ftake '%elliptic_kc
(sub 1 z
))))))
2469 ((float-numerical-eval-p x y z
)
2470 (calc ($float x
) ($float y
) ($float z
)))
2471 ((bigfloat-numerical-eval-p x y z
)
2472 (calc ($bfloat x
) ($bfloat y
) ($bfloat z
)))
2473 ((setf args
(complex-float-numerical-eval-p x y z
))
2474 (destructuring-bind (x y z
)
2476 (calc ($float x
) ($float y
) ($float z
))))
2477 ((setf args
(complex-bigfloat-numerical-eval-p x y z
))
2478 (destructuring-bind (x y z
)
2480 (calc ($bfloat x
) ($bfloat y
) ($bfloat z
))))
2484 (def-simplifier carlson_rf
(x y z
)
2486 (flet ((calc (x y z
)
2487 (to (bigfloat::bf-rf
(bigfloat:to x
)
2490 ;; See https://dlmf.nist.gov/19.20.i
2491 (cond ((and (alike1 x y
)
2493 ;; Rf(x,x,x) = x^(-1/2)
2497 ;; Rf(0,y,y) = 1/2*%pi*y^(-1/2)
2501 (ftake '%carlson_rc x y
))
2502 ((some #'(lambda (args)
2503 (destructuring-bind (x y z
)
2514 ;; Rf(0,1,2) = (gamma(1/4))^2/(4*sqrt(2*%pi))
2516 ;; And Rf is symmetric in all the args, so check every
2517 ;; permutation too. This could probably be simplified
2518 ;; without consing all the lists, but I'm lazy.
2519 (div (power (ftake '%gamma
(div 1 4)) 2)
2520 (mul 4 (power (mul 2 '$%pi
) 1//2))))
2521 ((some #'(lambda (args)
2522 (destructuring-bind (x y z
)
2524 (and (alike1 x
'$%i
)
2525 (alike1 y
(mul -
1 '$%i
))
2534 ;; = 1/2*integrate(1/sqrt(t^2+1)/sqrt(t),t,0,inf)
2535 ;; = beta(1/4,1/4)/4;
2537 ;; = gamma(1/4)^2/(4*sqrt(%pi))
2539 ;; Rf is symmetric, so check all the permutations too.
2540 (div (power (ftake '%gamma
(div 1 4)) 2)
2541 (mul 4 (power '$%pi
1//2))))
2543 (some #'(lambda (args)
2544 (destructuring-bind (x y z
)
2546 ;; Check that x = 0 and z = 1, and
2557 ;; Rf(0,1-m,1) = elliptic_kc(m).
2558 ;; See https://dlmf.nist.gov/19.25.E1
2559 (ftake '%elliptic_kc
(sub 1 args
)))
2560 ((some #'(lambda (args)
2561 (destructuring-bind (x y z
)
2563 (and (alike1 x
'$%i
)
2564 (alike1 y
(mul -
1 '$%i
))
2573 ;; = 1/2*integrate(1/sqrt(t^2+1)/sqrt(t),t,0,inf)
2574 ;; = beta(1/4,1/4)/4;
2576 ;; = gamma(1/4)^2/(4*sqrt(%pi))
2578 ;; Rf is symmetric, so check all the permutations too.
2579 (div (power (ftake '%gamma
(div 1 4)) 2)
2580 (mul 4 (power '$%pi
1//2))))
2581 ((float-numerical-eval-p x y z
)
2582 (calc ($float x
) ($float y
) ($float z
)))
2583 ((bigfloat-numerical-eval-p x y z
)
2584 (calc ($bfloat x
) ($bfloat y
) ($bfloat z
)))
2585 ((setf args
(complex-float-numerical-eval-p x y z
))
2586 (destructuring-bind (x y z
)
2588 (calc ($float x
) ($float y
) ($float z
))))
2589 ((setf args
(complex-bigfloat-numerical-eval-p x y z
))
2590 (destructuring-bind (x y z
)
2592 (calc ($bfloat x
) ($bfloat y
) ($bfloat z
))))
2596 (def-simplifier carlson_rj
(x y z p
)
2598 (flet ((calc (x y z p
)
2599 (to (bigfloat::bf-rj
(bigfloat:to x
)
2603 ;; See https://dlmf.nist.gov/19.20.iii
2604 (cond ((and (alike1 x y
)
2607 ;; Rj(x,x,x,x) = x^(-3/2)
2608 (power x
(div -
3 2)))
2610 ;; Rj(x,y,z,z) = Rd(x,y,z)
2611 (ftake '%carlson_rd x y z
))
2614 ;; Rj(0,y,y,p) = 3*%pi/(2*(y*sqrt(p)+p*sqrt(y)))
2617 (add (mul y
(power p
1//2))
2618 (mul p
(power y
1//2))))))
2620 ;; Rj(x,y,y,p) = 3/(p-y)*(Rc(x,y) - Rc(x,p))
2621 (mul (div 3 (sub p y
))
2622 (sub (ftake '%carlson_rc x y
)
2623 (ftake '%carlson_rc x p
))))
2626 ;; Rj(x,y,y,y) = Rd(x,y,y)
2627 (ftake '%carlson_rd x y y
))
2628 ((float-numerical-eval-p x y z p
)
2629 (calc ($float x
) ($float y
) ($float z
) ($float p
)))
2630 ((bigfloat-numerical-eval-p x y z p
)
2631 (calc ($bfloat x
) ($bfloat y
) ($bfloat z
) ($bfloat p
)))
2632 ((setf args
(complex-float-numerical-eval-p x y z p
))
2633 (destructuring-bind (x y z p
)
2635 (calc ($float x
) ($float y
) ($float z
) ($float p
))))
2636 ((setf args
(complex-bigfloat-numerical-eval-p x y z p
))
2637 (destructuring-bind (x y z p
)
2639 (calc ($bfloat x
) ($bfloat y
) ($bfloat z
) ($bfloat p
))))
2643 ;;; Other Jacobian elliptic functions
2645 ;; jacobi_ns(u,m) = 1/jacobi_sn(u,m)
2649 ((mtimes) -
1 ((%jacobi_cn
) u m
) ((%jacobi_dn
) u m
)
2650 ((mexpt) ((%jacobi_sn
) u m
) -
2))
2652 ((mtimes) -
1 ((mexpt) ((%jacobi_sn
) u m
) -
2)
2654 ((mtimes) ((rat) 1 2)
2655 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
2656 ((mexpt) ((%jacobi_cn
) u m
) 2)
2658 ((mtimes) ((rat) 1 2) ((mexpt) m -
1)
2659 ((%jacobi_cn
) u m
) ((%jacobi_dn
) u m
)
2662 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
2663 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
2667 (def-simplifier jacobi_ns
(u m
)
2670 ((float-numerical-eval-p u m
)
2671 (to (bigfloat:/ (bigfloat::sn
(bigfloat:to
($float u
))
2672 (bigfloat:to
($float m
))))))
2673 ((setf args
(complex-float-numerical-eval-p u m
))
2674 (destructuring-bind (u m
)
2676 (to (bigfloat:/ (bigfloat::sn
(bigfloat:to
($float u
))
2677 (bigfloat:to
($float m
)))))))
2678 ((bigfloat-numerical-eval-p u m
)
2679 (let ((uu (bigfloat:to
($bfloat u
)))
2680 (mm (bigfloat:to
($bfloat m
))))
2681 (to (bigfloat:/ (bigfloat::sn uu mm
)))))
2682 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
2683 (destructuring-bind (u m
)
2685 (let ((uu (bigfloat:to
($bfloat u
)))
2686 (mm (bigfloat:to
($bfloat m
))))
2687 (to (bigfloat:/ (bigfloat::sn uu mm
))))))
2695 (dbz-err1 'jacobi_ns
))
2696 ((and $trigsign
(mminusp* u
))
2698 (neg (ftake* '%jacobi_ns
(neg u
) m
)))
2701 (member (caar u
) '(%inverse_jacobi_sn
2712 %inverse_jacobi_dc
))
2713 (alike1 (third u
) m
))
2714 (cond ((eq (caar u
) '%inverse_jacobi_ns
)
2717 ;; Express in terms of sn:
2719 (div 1 (ftake '%jacobi_sn u m
)))))
2720 ;; A&S 16.20 (Jacobi's Imaginary transformation)
2721 ((and $%iargs
(multiplep u
'$%i
))
2722 ;; ns(i*u) = 1/sn(i*u) = -i/sc(u,m1) = -i*cs(u,m1)
2724 (ftake* '%jacobi_cs
(coeff u
'$%i
1) (add 1 (neg m
))))))
2725 ((setq coef
(kc-arg2 u m
))
2728 ;; ns(m*K+u) = 1/sn(m*K+u)
2730 (destructuring-bind (lin const
)
2732 (cond ((integerp lin
)
2735 ;; ns(4*m*K+u) = ns(u)
2738 (dbz-err1 'jacobi_ns
)
2739 (ftake '%jacobi_ns const m
)))
2741 ;; ns(4*m*K + K + u) = ns(K+u) = dc(u)
2745 (ftake '%jacobi_dc const m
)))
2747 ;; ns(4*m*K+2*K + u) = ns(2*K+u) = -ns(u)
2748 ;; ns(2*K) = infinity
2750 (dbz-err1 'jacobi_ns
)
2751 (neg (ftake '%jacobi_ns const m
))))
2753 ;; ns(4*m*K+3*K+u) = ns(2*K + K + u) = -ns(K+u) = -dc(u)
2757 (neg (ftake '%jacobi_dc const m
))))))
2758 ((and (alike1 lin
1//2)
2760 (div 1 (ftake '%jacobi_sn u m
)))
2767 ;; jacobi_nc(u,m) = 1/jacobi_cn(u,m)
2771 ((mtimes) ((mexpt) ((%jacobi_cn
) u m
) -
2)
2772 ((%jacobi_dn
) u m
) ((%jacobi_sn
) u m
))
2774 ((mtimes) -
1 ((mexpt) ((%jacobi_cn
) u m
) -
2)
2776 ((mtimes) ((rat) -
1 2)
2777 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
2778 ((%jacobi_cn
) u m
) ((mexpt) ((%jacobi_sn
) u m
) 2))
2779 ((mtimes) ((rat) -
1 2) ((mexpt) m -
1)
2780 ((%jacobi_dn
) u m
) ((%jacobi_sn
) u m
)
2782 ((mtimes) -
1 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
2783 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
)) m
)))))))
2786 (def-simplifier jacobi_nc
(u m
)
2789 ((float-numerical-eval-p u m
)
2790 (to (bigfloat:/ (bigfloat::cn
(bigfloat:to
($float u
))
2791 (bigfloat:to
($float m
))))))
2792 ((setf args
(complex-float-numerical-eval-p u m
))
2793 (destructuring-bind (u m
)
2795 (to (bigfloat:/ (bigfloat::cn
(bigfloat:to
($float u
))
2796 (bigfloat:to
($float m
)))))))
2797 ((bigfloat-numerical-eval-p u m
)
2798 (let ((uu (bigfloat:to
($bfloat u
)))
2799 (mm (bigfloat:to
($bfloat m
))))
2800 (to (bigfloat:/ (bigfloat::cn uu mm
)))))
2801 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
2802 (destructuring-bind (u m
)
2804 (let ((uu (bigfloat:to
($bfloat u
)))
2805 (mm (bigfloat:to
($bfloat m
))))
2806 (to (bigfloat:/ (bigfloat::cn uu mm
))))))
2815 ((and $trigsign
(mminusp* u
))
2817 (ftake* '%jacobi_nc
(neg u
) m
))
2820 (member (caar u
) '(%inverse_jacobi_sn
2831 %inverse_jacobi_dc
))
2832 (alike1 (third u
) m
))
2833 (cond ((eq (caar u
) '%inverse_jacobi_nc
)
2836 ;; Express in terms of cn:
2838 (div 1 (ftake '%jacobi_cn u m
)))))
2839 ;; A&S 16.20 (Jacobi's Imaginary transformation)
2840 ((and $%iargs
(multiplep u
'$%i
))
2841 ;; nc(i*u) = 1/cn(i*u) = 1/nc(u,1-m) = cn(u,1-m)
2842 (ftake* '%jacobi_cn
(coeff u
'$%i
1) (add 1 (neg m
))))
2843 ((setq coef
(kc-arg2 u m
))
2848 (destructuring-bind (lin const
)
2850 (cond ((integerp lin
)
2853 ;; nc(4*m*K+u) = nc(u)
2857 (ftake '%jacobi_nc const m
)))
2859 ;; nc(4*m*K+K+u) = nc(K+u) = -ds(u)/sqrt(1-m)
2862 (dbz-err1 'jacobi_nc
)
2863 (neg (div (ftake '%jacobi_ds const m
)
2864 (power (sub 1 m
) 1//2)))))
2866 ;; nc(4*m*K+2*K+u) = nc(2*K+u) = -nc(u)
2870 (neg (ftake '%jacobi_nc const m
))))
2872 ;; nc(4*m*K+3*K+u) = nc(3*K+u) = nc(2*K+K+u) =
2873 ;; -nc(K+u) = ds(u)/sqrt(1-m)
2875 ;; nc(3*K) = infinity
2877 (dbz-err1 'jacobi_nc
)
2878 (div (ftake '%jacobi_ds const m
)
2879 (power (sub 1 m
) 1//2))))))
2880 ((and (alike1 1//2 lin
)
2882 (div 1 (ftake '%jacobi_cn u m
)))
2889 ;; jacobi_nd(u,m) = 1/jacobi_dn(u,m)
2893 ((mtimes) m
((%jacobi_cn
) u m
)
2894 ((mexpt) ((%jacobi_dn
) u m
) -
2) ((%jacobi_sn
) u m
))
2896 ((mtimes) -
1 ((mexpt) ((%jacobi_dn
) u m
) -
2)
2898 ((mtimes) ((rat) -
1 2)
2899 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
2901 ((mexpt) ((%jacobi_sn
) u m
) 2))
2902 ((mtimes) ((rat) -
1 2) ((%jacobi_cn
) u m
)
2906 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
2907 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
2911 (def-simplifier jacobi_nd
(u m
)
2914 ((float-numerical-eval-p u m
)
2915 (to (bigfloat:/ (bigfloat::dn
(bigfloat:to
($float u
))
2916 (bigfloat:to
($float m
))))))
2917 ((setf args
(complex-float-numerical-eval-p u m
))
2918 (destructuring-bind (u m
)
2920 (to (bigfloat:/ (bigfloat::dn
(bigfloat:to
($float u
))
2921 (bigfloat:to
($float m
)))))))
2922 ((bigfloat-numerical-eval-p u m
)
2923 (let ((uu (bigfloat:to
($bfloat u
)))
2924 (mm (bigfloat:to
($bfloat m
))))
2925 (to (bigfloat:/ (bigfloat::dn uu mm
)))))
2926 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
2927 (destructuring-bind (u m
)
2929 (let ((uu (bigfloat:to
($bfloat u
)))
2930 (mm (bigfloat:to
($bfloat m
))))
2931 (to (bigfloat:/ (bigfloat::dn uu mm
))))))
2940 ((and $trigsign
(mminusp* u
))
2942 (ftake* '%jacobi_nd
(neg u
) m
))
2945 (member (caar u
) '(%inverse_jacobi_sn
2956 %inverse_jacobi_dc
))
2957 (alike1 (third u
) m
))
2958 (cond ((eq (caar u
) '%inverse_jacobi_nd
)
2961 ;; Express in terms of dn:
2963 (div 1 (ftake '%jacobi_dn u m
)))))
2964 ;; A&S 16.20 (Jacobi's Imaginary transformation)
2965 ((and $%iargs
(multiplep u
'$%i
))
2966 ;; nd(i*u) = 1/dn(i*u) = 1/dc(u,1-m) = cd(u,1-m)
2967 (ftake* '%jacobi_cd
(coeff u
'$%i
1) (add 1 (neg m
))))
2968 ((setq coef
(kc-arg2 u m
))
2971 (destructuring-bind (lin const
)
2973 (cond ((integerp lin
)
2977 ;; nd(2*m*K+u) = nd(u)
2981 (ftake '%jacobi_nd const m
)))
2983 ;; nd(2*m*K+K+u) = nd(K+u) = dn(u)/sqrt(1-m)
2984 ;; nd(K) = 1/sqrt(1-m)
2986 (power (sub 1 m
) -
1//2)
2987 (div (ftake '%jacobi_nd const m
)
2988 (power (sub 1 m
) 1//2))))))
2995 ;; jacobi_sc(u,m) = jacobi_sn/jacobi_cn
2999 ((mtimes) ((mexpt) ((%jacobi_cn
) u m
) -
2)
3003 ((mtimes) ((mexpt) ((%jacobi_cn
) u m
) -
1)
3005 ((mtimes) ((rat) 1 2)
3006 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3007 ((mexpt) ((%jacobi_cn
) u m
) 2)
3009 ((mtimes) ((rat) 1 2) ((mexpt) m -
1)
3010 ((%jacobi_cn
) u m
) ((%jacobi_dn
) u m
)
3013 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3014 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3016 ((mtimes) -
1 ((mexpt) ((%jacobi_cn
) u m
) -
2)
3019 ((mtimes) ((rat) -
1 2)
3020 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3022 ((mexpt) ((%jacobi_sn
) u m
) 2))
3023 ((mtimes) ((rat) -
1 2) ((mexpt) m -
1)
3024 ((%jacobi_dn
) u m
) ((%jacobi_sn
) u m
)
3027 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3028 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3032 (def-simplifier jacobi_sc
(u m
)
3035 ((float-numerical-eval-p u m
)
3036 (let ((fu (bigfloat:to
($float u
)))
3037 (fm (bigfloat:to
($float m
))))
3038 (to (bigfloat:/ (bigfloat::sn fu fm
) (bigfloat::cn fu fm
)))))
3039 ((setf args
(complex-float-numerical-eval-p u m
))
3040 (destructuring-bind (u m
)
3042 (let ((fu (bigfloat:to
($float u
)))
3043 (fm (bigfloat:to
($float m
))))
3044 (to (bigfloat:/ (bigfloat::sn fu fm
) (bigfloat::cn fu fm
))))))
3045 ((bigfloat-numerical-eval-p u m
)
3046 (let ((uu (bigfloat:to
($bfloat u
)))
3047 (mm (bigfloat:to
($bfloat m
))))
3048 (to (bigfloat:/ (bigfloat::sn uu mm
)
3049 (bigfloat::cn uu mm
)))))
3050 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
3051 (destructuring-bind (u m
)
3053 (let ((uu (bigfloat:to
($bfloat u
)))
3054 (mm (bigfloat:to
($bfloat m
))))
3055 (to (bigfloat:/ (bigfloat::sn uu mm
)
3056 (bigfloat::cn uu mm
))))))
3065 ((and $trigsign
(mminusp* u
))
3067 (neg (ftake* '%jacobi_sc
(neg u
) m
)))
3070 (member (caar u
) '(%inverse_jacobi_sn
3081 %inverse_jacobi_dc
))
3082 (alike1 (third u
) m
))
3083 (cond ((eq (caar u
) '%inverse_jacobi_sc
)
3086 ;; Express in terms of sn and cn
3087 ;; sc(x) = sn(x)/cn(x)
3088 (div (ftake '%jacobi_sn u m
)
3089 (ftake '%jacobi_cn u m
)))))
3090 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3091 ((and $%iargs
(multiplep u
'$%i
))
3092 ;; sc(i*u) = sn(i*u)/cn(i*u) = i*sc(u,m1)/nc(u,m1) = i*sn(u,m1)
3094 (ftake* '%jacobi_sn
(coeff u
'$%i
1) (add 1 (neg m
)))))
3095 ((setq coef
(kc-arg2 u m
))
3097 ;; sc(2*m*K+u) = sc(u)
3098 (destructuring-bind (lin const
)
3100 (cond ((integerp lin
)
3103 ;; sc(2*m*K+ u) = sc(u)
3107 (ftake '%jacobi_sc const m
)))
3109 ;; sc(2*m*K + K + u) = sc(K+u)= - cs(u)/sqrt(1-m)
3112 (dbz-err1 'jacobi_sc
)
3114 (div (ftake* '%jacobi_cs const m
)
3115 (power (sub 1 m
) 1//2)))))))
3116 ((and (alike1 lin
1//2)
3118 ;; From A&S 16.3.3 and 16.5.2:
3119 ;; sc(1/2*K) = 1/(1-m)^(1/4)
3120 (power (sub 1 m
) (div -
1 4)))
3127 ;; jacobi_sd(u,m) = jacobi_sn/jacobi_dn
3131 ((mtimes) ((%jacobi_cn
) u m
)
3132 ((mexpt) ((%jacobi_dn
) u m
) -
2))
3135 ((mtimes) ((mexpt) ((%jacobi_dn
) u m
) -
1)
3137 ((mtimes) ((rat) 1 2)
3138 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3139 ((mexpt) ((%jacobi_cn
) u m
) 2)
3141 ((mtimes) ((rat) 1 2) ((mexpt) m -
1)
3142 ((%jacobi_cn
) u m
) ((%jacobi_dn
) u m
)
3145 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3146 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3148 ((mtimes) -
1 ((mexpt) ((%jacobi_dn
) u m
) -
2)
3151 ((mtimes) ((rat) -
1 2)
3152 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3154 ((mexpt) ((%jacobi_sn
) u m
) 2))
3155 ((mtimes) ((rat) -
1 2) ((%jacobi_cn
) u m
)
3159 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3160 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3164 (def-simplifier jacobi_sd
(u m
)
3167 ((float-numerical-eval-p u m
)
3168 (let ((fu (bigfloat:to
($float u
)))
3169 (fm (bigfloat:to
($float m
))))
3170 (to (bigfloat:/ (bigfloat::sn fu fm
) (bigfloat::dn fu fm
)))))
3171 ((setf args
(complex-float-numerical-eval-p u m
))
3172 (destructuring-bind (u m
)
3174 (let ((fu (bigfloat:to
($float u
)))
3175 (fm (bigfloat:to
($float m
))))
3176 (to (bigfloat:/ (bigfloat::sn fu fm
) (bigfloat::dn fu fm
))))))
3177 ((bigfloat-numerical-eval-p u m
)
3178 (let ((uu (bigfloat:to
($bfloat u
)))
3179 (mm (bigfloat:to
($bfloat m
))))
3180 (to (bigfloat:/ (bigfloat::sn uu mm
)
3181 (bigfloat::dn uu mm
)))))
3182 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
3183 (destructuring-bind (u m
)
3185 (let ((uu (bigfloat:to
($bfloat u
)))
3186 (mm (bigfloat:to
($bfloat m
))))
3187 (to (bigfloat:/ (bigfloat::sn uu mm
)
3188 (bigfloat::dn uu mm
))))))
3197 ((and $trigsign
(mminusp* u
))
3199 (neg (ftake* '%jacobi_sd
(neg u
) m
)))
3202 (member (caar u
) '(%inverse_jacobi_sn
3213 %inverse_jacobi_dc
))
3214 (alike1 (third u
) m
))
3215 (cond ((eq (caar u
) '%inverse_jacobi_sd
)
3218 ;; Express in terms of sn and dn
3219 (div (ftake '%jacobi_sn u m
)
3220 (ftake '%jacobi_dn u m
)))))
3221 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3222 ((and $%iargs
(multiplep u
'$%i
))
3223 ;; sd(i*u) = sn(i*u)/dn(i*u) = i*sc(u,m1)/dc(u,m1) = i*sd(u,m1)
3225 (ftake* '%jacobi_sd
(coeff u
'$%i
1) (add 1 (neg m
)))))
3226 ((setq coef
(kc-arg2 u m
))
3228 ;; sd(4*m*K+u) = sd(u)
3229 (destructuring-bind (lin const
)
3231 (cond ((integerp lin
)
3234 ;; sd(4*m*K+u) = sd(u)
3238 (ftake '%jacobi_sd const m
)))
3240 ;; sd(4*m*K+K+u) = sd(K+u) = cn(u)/sqrt(1-m)
3241 ;; sd(K) = 1/sqrt(m1)
3243 (power (sub 1 m
) 1//2)
3244 (div (ftake '%jacobi_cn const m
)
3245 (power (sub 1 m
) 1//2))))
3247 ;; sd(4*m*K+2*K+u) = sd(2*K+u) = -sd(u)
3251 (neg (ftake '%jacobi_sd const m
))))
3253 ;; sd(4*m*K+3*K+u) = sd(3*K+u) = sd(2*K+K+u) =
3254 ;; -sd(K+u) = -cn(u)/sqrt(1-m)
3255 ;; sd(3*K) = -1/sqrt(m1)
3257 (neg (power (sub 1 m
) -
1//2))
3258 (neg (div (ftake '%jacobi_cn const m
)
3259 (power (sub 1 m
) 1//2)))))))
3260 ((and (alike1 lin
1//2)
3262 ;; jacobi_sn/jacobi_dn
3263 (div (ftake '%jacobi_sn
3265 (ftake '%elliptic_kc m
))
3269 (ftake '%elliptic_kc m
))
3277 ;; jacobi_cs(u,m) = jacobi_cn/jacobi_sn
3281 ((mtimes) -
1 ((%jacobi_dn
) u m
)
3282 ((mexpt) ((%jacobi_sn
) u m
) -
2))
3285 ((mtimes) -
1 ((%jacobi_cn
) u m
)
3286 ((mexpt) ((%jacobi_sn
) u m
) -
2)
3288 ((mtimes) ((rat) 1 2)
3289 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3290 ((mexpt) ((%jacobi_cn
) u m
) 2)
3292 ((mtimes) ((rat) 1 2) ((mexpt) m -
1)
3293 ((%jacobi_cn
) u m
) ((%jacobi_dn
) u m
)
3296 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3297 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3299 ((mtimes) ((mexpt) ((%jacobi_sn
) u m
) -
1)
3301 ((mtimes) ((rat) -
1 2)
3302 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3304 ((mexpt) ((%jacobi_sn
) u m
) 2))
3305 ((mtimes) ((rat) -
1 2) ((mexpt) m -
1)
3306 ((%jacobi_dn
) u m
) ((%jacobi_sn
) u m
)
3309 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3310 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3314 (def-simplifier jacobi_cs
(u m
)
3317 ((float-numerical-eval-p u m
)
3318 (let ((fu (bigfloat:to
($float u
)))
3319 (fm (bigfloat:to
($float m
))))
3320 (to (bigfloat:/ (bigfloat::cn fu fm
) (bigfloat::sn fu fm
)))))
3321 ((setf args
(complex-float-numerical-eval-p u m
))
3322 (destructuring-bind (u m
)
3324 (let ((fu (bigfloat:to
($float u
)))
3325 (fm (bigfloat:to
($float m
))))
3326 (to (bigfloat:/ (bigfloat::cn fu fm
) (bigfloat::sn fu fm
))))))
3327 ((bigfloat-numerical-eval-p u m
)
3328 (let ((uu (bigfloat:to
($bfloat u
)))
3329 (mm (bigfloat:to
($bfloat m
))))
3330 (to (bigfloat:/ (bigfloat::cn uu mm
)
3331 (bigfloat::sn uu mm
)))))
3332 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
3333 (destructuring-bind (u m
)
3335 (let ((uu (bigfloat:to
($bfloat u
)))
3336 (mm (bigfloat:to
($bfloat m
))))
3337 (to (bigfloat:/ (bigfloat::cn uu mm
)
3338 (bigfloat::sn uu mm
))))))
3346 (dbz-err1 'jacobi_cs
))
3347 ((and $trigsign
(mminusp* u
))
3349 (neg (ftake* '%jacobi_cs
(neg u
) m
)))
3352 (member (caar u
) '(%inverse_jacobi_sn
3363 %inverse_jacobi_dc
))
3364 (alike1 (third u
) m
))
3365 (cond ((eq (caar u
) '%inverse_jacobi_cs
)
3368 ;; Express in terms of cn an sn
3369 (div (ftake '%jacobi_cn u m
)
3370 (ftake '%jacobi_sn u m
)))))
3371 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3372 ((and $%iargs
(multiplep u
'$%i
))
3373 ;; cs(i*u) = cn(i*u)/sn(i*u) = -i*nc(u,m1)/sc(u,m1) = -i*ns(u,m1)
3375 (ftake* '%jacobi_ns
(coeff u
'$%i
1) (add 1 (neg m
))))))
3376 ((setq coef
(kc-arg2 u m
))
3379 ;; cs(2*m*K + u) = cs(u)
3380 (destructuring-bind (lin const
)
3382 (cond ((integerp lin
)
3385 ;; cs(2*m*K + u) = cs(u)
3388 (dbz-err1 'jacobi_cs
)
3389 (ftake '%jacobi_cs const m
)))
3391 ;; cs(K+u) = -sqrt(1-m)*sc(u)
3395 (neg (mul (power (sub 1 m
) 1//2)
3396 (ftake '%jacobi_sc const m
)))))))
3397 ((and (alike1 lin
1//2)
3401 (ftake '%jacobi_sc
(mul 1//2
3402 (ftake '%elliptic_kc m
))
3410 ;; jacobi_cd(u,m) = jacobi_cn/jacobi_dn
3414 ((mtimes) ((mplus) -
1 m
)
3415 ((mexpt) ((%jacobi_dn
) u m
) -
2)
3419 ((mtimes) -
1 ((%jacobi_cn
) u m
)
3420 ((mexpt) ((%jacobi_dn
) u m
) -
2)
3422 ((mtimes) ((rat) -
1 2)
3423 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3425 ((mexpt) ((%jacobi_sn
) u m
) 2))
3426 ((mtimes) ((rat) -
1 2) ((%jacobi_cn
) u m
)
3430 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3431 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3433 ((mtimes) ((mexpt) ((%jacobi_dn
) u m
) -
1)
3435 ((mtimes) ((rat) -
1 2)
3436 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3438 ((mexpt) ((%jacobi_sn
) u m
) 2))
3439 ((mtimes) ((rat) -
1 2) ((mexpt) m -
1)
3440 ((%jacobi_dn
) u m
) ((%jacobi_sn
) u m
)
3443 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3444 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3448 (def-simplifier jacobi_cd
(u m
)
3451 ((float-numerical-eval-p u m
)
3452 (let ((fu (bigfloat:to
($float u
)))
3453 (fm (bigfloat:to
($float m
))))
3454 (to (bigfloat:/ (bigfloat::cn fu fm
) (bigfloat::dn fu fm
)))))
3455 ((setf args
(complex-float-numerical-eval-p u m
))
3456 (destructuring-bind (u m
)
3458 (let ((fu (bigfloat:to
($float u
)))
3459 (fm (bigfloat:to
($float m
))))
3460 (to (bigfloat:/ (bigfloat::cn fu fm
) (bigfloat::dn fu fm
))))))
3461 ((bigfloat-numerical-eval-p u m
)
3462 (let ((uu (bigfloat:to
($bfloat u
)))
3463 (mm (bigfloat:to
($bfloat m
))))
3464 (to (bigfloat:/ (bigfloat::cn uu mm
) (bigfloat::dn uu mm
)))))
3465 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
3466 (destructuring-bind (u m
)
3468 (let ((uu (bigfloat:to
($bfloat u
)))
3469 (mm (bigfloat:to
($bfloat m
))))
3470 (to (bigfloat:/ (bigfloat::cn uu mm
) (bigfloat::dn uu mm
))))))
3479 ((and $trigsign
(mminusp* u
))
3481 (ftake* '%jacobi_cd
(neg u
) m
))
3484 (member (caar u
) '(%inverse_jacobi_sn
3495 %inverse_jacobi_dc
))
3496 (alike1 (third u
) m
))
3497 (cond ((eq (caar u
) '%inverse_jacobi_cd
)
3500 ;; Express in terms of cn and dn
3501 (div (ftake '%jacobi_cn u m
)
3502 (ftake '%jacobi_dn u m
)))))
3503 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3504 ((and $%iargs
(multiplep u
'$%i
))
3505 ;; cd(i*u) = cn(i*u)/dn(i*u) = nc(u,m1)/dc(u,m1) = nd(u,m1)
3506 (ftake* '%jacobi_nd
(coeff u
'$%i
1) (add 1 (neg m
))))
3507 ((setf coef
(kc-arg2 u m
))
3510 (destructuring-bind (lin const
)
3512 (cond ((integerp lin
)
3515 ;; cd(4*m*K + u) = cd(u)
3519 (ftake '%jacobi_cd const m
)))
3521 ;; cd(4*m*K + K + u) = cd(K+u) = -sn(u)
3525 (neg (ftake '%jacobi_sn const m
))))
3527 ;; cd(4*m*K + 2*K + u) = cd(2*K+u) = -cd(u)
3531 (neg (ftake '%jacobi_cd const m
))))
3533 ;; cd(4*m*K + 3*K + u) = cd(2*K + K + u) =
3538 (ftake '%jacobi_sn const m
)))))
3539 ((and (alike1 lin
1//2)
3541 ;; jacobi_cn/jacobi_dn
3542 (div (ftake '%jacobi_cn
3544 (ftake '%elliptic_kc m
))
3548 (ftake '%elliptic_kc m
))
3557 ;; jacobi_ds(u,m) = jacobi_dn/jacobi_sn
3561 ((mtimes) -
1 ((%jacobi_cn
) u m
)
3562 ((mexpt) ((%jacobi_sn
) u m
) -
2))
3565 ((mtimes) -
1 ((%jacobi_dn
) u m
)
3566 ((mexpt) ((%jacobi_sn
) u m
) -
2)
3568 ((mtimes) ((rat) 1 2)
3569 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3570 ((mexpt) ((%jacobi_cn
) u m
) 2)
3572 ((mtimes) ((rat) 1 2) ((mexpt) m -
1)
3573 ((%jacobi_cn
) u m
) ((%jacobi_dn
) u m
)
3576 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3577 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3579 ((mtimes) ((mexpt) ((%jacobi_sn
) u m
) -
1)
3581 ((mtimes) ((rat) -
1 2)
3582 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3584 ((mexpt) ((%jacobi_sn
) u m
) 2))
3585 ((mtimes) ((rat) -
1 2) ((%jacobi_cn
) u m
)
3589 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3590 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3594 (def-simplifier jacobi_ds
(u m
)
3597 ((float-numerical-eval-p u m
)
3598 (let ((fu (bigfloat:to
($float u
)))
3599 (fm (bigfloat:to
($float m
))))
3600 (to (bigfloat:/ (bigfloat::dn fu fm
) (bigfloat::sn fu fm
)))))
3601 ((setf args
(complex-float-numerical-eval-p u m
))
3602 (destructuring-bind (u m
)
3604 (let ((fu (bigfloat:to
($float u
)))
3605 (fm (bigfloat:to
($float m
))))
3606 (to (bigfloat:/ (bigfloat::dn fu fm
) (bigfloat::sn fu fm
))))))
3607 ((bigfloat-numerical-eval-p u m
)
3608 (let ((uu (bigfloat:to
($bfloat u
)))
3609 (mm (bigfloat:to
($bfloat m
))))
3610 (to (bigfloat:/ (bigfloat::dn uu mm
)
3611 (bigfloat::sn uu mm
)))))
3612 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
3613 (destructuring-bind (u m
)
3615 (let ((uu (bigfloat:to
($bfloat u
)))
3616 (mm (bigfloat:to
($bfloat m
))))
3617 (to (bigfloat:/ (bigfloat::dn uu mm
)
3618 (bigfloat::sn uu mm
))))))
3626 (dbz-err1 'jacobi_ds
))
3627 ((and $trigsign
(mminusp* u
))
3628 (neg (ftake* '%jacobi_ds
(neg u
) m
)))
3631 (member (caar u
) '(%inverse_jacobi_sn
3642 %inverse_jacobi_dc
))
3643 (alike1 (third u
) m
))
3644 (cond ((eq (caar u
) '%inverse_jacobi_ds
)
3647 ;; Express in terms of dn and sn
3648 (div (ftake '%jacobi_dn u m
)
3649 (ftake '%jacobi_sn u m
)))))
3650 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3651 ((and $%iargs
(multiplep u
'$%i
))
3652 ;; ds(i*u) = dn(i*u)/sn(i*u) = -i*dc(u,m1)/sc(u,m1) = -i*ds(u,m1)
3654 (ftake* '%jacobi_ds
(coeff u
'$%i
1) (add 1 (neg m
))))))
3655 ((setf coef
(kc-arg2 u m
))
3657 (destructuring-bind (lin const
)
3659 (cond ((integerp lin
)
3662 ;; ds(4*m*K + u) = ds(u)
3665 (dbz-err1 'jacobi_ds
)
3666 (ftake '%jacobi_ds const m
)))
3668 ;; ds(4*m*K + K + u) = ds(K+u) = sqrt(1-m)*nc(u)
3669 ;; ds(K) = sqrt(1-m)
3671 (power (sub 1 m
) 1//2)
3672 (mul (power (sub 1 m
) 1//2)
3673 (ftake '%jacobi_nc const m
))))
3675 ;; ds(4*m*K + 2*K + u) = ds(2*K+u) = -ds(u)
3678 (dbz-err1 'jacobi_ds
)
3679 (neg (ftake '%jacobi_ds const m
))))
3681 ;; ds(4*m*K + 3*K + u) = ds(2*K + K + u) =
3682 ;; -ds(K+u) = -sqrt(1-m)*nc(u)
3683 ;; ds(3*K) = -sqrt(1-m)
3685 (neg (power (sub 1 m
) 1//2))
3686 (neg (mul (power (sub 1 m
) 1//2)
3687 (ftake '%jacobi_nc u m
)))))))
3688 ((and (alike1 lin
1//2)
3690 ;; jacobi_dn/jacobi_sn
3693 (mul 1//2 (ftake '%elliptic_kc m
))
3696 (mul 1//2 (ftake '%elliptic_kc m
))
3705 ;; jacobi_dc(u,m) = jacobi_dn/jacobi_cn
3709 ((mtimes) ((mplus) 1 ((mtimes) -
1 m
))
3710 ((mexpt) ((%jacobi_cn
) u m
) -
2)
3714 ((mtimes) ((mexpt) ((%jacobi_cn
) u m
) -
1)
3716 ((mtimes) ((rat) -
1 2)
3717 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3719 ((mexpt) ((%jacobi_sn
) u m
) 2))
3720 ((mtimes) ((rat) -
1 2) ((%jacobi_cn
) u m
)
3724 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3725 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3727 ((mtimes) -
1 ((mexpt) ((%jacobi_cn
) u m
) -
2)
3730 ((mtimes) ((rat) -
1 2)
3731 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3733 ((mexpt) ((%jacobi_sn
) u m
) 2))
3734 ((mtimes) ((rat) -
1 2) ((mexpt) m -
1)
3735 ((%jacobi_dn
) u m
) ((%jacobi_sn
) u m
)
3738 ((mexpt) ((mplus) 1 ((mtimes) -
1 m
)) -
1)
3739 ((%elliptic_e
) ((%asin
) ((%jacobi_sn
) u m
))
3743 (def-simplifier jacobi_dc
(u m
)
3746 ((float-numerical-eval-p u m
)
3747 (let ((fu (bigfloat:to
($float u
)))
3748 (fm (bigfloat:to
($float m
))))
3749 (to (bigfloat:/ (bigfloat::dn fu fm
) (bigfloat::cn fu fm
)))))
3750 ((setf args
(complex-float-numerical-eval-p u m
))
3751 (destructuring-bind (u m
)
3753 (let ((fu (bigfloat:to
($float u
)))
3754 (fm (bigfloat:to
($float m
))))
3755 (to (bigfloat:/ (bigfloat::dn fu fm
) (bigfloat::cn fu fm
))))))
3756 ((bigfloat-numerical-eval-p u m
)
3757 (let ((uu (bigfloat:to
($bfloat u
)))
3758 (mm (bigfloat:to
($bfloat m
))))
3759 (to (bigfloat:/ (bigfloat::dn uu mm
)
3760 (bigfloat::cn uu mm
)))))
3761 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
3762 (destructuring-bind (u m
)
3764 (let ((uu (bigfloat:to
($bfloat u
)))
3765 (mm (bigfloat:to
($bfloat m
))))
3766 (to (bigfloat:/ (bigfloat::dn uu mm
)
3767 (bigfloat::cn uu mm
))))))
3776 ((and $trigsign
(mminusp* u
))
3777 (ftake* '%jacobi_dc
(neg u
) m
))
3780 (member (caar u
) '(%inverse_jacobi_sn
3791 %inverse_jacobi_dc
))
3792 (alike1 (third u
) m
))
3793 (cond ((eq (caar u
) '%inverse_jacobi_dc
)
3796 ;; Express in terms of dn and cn
3797 (div (ftake '%jacobi_dn u m
)
3798 (ftake '%jacobi_cn u m
)))))
3799 ;; A&S 16.20 (Jacobi's Imaginary transformation)
3800 ((and $%iargs
(multiplep u
'$%i
))
3801 ;; dc(i*u) = dn(i*u)/cn(i*u) = dc(u,m1)/nc(u,m1) = dn(u,m1)
3802 (ftake* '%jacobi_dn
(coeff u
'$%i
1) (add 1 (neg m
))))
3803 ((setf coef
(kc-arg2 u m
))
3805 (destructuring-bind (lin const
)
3807 (cond ((integerp lin
)
3810 ;; dc(4*m*K + u) = dc(u)
3814 (ftake '%jacobi_dc const m
)))
3816 ;; dc(4*m*K + K + u) = dc(K+u) = -ns(u)
3819 (dbz-err1 'jacobi_dc
)
3820 (neg (ftake '%jacobi_ns const m
))))
3822 ;; dc(4*m*K + 2*K + u) = dc(2*K+u) = -dc(u)
3826 (neg (ftake '%jacobi_dc const m
))))
3828 ;; dc(4*m*K + 3*K + u) = dc(2*K + K + u) =
3830 ;; dc(3*K) = ns(0) = inf
3832 (dbz-err1 'jacobi_dc
)
3833 (ftake '%jacobi_dc const m
)))))
3834 ((and (alike1 lin
1//2)
3836 ;; jacobi_dn/jacobi_cn
3839 (mul 1//2 (ftake '%elliptic_kc m
))
3842 (mul 1//2 (ftake '%elliptic_kc m
))
3851 ;;; Other inverse Jacobian functions
3853 ;; inverse_jacobi_ns(x)
3855 ;; Let u = inverse_jacobi_ns(x). Then jacobi_ns(u) = x or
3856 ;; 1/jacobi_sn(u) = x or
3858 ;; jacobi_sn(u) = 1/x
3860 ;; so u = inverse_jacobi_sn(1/x)
3861 (defprop %inverse_jacobi_ns
3863 ;; Whittaker and Watson, example in 22.122
3864 ;; inverse_jacobi_ns(u,m) = integrate(1/sqrt(t^2-1)/sqrt(t^2-m), t, u, inf)
3865 ;; -> -1/sqrt(x^2-1)/sqrt(x^2-m)
3867 ((mexpt) ((mplus) -
1 ((mexpt) x
2)) ((rat) -
1 2))
3869 ((mplus) ((mtimes simp ratsimp
) -
1 m
) ((mexpt) x
2))
3872 ; ((%derivative) ((%inverse_jacobi_ns) x m) m 1)
3876 (def-simplifier inverse_jacobi_ns
(u m
)
3879 ((float-numerical-eval-p u m
)
3880 ;; Numerically evaluate asn
3882 ;; ans(x,m) = asn(1/x,m) = F(asin(1/x),m)
3883 (to (elliptic-f (cl:asin
(/ ($float u
))) ($float m
))))
3884 ((complex-float-numerical-eval-p u m
)
3885 (to (elliptic-f (cl:asin
(/ (complex ($realpart
($float u
)) ($imagpart
($float u
)))))
3886 (complex ($realpart
($float m
)) ($imagpart
($float m
))))))
3887 ((bigfloat-numerical-eval-p u m
)
3888 (to (bigfloat::bf-elliptic-f
(bigfloat:asin
(bigfloat:/ (bigfloat:to
($bfloat u
))))
3889 (bigfloat:to
($bfloat m
)))))
3890 ((setf args
(complex-bigfloat-numerical-eval-p u m
))
3891 (destructuring-bind (u m
)
3893 (to (bigfloat::bf-elliptic-f
(bigfloat:asin
(bigfloat:/ (bigfloat:to
($bfloat u
))))
3894 (bigfloat:to
($bfloat m
))))))
3896 ;; ans(x,0) = F(asin(1/x),0) = asin(1/x)
3897 (ftake '%elliptic_f
(ftake '%asin
(div 1 u
)) 0))
3899 ;; ans(x,1) = F(asin(1/x),1) = log(tan(pi/2+asin(1/x)/2))
3900 (ftake '%elliptic_f
(ftake '%asin
(div 1 u
)) 1))
3902 (ftake '%elliptic_kc m
))
3904 (neg (ftake '%elliptic_kc m
)))
3905 ((and (eq $triginverses
'$all
)
3907 (eq (caar u
) '%jacobi_ns
)
3908 (alike1 (third u
) m
))
3909 ;; inverse_jacobi_ns(ns(u)) = u
3915 ;; inverse_jacobi_nc(x)
3917 ;; Let u = inverse_jacobi_nc(x). Then jacobi_nc(u) = x or
3918 ;; 1/jacobi_cn(u) = x or
3920 ;; jacobi_cn(u) = 1/x
3922 ;; so u = inverse_jacobi_cn(1/x)
3923 (defprop %inverse_jacobi_nc
3925 ;; Whittaker and Watson, example in 22.122
3926 ;; inverse_jacobi_nc(u,m) = integrate(1/sqrt(t^2-1)/sqrt((1-m)*t^2+m), t, 1, u)
3927 ;; -> 1/sqrt(x^2-1)/sqrt((1-m)*x^2+m)
3929 ((mexpt) ((mplus) -
1 ((mexpt) x
2)) ((rat) -
1 2))
3932 ((mtimes) -
1 ((mplus) -
1 m
) ((mexpt) x
2)))
3935 ; ((%derivative) ((%inverse_jacobi_nc) x m) m 1)
3939 (def-simplifier inverse_jacobi_nc
(u m
)
3940 (cond ((or (float-numerical-eval-p u m
)
3941 (complex-float-numerical-eval-p u m
)
3942 (bigfloat-numerical-eval-p u m
)
3943 (complex-bigfloat-numerical-eval-p u m
))
3945 (ftake '%inverse_jacobi_cn
($rectform
(div 1 u
)) m
))
3949 (mul 2 (ftake '%elliptic_kc m
)))
3950 ((and (eq $triginverses
'$all
)
3952 (eq (caar u
) '%jacobi_nc
)
3953 (alike1 (third u
) m
))
3954 ;; inverse_jacobi_nc(nc(u)) = u
3960 ;; inverse_jacobi_nd(x)
3962 ;; Let u = inverse_jacobi_nd(x). Then jacobi_nd(u) = x or
3963 ;; 1/jacobi_dn(u) = x or
3965 ;; jacobi_dn(u) = 1/x
3967 ;; so u = inverse_jacobi_dn(1/x)
3968 (defprop %inverse_jacobi_nd
3970 ;; Whittaker and Watson, example in 22.122
3971 ;; inverse_jacobi_nd(u,m) = integrate(1/sqrt(t^2-1)/sqrt(1-(1-m)*t^2), t, 1, u)
3972 ;; -> 1/sqrt(u^2-1)/sqrt(1-(1-m)*t^2)
3974 ((mexpt) ((mplus) -
1 ((mexpt simp ratsimp
) x
2))
3978 ((mtimes) ((mplus) -
1 m
) ((mexpt simp ratsimp
) x
2)))
3981 ; ((%derivative) ((%inverse_jacobi_nd) x m) m 1)
3985 (def-simplifier inverse_jacobi_nd
(u m
)
3986 (cond ((or (float-numerical-eval-p u m
)
3987 (complex-float-numerical-eval-p u m
)
3988 (bigfloat-numerical-eval-p u m
)
3989 (complex-bigfloat-numerical-eval-p u m
))
3990 (ftake '%inverse_jacobi_dn
($rectform
(div 1 u
)) m
))
3993 ((onep1 ($ratsimp
(mul (power (sub 1 m
) 1//2) u
)))
3994 ;; jacobi_nd(1/sqrt(1-m),m) = K(m). This follows from
3995 ;; jacobi_dn(sqrt(1-m),m) = K(m).
3996 (ftake '%elliptic_kc m
))
3997 ((and (eq $triginverses
'$all
)
3999 (eq (caar u
) '%jacobi_nd
)
4000 (alike1 (third u
) m
))
4001 ;; inverse_jacobi_nd(nd(u)) = u
4007 ;; inverse_jacobi_sc(x)
4009 ;; Let u = inverse_jacobi_sc(x). Then jacobi_sc(u) = x or
4010 ;; x = jacobi_sn(u)/jacobi_cn(u)
4017 ;; sn^2 = x^2/(1+x^2)
4019 ;; sn(u) = x/sqrt(1+x^2)
4021 ;; u = inverse_sn(x/sqrt(1+x^2))
4023 (defprop %inverse_jacobi_sc
4025 ;; Whittaker and Watson, example in 22.122
4026 ;; inverse_jacobi_sc(u,m) = integrate(1/sqrt(1+t^2)/sqrt(1+(1-m)*t^2), t, 0, u)
4027 ;; -> 1/sqrt(1+x^2)/sqrt(1+(1-m)*x^2)
4029 ((mexpt) ((mplus) 1 ((mexpt) x
2))
4033 ((mtimes) -
1 ((mplus) -
1 m
) ((mexpt) x
2)))
4036 ; ((%derivative) ((%inverse_jacobi_sc) x m) m 1)
4040 (def-simplifier inverse_jacobi_sc
(u m
)
4041 (cond ((or (float-numerical-eval-p u m
)
4042 (complex-float-numerical-eval-p u m
)
4043 (bigfloat-numerical-eval-p u m
)
4044 (complex-bigfloat-numerical-eval-p u m
))
4045 (ftake '%inverse_jacobi_sn
4046 ($rectform
(div u
(power (add 1 (mul u u
)) 1//2)))
4049 ;; jacobi_sc(0,m) = 0
4051 ((and (eq $triginverses
'$all
)
4053 (eq (caar u
) '%jacobi_sc
)
4054 (alike1 (third u
) m
))
4055 ;; inverse_jacobi_sc(sc(u)) = u
4061 ;; inverse_jacobi_sd(x)
4063 ;; Let u = inverse_jacobi_sd(x). Then jacobi_sd(u) = x or
4064 ;; x = jacobi_sn(u)/jacobi_dn(u)
4067 ;; = sn^2/(1-m*sn^2)
4071 ;; sn^2 = x^2/(1+m*x^2)
4073 ;; sn(u) = x/sqrt(1+m*x^2)
4075 ;; u = inverse_sn(x/sqrt(1+m*x^2))
4077 (defprop %inverse_jacobi_sd
4079 ;; Whittaker and Watson, example in 22.122
4080 ;; inverse_jacobi_sd(u,m) = integrate(1/sqrt(1-(1-m)*t^2)/sqrt(1+m*t^2), t, 0, u)
4081 ;; -> 1/sqrt(1-(1-m)*x^2)/sqrt(1+m*x^2)
4084 ((mplus) 1 ((mtimes) ((mplus) -
1 m
) ((mexpt) x
2)))
4086 ((mexpt) ((mplus) 1 ((mtimes) m
((mexpt) x
2)))
4089 ; ((%derivative) ((%inverse_jacobi_sd) x m) m 1)
4093 (def-simplifier inverse_jacobi_sd
(u m
)
4094 (cond ((or (float-numerical-eval-p u m
)
4095 (complex-float-numerical-eval-p u m
)
4096 (bigfloat-numerical-eval-p u m
)
4097 (complex-bigfloat-numerical-eval-p u m
))
4098 (ftake '%inverse_jacobi_sn
4099 ($rectform
(div u
(power (add 1 (mul m
(mul u u
))) 1//2)))
4103 ((eql 0 ($ratsimp
(sub u
(div 1 (power (sub 1 m
) 1//2)))))
4104 ;; inverse_jacobi_sd(1/sqrt(1-m), m) = elliptic_kc(m)
4106 ;; We can see this from inverse_jacobi_sd(x,m) =
4107 ;; inverse_jacobi_sn(x/sqrt(1+m*x^2), m). So
4108 ;; inverse_jacobi_sd(1/sqrt(1-m),m) = inverse_jacobi_sn(1,m)
4109 (ftake '%elliptic_kc m
))
4110 ((and (eq $triginverses
'$all
)
4112 (eq (caar u
) '%jacobi_sd
)
4113 (alike1 (third u
) m
))
4114 ;; inverse_jacobi_sd(sd(u)) = u
4120 ;; inverse_jacobi_cs(x)
4122 ;; Let u = inverse_jacobi_cs(x). Then jacobi_cs(u) = x or
4123 ;; 1/x = 1/jacobi_cs(u) = jacobi_sc(u)
4125 ;; u = inverse_sc(1/x)
4127 (defprop %inverse_jacobi_cs
4129 ;; Whittaker and Watson, example in 22.122
4130 ;; inverse_jacobi_cs(u,m) = integrate(1/sqrt(t^2+1)/sqrt(t^2+(1-m)), t, u, inf)
4131 ;; -> -1/sqrt(x^2+1)/sqrt(x^2+(1-m))
4133 ((mexpt) ((mplus) 1 ((mexpt simp ratsimp
) x
2))
4136 ((mtimes simp ratsimp
) -
1 m
)
4137 ((mexpt simp ratsimp
) x
2))
4140 ; ((%derivative) ((%inverse_jacobi_cs) x m) m 1)
4144 (def-simplifier inverse_jacobi_cs
(u m
)
4145 (cond ((or (float-numerical-eval-p u m
)
4146 (complex-float-numerical-eval-p u m
)
4147 (bigfloat-numerical-eval-p u m
)
4148 (complex-bigfloat-numerical-eval-p u m
))
4149 (ftake '%inverse_jacobi_sc
($rectform
(div 1 u
)) m
))
4151 (ftake '%elliptic_kc m
))
4156 ;; inverse_jacobi_cd(x)
4158 ;; Let u = inverse_jacobi_cd(x). Then jacobi_cd(u) = x or
4159 ;; x = jacobi_cn(u)/jacobi_dn(u)
4162 ;; = (1-sn^2)/(1-m*sn^2)
4166 ;; sn^2 = (1-x^2)/(1-m*x^2)
4168 ;; sn(u) = sqrt(1-x^2)/sqrt(1-m*x^2)
4170 ;; u = inverse_sn(sqrt(1-x^2)/sqrt(1-m*x^2))
4172 (defprop %inverse_jacobi_cd
4174 ;; Whittaker and Watson, example in 22.122
4175 ;; inverse_jacobi_cd(u,m) = integrate(1/sqrt(1-t^2)/sqrt(1-m*t^2), t, u, 1)
4176 ;; -> -1/sqrt(1-x^2)/sqrt(1-m*x^2)
4179 ((mplus) 1 ((mtimes) -
1 ((mexpt) x
2)))
4182 ((mplus) 1 ((mtimes) -
1 m
((mexpt) x
2)))
4185 ; ((%derivative) ((%inverse_jacobi_cd) x m) m 1)
4189 (def-simplifier inverse_jacobi_cd
(u m
)
4190 (cond ((or (complex-float-numerical-eval-p u m
)
4191 (complex-bigfloat-numerical-eval-p u m
))
4193 (ftake '%inverse_jacobi_sn
4194 ($rectform
(div (power (mul (sub 1 u
) (add 1 u
)) 1//2)
4195 (power (sub 1 (mul m
(mul u u
))) 1//2)))
4200 (ftake '%elliptic_kc m
))
4201 ((and (eq $triginverses
'$all
)
4203 (eq (caar u
) '%jacobi_cd
)
4204 (alike1 (third u
) m
))
4205 ;; inverse_jacobi_cd(cd(u)) = u
4211 ;; inverse_jacobi_ds(x)
4213 ;; Let u = inverse_jacobi_ds(x). Then jacobi_ds(u) = x or
4214 ;; 1/x = 1/jacobi_ds(u) = jacobi_sd(u)
4216 ;; u = inverse_sd(1/x)
4218 (defprop %inverse_jacobi_ds
4220 ;; Whittaker and Watson, example in 22.122
4221 ;; inverse_jacobi_ds(u,m) = integrate(1/sqrt(t^2-(1-m))/sqrt(t^2+m), t, u, inf)
4222 ;; -> -1/sqrt(x^2-(1-m))/sqrt(x^2+m)
4225 ((mplus) -
1 m
((mexpt simp ratsimp
) x
2))
4228 ((mplus) m
((mexpt simp ratsimp
) x
2))
4231 ; ((%derivative) ((%inverse_jacobi_ds) x m) m 1)
4235 (def-simplifier inverse_jacobi_ds
(u m
)
4236 (cond ((or (float-numerical-eval-p u m
)
4237 (complex-float-numerical-eval-p u m
)
4238 (bigfloat-numerical-eval-p u m
)
4239 (complex-bigfloat-numerical-eval-p u m
))
4240 (ftake '%inverse_jacobi_sd
($rectform
(div 1 u
)) m
))
4241 ((and $trigsign
(mminusp* u
))
4242 (neg (ftake* '%inverse_jacobi_ds
(neg u
) m
)))
4243 ((eql 0 ($ratsimp
(sub u
(power (sub 1 m
) 1//2))))
4244 ;; inverse_jacobi_ds(sqrt(1-m),m) = elliptic_kc(m)
4246 ;; Since inverse_jacobi_ds(sqrt(1-m), m) =
4247 ;; inverse_jacobi_sd(1/sqrt(1-m),m). And we know from
4248 ;; above that this is elliptic_kc(m)
4249 (ftake '%elliptic_kc m
))
4250 ((and (eq $triginverses
'$all
)
4252 (eq (caar u
) '%jacobi_ds
)
4253 (alike1 (third u
) m
))
4254 ;; inverse_jacobi_ds(ds(u)) = u
4261 ;; inverse_jacobi_dc(x)
4263 ;; Let u = inverse_jacobi_dc(x). Then jacobi_dc(u) = x or
4264 ;; 1/x = 1/jacobi_dc(u) = jacobi_cd(u)
4266 ;; u = inverse_cd(1/x)
4268 (defprop %inverse_jacobi_dc
4270 ;; Note: Whittaker and Watson, example in 22.122 says
4271 ;; inverse_jacobi_dc(u,m) = integrate(1/sqrt(t^2-1)/sqrt(t^2-m),
4272 ;; t, u, 1) but that seems wrong. A&S 17.4.47 says
4273 ;; integrate(1/sqrt(t^2-1)/sqrt(t^2-m), t, a, u) =
4274 ;; inverse_jacobi_cd(x,m). Lawden 3.2.8 says the same.
4275 ;; functions.wolfram.com says the derivative is
4276 ;; 1/sqrt(t^2-1)/sqrt(t^2-m).
4279 ((mplus) -
1 ((mexpt simp ratsimp
) x
2))
4283 ((mtimes simp ratsimp
) -
1 m
)
4284 ((mexpt simp ratsimp
) x
2))
4287 ; ((%derivative) ((%inverse_jacobi_dc) x m) m 1)
4291 (def-simplifier inverse_jacobi_dc
(u m
)
4292 (cond ((or (complex-float-numerical-eval-p u m
)
4293 (complex-bigfloat-numerical-eval-p u m
))
4294 (ftake '%inverse_jacobi_cd
($rectform
(div 1 u
)) m
))
4297 ((and (eq $triginverses
'$all
)
4299 (eq (caar u
) '%jacobi_dc
)
4300 (alike1 (third u
) m
))
4301 ;; inverse_jacobi_dc(dc(u)) = u
4307 ;; Convert an inverse Jacobian function into the equivalent elliptic
4310 ;; See A&S 17.4.41-17.4.52.
4311 (defun make-elliptic-f (e)
4314 ((member (caar e
) '(%inverse_jacobi_sc %inverse_jacobi_cs
4315 %inverse_jacobi_nd %inverse_jacobi_dn
4316 %inverse_jacobi_sn %inverse_jacobi_cd
4317 %inverse_jacobi_dc %inverse_jacobi_ns
4318 %inverse_jacobi_nc %inverse_jacobi_ds
4319 %inverse_jacobi_sd %inverse_jacobi_cn
))
4320 ;; We have some inverse Jacobi function. Convert it to the F form.
4321 (destructuring-bind ((fn &rest ops
) u m
)
4323 (declare (ignore ops
))
4327 (ftake '%elliptic_f
(ftake '%atan u
) m
))
4330 (ftake '%elliptic_f
(ftake '%atan
(div 1 u
)) m
))
4335 (mul (power m -
1//2)
4337 (power (add -
1 (mul u u
))
4345 (power (sub 1 (power u
2)) 1//2)))
4349 (ftake '%elliptic_f
(ftake '%asin u
) m
))
4354 (power (mul (sub 1 (mul u u
))
4355 (sub 1 (mul m u u
)))
4362 (power (mul (sub (mul u u
) 1)
4368 (ftake '%elliptic_f
(ftake '%asin
(div 1 u
)) m
))
4371 (ftake '%elliptic_f
(ftake '%acos
(div 1 u
)) m
))
4376 (power (add m
(mul u u
))
4384 (power (add 1 (mul m u u
))
4389 (ftake '%elliptic_f
(ftake '%acos u
) m
)))))
4391 (recur-apply #'make-elliptic-f e
))))
4393 (defmfun $make_elliptic_f
(e)
4396 (simplify (make-elliptic-f e
))))
4398 (defun make-elliptic-e (e)
4400 ((eq (caar e
) '$elliptic_eu
)
4401 (destructuring-bind ((ffun &rest ops
) u m
) e
4402 (declare (ignore ffun ops
))
4403 (ftake '%elliptic_e
(ftake '%asin
(ftake '%jacobi_sn u m
)) m
)))
4405 (recur-apply #'make-elliptic-e e
))))
4407 (defmfun $make_elliptic_e
(e)
4410 (simplify (make-elliptic-e e
))))
4413 ;; Eu(u,m) = integrate(jacobi_dn(v,m)^2,v,0,u)
4414 ;; = integrate(sqrt((1-m*t^2)/(1-t^2)),t,0,jacobi_sn(u,m))
4416 ;; Eu(u,m) = E(am(u),m)
4418 ;; where E(u,m) is elliptic-e above.
4421 ;; Lawden gives the following relationships
4423 ;; E(u+v) = E(u) + E(v) - m*sn(u)*sn(v)*sn(u+v)
4424 ;; E(u,0) = u, E(u,1) = tanh u
4426 ;; E(i*u,k) = i*sc(u,k')*dn(u,k') - i*E(u,k') + i*u
4428 ;; E(2*i*K') = 2*i*(K'-E')
4430 ;; E(u + 2*i*K') = E(u) + 2*i*(K' - E')
4432 ;; E(u+K) = E(u) + E - k^2*sn(u)*cd(u)
4433 (defun elliptic-eu (u m
)
4435 ;; E(u + 2*n*K) = E(u) + 2*n*E
4436 (let ((ell-k (to (elliptic-k m
)))
4437 (ell-e (elliptic-ec m
)))
4438 (multiple-value-bind (n u-rem
)
4439 (floor u
(* 2 ell-k
))
4442 (cond ((>= u-rem ell-k
)
4443 ;; 0 <= u-rem < K so
4444 ;; E(u + K) = E(u) + E - m*sn(u)*cd(u)
4445 (let ((u-k (- u ell-k
)))
4446 (- (+ (elliptic-e (cl:asin
(bigfloat::sn u-k m
)) m
)
4448 (/ (* m
(bigfloat::sn u-k m
) (bigfloat::cn u-k m
))
4449 (bigfloat::dn u-k m
)))))
4451 (elliptic-e (cl:asin
(bigfloat::sn u m
)) m
)))))))
4455 ;; E(u+i*v, m) = E(u,m) -i*E(v,m') + i*v + i*sc(v,m')*dn(v,m')
4456 ;; -i*m*sn(u,m)*sc(v,m')*sn(u+i*v,m)
4458 (let ((u-r (realpart u
))
4461 (+ (elliptic-eu u-r m
)
4464 (/ (* (bigfloat::sn u-i m1
) (bigfloat::dn u-i m1
))
4465 (bigfloat::cn u-i m1
)))
4466 (+ (elliptic-eu u-i m1
)
4467 (/ (* m
(bigfloat::sn u-r m
) (bigfloat::sn u-i m1
) (bigfloat::sn u m
))
4468 (bigfloat::cn u-i m1
))))))))))
4470 (defprop $elliptic_eu
4472 ((mexpt) ((%jacobi_dn
) u m
) 2)
4477 (def-simplifier elliptic_eu
(u m
)
4479 ;; as it stands, ELLIPTIC-EU can't handle bigfloats or complex bigfloats,
4480 ;; so handle only floats and complex floats here.
4481 ((float-numerical-eval-p u m
)
4482 (elliptic-eu ($float u
) ($float m
)))
4483 ((complex-float-numerical-eval-p u m
)
4484 (let ((u-r ($realpart u
))
4487 (complexify (elliptic-eu (complex u-r u-i
) m
))))
4491 (def-simplifier jacobi_am
(u m
)
4493 ;; as it stands, BIGFLOAT::SN can't handle bigfloats or complex bigfloats,
4494 ;; so handle only floats and complex floats here.
4495 ((float-numerical-eval-p u m
)
4496 (cl:asin
(bigfloat::sn
($float u
) ($float m
))))
4497 ((complex-float-numerical-eval-p u m
)
4498 (let ((u-r ($realpart
($float u
)))
4499 (u-i ($imagpart
($float u
)))
4501 (complexify (cl:asin
(bigfloat::sn
(complex u-r u-i
) m
)))))
4506 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4507 ;; Integrals. At present with respect to first argument only.
4508 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4510 ;; A&S 16.24.1: integrate(jacobi_sn(u,m),u)
4511 ;; = log(jacobi_dn(u,m)-sqrt(m)*jacobi_cn(u,m))/sqrt(m)
4514 ((mtimes simp
) ((mexpt simp
) m
((rat simp
) -
1 2))
4517 ((mtimes simp
) -
1 ((mexpt simp
) m
((rat simp
) 1 2))
4518 ((%jacobi_cn simp
) u m
))
4519 ((%jacobi_dn simp
) u m
))))
4523 ;; A&S 16.24.2: integrate(jacobi_cn(u,m),u) = acos(jacobi_dn(u,m))/sqrt(m)
4526 ((mtimes simp
) ((mexpt simp
) m
((rat simp
) -
1 2))
4527 ((%acos simp
) ((%jacobi_dn simp
) u m
)))
4531 ;; A&S 16.24.3: integrate(jacobi_dn(u,m),u) = asin(jacobi_sn(u,m))
4534 ((%asin simp
) ((%jacobi_sn simp
) u m
))
4538 ;; A&S 16.24.4: integrate(jacobi_cd(u,m),u)
4539 ;; = log(jacobi_nd(u,m)+sqrt(m)*jacobi_sd(u,m))/sqrt(m)
4542 ((mtimes simp
) ((mexpt simp
) m
((rat simp
) -
1 2))
4544 ((mplus simp
) ((%jacobi_nd simp
) u m
)
4545 ((mtimes simp
) ((mexpt simp
) m
((rat simp
) 1 2))
4546 ((%jacobi_sd simp
) u m
)))))
4550 ;; integrate(jacobi_sd(u,m),u)
4552 ;; A&S 16.24.5 gives
4553 ;; asin(-sqrt(m)*jacobi_cd(u,m))/sqrt(m*m_1), where m + m_1 = 1
4554 ;; but this does not pass some simple tests.
4556 ;; functions.wolfram.com 09.35.21.001.01 gives
4557 ;; -asin(sqrt(m)*jacobi_cd(u,m))*sqrt(1-m*jacobi_cd(u,m)^2)*jacobi_dn(u,m)/((1-m)*sqrt(m))
4558 ;; and this does pass.
4562 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
4563 ((mexpt simp
) m
((rat simp
) -
1 2))
4566 ((mtimes simp
) -
1 $m
((mexpt simp
) ((%jacobi_cd simp
) u m
) 2)))
4568 ((%jacobi_dn simp
) u m
)
4570 ((mtimes simp
) ((mexpt simp
) m
((rat simp
) 1 2))
4571 ((%jacobi_cd simp
) u m
))))
4575 ;; integrate(jacobi_nd(u,m),u)
4577 ;; A&S 16.24.6 gives
4578 ;; acos(jacobi_cd(u,m))/sqrt(m_1), where m + m_1 = 1
4579 ;; but this does not pass some simple tests.
4581 ;; functions.wolfram.com 09.32.21.0001.01 gives
4582 ;; sqrt(1-jacobi_cd(u,m)^2)*acos(jacobi_cd(u,m))/((1-m)*jacobi_sd(u,m))
4583 ;; and this does pass.
4586 ((mtimes simp
) ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
)) -
1)
4589 ((mtimes simp
) -
1 ((mexpt simp
) ((%jacobi_cd simp
) u m
) 2)))
4591 ((mexpt simp
) ((%jacobi_sd simp
) u m
) -
1)
4592 ((%acos simp
) ((%jacobi_cd simp
) u m
)))
4596 ;; A&S 16.24.7: integrate(jacobi_dc(u,m),u) = log(jacobi_nc(u,m)+jacobi_sc(u,m))
4599 ((%log simp
) ((mplus simp
) ((%jacobi_nc simp
) u m
) ((%jacobi_sc simp
) u m
)))
4603 ;; A&S 16.24.8: integrate(jacobi_nc(u,m),u)
4604 ;; = log(jacobi_dc(u,m)+sqrt(m_1)*jacobi_sc(u,m))/sqrt(m_1), where m + m_1 = 1
4608 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
))
4611 ((mplus simp
) ((%jacobi_dc simp
) u m
)
4613 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
))
4615 ((%jacobi_sc simp
) u m
)))))
4619 ;; A&S 16.24.9: integrate(jacobi_sc(u,m),u)
4620 ;; = log(jacobi_dc(u,m)+sqrt(m_1)*jacobi_nc(u,m))/sqrt(m_1), where m + m_1 = 1
4624 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
))
4627 ((mplus simp
) ((%jacobi_dc simp
) u m
)
4629 ((mexpt simp
) ((mplus simp
) 1 ((mtimes simp
) -
1 m
))
4631 ((%jacobi_nc simp
) u m
)))))
4635 ;; A&S 16.24.10: integrate(jacobi_ns(u,m),u)
4636 ;; = log(jacobi_ds(u,m)-jacobi_cs(u,m))
4640 ((mplus simp
) ((mtimes simp
) -
1 ((%jacobi_cs simp
) u m
))
4641 ((%jacobi_ds simp
) u m
)))
4645 ;; integrate(jacobi_ds(u,m),u)
4647 ;; A&S 16.24.11 gives
4648 ;; log(jacobi_ds(u,m)-jacobi_cs(u,m))
4649 ;; but this does not pass some simple tests.
4651 ;; functions.wolfram.com 09.30.21.0001.01 gives
4652 ;; log((1-jacobi_cn(u,m))/jacobi_sn(u,m))
4658 ((mplus simp
) 1 ((mtimes simp
) -
1 ((%jacobi_cn simp
) u m
)))
4659 ((mexpt simp
) ((%jacobi_sn simp
) u m
) -
1)))
4663 ;; A&S 16.24.12: integrate(jacobi_cs(u,m),u) = log(jacobi_ns(u,m)-jacobi_ds(u,m))
4667 ((mplus simp
) ((mtimes simp
) -
1 ((%jacobi_ds simp
) u m
))
4668 ((%jacobi_ns simp
) u m
)))
4672 ;; functions.wolfram.com 09.48.21.0001.01
4673 ;; integrate(inverse_jacobi_sn(u,m),u) =
4674 ;; inverse_jacobi_sn(u,m)*u
4675 ;; - log( jacobi_dn(inverse_jacobi_sn(u,m),m)
4676 ;; -sqrt(m)*jacobi_cn(inverse_jacobi_sn(u,m),m)) / sqrt(m)
4677 (defprop %inverse_jacobi_sn
4679 ((mplus simp
) ((mtimes simp
) u
((%inverse_jacobi_sn simp
) u m
))
4680 ((mtimes simp
) -
1 ((mexpt simp
) m
((rat simp
) -
1 2))
4683 ((mtimes simp
) -
1 ((mexpt simp
) m
((rat simp
) 1 2))
4684 ((%jacobi_cn simp
) ((%inverse_jacobi_sn simp
) u m
) m
))
4685 ((%jacobi_dn simp
) ((%inverse_jacobi_sn simp
) u m
) m
)))))
4689 ;; functions.wolfram.com 09.38.21.0001.01
4690 ;; integrate(inverse_jacobi_cn(u,m),u) =
4691 ;; u*inverse_jacobi_cn(u,m)
4692 ;; -%i*log(%i*jacobi_dn(inverse_jacobi_cn(u,m),m)/sqrt(m)
4693 ;; -jacobi_sn(inverse_jacobi_cn(u,m),m))
4695 (defprop %inverse_jacobi_cn
4697 ((mplus simp
) ((mtimes simp
) u
((%inverse_jacobi_cn simp
) u m
))
4698 ((mtimes simp
) -
1 $%i
((mexpt simp
) m
((rat simp
) -
1 2))
4701 ((mtimes simp
) $%i
((mexpt simp
) m
((rat simp
) -
1 2))
4702 ((%jacobi_dn simp
) ((%inverse_jacobi_cn simp
) u m
) m
))
4704 ((%jacobi_sn simp
) ((%inverse_jacobi_cn simp
) u m
) m
))))))
4708 ;; functions.wolfram.com 09.41.21.0001.01
4709 ;; integrate(inverse_jacobi_dn(u,m),u) =
4710 ;; u*inverse_jacobi_dn(u,m)
4711 ;; - %i*log(%i*jacobi_cn(inverse_jacobi_dn(u,m),m)
4712 ;; +jacobi_sn(inverse_jacobi_dn(u,m),m))
4713 (defprop %inverse_jacobi_dn
4715 ((mplus simp
) ((mtimes simp
) u
((%inverse_jacobi_dn simp
) u m
))
4716 ((mtimes simp
) -
1 $%i
4720 ((%jacobi_cn simp
) ((%inverse_jacobi_dn simp
) u m
) m
))
4721 ((%jacobi_sn simp
) ((%inverse_jacobi_dn simp
) u m
) m
)))))
4726 ;; Real and imaginary part for Jacobi elliptic functions.
4727 (defprop %jacobi_sn risplit-sn-cn-dn risplit-function
)
4728 (defprop %jacobi_cn risplit-sn-cn-dn risplit-function
)
4729 (defprop %jacobi_dn risplit-sn-cn-dn risplit-function
)
4731 (defun risplit-sn-cn-dn (expr)
4732 (let* ((arg (second expr
))
4733 (param (third expr
)))
4734 ;; We only split on the argument, not the order
4735 (destructuring-bind (arg-r . arg-i
)
4739 (cons (take (first expr
) arg-r param
)
4742 (let* ((s (ftake '%jacobi_sn arg-r param
))
4743 (c (ftake '%jacobi_cn arg-r param
))
4744 (d (ftake '%jacobi_dn arg-r param
))
4745 (s1 (ftake '%jacobi_sn arg-i
(sub 1 param
)))
4746 (c1 (ftake '%jacobi_cn arg-i
(sub 1 param
)))
4747 (d1 (ftake '%jacobi_dn arg-i
(sub 1 param
)))
4748 (den (add (mul c1 c1
)
4752 ;; Let s = jacobi_sn(x,m)
4753 ;; c = jacobi_cn(x,m)
4754 ;; d = jacobi_dn(x,m)
4755 ;; s1 = jacobi_sn(y,1-m)
4756 ;; c1 = jacobi_cn(y,1-m)
4757 ;; d1 = jacobi_dn(y,1-m)
4761 ;; jacobi_sn(x+%i*y,m) =
4763 ;; s*d1 + %i*c*d*s1*c1
4764 ;; -------------------
4767 (cons (div (mul s d1
) den
)
4768 (div (mul c
(mul d
(mul s1 c1
)))
4775 ;; c*c1 - %i*s*d*s1*d1
4776 ;; -------------------
4778 (cons (div (mul c c1
) den
)
4780 (mul s
(mul d
(mul s1 d1
))))
4787 ;; d*c1*d1 - %i*m*s*c*s1
4788 ;; ---------------------
4790 (cons (div (mul d
(mul c1 d1
))
4792 (div (mul -
1 (mul param
(mul s
(mul c s1
))))