Coverage for /usr/lib/python3/dist-packages/mpmath/functions/elliptic.py: 8%
455 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1r"""
2Elliptic functions historically comprise the elliptic integrals
3and their inverses, and originate from the problem of computing the
4arc length of an ellipse. From a more modern point of view,
5an elliptic function is defined as a doubly periodic function, i.e.
6a function which satisfies
8.. math ::
10 f(z + 2 \omega_1) = f(z + 2 \omega_2) = f(z)
12for some half-periods `\omega_1, \omega_2` with
13`\mathrm{Im}[\omega_1 / \omega_2] > 0`. The canonical elliptic
14functions are the Jacobi elliptic functions. More broadly, this section
15includes quasi-doubly periodic functions (such as the Jacobi theta
16functions) and other functions useful in the study of elliptic functions.
18Many different conventions for the arguments of
19elliptic functions are in use. It is even standard to use
20different parameterizations for different functions in the same
21text or software (and mpmath is no exception).
22The usual parameters are the elliptic nome `q`, which usually
23must satisfy `|q| < 1`; the elliptic parameter `m` (an arbitrary
24complex number); the elliptic modulus `k` (an arbitrary complex
25number); and the half-period ratio `\tau`, which usually must
26satisfy `\mathrm{Im}[\tau] > 0`.
27These quantities can be expressed in terms of each other
28using the following relations:
30.. math ::
32 m = k^2
34.. math ::
36 \tau = i \frac{K(1-m)}{K(m)}
38.. math ::
40 q = e^{i \pi \tau}
42.. math ::
44 k = \frac{\vartheta_2^2(q)}{\vartheta_3^2(q)}
46In addition, an alternative definition is used for the nome in
47number theory, which we here denote by q-bar:
49.. math ::
51 \bar{q} = q^2 = e^{2 i \pi \tau}
53For convenience, mpmath provides functions to convert
54between the various parameters (:func:`~mpmath.qfrom`, :func:`~mpmath.mfrom`,
55:func:`~mpmath.kfrom`, :func:`~mpmath.taufrom`, :func:`~mpmath.qbarfrom`).
57**References**
591. [AbramowitzStegun]_
612. [WhittakerWatson]_
63"""
65from .functions import defun, defun_wrapped
67@defun_wrapped
68def eta(ctx, tau):
69 r"""
70 Returns the Dedekind eta function of tau in the upper half-plane.
72 >>> from mpmath import *
73 >>> mp.dps = 25; mp.pretty = True
74 >>> eta(1j); gamma(0.25) / (2*pi**0.75)
75 (0.7682254223260566590025942 + 0.0j)
76 0.7682254223260566590025942
77 >>> tau = sqrt(2) + sqrt(5)*1j
78 >>> eta(-1/tau); sqrt(-1j*tau) * eta(tau)
79 (0.9022859908439376463573294 + 0.07985093673948098408048575j)
80 (0.9022859908439376463573295 + 0.07985093673948098408048575j)
81 >>> eta(tau+1); exp(pi*1j/12) * eta(tau)
82 (0.4493066139717553786223114 + 0.3290014793877986663915939j)
83 (0.4493066139717553786223114 + 0.3290014793877986663915939j)
84 >>> f = lambda z: diff(eta, z) / eta(z)
85 >>> chop(36*diff(f,tau)**2 - 24*diff(f,tau,2)*f(tau) + diff(f,tau,3))
86 0.0
88 """
89 if ctx.im(tau) <= 0.0:
90 raise ValueError("eta is only defined in the upper half-plane")
91 q = ctx.expjpi(tau/12)
92 return q * ctx.qp(q**24)
94def nome(ctx, m):
95 m = ctx.convert(m)
96 if not m:
97 return m
98 if m == ctx.one:
99 return m
100 if ctx.isnan(m):
101 return m
102 if ctx.isinf(m):
103 if m == ctx.ninf:
104 return type(m)(-1)
105 else:
106 return ctx.mpc(-1)
107 a = ctx.ellipk(ctx.one-m)
108 b = ctx.ellipk(m)
109 v = ctx.exp(-ctx.pi*a/b)
110 if not ctx._im(m) and ctx._re(m) < 1:
111 if ctx._is_real_type(m):
112 return v.real
113 else:
114 return v.real + 0j
115 elif m == 2:
116 v = ctx.mpc(0, v.imag)
117 return v
119@defun_wrapped
120def qfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
121 r"""
122 Returns the elliptic nome `q`, given any of `q, m, k, \tau, \bar{q}`::
124 >>> from mpmath import *
125 >>> mp.dps = 25; mp.pretty = True
126 >>> qfrom(q=0.25)
127 0.25
128 >>> qfrom(m=mfrom(q=0.25))
129 0.25
130 >>> qfrom(k=kfrom(q=0.25))
131 0.25
132 >>> qfrom(tau=taufrom(q=0.25))
133 (0.25 + 0.0j)
134 >>> qfrom(qbar=qbarfrom(q=0.25))
135 0.25
137 """
138 if q is not None:
139 return ctx.convert(q)
140 if m is not None:
141 return nome(ctx, m)
142 if k is not None:
143 return nome(ctx, ctx.convert(k)**2)
144 if tau is not None:
145 return ctx.expjpi(tau)
146 if qbar is not None:
147 return ctx.sqrt(qbar)
149@defun_wrapped
150def qbarfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
151 r"""
152 Returns the number-theoretic nome `\bar q`, given any of
153 `q, m, k, \tau, \bar{q}`::
155 >>> from mpmath import *
156 >>> mp.dps = 25; mp.pretty = True
157 >>> qbarfrom(qbar=0.25)
158 0.25
159 >>> qbarfrom(q=qfrom(qbar=0.25))
160 0.25
161 >>> qbarfrom(m=extraprec(20)(mfrom)(qbar=0.25)) # ill-conditioned
162 0.25
163 >>> qbarfrom(k=extraprec(20)(kfrom)(qbar=0.25)) # ill-conditioned
164 0.25
165 >>> qbarfrom(tau=taufrom(qbar=0.25))
166 (0.25 + 0.0j)
168 """
169 if qbar is not None:
170 return ctx.convert(qbar)
171 if q is not None:
172 return ctx.convert(q) ** 2
173 if m is not None:
174 return nome(ctx, m) ** 2
175 if k is not None:
176 return nome(ctx, ctx.convert(k)**2) ** 2
177 if tau is not None:
178 return ctx.expjpi(2*tau)
180@defun_wrapped
181def taufrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
182 r"""
183 Returns the elliptic half-period ratio `\tau`, given any of
184 `q, m, k, \tau, \bar{q}`::
186 >>> from mpmath import *
187 >>> mp.dps = 25; mp.pretty = True
188 >>> taufrom(tau=0.5j)
189 (0.0 + 0.5j)
190 >>> taufrom(q=qfrom(tau=0.5j))
191 (0.0 + 0.5j)
192 >>> taufrom(m=mfrom(tau=0.5j))
193 (0.0 + 0.5j)
194 >>> taufrom(k=kfrom(tau=0.5j))
195 (0.0 + 0.5j)
196 >>> taufrom(qbar=qbarfrom(tau=0.5j))
197 (0.0 + 0.5j)
199 """
200 if tau is not None:
201 return ctx.convert(tau)
202 if m is not None:
203 m = ctx.convert(m)
204 return ctx.j*ctx.ellipk(1-m)/ctx.ellipk(m)
205 if k is not None:
206 k = ctx.convert(k)
207 return ctx.j*ctx.ellipk(1-k**2)/ctx.ellipk(k**2)
208 if q is not None:
209 return ctx.log(q) / (ctx.pi*ctx.j)
210 if qbar is not None:
211 qbar = ctx.convert(qbar)
212 return ctx.log(qbar) / (2*ctx.pi*ctx.j)
214@defun_wrapped
215def kfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
216 r"""
217 Returns the elliptic modulus `k`, given any of
218 `q, m, k, \tau, \bar{q}`::
220 >>> from mpmath import *
221 >>> mp.dps = 25; mp.pretty = True
222 >>> kfrom(k=0.25)
223 0.25
224 >>> kfrom(m=mfrom(k=0.25))
225 0.25
226 >>> kfrom(q=qfrom(k=0.25))
227 0.25
228 >>> kfrom(tau=taufrom(k=0.25))
229 (0.25 + 0.0j)
230 >>> kfrom(qbar=qbarfrom(k=0.25))
231 0.25
233 As `q \to 1` and `q \to -1`, `k` rapidly approaches
234 `1` and `i \infty` respectively::
236 >>> kfrom(q=0.75)
237 0.9999999999999899166471767
238 >>> kfrom(q=-0.75)
239 (0.0 + 7041781.096692038332790615j)
240 >>> kfrom(q=1)
241 1
242 >>> kfrom(q=-1)
243 (0.0 + +infj)
244 """
245 if k is not None:
246 return ctx.convert(k)
247 if m is not None:
248 return ctx.sqrt(m)
249 if tau is not None:
250 q = ctx.expjpi(tau)
251 if qbar is not None:
252 q = ctx.sqrt(qbar)
253 if q == 1:
254 return q
255 if q == -1:
256 return ctx.mpc(0,'inf')
257 return (ctx.jtheta(2,0,q)/ctx.jtheta(3,0,q))**2
259@defun_wrapped
260def mfrom(ctx, q=None, m=None, k=None, tau=None, qbar=None):
261 r"""
262 Returns the elliptic parameter `m`, given any of
263 `q, m, k, \tau, \bar{q}`::
265 >>> from mpmath import *
266 >>> mp.dps = 25; mp.pretty = True
267 >>> mfrom(m=0.25)
268 0.25
269 >>> mfrom(q=qfrom(m=0.25))
270 0.25
271 >>> mfrom(k=kfrom(m=0.25))
272 0.25
273 >>> mfrom(tau=taufrom(m=0.25))
274 (0.25 + 0.0j)
275 >>> mfrom(qbar=qbarfrom(m=0.25))
276 0.25
278 As `q \to 1` and `q \to -1`, `m` rapidly approaches
279 `1` and `-\infty` respectively::
281 >>> mfrom(q=0.75)
282 0.9999999999999798332943533
283 >>> mfrom(q=-0.75)
284 -49586681013729.32611558353
285 >>> mfrom(q=1)
286 1.0
287 >>> mfrom(q=-1)
288 -inf
290 The inverse nome as a function of `q` has an integer
291 Taylor series expansion::
293 >>> taylor(lambda q: mfrom(q), 0, 7)
294 [0.0, 16.0, -128.0, 704.0, -3072.0, 11488.0, -38400.0, 117632.0]
296 """
297 if m is not None:
298 return m
299 if k is not None:
300 return k**2
301 if tau is not None:
302 q = ctx.expjpi(tau)
303 if qbar is not None:
304 q = ctx.sqrt(qbar)
305 if q == 1:
306 return ctx.convert(q)
307 if q == -1:
308 return q*ctx.inf
309 v = (ctx.jtheta(2,0,q)/ctx.jtheta(3,0,q))**4
310 if ctx._is_real_type(q) and q < 0:
311 v = v.real
312 return v
314jacobi_spec = {
315 'sn' : ([3],[2],[1],[4], 'sin', 'tanh'),
316 'cn' : ([4],[2],[2],[4], 'cos', 'sech'),
317 'dn' : ([4],[3],[3],[4], '1', 'sech'),
318 'ns' : ([2],[3],[4],[1], 'csc', 'coth'),
319 'nc' : ([2],[4],[4],[2], 'sec', 'cosh'),
320 'nd' : ([3],[4],[4],[3], '1', 'cosh'),
321 'sc' : ([3],[4],[1],[2], 'tan', 'sinh'),
322 'sd' : ([3,3],[2,4],[1],[3], 'sin', 'sinh'),
323 'cd' : ([3],[2],[2],[3], 'cos', '1'),
324 'cs' : ([4],[3],[2],[1], 'cot', 'csch'),
325 'dc' : ([2],[3],[3],[2], 'sec', '1'),
326 'ds' : ([2,4],[3,3],[3],[1], 'csc', 'csch'),
327 'cc' : None,
328 'ss' : None,
329 'nn' : None,
330 'dd' : None
331}
333@defun
334def ellipfun(ctx, kind, u=None, m=None, q=None, k=None, tau=None):
335 try:
336 S = jacobi_spec[kind]
337 except KeyError:
338 raise ValueError("First argument must be a two-character string "
339 "containing 's', 'c', 'd' or 'n', e.g.: 'sn'")
340 if u is None:
341 def f(*args, **kwargs):
342 return ctx.ellipfun(kind, *args, **kwargs)
343 f.__name__ = kind
344 return f
345 prec = ctx.prec
346 try:
347 ctx.prec += 10
348 u = ctx.convert(u)
349 q = ctx.qfrom(m=m, q=q, k=k, tau=tau)
350 if S is None:
351 v = ctx.one + 0*q*u
352 elif q == ctx.zero:
353 if S[4] == '1': v = ctx.one
354 else: v = getattr(ctx, S[4])(u)
355 v += 0*q*u
356 elif q == ctx.one:
357 if S[5] == '1': v = ctx.one
358 else: v = getattr(ctx, S[5])(u)
359 v += 0*q*u
360 else:
361 t = u / ctx.jtheta(3, 0, q)**2
362 v = ctx.one
363 for a in S[0]: v *= ctx.jtheta(a, 0, q)
364 for b in S[1]: v /= ctx.jtheta(b, 0, q)
365 for c in S[2]: v *= ctx.jtheta(c, t, q)
366 for d in S[3]: v /= ctx.jtheta(d, t, q)
367 finally:
368 ctx.prec = prec
369 return +v
371@defun_wrapped
372def kleinj(ctx, tau=None, **kwargs):
373 r"""
374 Evaluates the Klein j-invariant, which is a modular function defined for
375 `\tau` in the upper half-plane as
377 .. math ::
379 J(\tau) = \frac{g_2^3(\tau)}{g_2^3(\tau) - 27 g_3^2(\tau)}
381 where `g_2` and `g_3` are the modular invariants of the Weierstrass
382 elliptic function,
384 .. math ::
386 g_2(\tau) = 60 \sum_{(m,n) \in \mathbb{Z}^2 \setminus (0,0)} (m \tau+n)^{-4}
388 g_3(\tau) = 140 \sum_{(m,n) \in \mathbb{Z}^2 \setminus (0,0)} (m \tau+n)^{-6}.
390 An alternative, common notation is that of the j-function
391 `j(\tau) = 1728 J(\tau)`.
393 **Plots**
395 .. literalinclude :: /plots/kleinj.py
396 .. image :: /plots/kleinj.png
397 .. literalinclude :: /plots/kleinj2.py
398 .. image :: /plots/kleinj2.png
400 **Examples**
402 Verifying the functional equation `J(\tau) = J(\tau+1) = J(-\tau^{-1})`::
404 >>> from mpmath import *
405 >>> mp.dps = 25; mp.pretty = True
406 >>> tau = 0.625+0.75*j
407 >>> tau = 0.625+0.75*j
408 >>> kleinj(tau)
409 (-0.1507492166511182267125242 + 0.07595948379084571927228948j)
410 >>> kleinj(tau+1)
411 (-0.1507492166511182267125242 + 0.07595948379084571927228948j)
412 >>> kleinj(-1/tau)
413 (-0.1507492166511182267125242 + 0.07595948379084571927228946j)
415 The j-function has a famous Laurent series expansion in terms of the nome
416 `\bar{q}`, `j(\tau) = \bar{q}^{-1} + 744 + 196884\bar{q} + \ldots`::
418 >>> mp.dps = 15
419 >>> taylor(lambda q: 1728*q*kleinj(qbar=q), 0, 5, singular=True)
420 [1.0, 744.0, 196884.0, 21493760.0, 864299970.0, 20245856256.0]
422 The j-function admits exact evaluation at special algebraic points
423 related to the Heegner numbers 1, 2, 3, 7, 11, 19, 43, 67, 163::
425 >>> @extraprec(10)
426 ... def h(n):
427 ... v = (1+sqrt(n)*j)
428 ... if n > 2:
429 ... v *= 0.5
430 ... return v
431 ...
432 >>> mp.dps = 25
433 >>> for n in [1,2,3,7,11,19,43,67,163]:
434 ... n, chop(1728*kleinj(h(n)))
435 ...
436 (1, 1728.0)
437 (2, 8000.0)
438 (3, 0.0)
439 (7, -3375.0)
440 (11, -32768.0)
441 (19, -884736.0)
442 (43, -884736000.0)
443 (67, -147197952000.0)
444 (163, -262537412640768000.0)
446 Also at other special points, the j-function assumes explicit
447 algebraic values, e.g.::
449 >>> chop(1728*kleinj(j*sqrt(5)))
450 1264538.909475140509320227
451 >>> identify(cbrt(_)) # note: not simplified
452 '((100+sqrt(13520))/2)'
453 >>> (50+26*sqrt(5))**3
454 1264538.909475140509320227
456 """
457 q = ctx.qfrom(tau=tau, **kwargs)
458 t2 = ctx.jtheta(2,0,q)
459 t3 = ctx.jtheta(3,0,q)
460 t4 = ctx.jtheta(4,0,q)
461 P = (t2**8 + t3**8 + t4**8)**3
462 Q = 54*(t2*t3*t4)**8
463 return P/Q
466def RF_calc(ctx, x, y, z, r):
467 if y == z: return RC_calc(ctx, x, y, r)
468 if x == z: return RC_calc(ctx, y, x, r)
469 if x == y: return RC_calc(ctx, z, x, r)
470 if not (ctx.isnormal(x) and ctx.isnormal(y) and ctx.isnormal(z)):
471 if ctx.isnan(x) or ctx.isnan(y) or ctx.isnan(z):
472 return x*y*z
473 if ctx.isinf(x) or ctx.isinf(y) or ctx.isinf(z):
474 return ctx.zero
475 xm,ym,zm = x,y,z
476 A0 = Am = (x+y+z)/3
477 Q = ctx.root(3*r, -6) * max(abs(A0-x),abs(A0-y),abs(A0-z))
478 g = ctx.mpf(0.25)
479 pow4 = ctx.one
480 m = 0
481 while 1:
482 xs = ctx.sqrt(xm)
483 ys = ctx.sqrt(ym)
484 zs = ctx.sqrt(zm)
485 lm = xs*ys + xs*zs + ys*zs
486 Am1 = (Am+lm)*g
487 xm, ym, zm = (xm+lm)*g, (ym+lm)*g, (zm+lm)*g
488 if pow4 * Q < abs(Am):
489 break
490 Am = Am1
491 m += 1
492 pow4 *= g
493 t = pow4/Am
494 X = (A0-x)*t
495 Y = (A0-y)*t
496 Z = -X-Y
497 E2 = X*Y-Z**2
498 E3 = X*Y*Z
499 return ctx.power(Am,-0.5) * (9240-924*E2+385*E2**2+660*E3-630*E2*E3)/9240
501def RC_calc(ctx, x, y, r, pv=True):
502 if not (ctx.isnormal(x) and ctx.isnormal(y)):
503 if ctx.isinf(x) or ctx.isinf(y):
504 return 1/(x*y)
505 if y == 0:
506 return ctx.inf
507 if x == 0:
508 return ctx.pi / ctx.sqrt(y) / 2
509 raise ValueError
510 # Cauchy principal value
511 if pv and ctx._im(y) == 0 and ctx._re(y) < 0:
512 return ctx.sqrt(x/(x-y)) * RC_calc(ctx, x-y, -y, r)
513 if x == y:
514 return 1/ctx.sqrt(x)
515 extraprec = 2*max(0,-ctx.mag(x-y)+ctx.mag(x))
516 ctx.prec += extraprec
517 if ctx._is_real_type(x) and ctx._is_real_type(y):
518 x = ctx._re(x)
519 y = ctx._re(y)
520 a = ctx.sqrt(x/y)
521 if x < y:
522 b = ctx.sqrt(y-x)
523 v = ctx.acos(a)/b
524 else:
525 b = ctx.sqrt(x-y)
526 v = ctx.acosh(a)/b
527 else:
528 sx = ctx.sqrt(x)
529 sy = ctx.sqrt(y)
530 v = ctx.acos(sx/sy)/(ctx.sqrt((1-x/y))*sy)
531 ctx.prec -= extraprec
532 return v
534def RJ_calc(ctx, x, y, z, p, r, integration):
535 """
536 With integration == 0, computes RJ only using Carlson's algorithm
537 (may be wrong for some values).
538 With integration == 1, uses an initial integration to make sure
539 Carlson's algorithm is correct.
540 With integration == 2, uses only integration.
541 """
542 if not (ctx.isnormal(x) and ctx.isnormal(y) and \
543 ctx.isnormal(z) and ctx.isnormal(p)):
544 if ctx.isnan(x) or ctx.isnan(y) or ctx.isnan(z) or ctx.isnan(p):
545 return x*y*z
546 if ctx.isinf(x) or ctx.isinf(y) or ctx.isinf(z) or ctx.isinf(p):
547 return ctx.zero
548 if not p:
549 return ctx.inf
550 if (not x) + (not y) + (not z) > 1:
551 return ctx.inf
552 # Check conditions and fall back on integration for argument
553 # reduction if needed. The following conditions might be needlessly
554 # restrictive.
555 initial_integral = ctx.zero
556 if integration >= 1:
557 ok = (x.real >= 0 and y.real >= 0 and z.real >= 0 and p.real > 0)
558 if not ok:
559 if x == p or y == p or z == p:
560 ok = True
561 if not ok:
562 if p.imag != 0 or p.real >= 0:
563 if (x.imag == 0 and x.real >= 0 and ctx.conj(y) == z):
564 ok = True
565 if (y.imag == 0 and y.real >= 0 and ctx.conj(x) == z):
566 ok = True
567 if (z.imag == 0 and z.real >= 0 and ctx.conj(x) == y):
568 ok = True
569 if not ok or (integration == 2):
570 N = ctx.ceil(-min(x.real, y.real, z.real, p.real)) + 1
571 # Integrate around any singularities
572 if all((t.imag >= 0 or t.real > 0) for t in [x, y, z, p]):
573 margin = ctx.j
574 elif all((t.imag < 0 or t.real > 0) for t in [x, y, z, p]):
575 margin = -ctx.j
576 else:
577 margin = 1
578 # Go through the upper half-plane, but low enough that any
579 # parameter starting in the lower plane doesn't cross the
580 # branch cut
581 for t in [x, y, z, p]:
582 if t.imag >= 0 or t.real > 0:
583 continue
584 margin = min(margin, abs(t.imag) * 0.5)
585 margin *= ctx.j
586 N += margin
587 F = lambda t: 1/(ctx.sqrt(t+x)*ctx.sqrt(t+y)*ctx.sqrt(t+z)*(t+p))
588 if integration == 2:
589 return 1.5 * ctx.quad(F, [0, N, ctx.inf])
590 initial_integral = 1.5 * ctx.quad(F, [0, N])
591 x += N; y += N; z += N; p += N
592 xm,ym,zm,pm = x,y,z,p
593 A0 = Am = (x + y + z + 2*p)/5
594 delta = (p-x)*(p-y)*(p-z)
595 Q = ctx.root(0.25*r, -6) * max(abs(A0-x),abs(A0-y),abs(A0-z),abs(A0-p))
596 m = 0
597 g = ctx.mpf(0.25)
598 pow4 = ctx.one
599 S = 0
600 while 1:
601 sx = ctx.sqrt(xm)
602 sy = ctx.sqrt(ym)
603 sz = ctx.sqrt(zm)
604 sp = ctx.sqrt(pm)
605 lm = sx*sy + sx*sz + sy*sz
606 Am1 = (Am+lm)*g
607 xm = (xm+lm)*g; ym = (ym+lm)*g; zm = (zm+lm)*g; pm = (pm+lm)*g
608 dm = (sp+sx) * (sp+sy) * (sp+sz)
609 em = delta * ctx.power(4, -3*m) / dm**2
610 if pow4 * Q < abs(Am):
611 break
612 T = RC_calc(ctx, ctx.one, ctx.one+em, r) * pow4 / dm
613 S += T
614 pow4 *= g
615 m += 1
616 Am = Am1
617 t = ctx.ldexp(1,-2*m) / Am
618 X = (A0-x)*t
619 Y = (A0-y)*t
620 Z = (A0-z)*t
621 P = (-X-Y-Z)/2
622 E2 = X*Y + X*Z + Y*Z - 3*P**2
623 E3 = X*Y*Z + 2*E2*P + 4*P**3
624 E4 = (2*X*Y*Z + E2*P + 3*P**3)*P
625 E5 = X*Y*Z*P**2
626 P = 24024 - 5148*E2 + 2457*E2**2 + 4004*E3 - 4158*E2*E3 - 3276*E4 + 2772*E5
627 Q = 24024
628 v1 = g**m * ctx.power(Am, -1.5) * P/Q
629 v2 = 6*S
630 return initial_integral + v1 + v2
632@defun
633def elliprf(ctx, x, y, z):
634 r"""
635 Evaluates the Carlson symmetric elliptic integral of the first kind
637 .. math ::
639 R_F(x,y,z) = \frac{1}{2}
640 \int_0^{\infty} \frac{dt}{\sqrt{(t+x)(t+y)(t+z)}}
642 which is defined for `x,y,z \notin (-\infty,0)`, and with
643 at most one of `x,y,z` being zero.
645 For real `x,y,z \ge 0`, the principal square root is taken in the integrand.
646 For complex `x,y,z`, the principal square root is taken as `t \to \infty`
647 and as `t \to 0` non-principal branches are chosen as necessary so as to
648 make the integrand continuous.
650 **Examples**
652 Some basic values and limits::
654 >>> from mpmath import *
655 >>> mp.dps = 25; mp.pretty = True
656 >>> elliprf(0,1,1); pi/2
657 1.570796326794896619231322
658 1.570796326794896619231322
659 >>> elliprf(0,1,inf)
660 0.0
661 >>> elliprf(1,1,1)
662 1.0
663 >>> elliprf(2,2,2)**2
664 0.5
665 >>> elliprf(1,0,0); elliprf(0,0,1); elliprf(0,1,0); elliprf(0,0,0)
666 +inf
667 +inf
668 +inf
669 +inf
671 Representing complete elliptic integrals in terms of `R_F`::
673 >>> m = mpf(0.75)
674 >>> ellipk(m); elliprf(0,1-m,1)
675 2.156515647499643235438675
676 2.156515647499643235438675
677 >>> ellipe(m); elliprf(0,1-m,1)-m*elliprd(0,1-m,1)/3
678 1.211056027568459524803563
679 1.211056027568459524803563
681 Some symmetries and argument transformations::
683 >>> x,y,z = 2,3,4
684 >>> elliprf(x,y,z); elliprf(y,x,z); elliprf(z,y,x)
685 0.5840828416771517066928492
686 0.5840828416771517066928492
687 0.5840828416771517066928492
688 >>> k = mpf(100000)
689 >>> elliprf(k*x,k*y,k*z); k**(-0.5) * elliprf(x,y,z)
690 0.001847032121923321253219284
691 0.001847032121923321253219284
692 >>> l = sqrt(x*y) + sqrt(y*z) + sqrt(z*x)
693 >>> elliprf(x,y,z); 2*elliprf(x+l,y+l,z+l)
694 0.5840828416771517066928492
695 0.5840828416771517066928492
696 >>> elliprf((x+l)/4,(y+l)/4,(z+l)/4)
697 0.5840828416771517066928492
699 Comparing with numerical integration::
701 >>> x,y,z = 2,3,4
702 >>> elliprf(x,y,z)
703 0.5840828416771517066928492
704 >>> f = lambda t: 0.5*((t+x)*(t+y)*(t+z))**(-0.5)
705 >>> q = extradps(25)(quad)
706 >>> q(f, [0,inf])
707 0.5840828416771517066928492
709 With the following arguments, the square root in the integrand becomes
710 discontinuous at `t = 1/2` if the principal branch is used. To obtain
711 the right value, `-\sqrt{r}` must be taken instead of `\sqrt{r}`
712 on `t \in (0, 1/2)`::
714 >>> x,y,z = j-1,j,0
715 >>> elliprf(x,y,z)
716 (0.7961258658423391329305694 - 1.213856669836495986430094j)
717 >>> -q(f, [0,0.5]) + q(f, [0.5,inf])
718 (0.7961258658423391329305694 - 1.213856669836495986430094j)
720 The so-called *first lemniscate constant*, a transcendental number::
722 >>> elliprf(0,1,2)
723 1.31102877714605990523242
724 >>> extradps(25)(quad)(lambda t: 1/sqrt(1-t**4), [0,1])
725 1.31102877714605990523242
726 >>> gamma('1/4')**2/(4*sqrt(2*pi))
727 1.31102877714605990523242
729 **References**
731 1. [Carlson]_
732 2. [DLMF]_ Chapter 19. Elliptic Integrals
734 """
735 x = ctx.convert(x)
736 y = ctx.convert(y)
737 z = ctx.convert(z)
738 prec = ctx.prec
739 try:
740 ctx.prec += 20
741 tol = ctx.eps * 2**10
742 v = RF_calc(ctx, x, y, z, tol)
743 finally:
744 ctx.prec = prec
745 return +v
747@defun
748def elliprc(ctx, x, y, pv=True):
749 r"""
750 Evaluates the degenerate Carlson symmetric elliptic integral
751 of the first kind
753 .. math ::
755 R_C(x,y) = R_F(x,y,y) =
756 \frac{1}{2} \int_0^{\infty} \frac{dt}{(t+y) \sqrt{(t+x)}}.
758 If `y \in (-\infty,0)`, either a value defined by continuity,
759 or with *pv=True* the Cauchy principal value, can be computed.
761 If `x \ge 0, y > 0`, the value can be expressed in terms of
762 elementary functions as
764 .. math ::
766 R_C(x,y) =
767 \begin{cases}
768 \dfrac{1}{\sqrt{y-x}}
769 \cos^{-1}\left(\sqrt{\dfrac{x}{y}}\right), & x < y \\
770 \dfrac{1}{\sqrt{y}}, & x = y \\
771 \dfrac{1}{\sqrt{x-y}}
772 \cosh^{-1}\left(\sqrt{\dfrac{x}{y}}\right), & x > y \\
773 \end{cases}.
775 **Examples**
777 Some special values and limits::
779 >>> from mpmath import *
780 >>> mp.dps = 25; mp.pretty = True
781 >>> elliprc(1,2)*4; elliprc(0,1)*2; +pi
782 3.141592653589793238462643
783 3.141592653589793238462643
784 3.141592653589793238462643
785 >>> elliprc(1,0)
786 +inf
787 >>> elliprc(5,5)**2
788 0.2
789 >>> elliprc(1,inf); elliprc(inf,1); elliprc(inf,inf)
790 0.0
791 0.0
792 0.0
794 Comparing with the elementary closed-form solution::
796 >>> elliprc('1/3', '1/5'); sqrt(7.5)*acosh(sqrt('5/3'))
797 2.041630778983498390751238
798 2.041630778983498390751238
799 >>> elliprc('1/5', '1/3'); sqrt(7.5)*acos(sqrt('3/5'))
800 1.875180765206547065111085
801 1.875180765206547065111085
803 Comparing with numerical integration::
805 >>> q = extradps(25)(quad)
806 >>> elliprc(2, -3, pv=True)
807 0.3333969101113672670749334
808 >>> elliprc(2, -3, pv=False)
809 (0.3333969101113672670749334 + 0.7024814731040726393156375j)
810 >>> 0.5*q(lambda t: 1/(sqrt(t+2)*(t-3)), [0,3-j,6,inf])
811 (0.3333969101113672670749334 + 0.7024814731040726393156375j)
813 """
814 x = ctx.convert(x)
815 y = ctx.convert(y)
816 prec = ctx.prec
817 try:
818 ctx.prec += 20
819 tol = ctx.eps * 2**10
820 v = RC_calc(ctx, x, y, tol, pv)
821 finally:
822 ctx.prec = prec
823 return +v
825@defun
826def elliprj(ctx, x, y, z, p, integration=1):
827 r"""
828 Evaluates the Carlson symmetric elliptic integral of the third kind
830 .. math ::
832 R_J(x,y,z,p) = \frac{3}{2}
833 \int_0^{\infty} \frac{dt}{(t+p)\sqrt{(t+x)(t+y)(t+z)}}.
835 Like :func:`~mpmath.elliprf`, the branch of the square root in the integrand
836 is defined so as to be continuous along the path of integration for
837 complex values of the arguments.
839 **Examples**
841 Some values and limits::
843 >>> from mpmath import *
844 >>> mp.dps = 25; mp.pretty = True
845 >>> elliprj(1,1,1,1)
846 1.0
847 >>> elliprj(2,2,2,2); 1/(2*sqrt(2))
848 0.3535533905932737622004222
849 0.3535533905932737622004222
850 >>> elliprj(0,1,2,2)
851 1.067937989667395702268688
852 >>> 3*(2*gamma('5/4')**2-pi**2/gamma('1/4')**2)/(sqrt(2*pi))
853 1.067937989667395702268688
854 >>> elliprj(0,1,1,2); 3*pi*(2-sqrt(2))/4
855 1.380226776765915172432054
856 1.380226776765915172432054
857 >>> elliprj(1,3,2,0); elliprj(0,1,1,0); elliprj(0,0,0,0)
858 +inf
859 +inf
860 +inf
861 >>> elliprj(1,inf,1,0); elliprj(1,1,1,inf)
862 0.0
863 0.0
864 >>> chop(elliprj(1+j, 1-j, 1, 1))
865 0.8505007163686739432927844
867 Scale transformation::
869 >>> x,y,z,p = 2,3,4,5
870 >>> k = mpf(100000)
871 >>> elliprj(k*x,k*y,k*z,k*p); k**(-1.5)*elliprj(x,y,z,p)
872 4.521291677592745527851168e-9
873 4.521291677592745527851168e-9
875 Comparing with numerical integration::
877 >>> elliprj(1,2,3,4)
878 0.2398480997495677621758617
879 >>> f = lambda t: 1/((t+4)*sqrt((t+1)*(t+2)*(t+3)))
880 >>> 1.5*quad(f, [0,inf])
881 0.2398480997495677621758617
882 >>> elliprj(1,2+1j,3,4-2j)
883 (0.216888906014633498739952 + 0.04081912627366673332369512j)
884 >>> f = lambda t: 1/((t+4-2j)*sqrt((t+1)*(t+2+1j)*(t+3)))
885 >>> 1.5*quad(f, [0,inf])
886 (0.216888906014633498739952 + 0.04081912627366673332369511j)
888 """
889 x = ctx.convert(x)
890 y = ctx.convert(y)
891 z = ctx.convert(z)
892 p = ctx.convert(p)
893 prec = ctx.prec
894 try:
895 ctx.prec += 20
896 tol = ctx.eps * 2**10
897 v = RJ_calc(ctx, x, y, z, p, tol, integration)
898 finally:
899 ctx.prec = prec
900 return +v
902@defun
903def elliprd(ctx, x, y, z):
904 r"""
905 Evaluates the degenerate Carlson symmetric elliptic integral
906 of the third kind or Carlson elliptic integral of the
907 second kind `R_D(x,y,z) = R_J(x,y,z,z)`.
909 See :func:`~mpmath.elliprj` for additional information.
911 **Examples**
913 >>> from mpmath import *
914 >>> mp.dps = 25; mp.pretty = True
915 >>> elliprd(1,2,3)
916 0.2904602810289906442326534
917 >>> elliprj(1,2,3,3)
918 0.2904602810289906442326534
920 The so-called *second lemniscate constant*, a transcendental number::
922 >>> elliprd(0,2,1)/3
923 0.5990701173677961037199612
924 >>> extradps(25)(quad)(lambda t: t**2/sqrt(1-t**4), [0,1])
925 0.5990701173677961037199612
926 >>> gamma('3/4')**2/sqrt(2*pi)
927 0.5990701173677961037199612
929 """
930 return ctx.elliprj(x,y,z,z)
932@defun
933def elliprg(ctx, x, y, z):
934 r"""
935 Evaluates the Carlson completely symmetric elliptic integral
936 of the second kind
938 .. math ::
940 R_G(x,y,z) = \frac{1}{4} \int_0^{\infty}
941 \frac{t}{\sqrt{(t+x)(t+y)(t+z)}}
942 \left( \frac{x}{t+x} + \frac{y}{t+y} + \frac{z}{t+z}\right) dt.
944 **Examples**
946 Evaluation for real and complex arguments::
948 >>> from mpmath import *
949 >>> mp.dps = 25; mp.pretty = True
950 >>> elliprg(0,1,1)*4; +pi
951 3.141592653589793238462643
952 3.141592653589793238462643
953 >>> elliprg(0,0.5,1)
954 0.6753219405238377512600874
955 >>> chop(elliprg(1+j, 1-j, 2))
956 1.172431327676416604532822
958 A double integral that can be evaluated in terms of `R_G`::
960 >>> x,y,z = 2,3,4
961 >>> def f(t,u):
962 ... st = fp.sin(t); ct = fp.cos(t)
963 ... su = fp.sin(u); cu = fp.cos(u)
964 ... return (x*(st*cu)**2 + y*(st*su)**2 + z*ct**2)**0.5 * st
965 ...
966 >>> nprint(mpf(fp.quad(f, [0,fp.pi], [0,2*fp.pi])/(4*fp.pi)), 13)
967 1.725503028069
968 >>> nprint(elliprg(x,y,z), 13)
969 1.725503028069
971 """
972 x = ctx.convert(x)
973 y = ctx.convert(y)
974 z = ctx.convert(z)
975 zeros = (not x) + (not y) + (not z)
976 if zeros == 3:
977 return (x+y+z)*0
978 if zeros == 2:
979 if x: return 0.5*ctx.sqrt(x)
980 if y: return 0.5*ctx.sqrt(y)
981 return 0.5*ctx.sqrt(z)
982 if zeros == 1:
983 if not z:
984 x, z = z, x
985 def terms():
986 T1 = 0.5*z*ctx.elliprf(x,y,z)
987 T2 = -0.5*(x-z)*(y-z)*ctx.elliprd(x,y,z)/3
988 T3 = 0.5*ctx.sqrt(x)*ctx.sqrt(y)/ctx.sqrt(z)
989 return T1,T2,T3
990 return ctx.sum_accurately(terms)
993@defun_wrapped
994def ellipf(ctx, phi, m):
995 r"""
996 Evaluates the Legendre incomplete elliptic integral of the first kind
998 .. math ::
1000 F(\phi,m) = \int_0^{\phi} \frac{dt}{\sqrt{1-m \sin^2 t}}
1002 or equivalently
1004 .. math ::
1006 F(\phi,m) = \int_0^{\sin \phi}
1007 \frac{dt}{\left(\sqrt{1-t^2}\right)\left(\sqrt{1-mt^2}\right)}.
1009 The function reduces to a complete elliptic integral of the first kind
1010 (see :func:`~mpmath.ellipk`) when `\phi = \frac{\pi}{2}`; that is,
1012 .. math ::
1014 F\left(\frac{\pi}{2}, m\right) = K(m).
1016 In the defining integral, it is assumed that the principal branch
1017 of the square root is taken and that the path of integration avoids
1018 crossing any branch cuts. Outside `-\pi/2 \le \Re(\phi) \le \pi/2`,
1019 the function extends quasi-periodically as
1021 .. math ::
1023 F(\phi + n \pi, m) = 2 n K(m) + F(\phi,m), n \in \mathbb{Z}.
1025 **Plots**
1027 .. literalinclude :: /plots/ellipf.py
1028 .. image :: /plots/ellipf.png
1030 **Examples**
1032 Basic values and limits::
1034 >>> from mpmath import *
1035 >>> mp.dps = 25; mp.pretty = True
1036 >>> ellipf(0,1)
1037 0.0
1038 >>> ellipf(0,0)
1039 0.0
1040 >>> ellipf(1,0); ellipf(2+3j,0)
1041 1.0
1042 (2.0 + 3.0j)
1043 >>> ellipf(1,1); log(sec(1)+tan(1))
1044 1.226191170883517070813061
1045 1.226191170883517070813061
1046 >>> ellipf(pi/2, -0.5); ellipk(-0.5)
1047 1.415737208425956198892166
1048 1.415737208425956198892166
1049 >>> ellipf(pi/2+eps, 1); ellipf(-pi/2-eps, 1)
1050 +inf
1051 +inf
1052 >>> ellipf(1.5, 1)
1053 3.340677542798311003320813
1055 Comparing with numerical integration::
1057 >>> z,m = 0.5, 1.25
1058 >>> ellipf(z,m)
1059 0.5287219202206327872978255
1060 >>> quad(lambda t: (1-m*sin(t)**2)**(-0.5), [0,z])
1061 0.5287219202206327872978255
1063 The arguments may be complex numbers::
1065 >>> ellipf(3j, 0.5)
1066 (0.0 + 1.713602407841590234804143j)
1067 >>> ellipf(3+4j, 5-6j)
1068 (1.269131241950351323305741 - 0.3561052815014558335412538j)
1069 >>> z,m = 2+3j, 1.25
1070 >>> k = 1011
1071 >>> ellipf(z+pi*k,m); ellipf(z,m) + 2*k*ellipk(m)
1072 (4086.184383622179764082821 - 3003.003538923749396546871j)
1073 (4086.184383622179764082821 - 3003.003538923749396546871j)
1075 For `|\Re(z)| < \pi/2`, the function can be expressed as a
1076 hypergeometric series of two variables
1077 (see :func:`~mpmath.appellf1`)::
1079 >>> z,m = 0.5, 0.25
1080 >>> ellipf(z,m)
1081 0.5050887275786480788831083
1082 >>> sin(z)*appellf1(0.5,0.5,0.5,1.5,sin(z)**2,m*sin(z)**2)
1083 0.5050887275786480788831083
1085 """
1086 z = phi
1087 if not (ctx.isnormal(z) and ctx.isnormal(m)):
1088 if m == 0:
1089 return z + m
1090 if z == 0:
1091 return z * m
1092 if m == ctx.inf or m == ctx.ninf: return z/m
1093 raise ValueError
1094 x = z.real
1095 ctx.prec += max(0, ctx.mag(x))
1096 pi = +ctx.pi
1097 away = abs(x) > pi/2
1098 if m == 1:
1099 if away:
1100 return ctx.inf
1101 if away:
1102 d = ctx.nint(x/pi)
1103 z = z-pi*d
1104 P = 2*d*ctx.ellipk(m)
1105 else:
1106 P = 0
1107 c, s = ctx.cos_sin(z)
1108 return s * ctx.elliprf(c**2, 1-m*s**2, 1) + P
1110@defun_wrapped
1111def ellipe(ctx, *args):
1112 r"""
1113 Called with a single argument `m`, evaluates the Legendre complete
1114 elliptic integral of the second kind, `E(m)`, defined by
1116 .. math :: E(m) = \int_0^{\pi/2} \sqrt{1-m \sin^2 t} \, dt \,=\,
1117 \frac{\pi}{2}
1118 \,_2F_1\left(\frac{1}{2}, -\frac{1}{2}, 1, m\right).
1120 Called with two arguments `\phi, m`, evaluates the incomplete elliptic
1121 integral of the second kind
1123 .. math ::
1125 E(\phi,m) = \int_0^{\phi} \sqrt{1-m \sin^2 t} \, dt =
1126 \int_0^{\sin z}
1127 \frac{\sqrt{1-mt^2}}{\sqrt{1-t^2}} \, dt.
1129 The incomplete integral reduces to a complete integral when
1130 `\phi = \frac{\pi}{2}`; that is,
1132 .. math ::
1134 E\left(\frac{\pi}{2}, m\right) = E(m).
1136 In the defining integral, it is assumed that the principal branch
1137 of the square root is taken and that the path of integration avoids
1138 crossing any branch cuts. Outside `-\pi/2 \le \Re(z) \le \pi/2`,
1139 the function extends quasi-periodically as
1141 .. math ::
1143 E(\phi + n \pi, m) = 2 n E(m) + E(\phi,m), n \in \mathbb{Z}.
1145 **Plots**
1147 .. literalinclude :: /plots/ellipe.py
1148 .. image :: /plots/ellipe.png
1150 **Examples for the complete integral**
1152 Basic values and limits::
1154 >>> from mpmath import *
1155 >>> mp.dps = 25; mp.pretty = True
1156 >>> ellipe(0)
1157 1.570796326794896619231322
1158 >>> ellipe(1)
1159 1.0
1160 >>> ellipe(-1)
1161 1.910098894513856008952381
1162 >>> ellipe(2)
1163 (0.5990701173677961037199612 + 0.5990701173677961037199612j)
1164 >>> ellipe(inf)
1165 (0.0 + +infj)
1166 >>> ellipe(-inf)
1167 +inf
1169 Verifying the defining integral and hypergeometric
1170 representation::
1172 >>> ellipe(0.5)
1173 1.350643881047675502520175
1174 >>> quad(lambda t: sqrt(1-0.5*sin(t)**2), [0, pi/2])
1175 1.350643881047675502520175
1176 >>> pi/2*hyp2f1(0.5,-0.5,1,0.5)
1177 1.350643881047675502520175
1179 Evaluation is supported for arbitrary complex `m`::
1181 >>> ellipe(0.5+0.25j)
1182 (1.360868682163129682716687 - 0.1238733442561786843557315j)
1183 >>> ellipe(3+4j)
1184 (1.499553520933346954333612 - 1.577879007912758274533309j)
1186 A definite integral::
1188 >>> quad(ellipe, [0,1])
1189 1.333333333333333333333333
1191 **Examples for the incomplete integral**
1193 Basic values and limits::
1195 >>> ellipe(0,1)
1196 0.0
1197 >>> ellipe(0,0)
1198 0.0
1199 >>> ellipe(1,0)
1200 1.0
1201 >>> ellipe(2+3j,0)
1202 (2.0 + 3.0j)
1203 >>> ellipe(1,1); sin(1)
1204 0.8414709848078965066525023
1205 0.8414709848078965066525023
1206 >>> ellipe(pi/2, -0.5); ellipe(-0.5)
1207 1.751771275694817862026502
1208 1.751771275694817862026502
1209 >>> ellipe(pi/2, 1); ellipe(-pi/2, 1)
1210 1.0
1211 -1.0
1212 >>> ellipe(1.5, 1)
1213 0.9974949866040544309417234
1215 Comparing with numerical integration::
1217 >>> z,m = 0.5, 1.25
1218 >>> ellipe(z,m)
1219 0.4740152182652628394264449
1220 >>> quad(lambda t: sqrt(1-m*sin(t)**2), [0,z])
1221 0.4740152182652628394264449
1223 The arguments may be complex numbers::
1225 >>> ellipe(3j, 0.5)
1226 (0.0 + 7.551991234890371873502105j)
1227 >>> ellipe(3+4j, 5-6j)
1228 (24.15299022574220502424466 + 75.2503670480325997418156j)
1229 >>> k = 35
1230 >>> z,m = 2+3j, 1.25
1231 >>> ellipe(z+pi*k,m); ellipe(z,m) + 2*k*ellipe(m)
1232 (48.30138799412005235090766 + 17.47255216721987688224357j)
1233 (48.30138799412005235090766 + 17.47255216721987688224357j)
1235 For `|\Re(z)| < \pi/2`, the function can be expressed as a
1236 hypergeometric series of two variables
1237 (see :func:`~mpmath.appellf1`)::
1239 >>> z,m = 0.5, 0.25
1240 >>> ellipe(z,m)
1241 0.4950017030164151928870375
1242 >>> sin(z)*appellf1(0.5,0.5,-0.5,1.5,sin(z)**2,m*sin(z)**2)
1243 0.4950017030164151928870376
1245 """
1246 if len(args) == 1:
1247 return ctx._ellipe(args[0])
1248 else:
1249 phi, m = args
1250 z = phi
1251 if not (ctx.isnormal(z) and ctx.isnormal(m)):
1252 if m == 0:
1253 return z + m
1254 if z == 0:
1255 return z * m
1256 if m == ctx.inf or m == ctx.ninf:
1257 return ctx.inf
1258 raise ValueError
1259 x = z.real
1260 ctx.prec += max(0, ctx.mag(x))
1261 pi = +ctx.pi
1262 away = abs(x) > pi/2
1263 if away:
1264 d = ctx.nint(x/pi)
1265 z = z-pi*d
1266 P = 2*d*ctx.ellipe(m)
1267 else:
1268 P = 0
1269 def terms():
1270 c, s = ctx.cos_sin(z)
1271 x = c**2
1272 y = 1-m*s**2
1273 RF = ctx.elliprf(x, y, 1)
1274 RD = ctx.elliprd(x, y, 1)
1275 return s*RF, -m*s**3*RD/3
1276 return ctx.sum_accurately(terms) + P
1278@defun_wrapped
1279def ellippi(ctx, *args):
1280 r"""
1281 Called with three arguments `n, \phi, m`, evaluates the Legendre
1282 incomplete elliptic integral of the third kind
1284 .. math ::
1286 \Pi(n; \phi, m) = \int_0^{\phi}
1287 \frac{dt}{(1-n \sin^2 t) \sqrt{1-m \sin^2 t}} =
1288 \int_0^{\sin \phi}
1289 \frac{dt}{(1-nt^2) \sqrt{1-t^2} \sqrt{1-mt^2}}.
1291 Called with two arguments `n, m`, evaluates the complete
1292 elliptic integral of the third kind
1293 `\Pi(n,m) = \Pi(n; \frac{\pi}{2},m)`.
1295 In the defining integral, it is assumed that the principal branch
1296 of the square root is taken and that the path of integration avoids
1297 crossing any branch cuts. Outside `-\pi/2 \le \Re(\phi) \le \pi/2`,
1298 the function extends quasi-periodically as
1300 .. math ::
1302 \Pi(n,\phi+k\pi,m) = 2k\Pi(n,m) + \Pi(n,\phi,m), k \in \mathbb{Z}.
1304 **Plots**
1306 .. literalinclude :: /plots/ellippi.py
1307 .. image :: /plots/ellippi.png
1309 **Examples for the complete integral**
1311 Some basic values and limits::
1313 >>> from mpmath import *
1314 >>> mp.dps = 25; mp.pretty = True
1315 >>> ellippi(0,-5); ellipk(-5)
1316 0.9555039270640439337379334
1317 0.9555039270640439337379334
1318 >>> ellippi(inf,2)
1319 0.0
1320 >>> ellippi(2,inf)
1321 0.0
1322 >>> abs(ellippi(1,5))
1323 +inf
1324 >>> abs(ellippi(0.25,1))
1325 +inf
1327 Evaluation in terms of simpler functions::
1329 >>> ellippi(0.25,0.25); ellipe(0.25)/(1-0.25)
1330 1.956616279119236207279727
1331 1.956616279119236207279727
1332 >>> ellippi(3,0); pi/(2*sqrt(-2))
1333 (0.0 - 1.11072073453959156175397j)
1334 (0.0 - 1.11072073453959156175397j)
1335 >>> ellippi(-3,0); pi/(2*sqrt(4))
1336 0.7853981633974483096156609
1337 0.7853981633974483096156609
1339 **Examples for the incomplete integral**
1341 Basic values and limits::
1343 >>> ellippi(0.25,-0.5); ellippi(0.25,pi/2,-0.5)
1344 1.622944760954741603710555
1345 1.622944760954741603710555
1346 >>> ellippi(1,0,1)
1347 0.0
1348 >>> ellippi(inf,0,1)
1349 0.0
1350 >>> ellippi(0,0.25,0.5); ellipf(0.25,0.5)
1351 0.2513040086544925794134591
1352 0.2513040086544925794134591
1353 >>> ellippi(1,1,1); (log(sec(1)+tan(1))+sec(1)*tan(1))/2
1354 2.054332933256248668692452
1355 2.054332933256248668692452
1356 >>> ellippi(0.25, 53*pi/2, 0.75); 53*ellippi(0.25,0.75)
1357 135.240868757890840755058
1358 135.240868757890840755058
1359 >>> ellippi(0.5,pi/4,0.5); 2*ellipe(pi/4,0.5)-1/sqrt(3)
1360 0.9190227391656969903987269
1361 0.9190227391656969903987269
1363 Complex arguments are supported::
1365 >>> ellippi(0.5, 5+6j-2*pi, -7-8j)
1366 (-0.3612856620076747660410167 + 0.5217735339984807829755815j)
1368 Some degenerate cases::
1370 >>> ellippi(1,1)
1371 +inf
1372 >>> ellippi(1,0)
1373 +inf
1374 >>> ellippi(1,2,0)
1375 +inf
1376 >>> ellippi(1,2,1)
1377 +inf
1378 >>> ellippi(1,0,1)
1379 0.0
1381 """
1382 if len(args) == 2:
1383 n, m = args
1384 complete = True
1385 z = phi = ctx.pi/2
1386 else:
1387 n, phi, m = args
1388 complete = False
1389 z = phi
1390 if not (ctx.isnormal(n) and ctx.isnormal(z) and ctx.isnormal(m)):
1391 if ctx.isnan(n) or ctx.isnan(z) or ctx.isnan(m):
1392 raise ValueError
1393 if complete:
1394 if m == 0:
1395 if n == 1:
1396 return ctx.inf
1397 return ctx.pi/(2*ctx.sqrt(1-n))
1398 if n == 0: return ctx.ellipk(m)
1399 if ctx.isinf(n) or ctx.isinf(m): return ctx.zero
1400 else:
1401 if z == 0: return z
1402 if ctx.isinf(n): return ctx.zero
1403 if ctx.isinf(m): return ctx.zero
1404 if ctx.isinf(n) or ctx.isinf(z) or ctx.isinf(m):
1405 raise ValueError
1406 if complete:
1407 if m == 1:
1408 if n == 1:
1409 return ctx.inf
1410 return -ctx.inf/ctx.sign(n-1)
1411 away = False
1412 else:
1413 x = z.real
1414 ctx.prec += max(0, ctx.mag(x))
1415 pi = +ctx.pi
1416 away = abs(x) > pi/2
1417 if away:
1418 d = ctx.nint(x/pi)
1419 z = z-pi*d
1420 P = 2*d*ctx.ellippi(n,m)
1421 if ctx.isinf(P):
1422 return ctx.inf
1423 else:
1424 P = 0
1425 def terms():
1426 if complete:
1427 c, s = ctx.zero, ctx.one
1428 else:
1429 c, s = ctx.cos_sin(z)
1430 x = c**2
1431 y = 1-m*s**2
1432 RF = ctx.elliprf(x, y, 1)
1433 RJ = ctx.elliprj(x, y, 1, 1-n*s**2)
1434 return s*RF, n*s**3*RJ/3
1435 return ctx.sum_accurately(terms) + P