Ticket #14996: trac14996_4.patch

File trac14996_4.patch, 77.2 KB (added by eviatarbach, 8 years ago)
  • doc/en/reference/functions/index.rst

    # HG changeset patch
    # User Eviatar Bach <eviatarbach@gmail.com>
    # Date 1376452834 25200
    # Node ID 2e47a5b8215379b15fd17c1bfb51a0bda0962dae
    # Parent  8c4989a7e8d20a16ad1c464620ee616c4cde719a
    Trac #14996: Rewriting Jacobi elliptic function code, new numeric evaluation, Jacobi amplitude function
    
    diff --git a/doc/en/reference/functions/index.rst b/doc/en/reference/functions/index.rst
    a b  
    1111   sage/functions/orthogonal_polys
    1212   sage/functions/other
    1313   sage/functions/special
     14   sage/functions/jacobi
    1415   sage/functions/bessel
    1516   sage/functions/exp_integral
    1617   sage/functions/wigner
  • sage/functions/all.py

    diff --git a/sage/functions/all.py b/sage/functions/all.py
    a b  
    3131from special import (hypergeometric_U,
    3232                     spherical_bessel_J, spherical_bessel_Y,
    3333                     spherical_hankel1, spherical_hankel2,
    34                      spherical_harmonic, jacobi,
    35                      inverse_jacobi,
     34                     spherical_harmonic,
    3635                     lngamma, error_fcn, elliptic_e,
    3736                     elliptic_f, elliptic_ec, elliptic_eu,
    3837                     elliptic_kc, elliptic_pi, elliptic_j,
    3938                     airy_ai, airy_bi)
    4039
     40from jacobi import (jacobi, inverse_jacobi, jacobi_nd, jacobi_ns, jacobi_nc,
     41                    jacobi_dn, jacobi_ds, jacobi_dc, jacobi_sn, jacobi_sd,
     42                    jacobi_sc, jacobi_cn, jacobi_cd, jacobi_cs, jacobi_am,
     43                    inverse_jacobi_nd, inverse_jacobi_ns, inverse_jacobi_nc,
     44                    inverse_jacobi_dn, inverse_jacobi_ds, inverse_jacobi_dc,
     45                    inverse_jacobi_sn, inverse_jacobi_sd, inverse_jacobi_sc,
     46                    inverse_jacobi_cn, inverse_jacobi_cd, inverse_jacobi_cs)
     47
    4148from orthogonal_polys import (chebyshev_T,
    4249                              chebyshev_U,
    4350                              gen_laguerre,
  • new file sage/functions/jacobi.py

    diff --git a/sage/functions/jacobi.py b/sage/functions/jacobi.py
    new file mode 100644
    - +  
     1r"""
     2Jacobi elliptic functions
     3
     4This module implements the 12 Jacobi elliptic functions, along with their
     5inverses and the Jacobi amplitude function.
     6
     7Jacobi elliptic functions can be thought of as generalizations
     8of both ordinary and hyperbolic trig functions. There are twelve
     9Jacobian elliptic functions. Each of the twelve corresponds to an
     10arrow drawn from one corner of a rectangle to another.
     11
     12::
     13
     14                n ------------------- d
     15                |                     |
     16                |                     |
     17                |                     |
     18                s ------------------- c
     19
     20Each of the corners of the rectangle are labeled, by convention, ``s``,
     21``c``, ``d``, and ``n``. The rectangle is understood to be lying on the complex
     22plane, so that ``s`` is at the origin, ``c`` is on the real axis, and ``n`` is
     23on the imaginary axis. The twelve Jacobian elliptic functions are
     24then $\operatorname{pq}(x)$, where ``p`` and ``q`` are one of the letters
     25``s``, ``c``, ``d``, ``n``.
     26
     27The Jacobian elliptic functions are then the unique
     28doubly-periodic, meromorphic functions satisfying the following
     29three properties:
     30
     31#. There is a simple zero at the corner ``p``, and a simple pole at the
     32   corner ``q``.
     33#. The step from ``p`` to ``q`` is equal to half the period of the function
     34   `\operatorname{pq}(x)`; that is, the function `\operatorname{pq}(x)` is
     35   periodic in the direction ``pq``, with the period being twice the distance
     36   from ``p`` to ``q``. `\operatorname{pq}(x)` is periodic in the other two
     37   directions as well, with a period such that the distance from ``p`` to one
     38   of the other corners is a quarter period.
     39#. If the function $\operatorname{pq}(x)$ is expanded in terms of $x$ at one of
     40   the corners, the leading term in the expansion has a coefficient of 1.
     41   In other words, the leading term of the expansion of `\operatorname{pq}(x)`
     42   at the corner ``p`` is $x$; the leading term of the expansion at the corner
     43   ``q`` is $1/x$, and the leading term of an expansion at the other two
     44   corners is 1.
     45
     46We can write
     47
     48.. math::
     49
     50    \operatorname{pq}(x)=\frac{\operatorname{pr}(x)}{\operatorname{qr}(x)}
     51
     52where ``p``, ``q``, and ``r`` are any of the
     53letters ``s``, ``c``, ``d``, ``n``, with
     54the understanding that `\mathrm{ss}=\mathrm{cc}=\mathrm{dd}=\mathrm{nn}=1`.
     55
     56Let
     57
     58.. math::
     59
     60    u=\int_0^\phi \frac{d\theta} {\sqrt {1-m \sin^2 \theta}}
     61
     62
     63Then the *Jacobi elliptic function* `\operatorname{sn}(u)` is given by
     64
     65.. math::
     66
     67    \operatorname{sn}{u} = \sin{\phi}
     68
     69and `\operatorname{cn}(u)` is given by
     70
     71.. math::
     72
     73    \operatorname{cn}{u} = \cos{\phi}
     74
     75and
     76
     77.. math::
     78
     79    \operatorname{dn}{u} = \sqrt{1 - m\sin^2 \phi}.
     80
     81To emphasize the dependence on `m`, one can write
     82`\operatorname{sn}(u|m)` for example (and similarly for `\mathrm{cn}` and
     83`\mathrm{dn}`). This is the notation used below.
     84
     85For a given `k` with `0 < k < 1` they therefore are
     86solutions to the following nonlinear ordinary differential
     87equations:
     88
     89-  `\operatorname{sn}\,(x;k)` solves the differential equations
     90
     91  .. math::
     92
     93      \frac{d^2 y}{dx^2} + (1+k^2) y - 2 k^2 y^3 = 0,
     94
     95  and
     96
     97  .. math::
     98
     99     \left(\frac{dy}{dx}\right)^2 = (1-y^2) (1-k^2 y^2).
     100
     101-  `\operatorname{cn}(x;k)` solves the differential equations
     102
     103
     104  .. math::
     105
     106     \frac{d^2 y}{dx^2} + (1-2k^2) y + 2 k^2 y^3 = 0,
     107
     108  and `\left(\frac{dy}{dx}\right)^2 = (1-y^2)(1-k^2 + k^2 y^2)`.
     109
     110-  `\operatorname{dn}(x;k)` solves the differential equations
     111
     112  .. math::
     113
     114     \frac{d^2 y}{dx^2} - (2 - k^2) y + 2 y^3 = 0,
     115
     116  and `\left(\frac{dy}{dx}\right)^2= y^2 (1 - k^2 - y^2)`.
     117
     118  If `K(m)` denotes the complete elliptic integral of the
     119  first kind (named ``elliptic_kc`` in Sage), the elliptic functions
     120  `\operatorname{sn}(x|m)` and `\operatorname{cn}(x|m)` have real periods
     121  `4K(m)`, whereas `\operatorname{dn}(x|m)` has a period
     122  `2K(m)`. The limit `m\rightarrow 0` gives
     123  `K(0) = \pi/2` and trigonometric functions:
     124  `\operatorname{sn}(x|0) = \sin{x}`, `\operatorname{cn}(x|0) = \cos{x}`,
     125  `\operatorname{dn}(x|0) = 1`. The limit `m \rightarrow 1` gives
     126  `K(1) \rightarrow \infty` and hyperbolic functions:
     127  `\operatorname{sn}(x|1) = \tanh{x}`,
     128  `\operatorname{cn}(x|1) = \operatorname{sech}{x}`,
     129  `\operatorname{dn}(x|1) = \operatorname{sech}{x}`.
     130
     131REFERENCES:
     132
     133- http://en.wikipedia.org/wiki/Jacobi's_elliptic_functions
     134
     135- A. Khare, U. Sukhatme, "Cyclic Identities Involving
     136  Jacobi Elliptic Functions", Math ArXiv, math-ph/0201004
     137
     138AUTHORS:
     139
     140- David Joyner (2006): initial version
     141
     142- Eviatar Bach (2013): complete rewrite, new numerical evaluation, and addition
     143  of the Jacobi amplitude function
     144"""
     145
     146#*****************************************************************************
     147#       Copyright (C) 2006 David Joyner <wdj@usna.edu>
     148#       Copyright (C) 2013 Eviatar Bach <eviatarbach@gmail.com>
     149#
     150#  Distributed under the terms of the GNU General Public License (GPL)
     151#  as published by the Free Software Foundation; either version 2 of
     152#  the License, or (at your option) any later version.
     153#                  http://www.gnu.org/licenses/
     154#*****************************************************************************
     155
     156from sage.symbolic.function import BuiltinFunction, is_inexact
     157from sage.structure.coerce import parent
     158from sage.structure.element import get_coercion_model
     159from sage.symbolic.expression import Expression
     160from sage.functions.trig import (arctan, arcsin, arccos, arccot, arcsec,
     161                                 arccsc, csc, sec, sin, cos, tan, cot)
     162from sage.functions.hyperbolic import (arctanh, arccosh, arcsinh, arcsech,
     163                                       arccsch, arccoth, cosh, coth, sech,
     164                                       csch, tanh, sinh)
     165from sage.rings.rational_field import QQ
     166from sage.rings.integer import Integer
     167from sage.functions.special import elliptic_e, elliptic_kc
     168from sage.libs.mpmath import utils
     169from sage.misc.latex import latex
     170
     171HALF = QQ('1/2')
     172
     173
     174class Jacobi(BuiltinFunction):
     175    def __init__(self, kind):
     176        r"""
     177        Base class for the Jacobi elliptic functions
     178
     179        EXAMPLES::
     180
     181            sage: from sage.functions.jacobi import Jacobi
     182            sage: Jacobi('sn')
     183            jacobi_sn
     184        """
     185        if kind not in ['nd', 'ns', 'nc', 'dn', 'ds', 'dc', 'sn', 'sd',
     186                        'sc', 'cn', 'cd', 'cs']:
     187            raise ValueError("kind must be one of 'nd', 'ns', 'nc', 'dn', "
     188                             "'ds', 'dc', 'sn', 'sd', 'sc', 'cn', 'cd', 'cs'.")
     189        self.kind = kind
     190        BuiltinFunction.__init__(self,
     191                                 name='jacobi_{}'.format(kind),
     192                                 nargs=2, evalf_params_first=False,
     193                                 conversions=dict(maple=
     194                                                  ('Jacobi{}'
     195                                                   .format(kind.upper())),
     196                                                  mathematica=
     197                                                  ('Jacobi{}'
     198                                                   .format(kind.upper())),
     199                                                  maxima=
     200                                                  ('jacobi_{}'
     201                                                   .format(kind))))
     202
     203    def _eval_(self, x, m):
     204        r"""
     205        TESTS:
     206
     207        Check that the simplifications are correct::
     208
     209            sage: from mpmath import almosteq
     210            sage: almosteq(n(jacobi_nd(8, 0, hold=True)), n(jacobi_nd(8, 0)))
     211            True
     212            sage: almosteq(n(jacobi_nd(1, 1, hold=True)), n(jacobi_nd(1, 1)))
     213            True
     214            sage: almosteq(n(jacobi_nd(0, -5, hold=True)), n(jacobi_nd(0, -5)))
     215            True
     216            sage: almosteq(n(jacobi_ns(-4, 0, hold=True)), n(jacobi_ns(-4, 0)))
     217            True
     218            sage: almosteq(n(jacobi_ns(-2, 1, hold=True)), n(jacobi_ns(-2, 1)))
     219            True
     220            sage: almosteq(n(jacobi_nc(2, 0, hold=True)), n(jacobi_nc(2, 0)))
     221            True
     222            sage: almosteq(n(jacobi_nc(1, 1, hold=True)), n(jacobi_nc(1, 1)))
     223            True
     224            sage: almosteq(n(jacobi_nc(0, 0, hold=True)), n(jacobi_nc(0, 0)))
     225            True
     226            sage: almosteq(n(jacobi_dn(-10, 0, hold=True)), n(jacobi_dn(-10, 0)))
     227            True
     228            sage: almosteq(n(jacobi_dn(-1, 1, hold=True)), n(jacobi_dn(-1, 1)))
     229            True
     230            sage: almosteq(n(jacobi_dn(0, 3, hold=True)), n(jacobi_dn(0, 3)))
     231            True
     232            sage: almosteq(n(jacobi_ds(2, 0, hold=True)), n(jacobi_ds(2, 0)))
     233            True
     234            sage: almosteq(n(jacobi_dc(-1, 0, hold=True)), n(jacobi_dc(-1, 0)))
     235            True
     236            sage: almosteq(n(jacobi_dc(-8, 1, hold=True)), n(jacobi_dc(-8, 1)))
     237            True
     238            sage: almosteq(n(jacobi_dc(0, -10, hold=True)), n(jacobi_dc(0, -10)))
     239            True
     240            sage: almosteq(n(jacobi_sn(-7, 0, hold=True)), n(jacobi_sn(-7, 0)))
     241            True
     242            sage: almosteq(n(jacobi_sn(-3, 1, hold=True)), n(jacobi_sn(-3, 1)))
     243            True
     244            sage: almosteq(n(jacobi_sn(0, -6, hold=True)), n(jacobi_sn(0, -6)))
     245            True
     246            sage: almosteq(n(jacobi_sd(4, 0, hold=True)), n(jacobi_sd(4, 0)))
     247            True
     248            sage: almosteq(n(jacobi_sd(0, 1, hold=True)), n(jacobi_sd(0, 1)))
     249            True
     250            sage: almosteq(n(jacobi_sd(0, 3, hold=True)), n(jacobi_sd(0, 3)))
     251            True
     252            sage: almosteq(n(jacobi_sc(-9, 0, hold=True)), n(jacobi_sc(-9, 0)))
     253            True
     254            sage: almosteq(n(jacobi_sc(0, 1, hold=True)), n(jacobi_sc(0, 1)))
     255            True
     256            sage: almosteq(n(jacobi_sc(0, -10, hold=True)), n(jacobi_sc(0, -10)))
     257            True
     258            sage: almosteq(n(jacobi_cn(-2, 0, hold=True)), n(jacobi_cn(-2, 0)))
     259            True
     260            sage: almosteq(n(jacobi_cn(6, 1, hold=True)), n(jacobi_cn(6, 1)))
     261            True
     262            sage: almosteq(n(jacobi_cn(0, -10, hold=True)), n(jacobi_cn(0, -10)))
     263            True
     264            sage: almosteq(n(jacobi_cd(9, 0, hold=True)), n(jacobi_cd(9, 0)))
     265            True
     266            sage: almosteq(n(jacobi_cd(-8, 1, hold=True)), n(jacobi_cd(-8, 1)))
     267            True
     268            sage: almosteq(n(jacobi_cd(0, 1, hold=True)), n(jacobi_cd(0, 1)))
     269            True
     270            sage: almosteq(n(jacobi_cs(-9, 0, hold=True)), n(jacobi_cs(-9, 0)))
     271            True
     272            sage: almosteq(n(jacobi_cs(-6, 1, hold=True)), n(jacobi_cs(-6, 1)))
     273            True
     274        """
     275        coercion_model = get_coercion_model()
     276        co = coercion_model.canonical_coercion(x, m)[0]
     277        if is_inexact(co) and not isinstance(co, Expression):
     278            return self._evalf_(x, m, parent(co))
     279        if self.kind == 'nd':
     280            if m == 0:
     281                return Integer(1)
     282            elif m == 1:
     283                return cosh(x)
     284            elif x == 0:
     285                return Integer(1)
     286        elif self.kind == 'ns':
     287            if m == 0:
     288                return csc(x)
     289            elif m == 1:
     290                return coth(x)
     291        elif self.kind == 'nc':
     292            if m == 0:
     293                return sec(x)
     294            elif m == 1:
     295                return cosh(x)
     296            elif x == 0:
     297                return Integer(1)
     298        elif self.kind == 'dn':
     299            if m == 0:
     300                return Integer(1)
     301            elif m == 1:
     302                return sech(x)
     303            elif x == 0:
     304                return Integer(1)
     305        elif self.kind == 'ds':
     306            if m == 0:
     307                return csc(x)
     308        elif self.kind == 'dc':
     309            if m == 0:
     310                return sec(x)
     311            elif m == 1:
     312                return Integer(1)
     313            elif x == 0:
     314                return Integer(1)
     315        elif self.kind == 'sn':
     316            if m == 0:
     317                return sin(x)
     318            elif m == 1:
     319                return tanh(x)
     320            elif x == 0:
     321                return Integer(0)
     322        elif self.kind == 'sd':
     323            if m == 0:
     324                return sin(x)
     325            elif m == 1:
     326                return sinh(x)
     327            elif x == 0:
     328                return Integer(0)
     329        elif self.kind == 'sc':
     330            if m == 0:
     331                return tan(x)
     332            elif m == 1:
     333                return sinh(x)
     334            elif x == 0:
     335                return Integer(0)
     336        elif self.kind == 'cn':
     337            if m == 0:
     338                return cos(x)
     339            elif m == 1:
     340                return sech(x)
     341            elif x == 0:
     342                return Integer(1)
     343        elif self.kind == 'cd':
     344            if m == 0:
     345                return cos(x)
     346            elif m == 1:
     347                return Integer(1)
     348            elif x == 0:
     349                return Integer(1)
     350        elif self.kind == 'cs':
     351            if m == 0:
     352                return cot(x)
     353            elif m == 1:
     354                return csch(x)
     355        return
     356
     357    def _evalf_(self, x, m, parent):
     358        r"""
     359        TESTS::
     360
     361            sage: jacobi_sn(3, 4).n(100)
     362            -0.33260000892770027112809652714 + 1.7077912301715219199143891076e-33*I
     363            sage: jacobi_dn(I, I).n()
     364            0.874189950651018 + 0.667346865048825*I
     365            """
     366        from mpmath import ellipfun
     367        return utils.call(ellipfun, self.kind, x, m, parent=parent)
     368
     369    def _derivative_(self, x, m, diff_param):
     370        r"""
     371        TESTS:
     372
     373        sn, cn, and dn are analytic for all real x, so we can check that the
     374        derivatives are correct by computing the series::
     375
     376            sage: from mpmath import almosteq
     377            sage: a = 0.9327542442482303
     378            sage: b = 0.7402326293643771
     379            sage: almosteq(jacobi_sn(x, b).series(x, 10).subs(x=a),
     380            ....:          jacobi_sn(a, b), abs_eps=0.01)
     381            True
     382            sage: almosteq(jacobi_cn(x, b).series(x, 10).subs(x=a),
     383            ....:          jacobi_cn(a, b), abs_eps=0.01)
     384            True
     385            sage: almosteq(jacobi_dn(x, b).series(x, 10).subs(x=a),
     386            ....:          jacobi_dn(a, b), abs_eps=0.01)
     387            True
     388        """
     389        if diff_param == 0:
     390            # From Wolfram Functions Site
     391            if self.kind == 'cd':
     392                return (m - Integer(1)) * jacobi_nd(x, m) * jacobi_sd(x, m)
     393            elif self.kind == 'cn':
     394                return -jacobi_sn(x, m) * jacobi_dn(x, m)
     395            elif self.kind == 'cs':
     396                return -jacobi_ds(x, m) * jacobi_ns(x, m)
     397            elif self.kind == 'dc':
     398                return (Integer(1) - m) * jacobi_nc(x, m) * jacobi_sc(x, m)
     399            elif self.kind == 'dn':
     400                return -m * jacobi_sn(x, m) * jacobi_cn(x, m)
     401            elif self.kind == 'ds':
     402                return -jacobi_cs(x, m) * jacobi_ns(x, m)
     403            elif self.kind == 'nc':
     404                return jacobi_dc(x, m) * jacobi_sc(x, m)
     405            elif self.kind == 'nd':
     406                return m * jacobi_cd(x, m) * jacobi_sd(x, m)
     407            elif self.kind == 'ns':
     408                return -jacobi_cs(x, m) * jacobi_ds(x, m)
     409            elif self.kind == 'sc':
     410                return jacobi_dc(x, m) * jacobi_nc(x, m)
     411            elif self.kind == 'sd':
     412                return jacobi_cd(x, m) * jacobi_nd(x, m)
     413            elif self.kind == 'sn':
     414                return jacobi_cn(x, m) * jacobi_dn(x, m)
     415        elif diff_param == 1:
     416            # From Maxima
     417            if self.kind == 'nd':
     418                return (HALF*((x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     419                              (m - Integer(1)))*jacobi_sn(x, m)*jacobi_cn(x, m) -
     420                              jacobi_dn(x, m)*jacobi_sn(x, m)**Integer(2)/(m - Integer(1)))/
     421                        jacobi_dn(x, m)**Integer(2))
     422            elif self.kind == 'ns':
     423                return (HALF*(jacobi_sn(x, m)*jacobi_cn(x, m)**Integer(2)/(m - Integer(1)) -
     424                              (x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     425                               (m - Integer(1)))*jacobi_dn(x, m)*jacobi_cn(x, m)/m)/
     426                        jacobi_sn(x, m)**Integer(2))
     427            elif self.kind == 'nc':
     428                return (-HALF*(jacobi_sn(x, m)**Integer(2)*jacobi_cn(x, m)/(m - Integer(1)) -
     429                               (x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     430                                (m - Integer(1)))*jacobi_dn(x, m)*
     431                               jacobi_sn(x, m)/m)/jacobi_cn(x, m)**Integer(2))
     432            elif self.kind == 'dn':
     433                return (-HALF*(x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     434                               (m - Integer(1)))*jacobi_sn(x, m)*jacobi_cn(x, m) +
     435                        HALF*jacobi_dn(x, m)*jacobi_sn(x, m)**Integer(2)/(m - Integer(1)))
     436            elif self.kind == 'ds':
     437                return (HALF*(jacobi_sn(x, m)*jacobi_cn(x, m)**Integer(2)/(m - Integer(1)) -
     438                        (x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     439                         (m - Integer(1)))*jacobi_dn(x, m)*jacobi_cn(x, m)/m)*
     440                        jacobi_dn(x, m)/jacobi_sn(x, m)**Integer(2) -
     441                        HALF*((x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     442                              (m - Integer(1)))*jacobi_sn(x, m)*jacobi_cn(x, m) -
     443                              jacobi_dn(x, m)*jacobi_sn(x, m)**Integer(2)/(m - Integer(1)))/
     444                        jacobi_sn(x, m))
     445            elif self.kind == 'dc':
     446                return (-HALF*(jacobi_sn(x, m)**Integer(2)*jacobi_cn(x, m)/(m - Integer(1)) -
     447                               (x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     448                                (m - Integer(1)))*jacobi_dn(x, m)*
     449                               jacobi_sn(x, m)/m)*jacobi_dn(x, m)/
     450                        jacobi_cn(x, m)**Integer(2) -
     451                        HALF*((x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     452                               (m - Integer(1)))*jacobi_sn(x, m)*jacobi_cn(x, m) -
     453                              jacobi_dn(x, m)*jacobi_sn(x, m)**Integer(2)/(m - Integer(1)))/
     454                        jacobi_cn(x, m))
     455            elif self.kind == 'sn':
     456                return (-HALF*jacobi_sn(x, m)*jacobi_cn(x, m)**Integer(2)/(m - Integer(1)) +
     457                        HALF*(x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     458                              (m - Integer(1)))*jacobi_dn(x, m)*jacobi_cn(x, m)/m)
     459            elif self.kind == 'sd':
     460                return (-HALF*(jacobi_sn(x, m)*jacobi_cn(x, m)**Integer(2)/(m - Integer(1)) -
     461                        (x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     462                         (m - Integer(1)))*jacobi_dn(x, m)*jacobi_cn(x, m)/m)/
     463                        jacobi_dn(x, m) + HALF*
     464                        ((x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     465                          (m - Integer(1)))*jacobi_sn(x, m)*jacobi_cn(x, m) -
     466                         jacobi_dn(x, m)*jacobi_sn(x, m)**Integer(2)/(m - Integer(1)))*
     467                        jacobi_sn(x, m)/jacobi_dn(x, m)**Integer(2))
     468            elif self.kind == 'sc':
     469                return (-HALF*(jacobi_sn(x, m)*jacobi_cn(x, m)**Integer(2)/(m - Integer(1)) -
     470                               (x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     471                                (m - Integer(1)))*jacobi_dn(x, m)*
     472                               jacobi_cn(x, m)/m)/jacobi_cn(x, m) -
     473                        HALF*(jacobi_sn(x, m)**Integer(2)*jacobi_cn(x, m)/(m - Integer(1)) -
     474                              (x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     475                               (m - Integer(1)))*jacobi_dn(x, m)*jacobi_sn(x, m)/m)*
     476                        jacobi_sn(x, m)/jacobi_cn(x, m)**Integer(2))
     477            elif self.kind == 'cn':
     478                return (HALF*jacobi_sn(x, m)**Integer(2)*jacobi_cn(x, m)/(m - Integer(1)) -
     479                        HALF*(x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     480                              (m - Integer(1)))*jacobi_dn(x, m)*jacobi_sn(x, m)/m)
     481            elif self.kind == 'cd':
     482                return (HALF*(jacobi_sn(x, m)**Integer(2)*jacobi_cn(x, m)/(m - Integer(1)) -
     483                        (x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     484                         (m - Integer(1)))*jacobi_dn(x, m)*jacobi_sn(x, m)/m)/
     485                        jacobi_dn(x, m) +
     486                        HALF*((x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     487                               (m - Integer(1)))*jacobi_sn(x, m)*jacobi_cn(x, m) -
     488                              jacobi_dn(x, m)*jacobi_sn(x, m)**Integer(2)/(m - Integer(1)))*
     489                        jacobi_cn(x, m)/jacobi_dn(x, m)**Integer(2))
     490            elif self.kind == 'cs':
     491                return (HALF*(jacobi_sn(x, m)*jacobi_cn(x, m)**Integer(2)/(m - Integer(1)) -
     492                        (x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     493                         (m - Integer(1)))*jacobi_dn(x, m)*jacobi_cn(x, m)/m)*
     494                        jacobi_cn(x, m)/jacobi_sn(x, m)**Integer(2) +
     495                        HALF*(jacobi_sn(x, m)**Integer(2)*jacobi_cn(x, m)/(m - Integer(1)) -
     496                              (x + elliptic_e(arcsin(jacobi_sn(x, m)), m)/
     497                               (m - Integer(1)))*jacobi_dn(x, m)*jacobi_sn(x, m)/m)/
     498                        jacobi_sn(x, m))
     499
     500    def _latex_(self):
     501        r"""
     502        TESTS::
     503
     504            sage: latex(jacobi_sn)
     505            \operatorname{sn}
     506        """
     507        return r"\operatorname{{{}}}".format(self.kind)
     508
     509    def _print_latex_(self, x, m):
     510        r"""
     511        TESTS::
     512
     513            sage: latex(jacobi_sn(x, 3))
     514            \operatorname{sn}\left(x\middle|3\right)
     515        """
     516        return r"\operatorname{{{}}}\left({}\middle|{}\right)".format(self.kind,
     517                                                                      latex(x),
     518                                                                      latex(m))
     519
     520jacobi_nd = Jacobi('nd')
     521jacobi_ns = Jacobi('ns')
     522jacobi_nc = Jacobi('nc')
     523jacobi_dn = Jacobi('dn')
     524jacobi_ds = Jacobi('ds')
     525jacobi_dc = Jacobi('dc')
     526jacobi_sn = Jacobi('sn')
     527jacobi_sd = Jacobi('sd')
     528jacobi_sc = Jacobi('sc')
     529jacobi_cn = Jacobi('cn')
     530jacobi_cd = Jacobi('cd')
     531jacobi_cs = Jacobi('cs')
     532
     533
     534class InverseJacobi(BuiltinFunction):
     535    def __init__(self, kind):
     536        r"""
     537        Base class for the inverse Jacobi elliptic functions
     538
     539        EXAMPLES::
     540
     541            sage: from sage.functions.jacobi import InverseJacobi
     542            sage: InverseJacobi('sn')
     543            inverse_jacobi_sn
     544        """
     545        if kind not in ['nd', 'ns', 'nc', 'dn', 'ds', 'dc', 'sn', 'sd',
     546                        'sc', 'cn', 'cd', 'cs']:
     547            raise ValueError("kind must be one of 'nd', 'ns', 'nc', 'dn', "
     548                             "'ds', 'dc', 'sn', 'sd', 'sc', 'cn', 'cd', 'cs'.")
     549        self.kind = kind
     550        BuiltinFunction.__init__(self,
     551                                 name='inverse_jacobi_{}'.format(kind),
     552                                 nargs=2, evalf_params_first=False,
     553                                 conversions=dict(maple=
     554                                                  ('InverseJacobi{}'
     555                                                   .format(kind.upper())),
     556                                                  mathematica=
     557                                                  ('InverseJacobi{}'
     558                                                   .format(kind.upper())),
     559                                                  maxima=
     560                                                  ('inverse_jacobi_{}'
     561                                                   .format(kind))))
     562
     563    def _eval_(self, x, m):
     564        r"""
     565        TESTS:
     566
     567        Check that the simplifications are correct::
     568
     569            sage: from mpmath import almosteq
     570            sage: almosteq(n(inverse_jacobi_cd(1, -8, hold=True)),
     571            ....:          n(inverse_jacobi_cd(1, -8)))
     572            True
     573            sage: almosteq(n(inverse_jacobi_cn(0, -5, hold=True)),
     574            ....:          n(inverse_jacobi_cn(0, -5)))
     575            True
     576            sage: almosteq(n(inverse_jacobi_cn(1, -8, hold=True)),
     577            ....:          n(inverse_jacobi_cn(1, -8)))
     578            True
     579            sage: almosteq(n(inverse_jacobi_cs(7, 1, hold=True)),
     580            ....:          n(inverse_jacobi_cs(7, 1)))
     581            True
     582            sage: almosteq(n(inverse_jacobi_dc(3, 0, hold=True)),
     583            ....:          n(inverse_jacobi_dc(3, 0)))
     584            True
     585            sage: almosteq(n(inverse_jacobi_dc(1, 7, hold=True)),
     586            ....:          n(inverse_jacobi_dc(1, 7)))
     587            True
     588            sage: almosteq(n(inverse_jacobi_dn(1, -1, hold=True)),
     589            ....:          n(inverse_jacobi_dn(1, -1)))
     590            True
     591            sage: almosteq(n(inverse_jacobi_ds(7, 0, hold=True)),
     592            ....:          n(inverse_jacobi_ds(7, 0)))
     593            True
     594            sage: almosteq(n(inverse_jacobi_ds(5, 1, hold=True)),
     595            ....:          n(inverse_jacobi_ds(5, 1)))
     596            True
     597            sage: almosteq(n(inverse_jacobi_nc(-2, 0, hold=True)),
     598            ....:          n(inverse_jacobi_nc(-2, 0)))
     599            True
     600            sage: almosteq(n(inverse_jacobi_nc(-1, 1, hold=True)),
     601            ....:          n(inverse_jacobi_nc(-1, 1)))
     602            True
     603            sage: almosteq(n(inverse_jacobi_nc(1, 4, hold=True)),
     604            ....:          n(inverse_jacobi_nc(1, 4)))
     605            True
     606            sage: almosteq(n(inverse_jacobi_nd(9, 1, hold=True)),
     607            ....:          n(inverse_jacobi_nd(9, 1)))
     608            True
     609            sage: almosteq(n(inverse_jacobi_nd(1, -9, hold=True)),
     610            ....:          n(inverse_jacobi_nd(1, -9)))
     611            True
     612            sage: almosteq(n(inverse_jacobi_ns(-6, 0, hold=True)),
     613            ....:          n(inverse_jacobi_ns(-6, 0)))
     614            True
     615            sage: almosteq(n(inverse_jacobi_ns(6, 1, hold=True)),
     616            ....:          n(inverse_jacobi_ns(6, 1)))
     617            True
     618            sage: almosteq(n(inverse_jacobi_sc(9, 0, hold=True)),
     619            ....:          n(inverse_jacobi_sc(9, 0)))
     620            True
     621            sage: almosteq(n(inverse_jacobi_sc(8, 1, hold=True)),
     622            ....:          n(inverse_jacobi_sc(8, 1)))
     623            True
     624            sage: almosteq(n(inverse_jacobi_sc(0, -8, hold=True)),
     625            ....:          n(inverse_jacobi_sc(0, -8)))
     626            True
     627            sage: almosteq(n(inverse_jacobi_sd(-1, 0, hold=True)),
     628            ....:          n(inverse_jacobi_sd(-1, 0)))
     629            True
     630            sage: almosteq(n(inverse_jacobi_sd(-2, 1, hold=True)),
     631            ....:          n(inverse_jacobi_sd(-2, 1)))
     632            True
     633            sage: almosteq(n(inverse_jacobi_sd(0, -2, hold=True)),
     634            ....:          n(inverse_jacobi_sd(0, -2)))
     635            True
     636            sage: almosteq(n(inverse_jacobi_sn(0, 0, hold=True)),
     637            ....:          n(inverse_jacobi_sn(0, 0)))
     638            True
     639            sage: almosteq(n(inverse_jacobi_sn(0, 6, hold=True)),
     640            ....:          n(inverse_jacobi_sn(0, 6)))
     641            True
     642        """
     643        coercion_model = get_coercion_model()
     644        x, m = coercion_model.canonical_coercion(x, m)
     645        if is_inexact(x) and not isinstance(x, Expression):
     646            return self._evalf_(x, m, parent(x))
     647        if self.kind == 'cd':
     648            if m == 0:
     649                return arccos(x)
     650            elif x == 1:
     651                return Integer(0)
     652        elif self.kind == 'cn':
     653            if m == 0:
     654                return arccos(x)
     655            elif m == 1:
     656                return arcsech(x)
     657            elif x == 0:
     658                return elliptic_kc(m)
     659            elif x == 1:
     660                return Integer(0)
     661        elif self.kind == 'cs':
     662            if m == 0:
     663                return arccot(x)
     664            elif m == 1:
     665                return arccsch(x)
     666        elif self.kind == 'dc':
     667            if m == 0:
     668                return arcsec(x)
     669            elif x == 1:
     670                return Integer(0)
     671        elif self.kind == 'dn':
     672            if m == 1:
     673                return arcsech(x)
     674            elif x == 1:
     675                return Integer(0)
     676        elif self.kind == 'ds':
     677            if m == 0:
     678                return arccsc(x)
     679            elif m == 1:
     680                return arccsch(x)
     681        elif self.kind == 'nc':
     682            if m == 0:
     683                return arcsec(x)
     684            elif m == 1:
     685                return arccosh(x)
     686            elif x == 1:
     687                return Integer(0)
     688        elif self.kind == 'nd':
     689            if m == 1:
     690                return arccosh(x)
     691            elif x == 1:
     692                return Integer(0)
     693        elif self.kind == 'ns':
     694            if m == 0:
     695                return arccsc(x)
     696            elif m == 1:
     697                return arccoth(x)
     698        elif self.kind == 'sc':
     699            if m == 0:
     700                return arctan(x)
     701            elif m == 1:
     702                return arcsinh(x)
     703            elif x == 0:
     704                return Integer(0)
     705        elif self.kind == 'sd':
     706            if m == 0:
     707                return arcsin(x)
     708            elif m == 1:
     709                return arcsinh(x)
     710            elif x == 0:
     711                return Integer(0)
     712        elif self.kind == 'sn':
     713            if m == 0:
     714                return arcsin(x)
     715            elif m == 1:
     716                return arctanh(x)
     717            elif x == 0:
     718                return Integer(0)
     719        return
     720
     721    def _evalf_(self, x, m, parent):
     722        r"""
     723        TESTS::
     724
     725            sage: inverse_jacobi_cn(2, 3).n()
     726            0.859663746362987*I
     727            sage: inverse_jacobi_cd(3, 4).n(100)
     728            -0.67214752201235862490069823239 + 2.1565156474996432354386749988*I
     729        """
     730        return utils.call(inverse_jacobi_f, self.kind, x, m,
     731                          parent=parent)
     732
     733    def _derivative_(self, x, m, diff_param):
     734        r"""
     735        TESTS:
     736
     737        Check that dy/dx * dx/dy == 1, where y = jacobi_pq(x, m) and
     738        x = inverse_jacobi_pq(y, m)::
     739
     740            sage: from mpmath import almosteq
     741            sage: a = 0.130103220857094
     742            sage: b = 0.437176765041986
     743            sage: m = var('m')
     744            sage: almosteq(abs((diff(jacobi_cd(x, m), x) *
     745            ....:               diff(inverse_jacobi_cd(x, m), x).subs(x=jacobi_cd(x, m))).subs(x=a, m=b)),
     746            ....:          1, abs_eps=1e-14)
     747            True
     748            sage: almosteq(abs((diff(jacobi_cn(x, m), x) *
     749            ....:               diff(inverse_jacobi_cn(x, m), x).subs(x=jacobi_cn(x, m))).subs(x=a, m=b)),
     750            ....:          1, abs_eps=1e-14)
     751            True
     752            sage: almosteq(abs((diff(jacobi_cs(x, m), x) *
     753            ....:               diff(inverse_jacobi_cs(x, m), x).subs(x=jacobi_cs(x, m))).subs(x=a, m=b)),
     754            ....:          1, abs_eps=1e-14)
     755            True
     756            sage: almosteq(abs((diff(jacobi_dc(x, m), x) *
     757            ....:               diff(inverse_jacobi_dc(x, m), x).subs(x=jacobi_dc(x, m))).subs(x=a, m=b)),
     758            ....:          1, abs_eps=1e-14)
     759            True
     760            sage: almosteq(abs((diff(jacobi_dn(x, m), x) *
     761            ....:               diff(inverse_jacobi_dn(x, m), x).subs(x=jacobi_dn(x, m))).subs(x=a, m=b)),
     762            ....:          1, abs_eps=1e-14)
     763            True
     764            sage: almosteq(abs((diff(jacobi_ds(x, m), x) *
     765            ....:               diff(inverse_jacobi_ds(x, m), x).subs(x=jacobi_ds(x, m))).subs(x=a, m=b)),
     766            ....:          1, abs_eps=1e-14)
     767            True
     768            sage: almosteq(abs((diff(jacobi_nc(x, m), x) *
     769            ....:               diff(inverse_jacobi_nc(x, m), x).subs(x=jacobi_nc(x, m))).subs(x=a, m=b)),
     770            ....:          1, abs_eps=1e-14)
     771            True
     772            sage: almosteq(abs((diff(jacobi_nd(x, m), x) *
     773            ....:               diff(inverse_jacobi_nd(x, m), x).subs(x=jacobi_nd(x, m))).subs(x=a, m=b)),
     774            ....:          1, abs_eps=1e-14)
     775            True
     776            sage: almosteq(abs((diff(jacobi_ns(x, m), x) *
     777            ....:               diff(inverse_jacobi_ns(x, m), x).subs(x=jacobi_ns(x, m))).subs(x=a, m=b)),
     778            ....:          1, abs_eps=1e-14)
     779            True
     780            sage: almosteq(abs((diff(jacobi_sc(x, m), x) *
     781            ....:               diff(inverse_jacobi_sc(x, m), x).subs(x=jacobi_sc(x, m))).subs(x=a, m=b)),
     782            ....:          1, abs_eps=1e-14)
     783            True
     784            sage: almosteq(abs((diff(jacobi_sd(x, m), x) *
     785            ....:               diff(inverse_jacobi_sd(x, m), x).subs(x=jacobi_sd(x, m))).subs(x=a, m=b)),
     786            ....:          1, abs_eps=1e-14)
     787            True
     788            sage: almosteq(abs((diff(jacobi_sn(x, m), x) *
     789            ....:               diff(inverse_jacobi_sn(x, m), x).subs(x=jacobi_sn(x, m))).subs(x=a, m=b)),
     790            ....:          1, abs_eps=1e-14)
     791            True
     792        """
     793        # From Wolfram Functions Site
     794        if diff_param == 0:
     795            if self.kind == 'cd':
     796                return (jacobi_sn(inverse_jacobi_cd(x, m), m) /
     797                        (x ** Integer(2) - Integer(1)))
     798            elif self.kind == 'cn':
     799                return (jacobi_ds(inverse_jacobi_cn(x, m), m) /
     800                        (m * x ** Integer(2) - m + Integer(1)))
     801            elif self.kind == 'cs':
     802                return (jacobi_nd(inverse_jacobi_cs(x, m), m) /
     803                        (x ** Integer(2) + Integer(1)))
     804            elif self.kind == 'dc':
     805                return (jacobi_sn(inverse_jacobi_dc(x, m), m) /
     806                        (x ** Integer(2) - Integer(1)))
     807            elif self.kind == 'dn':
     808                return -(jacobi_cs(inverse_jacobi_dn(x, m), m) /
     809                         (x ** Integer(2) + m - Integer(1)))
     810            elif self.kind == 'ds':
     811                return (jacobi_nc(inverse_jacobi_ds(x, m), m) /
     812                        (x ** Integer(2) + m))
     813            elif self.kind == 'nc':
     814                return (jacobi_ds(inverse_jacobi_nc(x, m), m) /
     815                        (-m * x ** Integer(2) + x ** Integer(2) + m))
     816            elif self.kind == 'nd':
     817                return (jacobi_sc(inverse_jacobi_nd(x, m), m) /
     818                        (x ** Integer(2) - Integer(1)))
     819            elif self.kind == 'ns':
     820                return Integer(1) / (jacobi_cs(inverse_jacobi_ns(x, m), m) *
     821                                     jacobi_ds(inverse_jacobi_ns(x, m), m))
     822            elif self.kind == 'sc':
     823                return (jacobi_nd(inverse_jacobi_sc(x, m), m) /
     824                        (x ** Integer(2) + Integer(1)))
     825            elif self.kind == 'sd':
     826                return (jacobi_cn(inverse_jacobi_sd(x, m), m) /
     827                        ((m - Integer(1)) * x ** Integer(2) + Integer(1)))
     828            elif self.kind == 'sn':
     829                return (jacobi_cd(inverse_jacobi_sn(x, m), m) /
     830                        (Integer(1) - x ** Integer(2)))
     831        elif diff_param == 1:
     832            if self.kind == 'cd':
     833                return ((Integer(1) / (Integer(2) * (Integer(1) - m) * m)) *
     834                        ((m - Integer(1)) * inverse_jacobi_cd(x, m) +
     835                         elliptic_e(jacobi_am(inverse_jacobi_cd(x, m), m),
     836                                    m)))
     837            elif self.kind == 'cn':
     838                return ((-(Integer(1) / (Integer(2) * (-Integer(1) + m) * m))) *
     839                        (elliptic_e(jacobi_am(inverse_jacobi_cn(x, m), m),
     840                                    m) + (-Integer(1) + m) *
     841                         inverse_jacobi_cn(x, m) - m * x *
     842                         jacobi_sd(inverse_jacobi_cn(x, m), m)))
     843            elif self.kind == 'cs':
     844                return ((-(Integer(1) / (Integer(2) * (-Integer(1) + m) * m * (Integer(1) + x ** Integer(2))))) *
     845                        ((Integer(1) + x ** Integer(2)) *
     846                         elliptic_e(jacobi_am(inverse_jacobi_cs(x, m), m),
     847                                    m) + (-Integer(1) + m) * (Integer(1) + x ** Integer(2)) *
     848                         inverse_jacobi_cs(x, m) - m * x *
     849                         jacobi_nd(inverse_jacobi_cs(x, m), m)))
     850            elif self.kind == 'dc':
     851                return ((Integer(1) / (Integer(2) * (Integer(1) - m) * m)) *
     852                        (elliptic_e(jacobi_am(inverse_jacobi_dc(x, m), m),
     853                                    m) - (Integer(1) - m) *
     854                         inverse_jacobi_dc(x, m)))
     855            elif self.kind == 'dn':
     856                return ((Integer(1) / (Integer(2) * (Integer(1) - m) * m)) * ((m - Integer(1)) *
     857                        inverse_jacobi_dn(x, m) +
     858                        elliptic_e(jacobi_am(inverse_jacobi_dn(x, m), m), m) -
     859                        x * jacobi_sc(inverse_jacobi_dn(x, m), m)))
     860            elif self.kind == 'ds':
     861                return ((-(Integer(1) / (Integer(2) * (-Integer(1) + m) * m))) *
     862                        (elliptic_e(jacobi_am(inverse_jacobi_ds(x, m), m), m) +
     863                         (-Integer(1) + m) * inverse_jacobi_ds(x, m) -
     864                         (m * x * jacobi_nc(inverse_jacobi_ds(x, m), m)) /
     865                         (m + x ** Integer(2))))
     866            elif self.kind == 'nc':
     867                return ((Integer(1) / (Integer(2) * (-Integer(1) + m) * m * x)) * ((-x) *
     868                        (elliptic_e(jacobi_am(inverse_jacobi_nc(x, m), m), m) +
     869                         (-Integer(1) + m) * inverse_jacobi_nc(x, m)) + m *
     870                        jacobi_sd(inverse_jacobi_nc(x, m), m)))
     871            elif self.kind == 'nd':
     872                return ((Integer(1) / (Integer(2) * (m - Integer(1)) * m)) *
     873                        ((Integer(1) - m) * inverse_jacobi_nd(x, m) -
     874                         elliptic_e(jacobi_am(inverse_jacobi_nd(x, m), m), m) +
     875                         (Integer(1) / x) * jacobi_sc(inverse_jacobi_nd(x, m), m)))
     876            elif self.kind == 'ns':
     877                return ((Integer(1)/(Integer(2) * (m - Integer(1)) * m)) *
     878                        ((Integer(1) - m) * inverse_jacobi_ns(x, m) -
     879                         elliptic_e(jacobi_am(inverse_jacobi_ns(x, m), m), m) +
     880                         (m / x) * jacobi_cd(inverse_jacobi_ns(x, m), m)))
     881            elif self.kind == 'sc':
     882                return ((-(Integer(1) / (Integer(2) * (-Integer(1) + m) * m * (Integer(1) + x ** Integer(2))))) *
     883                        ((Integer(1) + x ** Integer(2)) *
     884                         elliptic_e(jacobi_am(inverse_jacobi_sc(x, m), m), m) +
     885                         (-Integer(1) + m) * (Integer(1) + x ** Integer(2)) * inverse_jacobi_sc(x, m) -
     886                         m * x * jacobi_nd(inverse_jacobi_sc(x, m), m)))
     887            elif self.kind == 'sd':
     888                return ((-(Integer(1) / (Integer(2) * (-Integer(1) + m) * m))) *
     889                        (elliptic_e(jacobi_am(inverse_jacobi_sd(x, m), m), m) +
     890                         (-Integer(1) + m) * inverse_jacobi_sd(x, m) -
     891                         (m * x * jacobi_nc(inverse_jacobi_sd(x, m), m)) /
     892                         (Integer(1) + m * x ** Integer(2))))
     893            elif self.kind == 'sn':
     894                return ((Integer(1) / (Integer(2) * (Integer(1) - m) * m)) *
     895                        (elliptic_e(jacobi_am(inverse_jacobi_sn(x, m), m), m) +
     896                         (-Integer(1) + m) * inverse_jacobi_sn(x, m) - m * x *
     897                         jacobi_cd(inverse_jacobi_sn(x, m), m)))
     898
     899    def _latex_(self):
     900        r"""
     901        TESTS::
     902
     903            sage: latex(inverse_jacobi_dn)
     904            \operatorname{arcdn}
     905        """
     906        return r"\operatorname{{arc{}}}".format(self.kind)
     907
     908    def _print_latex_(self, x, m):
     909        r"""
     910        TESTS::
     911
     912            sage: latex(inverse_jacobi_dn(x, 3))
     913            \operatorname{arcdn}\left(x\middle|3\right)
     914        """
     915        return r"\operatorname{{arc{}}}\left({}\middle|{}\right)".format(self.kind,
     916                                                                         latex(x),
     917                                                                         latex(m))
     918
     919inverse_jacobi_nd = InverseJacobi('nd')
     920inverse_jacobi_ns = InverseJacobi('ns')
     921inverse_jacobi_nc = InverseJacobi('nc')
     922inverse_jacobi_dn = InverseJacobi('dn')
     923inverse_jacobi_ds = InverseJacobi('ds')
     924inverse_jacobi_dc = InverseJacobi('dc')
     925inverse_jacobi_sn = InverseJacobi('sn')
     926inverse_jacobi_sd = InverseJacobi('sd')
     927inverse_jacobi_sc = InverseJacobi('sc')
     928inverse_jacobi_cn = InverseJacobi('cn')
     929inverse_jacobi_cd = InverseJacobi('cd')
     930inverse_jacobi_cs = InverseJacobi('cs')
     931
     932
     933def jacobi(kind, z, m, **kwargs):
     934    r"""
     935    The 12 Jacobi elliptic functions
     936
     937    INPUT:
     938
     939    - ``kind`` -- a string of the form ``'pq'``, where ``p``, ``q`` are in
     940      ``c``, ``d``, ``n``, ``s``
     941    - ``z`` -- a complex number
     942    - ``m`` -- a complex number. Note that `m=k^2`, where `k` is the elliptic
     943      modulus.
     944
     945    EXAMPLES::
     946
     947        sage: jacobi('sn', 1, 1)
     948        tanh(1)
     949        sage: jacobi('cd', 1, 1/2)
     950        jacobi_cd(1, 1/2)
     951        sage: RDF(jacobi('cd', 1, 1/2))
     952        0.724009721659
     953        sage: (RDF(jacobi('cn', 1, 1/2)), RDF(jacobi('dn', 1, 1/2)),
     954        ....:  RDF(jacobi('cn', 1, 1/2) / jacobi('dn', 1, 1/2)))
     955        (0.595976567672, 0.823161001632, 0.724009721659)
     956        sage: jsn = jacobi('sn', x, 1)
     957        sage: P = plot(jsn, 0, 1)
     958    """
     959    if kind == 'nd':
     960        return jacobi_nd(z, m, **kwargs)
     961    elif kind == 'ns':
     962        return jacobi_ns(z, m, **kwargs)
     963    elif kind == 'nc':
     964        return jacobi_nc(z, m, **kwargs)
     965    elif kind == 'dn':
     966        return jacobi_dn(z, m, **kwargs)
     967    elif kind == 'ds':
     968        return jacobi_ds(z, m, **kwargs)
     969    elif kind == 'dc':
     970        return jacobi_dc(z, m, **kwargs)
     971    elif kind == 'sn':
     972        return jacobi_sn(z, m, **kwargs)
     973    elif kind == 'sd':
     974        return jacobi_sd(z, m, **kwargs)
     975    elif kind == 'sc':
     976        return jacobi_sc(z, m, **kwargs)
     977    elif kind == 'cn':
     978        return jacobi_cn(z, m, **kwargs)
     979    elif kind == 'cd':
     980        return jacobi_cd(z, m, **kwargs)
     981    elif kind == 'cs':
     982        return jacobi_cs(z, m, **kwargs)
     983    else:
     984        raise ValueError("kind must be one of 'nd', 'ns', 'nc', 'dn', "
     985                         "'ds', 'dc', 'sn', 'sd', 'sc', 'cn', 'cd', 'cs'.")
     986
     987def inverse_jacobi(kind, x, m, **kwargs):
     988    """
     989    The inverses of the 12 Jacobi elliptic functions. They have the property
     990    that
     991
     992    .. math::
     993        \operatorname{pq}(\operatorname{arcpq}(x|m)|m) =
     994        \operatorname{pq}(\operatorname{pq}^{-1}(x|m)|m) = x.
     995
     996    INPUT:
     997
     998    - ``kind`` -- a string of the form ``'pq'``, where ``p``, ``q`` are in
     999      ``c``, ``d``, ``n``, ``s``
     1000    - ``x`` -- a real number
     1001    - ``m`` -- a real number. Note that `m=k^2`, where `k` is the elliptic
     1002      modulus.
     1003
     1004    EXAMPLES::
     1005
     1006        sage: jacobi('dn', inverse_jacobi('dn', 3, 0.4), 0.4)
     1007        3.00000000000000
     1008        sage: inverse_jacobi('dn', 10, 1/10).n(digits=50)
     1009        2.4777736267904273296523691232988240759001423661683*I
     1010        sage: inverse_jacobi_dn(x, 1)
     1011        arcsech(x)
     1012        sage: inverse_jacobi_dn(1, 3)
     1013        0
     1014        sage: m = var('m')
     1015        sage: z = inverse_jacobi_dn(x, m).series(x, 4).subs(x=0.1, m=0.7)
     1016        sage: jacobi_dn(z, 0.7)
     1017        0.0999892750039819...
     1018        sage: inverse_jacobi_nd(x, 1)
     1019        arccosh(x)
     1020        sage: inverse_jacobi_nd(1, 2)
     1021        0
     1022        sage: inverse_jacobi_ns(10^-5, 3).n()
     1023        5.77350269202456e-6 + 1.17142008414677*I
     1024        sage: jacobi('sn', 1/2, 1/2)
     1025        jacobi_sn(1/2, 1/2)
     1026        sage: jacobi('sn', 1/2, 1/2).n()
     1027        0.470750473655657
     1028        sage: inverse_jacobi('sn', 0.47, 1/2)
     1029        0.499098231322220
     1030        sage: inverse_jacobi('sn', 0.4707504, 0.5)
     1031        0.499999911466555
     1032        sage: P = plot(inverse_jacobi('sn', x, 0.5), 0, 1)
     1033    """
     1034    if kind == 'nd':
     1035        return inverse_jacobi_nd(x, m, **kwargs)
     1036    elif kind == 'ns':
     1037        return inverse_jacobi_ns(x, m, **kwargs)
     1038    elif kind == 'nc':
     1039        return inverse_jacobi_nc(x, m, **kwargs)
     1040    elif kind == 'dn':
     1041        return inverse_jacobi_dn(x, m, **kwargs)
     1042    elif kind == 'ds':
     1043        return inverse_jacobi_ds(x, m, **kwargs)
     1044    elif kind == 'dc':
     1045        return inverse_jacobi_dc(x, m, **kwargs)
     1046    elif kind == 'sn':
     1047        return inverse_jacobi_sn(x, m, **kwargs)
     1048    elif kind == 'sd':
     1049        return inverse_jacobi_sd(x, m, **kwargs)
     1050    elif kind == 'sc':
     1051        return inverse_jacobi_sc(x, m, **kwargs)
     1052    elif kind == 'cn':
     1053        return inverse_jacobi_cn(x, m, **kwargs)
     1054    elif kind == 'cd':
     1055        return inverse_jacobi_cd(x, m, **kwargs)
     1056    elif kind == 'cs':
     1057        return inverse_jacobi_cs(x, m, **kwargs)
     1058    else:
     1059        raise ValueError("kind must be one of 'nd', 'ns', 'nc', 'dn', "
     1060                         "'ds', 'dc', 'sn', 'sd', 'sc', 'cn', 'cd', 'cs'.")
     1061
     1062class JacobiAmplitude(BuiltinFunction):
     1063    r"""
     1064    The Jacobi amplitude function
     1065    `\operatorname{am}(x|m)=\int_0^x \operatorname{dn}(t|m) dt`. For
     1066    `-K(m) \leq x \leq K(m)`, `F(\operatorname{am}(x|m)|m)=x.`
     1067    """
     1068    def __init__(self):
     1069        r"""
     1070        TESTS::
     1071
     1072            sage: from sage.functions.jacobi import JacobiAmplitude
     1073            sage: JacobiAmplitude()
     1074            jacobi_am
     1075        """
     1076        BuiltinFunction.__init__(self, name='jacobi_am', nargs=2,
     1077                                 conversions=dict(maple='JacobiAM',
     1078                                                  mathematica=
     1079                                                  'JacobiAmplitude'),
     1080                                 evalf_params_first=False)
     1081
     1082    def _eval_(self, x, m):
     1083        r"""
     1084        TESTS::
     1085
     1086            sage: jacobi_am(x, 0)
     1087            x
     1088            sage: jacobi_am(0, x)
     1089            0
     1090            sage: jacobi_am(3, 4.)
     1091            -0.339059208303591
     1092        """
     1093        coercion_model = get_coercion_model()
     1094        x, m = coercion_model.canonical_coercion(x, m)
     1095        if is_inexact(x) and not isinstance(x, Expression):
     1096            return self._evalf_(x, m, parent(x))
     1097        if m == 0:
     1098            return x
     1099        elif x == 0:
     1100            return Integer(0)
     1101        return
     1102
     1103    def _evalf_(self, x, m, parent):
     1104        r"""
     1105        TESTS::
     1106
     1107            sage: jacobi_am(1, 2).n(100)
     1108            0.73704379494724574105101929735
     1109        """
     1110        return utils.call(jacobi_am_f, x, m, parent=parent)
     1111
     1112    def _derivative_(self, x, m, diff_param):
     1113        r"""
     1114        TESTS::
     1115
     1116            sage: diff(jacobi_am(x, 3), x)
     1117            jacobi_dn(x, 3)
     1118            sage: diff(jacobi_am(3, x), x)
     1119            -1/2*(x*jacobi_cn(3, x)*jacobi_sn(3, x) -...
     1120            (3*x + elliptic_e(jacobi_am(3, x), x) - 3)*jacobi_dn(3, x))/((x - 1)*x)
     1121        """
     1122        if diff_param == 0:
     1123            return jacobi_dn(x, m)
     1124        elif diff_param == 1:
     1125            return (((Integer(-1) + m) * x + elliptic_e(jacobi_am(x, m), m)) *
     1126                    jacobi('dn', x, m) - m * jacobi('cn', x, m) *
     1127                    jacobi('sn', x, m)) / (Integer(2) * (Integer(-1) + m) * m)
     1128
     1129    def _latex_(self):
     1130        r"""
     1131        TESTS::
     1132
     1133            sage: latex(jacobi_am)
     1134            \operatorname{am}
     1135        """
     1136        return r"\operatorname{am}"
     1137
     1138    def _print_latex_(self, x, m):
     1139        r"""
     1140        TESTS::
     1141
     1142            sage: latex(jacobi_am(3,x))
     1143            \operatorname{am}\left(3\middle|x\right)
     1144        """
     1145        return r"\operatorname{{am}}\left({}\middle|{}\right)".format(latex(x),
     1146                                                                      latex(m))
     1147
     1148jacobi_am = JacobiAmplitude()
     1149
     1150
     1151def inverse_jacobi_f(kind, x, m):
     1152    """
     1153    Internal function for numerical evaluation of a continous complex branch
     1154    of each inverse Jacobi function, as described in [Tee]_. Only accepts real
     1155    arguments.
     1156
     1157    REFERENCES:
     1158
     1159    .. [Tee] Tee, Garry J. "Continuous branches of inverses of the 12 Jacobi
     1160             elliptic functions for real argument". 1997.
     1161             https://researchspace.auckland.ac.nz/bitstream/handle/2292/5042/\
     1162390.pdf.
     1163
     1164    TESTS::
     1165
     1166        sage: from mpmath import ellipfun, chop
     1167        sage: from sage.functions.jacobi import inverse_jacobi_f
     1168
     1169        sage: chop(ellipfun('sn', inverse_jacobi_f('sn', 0.6, 0), 0))
     1170        mpf('0.59999999999999998')
     1171        sage: chop(ellipfun('sn', inverse_jacobi_f('sn', 0.6, 1), 1))
     1172        mpf('0.59999999999999998')
     1173        sage: chop(ellipfun('sn', inverse_jacobi_f('sn', 0, -3), -3))
     1174        mpf('0.0')
     1175        sage: chop(ellipfun('sn', inverse_jacobi_f('sn', -1, 4), 4))
     1176        mpf('-1.0')
     1177        sage: chop(ellipfun('sn', inverse_jacobi_f('sn', 0.3, 4), 4))
     1178        mpf('0.29999999999999999')
     1179        sage: chop(ellipfun('sn', inverse_jacobi_f('sn', 0.8, 4), 4))
     1180        mpf('0.80000000000000004')
     1181
     1182        sage: chop(ellipfun('ns', inverse_jacobi_f('ns', 0.8, 0), 0))
     1183        mpf('0.80000000000000004')
     1184        sage: chop(ellipfun('ns', inverse_jacobi_f('ns', -0.7, 1), 1))
     1185        mpf('-0.69999999999999996')
     1186        sage: chop(ellipfun('ns', inverse_jacobi_f('ns', 0.01, 2), 2))
     1187        mpf('0.01')
     1188        sage: chop(ellipfun('ns', inverse_jacobi_f('ns', 0, 2), 2))
     1189        mpf('0.0')
     1190        sage: chop(ellipfun('ns', inverse_jacobi_f('ns', -10, 6), 6))
     1191        mpf('-10.0')
     1192
     1193        sage: chop(ellipfun('cn', inverse_jacobi_f('cn', -10, 0), 0))
     1194        mpf('-9.9999999999999982')
     1195        sage: chop(ellipfun('cn', inverse_jacobi_f('cn', 50, 1), 1))
     1196        mpf('50.000000000000071')
     1197        sage: chop(ellipfun('cn', inverse_jacobi_f('cn', 1, 5), 5))
     1198        mpf('1.0')
     1199        sage: chop(ellipfun('cn', inverse_jacobi_f('cn', 0.5, -5), -5))
     1200        mpf('0.5')
     1201        sage: chop(ellipfun('cn', inverse_jacobi_f('cn', -0.75, -15), -15))
     1202        mpf('-0.75000000000000022')
     1203        sage: chop(ellipfun('cn', inverse_jacobi_f('cn', 10, 0.8), 0.8))
     1204        mpf('9.9999999999999982')
     1205        sage: chop(ellipfun('cn', inverse_jacobi_f('cn', -2, 0.9), 0.9))
     1206        mpf('-2.0')
     1207
     1208        sage: chop(ellipfun('nc', inverse_jacobi_f('nc', -4, 0), 0))
     1209        mpf('-3.9999999999999987')
     1210        sage: chop(ellipfun('nc', inverse_jacobi_f('nc', 7, 1), 1))
     1211        mpf('7.0000000000000009')
     1212        sage: chop(ellipfun('nc', inverse_jacobi_f('nc', 7, 3), 3))
     1213        mpf('7.0')
     1214        sage: chop(ellipfun('nc', inverse_jacobi_f('nc', 0, 2), 2))
     1215        mpf('0.0')
     1216        sage: chop(ellipfun('nc', inverse_jacobi_f('nc', -18, -4), -4))
     1217        mpf('-17.999999999999925')
     1218
     1219        sage: chop(ellipfun('dn', inverse_jacobi_f('dn', -0.3, 1), 1))
     1220        mpf('-0.29999999999999999')
     1221        sage: chop(ellipfun('dn', inverse_jacobi_f('dn', 1, -1), -1))
     1222        mpf('1.0')
     1223        sage: chop(ellipfun('dn', inverse_jacobi_f('dn', 0.8, 0.5), 0.5))
     1224        mpf('0.80000000000000004')
     1225        sage: chop(ellipfun('dn', inverse_jacobi_f('dn', 5, -4), -4))
     1226        mpf('5.0')
     1227        sage: chop(ellipfun('dn', inverse_jacobi_f('dn', 0.4, 0.5), 0.5))
     1228        mpf('0.40000000000000002')
     1229        sage: chop(ellipfun('dn', inverse_jacobi_f('dn', -0.4, 0.5), 0.5))
     1230        mpf('-0.40000000000000002')
     1231        sage: chop(ellipfun('dn', inverse_jacobi_f('dn', -0.9, 0.5), 0.5))
     1232        mpf('-0.90000000000000002')
     1233        sage: chop(ellipfun('dn', inverse_jacobi_f('dn', -1.9, 0.2), 0.2))
     1234        mpf('-1.8999999999999999')
     1235
     1236        sage: chop(ellipfun('nd', inverse_jacobi_f('nd', -1.9, 1), 1))
     1237        mpf('-1.8999999999999999')
     1238        sage: chop(ellipfun('nd', inverse_jacobi_f('nd', 1, -1), -1))
     1239        mpf('1.0')
     1240        sage: chop(ellipfun('nd', inverse_jacobi_f('nd', 11, -6), -6))
     1241        mpf('11.0')
     1242        sage: chop(ellipfun('nd', inverse_jacobi_f('nd', 0, 8), 8))
     1243        mpf('0.0')
     1244        sage: chop(ellipfun('nd', inverse_jacobi_f('nd', -3, 0.8), 0.8))
     1245        mpf('-2.9999999999999996')
     1246
     1247        sage: chop(ellipfun('sc', inverse_jacobi_f('sc', -3, 0), 0))
     1248        mpf('-3.0')
     1249        sage: chop(ellipfun('sc', inverse_jacobi_f('sc', 2, 1), 1))
     1250        mpf('2.0')
     1251        sage: chop(ellipfun('sc', inverse_jacobi_f('sc', 0, 9), 9))
     1252        mpf('0.0')
     1253        sage: chop(ellipfun('sc', inverse_jacobi_f('sc', -7, 3), 3))
     1254        mpf('-7.0')
     1255
     1256        sage: chop(ellipfun('cs', inverse_jacobi_f('cs', -7, 0), 0))
     1257        mpf('-6.9999999999999991')
     1258        sage: chop(ellipfun('cs', inverse_jacobi_f('cs', 8, 1), 1))
     1259        mpf('8.0')
     1260        sage: chop(ellipfun('cs', inverse_jacobi_f('cs', 2, 6), 6))
     1261        mpf('2.0')
     1262        sage: chop(ellipfun('cs', inverse_jacobi_f('cs', 0, 4), 4))
     1263        mpf('0.0')
     1264        sage: chop(ellipfun('cs', inverse_jacobi_f('cs', -6, 8), 8))
     1265        mpf('-6.0000000000000018')
     1266
     1267        sage: chop(ellipfun('cd', inverse_jacobi_f('cd', -6, 0), 0))
     1268        mpf('-6.0000000000000009')
     1269        sage: chop(ellipfun('cd', inverse_jacobi_f('cd', 1, 3), 3))
     1270        mpf('1.0')
     1271        sage: chop(ellipfun('cd', inverse_jacobi_f('cd', 6, 8), 8))
     1272        mpf('6.0000000000000027')
     1273
     1274        sage: chop(ellipfun('dc', inverse_jacobi_f('dc', 5, 0), 0))
     1275        mpf('5.0000000000000018')
     1276        sage: chop(ellipfun('dc', inverse_jacobi_f('dc', -4, 2), 2))
     1277        mpf('-4.0000000000000018')
     1278
     1279        sage: chop(ellipfun('sd', inverse_jacobi_f('sd', -4, 0), 0))
     1280        mpf('-3.9999999999999991')
     1281        sage: chop(ellipfun('sd', inverse_jacobi_f('sd', 7, 1), 1))
     1282        mpf('7.0')
     1283        sage: chop(ellipfun('sd', inverse_jacobi_f('sd', 0, 9), 9))
     1284        mpf('0.0')
     1285        sage: chop(ellipfun('sd', inverse_jacobi_f('sd', 8, 0.8), 0.8))
     1286        mpf('7.9999999999999991')
     1287
     1288        sage: chop(ellipfun('ds', inverse_jacobi_f('ds', 4, 0.25), 0.25))
     1289        mpf('4.0')
     1290
     1291    """
     1292    from mpmath import mp
     1293
     1294    ctx = mp
     1295    prec = ctx.prec
     1296    try:
     1297        x = ctx.convert(x)
     1298        m = ctx.convert(m)
     1299        if not isinstance(x, ctx.mpf) or not isinstance(x, ctx.mpf):
     1300            raise ValueError('arguments must be real')
     1301        if kind == 'sn':
     1302            if m == 0:
     1303                return ctx.asin(x)
     1304            elif m == 1:
     1305                return ctx.atanh(x)
     1306            elif x == 0:
     1307                return ctx.zero
     1308            sign = ctx.sign(x)  # sn is odd in x, so operate with abs(x) and
     1309            x = abs(x)          # include the sign at the end
     1310            if x <= 1:
     1311                ctx.prec += 10
     1312                phi = ctx.asin(x)
     1313                return sign * ctx.ellipf(phi, m)
     1314            elif 1 < x <= 1 / ctx.sqrt(m):
     1315                K = ctx.ellipk(m)
     1316                ctx.prec += 10
     1317                xpn2 = x ** (-2)
     1318                m1 = 1 - m
     1319                ctx.prec += 10
     1320                omxpn2 = 1 - xpn2
     1321                ctx.prec += 10
     1322                omxpn2dm1 = omxpn2 / m1
     1323                ctx.prec += 10
     1324                phi = ctx.asin(omxpn2dm1.sqrt())
     1325                return sign * ctx.mpc(K, ctx.ellipf(phi, m1))
     1326            else:
     1327                ctx.prec += 10
     1328                m1 = 1 - m
     1329                K_prime = ctx.ellipk(m1)
     1330                sqrtm = ctx.sqrt(m)
     1331                ctx.prec += 10
     1332                xsqrtm = x * sqrtm
     1333                ctx.prec += 10
     1334                phi = ctx.asin(1 / xsqrtm)
     1335                ctx.prec += 10
     1336                return sign * ctx.mpc(ctx.ellipf(phi, m), K_prime)
     1337        if kind == 'ns':
     1338            if m == 0:
     1339                return ctx.acsc(x)
     1340            elif m == 1:
     1341                return ctx.acoth(x)
     1342            elif x > 0:
     1343                ctx.prec += 10
     1344                return inverse_jacobi_f('sn', 1 / x, m)
     1345            elif x == 0:
     1346                ctx.prec += 10
     1347                return ctx.j * ctx.ellipk(1 - m)
     1348            else:
     1349                ctx.prec += 10
     1350                K_prime = ctx.ellipk(1 - m)
     1351                odx = 1 / x
     1352                ctx.prec += 10
     1353                arcsnodx = inverse_jacobi_f('sn', odx, m)
     1354                itK_prime = ctx.j * 2 * K_prime
     1355                ctx.prec += 10
     1356                return arcsnodx + itK_prime
     1357        if kind == 'cn':
     1358            if m == 0:
     1359                return ctx.acos(x)
     1360            elif m == 1:
     1361                return ctx.asech(x)
     1362            elif x == 1:
     1363                return ctx.zero
     1364            elif 0 <= x < 1:
     1365                ctx.prec += 10
     1366                x2 = x ** 2
     1367                ctx.prec += 10
     1368                osx2 = 1 - x2
     1369                ctx.prec += 10
     1370                return ctx.ellipf(ctx.asin(ctx.sqrt(osx2)), m)
     1371            elif -1 <= x < 0:
     1372                K = ctx.ellipk(m)
     1373                ctx.prec += 10
     1374                x2 = x ** 2
     1375                ctx.prec += 10
     1376                osx2 = 1 - x2
     1377                ctx.prec += 10
     1378                return (2 * K) - ctx.ellipf(ctx.asin(ctx.sqrt(osx2)), m)
     1379            elif x > 1:
     1380                ctx.prec += 10
     1381                m1 = 1 - m
     1382                xn2 = x ** (-2)
     1383                ctx.prec += 10
     1384                osx2 = 1 - xn2
     1385                ctx.prec += 10
     1386                return ctx.j * ctx.ellipf(ctx.asin(ctx.sqrt(osx2)), m1)
     1387            elif x < -1:
     1388                K = ctx.ellipk(m)
     1389                ctx.prec += 10
     1390                m1 = 1 - m
     1391                xn2 = x ** (-2)
     1392                tK = 2 * K
     1393                ctx.prec += 10
     1394                osx2 = 1 - xn2
     1395                ctx.prec += 10
     1396                phi = ctx.asin(ctx.sqrt(osx2))
     1397                ctx.prec += 10
     1398                return tK - ctx.j * ctx.ellipf(phi, m1)
     1399        if kind == 'nc':
     1400            if m == 0:
     1401                return ctx.asec(x)
     1402            elif m == 1:
     1403                return ctx.acosh(x)
     1404            elif x == 1:
     1405                return ctx.zero
     1406            elif x > 0:
     1407                ctx.prec += 10
     1408                return inverse_jacobi_f('cn', 1 / x, m)
     1409            elif x == 0:
     1410                ctx.prec += 10
     1411                return ctx.j * ctx.ellipk(1 - m)
     1412            else:
     1413                K = ctx.ellipk(m)
     1414                ctx.prec += 10
     1415                K_prime = ctx.ellipk(1 - m)
     1416                odx = 1 / x
     1417                ctx.prec += 10
     1418                arccnodx = inverse_jacobi_f('cn', odx, m)
     1419                tK = 2 * K
     1420                ctx.prec += 10
     1421                return arccnodx - tK + ctx.j * 2 * K_prime
     1422        if kind == 'dn':
     1423            if x == 1:
     1424                return ctx.zero
     1425            if not m <= 1:
     1426                raise ValueError('m must be <= 1')
     1427            if m == 1:
     1428                return ctx.asech(x)
     1429            ctx.prec += 10
     1430            m1 = 1 - m
     1431            sqrtm1 = ctx.sqrt(m1)
     1432            if sqrtm1 <= x < 1:
     1433                ctx.prec += 10
     1434                x2 = x ** 2
     1435                ctx.prec += 10
     1436                osx2 = 1 - x2
     1437                ctx.prec += 10
     1438                osx2dm = osx2 / m
     1439                ctx.prec += 10
     1440                return ctx.ellipf(ctx.asin(ctx.sqrt(osx2dm)), m)
     1441            elif x > 1:
     1442                ctx.prec += 10
     1443                xn2 = x ** (-2)
     1444                ctx.prec += 10
     1445                osxn2 = 1 - xn2
     1446                m1xn2 = m1 * xn2
     1447                ctx.prec += 10
     1448                osm1xn2 = 1 - m1xn2
     1449                ctx.prec += 10
     1450                sqrtosxn2dosm1xn2 = ctx.sqrt(osxn2 / osm1xn2)
     1451                ctx.prec += 10
     1452                return ctx.j * ctx.ellipf(ctx.asin(sqrtosxn2dosm1xn2), m1)
     1453            elif 0 <= x < sqrtm1:
     1454                K = ctx.ellipk(m)
     1455                ctx.prec += 10
     1456                x2 = x ** 2
     1457                ctx.prec += 10
     1458                x2dm1 = x2 / m1
     1459                osx2 = 1 - x2
     1460                ctx.prec += 10
     1461                osx2dm1 = 1 - x2dm1
     1462                ctx.prec += 10
     1463                osx2dm1dosx2 = osx2dm1 / osx2
     1464                ctx.prec += 10
     1465                sqrtall = ctx.sqrt(osx2dm1dosx2)
     1466                ctx.prec += 10
     1467                phi = ctx.asin(sqrtall)
     1468                ctx.prec += 10
     1469                return K + ctx.j * ctx.ellipf(phi, m1)
     1470            elif -sqrtm1 <= x < 0:
     1471                K = ctx.ellipk(m)
     1472                K_prime = ctx.ellipk(m1)
     1473                ctx.prec += 10
     1474                tK_prime = 2 * K_prime
     1475                x2 = x ** 2
     1476                ctx.prec += 10
     1477                x2dm1 = x2 / m1
     1478                osx2 = 1 - x2
     1479                ctx.prec += 10
     1480                osx2dm1 = 1 - x2dm1
     1481                ctx.prec += 10
     1482                osx2dm1dosx2 = osx2dm1 / osx2
     1483                ctx.prec += 10
     1484                sqrtall = ctx.sqrt(osx2dm1dosx2)
     1485                ctx.prec += 10
     1486                phi = ctx.asin(sqrtall)
     1487                ctx.prec += 10
     1488                return K + ctx.j * (tK_prime - ctx.ellipf(phi, m1))
     1489            elif -1 <= x < -sqrtm1:
     1490                K = ctx.ellipk(m)
     1491                K_prime = ctx.ellipk(m1)
     1492                ctx.prec += 10
     1493                x2 = x ** 2
     1494                tK = 2 * K
     1495                # Note that the factor of 2 is missing in the reference
     1496                # (formula (81)), probably mistakenly so
     1497                tK_prime = 2 * K_prime
     1498                ctx.prec += 10
     1499                osx2 = 1 - x2
     1500                ctx.prec += 10
     1501                osx2dm = osx2 / m
     1502                sqrtall = ctx.sqrt(osx2dm)
     1503                ctx.prec += 10
     1504                phi = ctx.asin(sqrtall)
     1505                ctx.prec += 10
     1506                return (tK - ctx.ellipf(phi, m)) + (ctx.j * tK_prime)
     1507            elif x < -1:
     1508                K = ctx.ellipk(m)
     1509                K_prime = ctx.ellipk(m1)
     1510                ctx.prec += 10
     1511                tK = 2 * K
     1512                tK_prime = 2 * K_prime
     1513                xn2 = x ** (-2)
     1514                ctx.prec += 10
     1515                osxn2 = 1 - xn2
     1516                m1xn2 = m1 * xn2
     1517                ctx.prec += 10
     1518                osm1xn2 = 1 - m1xn2
     1519                ctx.prec += 10
     1520                sqrtosxn2dosm1xn2 = ctx.sqrt(osxn2 / osm1xn2)
     1521                ctx.prec += 10
     1522                phi = ctx.asin(sqrtosxn2dosm1xn2)
     1523                ctx.prec += 10
     1524                return tK + ctx.j * (tK_prime - ctx.ellipf(phi, m1))
     1525        if kind == 'nd':
     1526            if m == 1:
     1527                return ctx.acosh(x)
     1528            elif x == 1:
     1529                return ctx.zero
     1530            elif x > 0:
     1531                ctx.prec += 10
     1532                return inverse_jacobi_f('dn', 1 / x, m)
     1533            elif x == 0:
     1534                ctx.prec += 10
     1535                return ctx.j * ctx.ellipk(1 - m)
     1536            else:
     1537                K = ctx.ellipk(m)
     1538                ctx.prec += 10
     1539                tK = 2 * K
     1540                ctx.prec += 10
     1541                return inverse_jacobi_f('dn', 1 / x, m) - tK
     1542        if kind == 'sc':
     1543            if m == 0:
     1544                return ctx.atan(x)
     1545            elif m == 1:
     1546                return ctx.asinh(x)
     1547            elif x == 0:
     1548                return ctx.zero
     1549            else:
     1550                ctx.prec += 10
     1551                atanx = ctx.atan(x)
     1552                return ctx.ellipf(atanx, m)
     1553        if kind == 'cs':
     1554            if m == 0:
     1555                return ctx.acot(x)
     1556            elif m == 1:
     1557                return ctx.acsch(x)
     1558            elif x > 0:
     1559                ctx.prec += 10
     1560                odx = 1 / x
     1561                ctx.prec += 10
     1562                return ctx.ellipf(ctx.atan(odx), m)
     1563            elif x == 0:
     1564                return ctx.ellipk(m)
     1565            else:
     1566                K = ctx.ellipk(m)
     1567                ctx.prec += 10
     1568                odx = 1 / x
     1569                ctx.prec += 10
     1570                phi = ctx.atan(odx)
     1571                ctx.prec += 10
     1572                return ctx.ellipf(phi, m) + (2 * K)
     1573        if kind == 'cd':
     1574            if m == 0:
     1575                return ctx.acos(x)
     1576            elif x == 1:
     1577                return ctx.zero
     1578            else:
     1579                K = ctx.ellipk(m)
     1580                ctx.prec += 10
     1581                return inverse_jacobi_f('sn', x, m) - K
     1582        if kind == 'dc':
     1583            if m == 0:
     1584                return ctx.asec(x)
     1585            K = ctx.ellipk(m)
     1586            ctx.prec += 10
     1587            return inverse_jacobi_f('ns', x, m) - K
     1588        if kind == 'sd':
     1589            if m == 0:
     1590                return ctx.asin(x)
     1591            elif m == 1:
     1592                return ctx.asinh(x)
     1593            elif x == 0:
     1594                return ctx.zero
     1595            else:
     1596                if m > 1:
     1597                    raise ValueError('m must be <= 1')
     1598                K = ctx.ellipk(m)
     1599                ctx.prec += 10
     1600                m1 = 1 - m
     1601                ctx.prec += 10
     1602                sqrtm1 = ctx.sqrt(m1)
     1603                ctx.prec += 10
     1604                xsqrtm1 = x * sqrtm1
     1605                ctx.prec += 10
     1606                return inverse_jacobi_f('cn', xsqrtm1, m) + K
     1607        if kind == 'ds':
     1608            if m == 0:
     1609                return ctx.acsc(x)
     1610            elif m == 1:
     1611                return ctx.acsch(x)
     1612            else:
     1613                if m > 1:
     1614                    raise ValueError('m must be <= 1')
     1615                K = ctx.ellipk(m)
     1616                ctx.prec += 10
     1617                m1 = 1 - m
     1618                ctx.prec += 10
     1619                sqrtm1 = ctx.sqrt(m1)
     1620                ctx.prec += 10
     1621                xdsqrtm1 = x / sqrtm1
     1622                ctx.prec += 10
     1623                return inverse_jacobi_f('nc', xdsqrtm1, m) + K
     1624    finally:
     1625        ctx.prec = prec
     1626
     1627
     1628def jacobi_am_f(x, m):
     1629    """
     1630    Internal function for numeric evaluation of the Jacobi amplitude function
     1631    for real arguments. Procedure described in [Ehrhardt]_.
     1632
     1633    REFERENCES:
     1634
     1635    .. [Ehrhardt] Ehrhardt, Wolfgang. "The AMath and DAMath Special Functions:
     1636                  Reference Manual and Implementation Notes, Version 1.3".
     1637                  2013. http://www.wolfgang-ehrhardt.de/specialfunctions.pdf.
     1638
     1639
     1640    TESTS::
     1641
     1642        sage: from mpmath import ellipf
     1643        sage: from sage.functions.jacobi import jacobi_am_f
     1644        sage: ellipf(jacobi_am_f(0.5, 1), 1)
     1645        mpf('0.5')
     1646        sage: ellipf(jacobi_am(3, 0.3), 0.3)
     1647        mpf('3.0')
     1648        sage: ellipf(jacobi_am_f(2, -0.5), -0.5)
     1649        mpf('2.0')
     1650        sage: jacobi_am_f(2, -0.5)
     1651        mpf('2.2680930777934176')
     1652        sage: jacobi_am_f(-2, -0.5)
     1653        mpf('-2.2680930777934176')
     1654        sage: jacobi_am_f(-3, 2)
     1655        mpf('0.36067407399586108')
     1656    """
     1657    from mpmath import mp
     1658
     1659    ctx = mp
     1660    prec = ctx.prec
     1661    try:
     1662        x = ctx.convert(x)
     1663        m = ctx.convert(m)
     1664        if not isinstance(x, ctx.mpf) or not isinstance(m, ctx.mpf):
     1665            raise ValueError('arguments must be real')
     1666        if abs(m) == 1:
     1667            # gd(x)
     1668            ctx.prec += 10
     1669            tanhx = ctx.tanh(x)
     1670            ctx.prec += 10
     1671            return ctx.asin(tanhx)
     1672        elif abs(m) > 1:
     1673            ctx.prec += 10
     1674            # Real values needed for atan2; as per "Handbook of Elliptic
     1675            # Integrals for Engineers and Scientists" 121.02, sn is real for
     1676            # real x. The imaginary components can thus be safely discarded.
     1677            snx = ctx.ellipfun('sn', x, m).real
     1678            cnx = ctx.ellipfun('cn', x, m).real
     1679            ctx.prec += 10
     1680            return ctx.atan2(snx, cnx)
     1681        else:
     1682            ctx.prec += 10
     1683            K = ctx.ellipk(m)
     1684            if abs(x) <= K:
     1685                snx = ctx.ellipfun('sn', x, m).real
     1686                cnx = ctx.ellipfun('cn', x, m).real
     1687                ctx.prec += 10
     1688                return ctx.atan2(snx, cnx)
     1689            else:
     1690                # Do argument reduction on x to end up with z = x - 2nK, with
     1691                # abs(z) <= K
     1692                ctx.prec += 10
     1693                tK = 2 * K
     1694                ctx.prec += 10
     1695                n = ctx.floor(x / tK)
     1696                ctx.prec += 10
     1697                tnK = n * tK
     1698                npi = n * ctx.pi()
     1699                ctx.prec += 10
     1700                z = x - tnK
     1701                ctx.prec += 10
     1702                # z (and therefore sn(z, m) and cn(z, m)) is real because K(m)
     1703                # is real for abs(m) <= 1.
     1704                snz = ctx.ellipfun('sn', z, m).real
     1705                cnz = ctx.ellipfun('cn', z, m).real
     1706                ctx.prec += 10
     1707                return ctx.atan2(snz, cnz) + npi
     1708    finally:
     1709        ctx.prec = prec
  • sage/functions/special.py

    diff --git a/sage/functions/special.py b/sage/functions/special.py
    a b  
    148148   solution, called Kummer's function `M(a,b,x)`, which is not
    149149   implemented.)
    150150
    151 -  Jacobi elliptic functions can be thought of as generalizations
    152    of both ordinary and hyperbolic trig functions. There are twelve
    153    Jacobian elliptic functions. Each of the twelve corresponds to an
    154    arrow drawn from one corner of a rectangle to another.
    155 
    156    ::
    157 
    158                     n ------------------- d
    159                     |                     |
    160                     |                     |
    161                     |                     |
    162                     s ------------------- c
    163 
    164    Each of the corners of the rectangle are labeled, by convention, s,
    165    c, d and n. The rectangle is understood to be lying on the complex
    166    plane, so that s is at the origin, c is on the real axis, and n is
    167    on the imaginary axis. The twelve Jacobian elliptic functions are
    168    then pq(x), where p and q are one of the letters s,c,d,n.
    169 
    170    The Jacobian elliptic functions are then the unique
    171    doubly-periodic, meromorphic functions satisfying the following
    172    three properties:
    173 
    174    
    175    #. There is a simple zero at the corner p, and a simple pole at the
    176       corner q.
    177 
    178    #. The step from p to q is equal to half the period of the function
    179       pq(x); that is, the function pq(x) is periodic in the direction pq,
    180       with the period being twice the distance from p to q. Also, pq(x)
    181       is also periodic in the other two directions as well, with a period
    182       such that the distance from p to one of the other corners is a
    183       quarter period.
    184 
    185    #. If the function pq(x) is expanded in terms of x at one of the
    186       corners, the leading term in the expansion has a coefficient of 1.
    187       In other words, the leading term of the expansion of pq(x) at the
    188       corner p is x; the leading term of the expansion at the corner q is
    189       1/x, and the leading term of an expansion at the other two corners
    190       is 1.
    191 
    192 
    193    We can write
    194    
    195    .. math::
    196 
    197       pq(x)=\frac{pr(x)}{qr(x)}
    198 
    199 
    200    where `p`, `q`, and `r` are any of the
    201    letters `s`, `c`, `d`, `n`, with
    202    the understanding that `ss=cc=dd=nn=1`.
    203 
    204    Let
    205    
    206    .. math::
    207 
    208          u=\int_0^\phi \frac{d\theta} {\sqrt {1-m \sin^2 \theta}}
    209 
    210 
    211    Then the *Jacobi elliptic function* `sn(u)` is given by
    212    
    213    .. math::
    214 
    215         {sn}\; u = \sin \phi
    216 
    217    and `cn(u)` is given by
    218    
    219    .. math::
    220 
    221        {cn}\; u = \cos \phi
    222 
    223    and
    224    
    225    .. math::
    226 
    227      {dn}\; u = \sqrt {1-m\sin^2 \phi}. 
    228 
    229    To emphasize the dependence on `m`, one can write
    230    `sn(u,m)` for example (and similarly for `cn` and
    231    `dn`). This is the notation used below.
    232 
    233    For a given `k` with `0 < k < 1` they therefore are
    234    solutions to the following nonlinear ordinary differential
    235    equations:
    236 
    237    
    238    -  `\mathrm{sn}\,(x;k)` solves the differential equations
    239      
    240       .. math::
    241 
    242           \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} + (1+k^2) y - 2 k^2 y^3 = 0,
    243          
    244       and
    245 
    246       .. math::
    247 
    248          \left(\frac{\mathrm{d} y}{\mathrm{d}x}\right)^2 = (1-y^2) (1-k^2 y^2).
    249 
    250    -  `\mathrm{cn}\,(x;k)` solves the differential equations
    251 
    252      
    253       .. math::
    254      
    255          \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} + (1-2k^2) y + 2 k^2 y^3 = 0,
    256 
    257       and `\left(\frac{\mathrm{d} y}{\mathrm{d}x}\right)^2 = (1-y^2) (1-k^2 + k^2 y^2)`.
    258 
    259    -  `\mathrm{dn}\,(x;k)` solves the differential equations
    260      
    261       .. math::
    262 
    263          \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} - (2 - k^2) y + 2 y^3 = 0,
    264 
    265       and `\left(\frac{\mathrm{d} y}{\mathrm{d}x}\right)^2= y^2 (1 - k^2 - y^2)`.
    266 
    267       If `K(m)` denotes the complete elliptic integral of the
    268       first kind (denoted ``elliptic_kc``), the elliptic functions
    269       `sn (x,m)` and `cn (x,m)` have real periods
    270       `4K(m)`, whereas `dn (x,m)` has a period
    271       `2K(m)`. The limit `m\rightarrow 0` gives
    272       `K(0) = \pi/2` and trigonometric functions:
    273       `sn(x, 0) = \sin x`, `cn(x, 0) = \cos x`,
    274       `dn(x, 0) = 1`. The limit `m \rightarrow 1` gives
    275       `K(1) \rightarrow \infty` and hyperbolic functions:
    276       `sn(x, 1) = \tanh x`,
    277       `cn(x, 1) = \mbox{\rm sech} x`,
    278       `dn(x, 1) = \mbox{\rm sech} x`.
    279 
    280151   -  The incomplete elliptic integrals (of the first kind, etc.) are:
    281152     
    282153      .. math::
     
    297168
    298169- http://en.wikipedia.org/wiki/Helmholtz_equation
    299170
    300 - http://en.wikipedia.org/wiki/Jacobi's_elliptic_functions
    301 
    302 - A. Khare, U. Sukhatme, 'Cyclic Identities Involving
    303   Jacobi Elliptic Functions', Math ArXiv, math-ph/0201004
    304 
    305171- Online Encyclopedia of Special Function
    306172  http://algo.inria.fr/esf/index.html
    307173
     
    433299        Check if complex numbers in the arguments are converted to maxima
    434300        correctly :trac:`7557`::
    435301
    436             sage: t = jacobi('sn',1.2+2*I*elliptic_kc(1-.5),.5)
     302            sage: t = f(1.2+2*I*elliptic_kc(1-.5),.5)
    437303            sage: t._maxima_init_(maxima)
    438304            '0.88771548861927...*%i'
    439305            sage: t.n() # abs tol 1e-13
     
    808674   from sage.libs.all import pari
    809675   return CC(pari(z).ellj())
    810676
    811 def jacobi(sym,x,m):
    812     r"""
    813     Here sym = "pq", where p,q in c,d,n,s. This returns the value of
    814     the Jacobi function pq(x,m), as described in the documentation for
    815     Sage's "special" module. There are a total of 12 functions
    816     described by this.
    817    
    818     EXAMPLES::
    819    
    820         sage: jacobi("sn",1,1)
    821         tanh(1)
    822         sage: jacobi("cd",1,1/2)
    823         jacobi_cd(1, 1/2)
    824         sage: RDF(jacobi("cd",1,1/2))
    825         0.724009721659
    826         sage: RDF(jacobi("cn",1,1/2)); RDF(jacobi("dn",1,1/2)); RDF(jacobi("cn",1,1/2)/jacobi("dn",1,1/2))
    827         0.595976567672
    828         0.823161001632
    829         0.724009721659
    830         sage: jsn = jacobi("sn",x,1)
    831         sage: P = plot(jsn,0,1)
    832    
    833     To view this, type P.show().
    834     """
    835     return maxima_function("jacobi_%s"%sym)(x,m)
    836 
    837 def inverse_jacobi(sym,x,m):
    838     """
    839     Here sym = "pq", where p,q in c,d,n,s. This returns the value of
    840     the inverse Jacobi function `pq^{-1}(x,m)`. There are a
    841     total of 12 functions described by this.
    842    
    843     EXAMPLES::
    844    
    845         sage: jacobi("sn",1/2,1/2)
    846         jacobi_sn(1/2, 1/2)
    847         sage: float(jacobi("sn",1/2,1/2))
    848         0.4707504736556572
    849         sage: float(inverse_jacobi("sn",0.47,1/2))
    850         0.4990982313222196
    851         sage: float(inverse_jacobi("sn",0.4707504,0.5))
    852         0.4999999114665548
    853         sage: P = plot(inverse_jacobi('sn', x, 0.5), 0, 1, plot_points=20)
    854    
    855     Now to view this, just type show(P).
    856     """
    857     return maxima_function("inverse_jacobi_%s"%sym)(x,m)
    858 
    859677#### elliptic integrals
    860678
    861679class EllipticE(MaximaFunction):
  • sage/plot/line.py

    diff --git a/sage/plot/line.py b/sage/plot/line.py
    a b  
    409409 
    410410    A red plot of the Jacobi elliptic function `\text{sn}(x,2)`, `-3 < x < 3`::
    411411 
    412         sage: L = [(i/100.0, jacobi('sn', i/100.0 ,2.0)) for i in range(-300,300,30)]   
    413         sage: line(L, rgbcolor=(3/4,1/4,1/8))
    414        
     412        sage: L = [(i/100.0, real_part(jacobi('sn', i/100.0, 2.0))) for i in
     413        ....:      range(-300, 300, 30)]
     414        sage: line(L, rgbcolor=(3/4, 1/4, 1/8))
     415
    415416    A red plot of `J`-Bessel function `J_2(x)`, `0 < x < 10`::
    416417
    417418        sage: L = [(i/10.0, bessel_J(2,i/10.0)) for i in range(100)]