Ticket #14996: trac14996_3.patch

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