Ticket #14996: trac14996_2.patch

File trac14996_2.patch, 61.6 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 c33bbd792c85c262a994dbe7f33a527615d6e746
    # Parent  0f8fd922eaed351e39f913f1317d319dcceb4c01
    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/exp_integral
    1516   sage/functions/wigner
    1617   sage/functions/generalized
  • sage/functions/all.py

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

    diff --git a/sage/functions/special.py b/sage/functions/special.py
    a b  
    190190   solution, called Kummer's function `M(a,b,x)`, which is not
    191191   implemented.)
    192192
    193 -  Jacobi elliptic functions can be thought of as generalizations
    194    of both ordinary and hyperbolic trig functions. There are twelve
    195    Jacobian elliptic functions. Each of the twelve corresponds to an
    196    arrow drawn from one corner of a rectangle to another.
    197 
    198    ::
    199 
    200                     n ------------------- d
    201                     |                     |
    202                     |                     |
    203                     |                     |
    204                     s ------------------- c
    205 
    206    Each of the corners of the rectangle are labeled, by convention, s,
    207    c, d and n. The rectangle is understood to be lying on the complex
    208    plane, so that s is at the origin, c is on the real axis, and n is
    209    on the imaginary axis. The twelve Jacobian elliptic functions are
    210    then pq(x), where p and q are one of the letters s,c,d,n.
    211 
    212    The Jacobian elliptic functions are then the unique
    213    doubly-periodic, meromorphic functions satisfying the following
    214    three properties:
    215 
    216    
    217    #. There is a simple zero at the corner p, and a simple pole at the
    218       corner q.
    219 
    220    #. The step from p to q is equal to half the period of the function
    221       pq(x); that is, the function pq(x) is periodic in the direction pq,
    222       with the period being twice the distance from p to q. Also, pq(x)
    223       is also periodic in the other two directions as well, with a period
    224       such that the distance from p to one of the other corners is a
    225       quarter period.
    226 
    227    #. If the function pq(x) is expanded in terms of x at one of the
    228       corners, the leading term in the expansion has a coefficient of 1.
    229       In other words, the leading term of the expansion of pq(x) at the
    230       corner p is x; the leading term of the expansion at the corner q is
    231       1/x, and the leading term of an expansion at the other two corners
    232       is 1.
    233 
    234 
    235    We can write
    236    
    237    .. math::
    238 
    239       pq(x)=\frac{pr(x)}{qr(x)}
    240 
    241 
    242    where `p`, `q`, and `r` are any of the
    243    letters `s`, `c`, `d`, `n`, with
    244    the understanding that `ss=cc=dd=nn=1`.
    245 
    246    Let
    247    
    248    .. math::
    249 
    250          u=\int_0^\phi \frac{d\theta} {\sqrt {1-m \sin^2 \theta}}
    251 
    252 
    253    Then the *Jacobi elliptic function* `sn(u)` is given by
    254    
    255    .. math::
    256 
    257         {sn}\; u = \sin \phi
    258 
    259    and `cn(u)` is given by
    260    
    261    .. math::
    262 
    263        {cn}\; u = \cos \phi
    264 
    265    and
    266    
    267    .. math::
    268 
    269      {dn}\; u = \sqrt {1-m\sin^2 \phi}. 
    270 
    271    To emphasize the dependence on `m`, one can write
    272    `sn(u,m)` for example (and similarly for `cn` and
    273    `dn`). This is the notation used below.
    274 
    275    For a given `k` with `0 < k < 1` they therefore are
    276    solutions to the following nonlinear ordinary differential
    277    equations:
    278 
    279    
    280    -  `\mathrm{sn}\,(x;k)` solves the differential equations
    281      
    282       .. math::
    283 
    284           \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} + (1+k^2) y - 2 k^2 y^3 = 0,
    285          
    286       and
    287 
    288       .. math::
    289 
    290          \left(\frac{\mathrm{d} y}{\mathrm{d}x}\right)^2 = (1-y^2) (1-k^2 y^2).
    291 
    292    -  `\mathrm{cn}\,(x;k)` solves the differential equations
    293 
    294      
    295       .. math::
    296      
    297          \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} + (1-2k^2) y + 2 k^2 y^3 = 0,
    298 
    299       and `\left(\frac{\mathrm{d} y}{\mathrm{d}x}\right)^2 = (1-y^2) (1-k^2 + k^2 y^2)`.
    300 
    301    -  `\mathrm{dn}\,(x;k)` solves the differential equations
    302      
    303       .. math::
    304 
    305          \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} - (2 - k^2) y + 2 y^3 = 0,
    306 
    307       and `\left(\frac{\mathrm{d} y}{\mathrm{d}x}\right)^2= y^2 (1 - k^2 - y^2)`.
    308 
    309       If `K(m)` denotes the complete elliptic integral of the
    310       first kind (denoted ``elliptic_kc``), the elliptic functions
    311       `sn (x,m)` and `cn (x,m)` have real periods
    312       `4K(m)`, whereas `dn (x,m)` has a period
    313       `2K(m)`. The limit `m\rightarrow 0` gives
    314       `K(0) = \pi/2` and trigonometric functions:
    315       `sn(x, 0) = \sin x`, `cn(x, 0) = \cos x`,
    316       `dn(x, 0) = 1`. The limit `m \rightarrow 1` gives
    317       `K(1) \rightarrow \infty` and hyperbolic functions:
    318       `sn(x, 1) = \tanh x`,
    319       `cn(x, 1) = \mbox{\rm sech} x`,
    320       `dn(x, 1) = \mbox{\rm sech} x`.
    321 
    322193   -  The incomplete elliptic integrals (of the first kind, etc.) are:
    323194     
    324195      .. math::
     
    341212
    342213- http://en.wikipedia.org/wiki/Helmholtz_equation
    343214
    344 - http://en.wikipedia.org/wiki/Jacobi's_elliptic_functions
    345 
    346 - A. Khare, U. Sukhatme, 'Cyclic Identities Involving
    347   Jacobi Elliptic Functions', Math ArXiv, math-ph/0201004
    348 
    349215- Online Encyclopedia of Special Function
    350216  http://algo.inria.fr/esf/index.html
    351217
     
    477343        Check if complex numbers in the arguments are converted to maxima
    478344        correctly :trac:`7557`::
    479345
    480             sage: t = jacobi('sn',1.2+2*I*elliptic_kc(1-.5),.5)
     346            sage: t = f(1.2+2*I*elliptic_kc(1-.5),.5)
    481347            sage: t._maxima_init_(maxima)
    482348            '0.88771548861927...*%i'
    483349            sage: t.n() # abs tol 1e-13
     
    14081274   from sage.libs.all import pari
    14091275   return CC(pari(z).ellj())
    14101276
    1411 def jacobi(sym,x,m):
    1412     r"""
    1413     Here sym = "pq", where p,q in c,d,n,s. This returns the value of
    1414     the Jacobi function pq(x,m), as described in the documentation for
    1415     Sage's "special" module. There are a total of 12 functions
    1416     described by this.
    1417    
    1418     EXAMPLES::
    1419    
    1420         sage: jacobi("sn",1,1)
    1421         tanh(1)
    1422         sage: jacobi("cd",1,1/2)
    1423         jacobi_cd(1, 1/2)
    1424         sage: RDF(jacobi("cd",1,1/2))
    1425         0.724009721659
    1426         sage: RDF(jacobi("cn",1,1/2)); RDF(jacobi("dn",1,1/2)); RDF(jacobi("cn",1,1/2)/jacobi("dn",1,1/2))
    1427         0.595976567672
    1428         0.823161001632
    1429         0.724009721659
    1430         sage: jsn = jacobi("sn",x,1)
    1431         sage: P = plot(jsn,0,1)
    1432    
    1433     To view this, type P.show().
    1434     """
    1435     return maxima_function("jacobi_%s"%sym)(x,m)
    1436 
    1437 def inverse_jacobi(sym,x,m):
    1438     """
    1439     Here sym = "pq", where p,q in c,d,n,s. This returns the value of
    1440     the inverse Jacobi function `pq^{-1}(x,m)`. There are a
    1441     total of 12 functions described by this.
    1442    
    1443     EXAMPLES::
    1444    
    1445         sage: jacobi("sn",1/2,1/2)
    1446         jacobi_sn(1/2, 1/2)
    1447         sage: float(jacobi("sn",1/2,1/2))
    1448         0.4707504736556572
    1449         sage: float(inverse_jacobi("sn",0.47,1/2))
    1450         0.4990982313222196
    1451         sage: float(inverse_jacobi("sn",0.4707504,0.5))
    1452         0.4999999114665548
    1453         sage: P = plot(inverse_jacobi('sn', x, 0.5), 0, 1, plot_points=20)
    1454    
    1455     Now to view this, just type show(P).
    1456     """
    1457     return maxima_function("inverse_jacobi_%s"%sym)(x,m)
    1458 
    14591277#### elliptic integrals
    14601278
    14611279class EllipticE(MaximaFunction):