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
7
8.. math ::
9
10    f(z + 2 \omega_1) = f(z + 2 \omega_2) = f(z)
11
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.
17
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:
29
30.. math ::
31
32    m = k^2
33
34.. math ::
35
36    \tau = i \frac{K(1-m)}{K(m)}
37
38.. math ::
39
40    q = e^{i \pi \tau}
41
42.. math ::
43
44    k = \frac{\vartheta_2^2(q)}{\vartheta_3^2(q)}
45
46In addition, an alternative definition is used for the nome in
47number theory, which we here denote by q-bar:
48
49.. math ::
50
51    \bar{q} = q^2 = e^{2 i \pi \tau}
52
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`).
56
57**References**
58
591. [AbramowitzStegun]_
60
612. [WhittakerWatson]_
62
63"""
64
65from .functions import defun, defun_wrapped
66
67@defun_wrapped
68def eta(ctx, tau):
69    r"""
70    Returns the Dedekind eta function of tau in the upper half-plane.
71
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
87
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)
93
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
118
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}`::
123
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
136
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)
148
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}`::
154
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)
167
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)
179
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}`::
185
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)
198
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)
213
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}`::
219
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
232
233    As `q \to 1` and `q \to -1`, `k` rapidly approaches
234    `1` and `i \infty` respectively::
235
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
258
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}`::
264
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
277
278    As `q \to 1` and `q \to -1`, `m` rapidly approaches
279    `1` and `-\infty` respectively::
280
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
289
290    The inverse nome as a function of `q` has an integer
291    Taylor series expansion::
292
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]
295
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
313
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}
332
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
370
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
376
377    .. math ::
378
379        J(\tau) = \frac{g_2^3(\tau)}{g_2^3(\tau) - 27 g_3^2(\tau)}
380
381    where `g_2` and `g_3` are the modular invariants of the Weierstrass
382    elliptic function,
383
384    .. math ::
385
386        g_2(\tau) = 60 \sum_{(m,n) \in \mathbb{Z}^2 \setminus (0,0)} (m \tau+n)^{-4}
387
388        g_3(\tau) = 140 \sum_{(m,n) \in \mathbb{Z}^2 \setminus (0,0)} (m \tau+n)^{-6}.
389
390    An alternative, common notation is that of the j-function
391    `j(\tau) = 1728 J(\tau)`.
392
393    **Plots**
394
395    .. literalinclude :: /plots/kleinj.py
396    .. image :: /plots/kleinj.png
397    .. literalinclude :: /plots/kleinj2.py
398    .. image :: /plots/kleinj2.png
399
400    **Examples**
401
402    Verifying the functional equation `J(\tau) = J(\tau+1) = J(-\tau^{-1})`::
403
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)
414
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`::
417
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]
421
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::
424
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)
445
446    Also at other special points, the j-function assumes explicit
447    algebraic values, e.g.::
448
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
455
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
464
465
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
500
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
533
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
631
632@defun
633def elliprf(ctx, x, y, z):
634    r"""
635    Evaluates the Carlson symmetric elliptic integral of the first kind
636
637    .. math ::
638
639        R_F(x,y,z) = \frac{1}{2}
640            \int_0^{\infty} \frac{dt}{\sqrt{(t+x)(t+y)(t+z)}}
641
642    which is defined for `x,y,z \notin (-\infty,0)`, and with
643    at most one of `x,y,z` being zero.
644
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.
649
650    **Examples**
651
652    Some basic values and limits::
653
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
670
671    Representing complete elliptic integrals in terms of `R_F`::
672
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
680
681    Some symmetries and argument transformations::
682
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
698
699    Comparing with numerical integration::
700
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
708
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)`::
713
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)
719
720    The so-called *first lemniscate constant*, a transcendental number::
721
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
728
729    **References**
730
731    1. [Carlson]_
732    2. [DLMF]_ Chapter 19. Elliptic Integrals
733
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
746
747@defun
748def elliprc(ctx, x, y, pv=True):
749    r"""
750    Evaluates the degenerate Carlson symmetric elliptic integral
751    of the first kind
752
753    .. math ::
754
755        R_C(x,y) = R_F(x,y,y) =
756            \frac{1}{2} \int_0^{\infty} \frac{dt}{(t+y) \sqrt{(t+x)}}.
757
758    If `y \in (-\infty,0)`, either a value defined by continuity,
759    or with *pv=True* the Cauchy principal value, can be computed.
760
761    If `x \ge 0, y > 0`, the value can be expressed in terms of
762    elementary functions as
763
764    .. math ::
765
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}.
774
775    **Examples**
776
777    Some special values and limits::
778
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
793
794    Comparing with the elementary closed-form solution::
795
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
802
803    Comparing with numerical integration::
804
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)
812
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
824
825@defun
826def elliprj(ctx, x, y, z, p, integration=1):
827    r"""
828    Evaluates the Carlson symmetric elliptic integral of the third kind
829
830    .. math ::
831
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)}}.
834
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.
838
839    **Examples**
840
841    Some values and limits::
842
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
866
867    Scale transformation::
868
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
874
875    Comparing with numerical integration::
876
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)
887
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
901
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)`.
908
909    See :func:`~mpmath.elliprj` for additional information.
910
911    **Examples**
912
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
919
920    The so-called *second lemniscate constant*, a transcendental number::
921
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
928
929    """
930    return ctx.elliprj(x,y,z,z)
931
932@defun
933def elliprg(ctx, x, y, z):
934    r"""
935    Evaluates the Carlson completely symmetric elliptic integral
936    of the second kind
937
938    .. math ::
939
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.
943
944    **Examples**
945
946    Evaluation for real and complex arguments::
947
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
957
958    A double integral that can be evaluated in terms of `R_G`::
959
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
970
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)
991
992
993@defun_wrapped
994def ellipf(ctx, phi, m):
995    r"""
996    Evaluates the Legendre incomplete elliptic integral of the first kind
997
998     .. math ::
999
1000        F(\phi,m) = \int_0^{\phi} \frac{dt}{\sqrt{1-m \sin^2 t}}
1001
1002    or equivalently
1003
1004    .. math ::
1005
1006        F(\phi,m) = \int_0^{\sin \phi}
1007        \frac{dt}{\left(\sqrt{1-t^2}\right)\left(\sqrt{1-mt^2}\right)}.
1008
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,
1011
1012    .. math ::
1013
1014        F\left(\frac{\pi}{2}, m\right) = K(m).
1015
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
1020
1021    .. math ::
1022
1023        F(\phi + n \pi, m) = 2 n K(m) + F(\phi,m), n \in \mathbb{Z}.
1024
1025    **Plots**
1026
1027    .. literalinclude :: /plots/ellipf.py
1028    .. image :: /plots/ellipf.png
1029
1030    **Examples**
1031
1032    Basic values and limits::
1033
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
1054
1055    Comparing with numerical integration::
1056
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
1062
1063    The arguments may be complex numbers::
1064
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)
1074
1075    For `|\Re(z)| < \pi/2`, the function can be expressed as a
1076    hypergeometric series of two variables
1077    (see :func:`~mpmath.appellf1`)::
1078
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
1084
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
1109
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
1115
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).
1119
1120    Called with two arguments `\phi, m`, evaluates the incomplete elliptic
1121    integral of the second kind
1122
1123     .. math ::
1124
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.
1128
1129    The incomplete integral reduces to a complete integral when
1130    `\phi = \frac{\pi}{2}`; that is,
1131
1132    .. math ::
1133
1134        E\left(\frac{\pi}{2}, m\right) = E(m).
1135
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
1140
1141    .. math ::
1142
1143        E(\phi + n \pi, m) = 2 n E(m) + E(\phi,m), n \in \mathbb{Z}.
1144
1145    **Plots**
1146
1147    .. literalinclude :: /plots/ellipe.py
1148    .. image :: /plots/ellipe.png
1149
1150    **Examples for the complete integral**
1151
1152    Basic values and limits::
1153
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
1168
1169    Verifying the defining integral and hypergeometric
1170    representation::
1171
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
1178
1179    Evaluation is supported for arbitrary complex `m`::
1180
1181        >>> ellipe(0.5+0.25j)
1182        (1.360868682163129682716687 - 0.1238733442561786843557315j)
1183        >>> ellipe(3+4j)
1184        (1.499553520933346954333612 - 1.577879007912758274533309j)
1185
1186    A definite integral::
1187
1188        >>> quad(ellipe, [0,1])
1189        1.333333333333333333333333
1190
1191    **Examples for the incomplete integral**
1192
1193    Basic values and limits::
1194
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
1214
1215    Comparing with numerical integration::
1216
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
1222
1223    The arguments may be complex numbers::
1224
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)
1234
1235    For `|\Re(z)| < \pi/2`, the function can be expressed as a
1236    hypergeometric series of two variables
1237    (see :func:`~mpmath.appellf1`)::
1238
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
1244
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
1277
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
1283
1284    .. math ::
1285
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}}.
1290
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)`.
1294
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
1299
1300    .. math ::
1301
1302        \Pi(n,\phi+k\pi,m) = 2k\Pi(n,m) + \Pi(n,\phi,m), k \in \mathbb{Z}.
1303
1304    **Plots**
1305
1306    .. literalinclude :: /plots/ellippi.py
1307    .. image :: /plots/ellippi.png
1308
1309    **Examples for the complete integral**
1310
1311    Some basic values and limits::
1312
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
1326
1327    Evaluation in terms of simpler functions::
1328
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
1338
1339    **Examples for the incomplete integral**
1340
1341    Basic values and limits::
1342
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
1362
1363    Complex arguments are supported::
1364
1365        >>> ellippi(0.5, 5+6j-2*pi, -7-8j)
1366        (-0.3612856620076747660410167 + 0.5217735339984807829755815j)
1367
1368    Some degenerate cases::
1369
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
1380
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
1436