20.19 ELLIPFN: A Package for Elliptic Functions and Integrals

Additional documentation and examples can be found in the files efellip.tex, eftheta.tst and efweier.tst in the packages/ellipfn directory.

Author: Lisa Temme and Alan Barnes with contributions from Winfried Neun and several others.

20.19.1 Elliptic Functions: Introduction

\(\newcommand {\f }[1]{\texttt {#1}}\)The package ELLIPFN is designed to provide algebraic and numeric manipulations of many elliptic functions, namely:

• Jacobi’s Elliptic Functions;

• Elliptic Integrals;

• Nome and Related Functions;

• Jacobi’s Theta Functions and their derivatives;

• Weierstrass Elliptic Functions and the Sigma Function;

• Other Sigma Functions;

• Period Lattice and Related Functions;

• Inverse Elliptic Functions.

The implementation of the functions in this and the next two subsections have been substantially revised by Alan Barnes in 2019. This is to bring the notation more into line with standard (British) texts such as Whittaker & Watson [WW69] and Lawden [Law89] and also to correct a number of errors and omissions. These changes and additions will be itemised in the relevant sections below. New subsections has been added starting in 2021 to support Weierstrassian Elliptic functions, Sigma functions, inverse Jacobi elliptic functions and finally symmetric elliptic integrals.

The functions in this subsection are for the most part autoloading; the exceptions being the subsidiary utility functions such as the AGM function, the quasi-period factors, lattice functions, derivatives of the theta functions and symmetric elliptic integrals.

20.19.2 Jacobi Elliptic Functions

The following functions have been implemented:

• The 12 Jacobi Elliptic Functions

• The Jacobi Amplitude Function

• Arithmetic Geometric Mean

• Descending Landen Transformation

The following Jacobi elliptic functions are available:-

• jacobisn(u,k)

• jacobidn(u,k)

• jacobicn(u,k)

• jacobicd(u,k)

• jacobisd(u,k)

• jacobind(u,k)

• jacobidc(u,k)

• jacobinc(u,k)

• jacobisc(u,k)

• jacobins(u,k)

• jacobids(u,k)

• jacobics(u,k)

These differ somewhat from the originals implemented by Lisa Temme in that the second argument is now the modulus (usually denoted by \(k\) in most texts rather than its square \(m\)). The notation for the most part follows Lawden [Law89]. The last nine Jacobi functions are related to the three basic ones: jacobisn(u,k), jacobicn(u,k) and jacobidn(u,k) and use Glaisher’s notation. For example \[ \mathrm {ns}(x,k) = \frac {1}{\mathrm {sn}(u,k)}, \qquad \mathrm {cs}(x,k) = \frac {\mathrm {cn}(u,k)}{\mathrm {sn}(u,k)}, \qquad \mathrm {cd}(x,k) = \frac {\mathrm {cn}(u,k)}{\mathrm {dn}(u,k)}. \]

All twelve functions are doubly periodic in the complex plane. The primitive periods and the positions of the zeros and poles of the functions are conveniently expressed in terms of the so-called quarter periods \(\mathrm {K}(k)\) and \(i\mathrm {K}'(k)\) which are complete elliptic integrals of the first kind. The details are displayed in the table below; in all cases the pole is single and the corresponding residue is given in the last column of the table. As the functions are doubly-periodic, there are, of course, infinitely many other zeros and poles. These occur at points ‘congruent’ to the points given in the table obtained by translating by \(2m\mathrm {K} +2i n\mathrm {K}'\) where \(m\) and \(n\) are arbitrary integers.

All other periods of the Jacobi functions can be expressed as linear combinations of the two primitive periods with integer coefficients. Thus, for example, \(4i\mathrm {K}'\) is a period of \(\mathrm {cn}\) as \[4i\mathrm {K}' = 2(\mathrm {K}+\mathrm {K}') -4\mathrm {K}\]

Extensive rule lists are provided for differentiation of these functions with respect to either argument, for argument shifts by integer multiples of the two quarter periods \(\mathrm {K}(k)\) and \(i\mathrm {K}'(k)\) and finally Jacobi’s transformations for a purely imaginary first argument.

Rules are also provided for the values of the twelve Jacobi functions at the ‘eighth’ period values: \(\mathrm {K}/2\), \(i\mathrm {K}'/2\) and \((\mathrm {K}+i\mathrm {K}')/2\). For these rules to yield correct values it is essential that the switch precise_complex is ON, otherwise Reduce will often incorrectly simplify the results if they involve complex values of \(k\) or \(k'\). It should also be noted that the rule for \(\mathrm {dn}((\mathrm {K}(k)+i\mathrm {K}'(k))/2, k)\) differs from that on DLMF:NIST which appears to give incorrect values for some values of \(k\). The rule used in REDUCE is not only simpler, but appears to give correct values in all cases as do the derived rules for \(\mathrm {cd}\), \(\mathrm {sd}\), \(\mathrm {nd}\), \(\mathrm {ds}\) and \(\mathrm {dc}\). It is derived from the third of the identities (2.2.19) of Lawden [Law89] with \(u =-(\mathrm {K}(k)+i\mathrm {K}'(k))/2\) and the corresponding identities for \(\mathrm {sn}\) and \(\mathrm {cn}\) on the DLMF NIST site.

Some care must be exercised when applying the ‘eighth’ period rules when the modulus \(k\) belongs to the negative real or imaginary axes as \(\mathrm {K}(k)\) and all twelve Jacobi functions \(\mathrm {sn}(x,k)\) etc. are even functions of the modulus \(k\). The rule should be applied when the modulus has no assigned value and then resimplifying the result after assigning the modulus its required value.

Complex split functions have recently been implemented for the 12 Jacobi functions for the cases when the modulus \(k\) is either real or imaginary. Thus REPART and IMPART will now return rational expressions involving the twelve Jacobi functions with real arguments and a real modulus \(|k|<1\). As far as I (AB) am aware there are no known methods of implementing split functions for general complex values of the modulus.

Four useful rule lists TO_SN, TO_CN, TO_DN and NO_GLAISHER are provided for user convenience. The first TO_SN replaces occurences of the squares of \(\mathrm {cn}\) and \(\mathrm {dn}\) in favour of the square of \(\mathrm {sn}\) using the ‘Pythagorean’ identities. TO_CN and TO_DN apply these identities to replace squares in favour of \(\mathrm {cn}\) and \(\mathrm {dn}\) respectively. At most one of these rule sets should be active at any one time to avoid a potential infinite recursion in the simplification. The fourth rule list NO_GLAISHER replaces all occurences of the nine subsidiary Jacobi elliptic functions \(\mathrm {ns}\), \(\mathrm {sc}\), \(\mathrm {sd}\), …by reciprocals and quotients of the ‘basic’ Jacobi elliptic functions \(\mathrm {sn}\), \(\mathrm {cn}\) and \(\mathrm {dn}\).

A fifth rule list JACOBIADDITIONRULES applies addition rules when the first argument of any Jacobi function is a sum. Note that this rule list is NO LONGER active by default. There are many equivalent ways of expressing the right-hand sides of these rules. Previously all the rules were expressed using only the ‘basic’ Jacobi functions: \(\mathrm {sn}\), \(\mathrm {cn}\) and \(\mathrm {dn}\). Currently, however, the right-hand sides of the rules for the reciprocal and quotient Glaisher functions \(\mathrm {ns}\), \(\mathrm {cs}\), \(\mathrm {sd}\) etc are expressed using the three Glaisher functions with the same ‘denominator’ as on the left-hand side. Thus the rule for \(\mathrm {ns}\) is expressed in terms of \(\mathrm {ns}\), \(\mathrm {cs}\) and \(\mathrm {ds}\) whilst that for \(\mathrm {cd}\) is expressed in terms of \(\mathrm {nd}\), \(\mathrm {cd}\) and \(\mathrm {sd}\) whereas the rules for the three ‘basic’ functions are unchanged and use only \(\mathrm {sn}\), \(\mathrm {cn}\) and \(\mathrm {dn}\). Note, however, that the previous form of these rules will be obtained if the rule-list NO_GLAISHER is also active.

If both the switches rounded and complex are ON and both arguments of a Jacobi function are purely numerical, it will be evaluated numerically The numerical evaluation of the Jacobi functions now uses their definitions in terms of theta functions as these are valid for all complex values of the argument and modulus. Note that as a consequence it is now necessary that the switch complex is ON even if both the argument and the modulus are real.

The traditional AGM and \(\phi \)-function based algorithm which was previously used when the argument and modulus are both real and \(|k| < 1\) as been superceded. It appears that there was a problem with the numerical evaluation of dn, nd, sd, ds, cd, dc resulting in largish rounding errors. There was no problem with previous evaluation of sn, ns, cn, nc, sc, cs and these produced the same results (modulo acceptable rounding errors) as the theta function based method now used.

Jacobi Amplitude Function

This function is defined as \[ \mathrm {am}(u,k) = \int _0^u\mathrm {dn}(z,k)\mathrm {d}z.\]

As \(\mathrm {dn}(z,k)\) has poles with residue \(-i\) at points \(z=(4m+1)i\mathrm {K}'(k)+2n\mathrm {K}(k)\) and others with residue \(+i\) when \(z=(4m+3)i\mathrm {K}'(k)+2n\mathrm {K}(k)\) where \(m\) and \(n\) are arbitrary integers, \(\mathrm {am}(z,k)\) has logarithmic singularities at these points. In fact the amplitude function is multivalued with its principal value \(v\), say, being given by choosing the contour in the defining integral to be the straight line segment between \(0\) and \(u\). Other values are given by \(v+2n\pi \) where \(n\) is an arbitrary integer and depend on how the chosen contour winds around the logarithmic branch points. To obtain a single valued function it is necessary to introduce branch cuts between the points \(i\mathrm {K}'(k)\) and \(3i\mathrm {K}'(k)\) and between corresponding congruent branch points.

A rule list is provided for to simplify this function for special values of the arguments and for differentiation with respect to either argument. When the switches ROUNDED and COMPLEX are both ON and both arguments are numeric, REDUCE will attempt to evaluate jacobiam(z,k) numerically. Several possible methods are available to evalaute the function, for example summation of its Fourier series and an AGM-based method due to Sala (see DLMF:NIST for more details. Currently the AGM-based method is used and this certainly produces reliable results if the modulus \(k\) is real and \(|\Im (z)| < \mathrm {K}'(k)\). These agree with those produced by the Fourier series method although rounding errors become significant as \(|\Im (z)|\) approaches \(\mathrm {K}'(k)\). If \(k\) is complex both methods produce the same results when the imaginary parts of \(k\) and \(z\) are not too large, but it is difficult to give precise bounds on the size of these.

For general values of \(z\) and \(k\) the AGM-based method always produces a value except perhaps at the logarithmic branch points but these, I believe, are not correct as they fail to satisfy the identity \[\mathrm {sn}(z,k) = \sin (\mathrm {am}(z,k)) \] where \(\mathrm {sn}\) is calculated using the theta function method. The Fourier series method fails to converge in these cases and so is not used except for comparison purposes. Currently an alternative method of evaluation due to Sala is being investigated; it uses the Poisson summation formula and is claimed to be valid for all \(z\) and \(k\) except at branch points.

20.19.3 Some Numerical Procedures

This section briefly describes several procedures which are primarily intended for use in the numerical evaluation of the various elliptic functions and integrals rather than for direct use by users.

Arithmetic Geometric Mean (AGM)

A procedure to evaluate the AGM of initial values \(a_0,b_0,c_0\) exists as
AGM_function(\(a_0,b_0,c_0\)) and will return
\(\{ N, AGM, \{ a_N, \ldots ,a_0\}, \{ b_N, \ldots ,b_0\}, \{c_N, \ldots ,c_0\}\}\), where \(N\) is the number of steps to compute the AGM to the desired accuracy.

To determine the Elliptic Integrals K(\(m\)), E(\(m\)) we use initial values \(a_0 = 1\); \(b_0 = \sqrt {1-k^2}\) ; \(c_0 = k\).

Descending Landen Transformation

The procedure to evaluate the Descending Landen Transformation of \(\phi \) and \(\alpha \) uses the following equations: \begin {align*} (1+\sin \alpha _{n+1})(1+\cos \alpha _n)=2 &\text { where } \alpha _{n+1}<\alpha _n, \\ \tan (\phi _{n+1}-\phi _n)=\cos \alpha _n \tan \phi _n & \text { where } \phi _{n+1}>\phi _n. \end {align*}

It can be called using landentrans(\(\phi _0\), \(\alpha _0\)) and will return
\(\{\{\phi _0, \ldots ,\phi _n\},\{\alpha _0, \ldots ,\alpha _n\}\}\).

20.19.4 Legendre’s Elliptic Integrals

The following functions have been implemented:

• Complete & Incomplete Elliptic Integrals of the First Kind

• Complete & Incomplete Elliptic Integrals of the Second Kind

• Jacobi’s Epsilon & Zeta Functions

• Complete & Incomplete Ellpitic Integrals of the Third Kind

• Some Utility Functions

These again differ somewhat from the originals implemented by Lisa Temme as the second argument is now the modulus \(k\) rather that its square. Also in the original implementation there was some confusion between Legendre’s form and Jacobi’s form of the incomplete elliptic integrals of the second kind; \(E(u,k)\) denoted the first in numerical evaluations and the second in the derivative formulae for the Jacobi elliptic functions with respect to their second argument. This confusion was perhaps understandable as in the literature some authors use the notation \(\mathrm {E}(u, k)\) for the Legendre form and others for Jacobi’s form.

To bring the notation more into line with that in the NIST Digital Library of Mathematical Functions and avoid any possible confusion, \(\mathrm {E}(u, k)\) is used for the Legendre form and \(\mathcal {E}(u, k)\) for Jacobi’s form. This differs from the 2019 version of this section which followed Lawden [Law89] chapter 3, where the notation \(\mathrm {D}(\phi , k)\) and \(\mathrm {E}(u, k)\) were used for the Legendre and Jacobi forms respectively.

A number of rule lists have been provided to implement, where appropriate, derivatives of these functions, addition rules and periodicity and quasi-periodicity properties and to provide simplifications for special values of the arguments.

Elliptic F

The Elliptic F function can be used as EllipticF(phi,k) and will return the value of the Incomplete Elliptic Integral of the First Kind: \[\mathrm {F}(\phi , k)=\int _0^\phi (1-k^2 \sin ^2 \theta )^{-1/2} \mathrm {d}\theta .\] This is actually closely related to the inverse Jacobi function \(\mathrm {arcsn}\); in fact for \(-\pi /2 <= \Re \phi <=\pi /2\): \[ \mathrm {F}(\phi , k) = \mathrm {arcsn}(\sin \phi , k) \].

Elliptic K

The Elliptic K function can be used as EllipticK(k) and will return the value of the Complete Elliptic Integral of the First Kind: \[\mathrm {K}(k)=\mathrm {F}(\pi /2, k) =\int _0^{\pi /2}(1-k^2 \sin ^2 \theta )^{-1/2}\mathrm {d}\theta .\] This is one of the quarter periods of the Jacobi elliptic functions and is often used in the calculation of other elliptic functions. The complementary Elliptic K\('\) function can be used as EllipticK!\('\)(k). Note that \(i\mathrm {K}'(k)\) is the other quarter period of the Jacobi elliptic functions.

Elliptic E

The Elliptic E function comes with either one or two arguments; used with two arguments as EllipticE(u,k) it will return the value of Legendre’s form of the Incomplete Elliptic Integral of the Second Kind: \[\mathrm {E}(\phi , k)=\int _0^\phi \sqrt {1-k^2 \sin ^2 \theta } \,\mathrm {d}\theta .\] When called with one argument EllipticE(k) will return the value of the Complete Elliptic Integral of the Second Kind: \[\mathrm {E}(k)=\mathrm {E}(\pi /2, k) = \int _0^{\pi /2} \sqrt {1-k^2 \sin ^2 \theta } \,\mathrm {d}\theta .\] The complementary Elliptic E\('\) function can be used as EllipticE!\('\)(k).

The complete integrals are actually multi-valued; to obtain single valued functions it is necessary to introduce branch cuts. For \(\mathrm {K}(k)\) and \(\mathrm {E}(k)\) these are \((-\infty , -1] \cup [1, +\infty )\). The functions are continuous if the two components of the branch cut are approached from the second and fourth quadrants respectively. For \(\mathrm {K}'(k)\) and \(\mathrm {E}'(k)\) the branch cut is \((-\infty , 0]\) with continuity if the cut is approached from the second quadrant. For more details, see Lawden [Law89] sections 8.12 to 8.14. The principal values of \(\mathrm {K}(k)\) and \(\mathrm {E}(k)\) are even functions of \(k\).

The numerical evaluation of the complete integrals is more robust and uses symmetric elliptic integrals. They should now work for all complex values of the modulus and return the principal value of the integral concerned. Note that for all complex values of \(k\), the well-known identities: \[\mathrm {K}'(k)=\mathrm {K}(\sqrt {1-k^2})\qquad \mathrm {E}'(k) = \mathrm {E}(\sqrt {1-k^2}).\] are not actually valid for the principal values of the functions concerned. It is necessary that \(\Re k >=0\) with in addition if \(\Re k=0\), \(\Im k>=0\). One of the consequences of this is that the principal values of \(\mathrm {K}'(k)\) and \(\mathrm {E}'(k)\) are not even functions.

General values for the four complete integrals together with the principal value when \(\Re {k} <0\) are given in the table below, where \(n\) is an arbitrary integer, \(k' = \sqrt {1-k^2}\) and the expressions in the second and third columns refer always to principal values:

In the third column the upper and lower alternative signs are taken when \(k\) lies in the second and third quadrants respectively.

One quite subtle point arises: although the twelve Jacobi elliptic functions and \(\mathrm {K}(k)\) are even functions of the modulus \(k\), the quarter period \(i\mathrm {K}'(k)\) is not (see the second row of the table above). If \(\Re (k)>0\), then for example \(4\mathrm {K}(k)\) and \(2i\mathrm {K}'(k)=2i\mathrm {K}(k')\) are primitive periods of \(\mathrm {sn}(x,k)\). However, if we use \(-k\) as the modulus, the first primitive period obtained is the same (since \(\mathrm {K}(k)\) is an even function) but the second is different namely: \(2i\mathrm {K}'(-k)=2i\mathrm {K}(k') \pm 4\mathrm {K}(k)\). This explains why some of the values returned by the ‘eighth’ period rules for \(\mathrm {sn}\) etc. are not even functions of k.

Elliptic D

The Elliptic D function also comes with either one or two arguments; used with two arguments as EllipticD(u,k) it will return the value of an alternative form of Legendre’s Incomplete Elliptic Integral of the Second Kind: \[\mathrm {D}(\phi , k)=\int _0^\phi \frac {\sin ^2 \theta }{\sqrt {1-k^2\sin ^2\theta }} \,\mathrm {d}\theta .\] When called with one argument EllipticD(k) will return the value of the Complete Elliptic Integral of the Second Kind: \[\mathrm {D}(k)=\mathrm {D}(\pi /2, k) = \int _0^{\pi /2} \frac {\sin ^2 theta}{\sqrt {1-k^2 \sin ^2 \theta }} \,\mathrm {d}\theta .\] The integrals of the first and second kind are related: \[\mathrm {E}(k) = \mathrm {K}(k)-k^2\mathrm {D}(k), \qquad \mathrm {E}(\phi , k) = \mathrm {F}(\phi , k)-k^2\mathrm {D}(\phi , k).\]

For all the elliptic integrals, rule lists are provided for simplification for special values of the argument(s), differentiation with respect to either argument and shift rules for the incomplete integrals so that their first argument lies in the range \(0 <= \phi <= \pi /2\).

If the arguments of the elliptic integral function are numeric and the switches ROUNDED and COMPLEX are both ON, the numerical value will be returned. The numerical evaluation routines should now succeed for all complex values of the arguments provided, of course, the corresponding integral exists. The incomplete integrals \(\mathrm {F}(\phi , k)\), \(\mathrm {E}(\phi , k)\), \(\mathrm {D}(\phi , k)\) and \(\Pi (\phi , a, k)\) are all multi-valued as there are branch points in the defining integrands when \(\sin (\phi ) =\pm 1\) and \(k\sin (\phi ) = \pm 1\).

For mumerical evaluation when the first argument \(\phi \) is complex there are issues that arise which are absent when \(\phi \) is real. One might use the the straight line contour from \(\theta =0\) to \(\theta =\phi \), but this is leads to integrals that are difficult to evaluate. When \(|\Re {\phi }| \leq \pi /2\), the change of variable \(x=\sin \theta \) is applied and the value returned uses the straight line contour between \(x=0\) and \(x=\sin (\phi )\). This will produce a different value if the loop formed by the two contours encloses a branch point. When \(|\Re {\phi }| > \pi /2\) the contour used is the straight line segment along the real axis from 0 to the multiple of \(\pi \) nearest to \(\phi \) followed by the straight line segment from between \(x=0\) and \(x=\sin (\phi )\). When \(\phi \) is real the contours in the \(\theta \) and \(x\) planes coincide. Other contours between the two endpoints will yield different values depending on how the contour winds around the two branch points.

Note also that if the contour contains a branch point, there is sign ambiguity as the branch point is crossed depending on which branch of the square root in the integrand is chosen. It seems logical to choose the principal branch.

For example when \(k>1\) & \(1/k < y <1\) \[ \int _0^y \frac {\mathrm {d}x}{\sqrt {1-x^2}\sqrt {1-k^2x^2}}\] is taken to be \[\int _0^{1/k} \frac {\mathrm {d}x}{\sqrt {1-x^2}\sqrt {1-k^2x^2}} -i\int _{1/k}^y \frac {\mathrm {d}x}{\sqrt {1-x^2}\sqrt {k^2x^2-1}}.\]

Elliptic \(\Pi \)

The Elliptic \(\Pi \) function can be used as EllipticPi( ) and will return the value of the Elliptic Integral of the Third Kind.

The Elliptic \(\Pi \) function comes with either two or three arguments; when called with three arguments EllipticPi(phi,a,k) will return the value of the Incomplete Elliptic Integral of the Third Kind: \[\Pi (\phi ,a,k)=\int _0^{\phi } \left ((1-a^2\sin ^2 \theta )\sqrt {1-k^2 \sin ^2 \theta }\right )^{-1} \,\mathrm {d}\theta .\] For \(-\pi /2 \leq \Re {\phi } \leq \pi /2\), the incomplete integral may be written in the form: \[\Pi (\phi ,a,k)=\int _0^{sin \phi } \left ((1-a^2t^2)\sqrt {(1-t^2)(1-k^2t^2)}\right )^{-1}\,\mathrm {d}t.\] When called with two arguments as EllipticPi(a,k) it will return the value of the Complete Elliptic Integral of the Third Kind:

For certain values of \(a\) namely \(0\), \(\pm 1\) and \(\pm k\) the integrals reduce to elliptic integrals of the first and second kinds.

20.19.5 Jacobi’s Elliptic Integrals

Jacobi E

The Jacobi E function can be used as jacobiE(u,k); it will return the value of Jacobi’s form of the Incomplete Elliptic Integral of the Second Kind: \[\mathcal {E}(u, k)=\int _0^u \mathrm {dn}^2 (v, k) \,\mathrm {d}v.\]

The relationship between the two forms of incomplete elliptic integrals can be expressed as \[\mathcal {E}(u, k) = \mathrm {E}(\mathrm {am}(u), k).\] Note that \[\mathrm {E}(k)=\mathcal {E}(\mathrm {K}(k), k) =\int _0^{\mathrm {K}(k)} \mathrm {dn}^2(v, k) \,\mathrm {d}v.\]

On a GUI that supports calligraphic characters (NB. this is now the case with the CSL GUI), there is no problem and it is rendered as \(\mathcal {E}(u,k)\) in accordance with NIST usage. On non-GUI interfaces the Jacobi E function is rendered as E_j.

Jacobi’s Zeta Function

This can be called as jacobiZeta(u,k) and refers to Jacobi’s (elliptic) Zeta function \(\mathrm {Z}(u,k)\) whereas the operator Zeta will invoke Riemann’s \(\zeta \) function. It is closely related to Jacobi’s Epsilon function jacobie; in fact \[\mathrm {Z}(u,k) = \mathcal {E}(u,k)-u\mathrm {E}(k)/\mathrm {K}(k).\]

20.19.6 Symmetric Elliptic Integrals

The functions in this section are currently only intended for use for in the numerical evaluation of the ‘classical’ elliptic integrals discussed above and in certain other ‘basic’ elliptic integrals. The switches ROUNDED and COMPLEX should both be ON.

Symmetric elliptic integrals are a relatively new development in the theory of elliptic functions mainly due to Carlson and his collaborators; an extensive bibliography is available on the NIST Digital Library of Mathematical Functions starting at the link DLMF:NIST, Bibliography. They have a number of advantages over more traditional approaches; see for example Advantages of Symmetry.

The fundamental symmetric integrals are all integrals over the positive real line and involve the function \(s(t)=\sqrt {t+x}\sqrt {t+y}\sqrt {t+z}\) where \(x,y,z\in \mathbb {C} \setminus (-\infty , 0]\) but, except where otherwise stated, at most one of \(x,y,z\) may be zero.

The first three integrals defined above are symmetric in \(x,y,z\) and are termed symmetric elliptic integrals of the first, second and third kinds respectively. \(\mathrm {RD}\) and \(\mathrm {RC}\) are degnerate versions of \(\mathrm {RJ}\) and \(\mathrm {RF}\) respectively as \[\mathrm {RD}(x,y,z) =\mathrm {RJ}(x,y,z,z) \qquad \mathrm {RC}(x,y) =\mathrm {RF}(x,y,y) \] \(\mathrm {RD}\) is an elliptic integral of the second kind and is only symmetric in \(x,y\) whilst \(\mathrm {RC}\) is an elementary integral, but can be numerically evaluated efficiently without needing to distinguish between the trigonometric and hyperbolic cases. For certain special values of \(p\) the integral \(\mathrm {RJ}\) of the third kind degenerates into integrals of the first and second kinds. For numeric evaluations it turns out that the use of \(\mathrm {RD}\) is more convenient than that of \(\mathrm {RG}\) (which is not currently implemented in REDUCE). For more information see §19.15-19.29 of the DLMF:NIST.

The Legendre elliptic integrals may all be expressed in terms of the symmetric integrals; the complete integrals satisfy

where here and below \(k^\prime =\sqrt {1-k^2}\). For the incomplete integrals, defining \(s=\sin \phi \) with \(-\pi /2 \leq \Re {\phi } \leq \pi /2\), we have

The above formulae are only valid when the range of integration does not include one or more branch points of the integrand (when \(s\) is real and \(s^2>1\) or when \(k s\) is real and \((k s)^2 >1\)). In these cases it is necessary to split the range of integration at the branch points resulting in two or three elliptic integrals.

Currently the Carlson’s duplication method is used to evaluate the symmetric elliptic integrals of the first, second and third kinds \(\mathrm {RF}\), \(\mathrm {RD}\) and \(\mathrm {RJ}\) (and also the related elementary integral \(\mathrm {RC}\)). For more details see the DLMF website: Duplication Method. In particular the REDUCE code is essentially a translation from Fortran of the code of Carlson & Notis [CN81] generalised for complex arguments and structured to avoid GO TO.

Alternatively the symmetric integral \(\mathrm {RF}\) may be evaluated using a sequence of quadratic transformations which converge rapidly to the elementary hyperbolic integral: \[\mathrm {RC}(X^2+Y^2,X^2) = \mathrm {RF}(X^2,X^2,X^2+Y^2) = \mathrm {arctanh} (Y/(X^2+Y^2))/Y.\] More information on this method due to Carlson in the 1990’s may be found on the DLMF website: Quadratic Transformations.

Basic Elliptic Integrals

Elliptic integrals involving the square root \(s(t)\) of a cubic or quartic polynomial in which the range of integration is an arbitrary interval \([y,x]\) of the real line are considered. Again there are three basic kinds: first, second and third. Carlson [Car88] & [Car99] showed how each of these can be expressed as a fundamental symmetric elliptic integral of the corresponding kind over the range \([0,\infty )\). Let

and where the four line segments with end points \(a_\alpha +b_\alpha y\) and \(a_\alpha +b_\alpha x\) for \(1\leq \alpha \leq 4\) lie in \( \mathbb {C} \setminus (-\infty , 0)\). Note that if \(s(t)\) is a cubic then one simply chooses supplies unity i.e. \(a_\alpha =1, b_\alpha =0\) as a fourth factor of \(s(t)\).

Then integrals of the first and second kinds take the form:

where in the second equation we have chosen (wlog) the distinguished factors to have indices 1 & 4.

For any permutation \(\alpha ,\beta ,\gamma ,\delta \) of 1,2,3,4, let

Clearly \(U_{\alpha \beta }=U_{\beta \alpha }=U_{\gamma \delta }=U_{\delta \gamma }\) and at most one of \(U_{12}, U_{13}, U_{23}\) is zero since \(U_{\alpha \beta }^2 - U_{\alpha \gamma }^2 = d_{\alpha \delta }d_{\beta \gamma } \neq 0\), thus at most one of the parameters of \(\mathrm {RF}\) and \(\mathrm {RD}\) is zero. If \(X_4\) or \(Y_4\) is zero,that is if \(a_4+b_4t\) vanishes at one end of the range of integration, the integral of the second kind diverges.

In the ‘awkward’ case when \(U_{23}=0\) or \(U_{12}^2+U_{13}^2=0\), the above method breaks down. However it is still possible to calculate the required integral; choose \(\beta =2\) or \(\beta =3\) so that \(b_\beta \neq 0\) and note that

The second term on the lhs of the second equation reduces to \(b_\beta \int 1/s(t)\mathrm {d}t\) and so is an integral of the first kind. Thus the required integral can be expressed in terms of an integral of the second kind and one of the first kind.

Integrals of the third kind take the form \[\int _y^x\frac {(a_1+b_1t)\mathrm {d}t}{(a_5+b_5t)s(t)} = \frac {2d_{12}d_{13} d_{14}}{3d_{15}}\mathrm {RJ}(U_{12}^2, U_{13}^2, U_{23}^2,U_{15}^2) + \mathrm {RC}(S_{15}^2,Q_{15}^2)\] where

where \(\beta ,\gamma ,\delta \) is a permutation of 2,3,4. Again at most one of the parameters of \(\mathrm {RJ}\) is zero. This method will break down when calculating \(S_{15}\) if \(a_1+b_1t\) vanishes at either the upper or lower limit of integration as either \(X_1\) or \(Y_1\) is zero. However the situation can be remedied by choosing \(\beta \) so that \(X_\beta \neq 0\), \(Y_\beta \neq 0\) and \(b_\beta \neq 0\) and using similar methods to that used for the awkward case for integrals of the second kind.

Note if the integration range is \((-\infty , +\infty )\), the above formulae for integrals of all three kinds are not valid and the integral needs to be split into two at, say, zero. Similarly if the integration range is such that the integrand has branch points, the integration range will need to be split into two at each of these branch points.

In REDUCE elliptic integrals of the three kinds may be evaluated numerically when the switches ROUNDED & COMPLEX are both ON by calling the functions ellint_1st, ellint_2nd and ellint_3rd. The first two parameters are the lower and upper limits of the integration range; these are followed by 4 (or 5 for ellint_3rd) two-element lists \(\{a_1,b_1\}, \{a_2,b_2\} \ldots \).

Known bug: As pointed out by Carlson & FitzSimmons [CF00], the algorithms for \(\mathrm {RF}\) & \(\mathrm {RJ}\) may break down when \(s(t)\) is the product of two quadratic factors each with complex conjugate zeros \(x_1\pm i y_1\) & \(x_2\pm i y_2\). One of the three quantities \(U_{12}, U_{13}, U_{23}\) may be negative and this causes the algorithm to return an incorrect result. This error only arises when the crossing point (neccessarily real) of the diagonals of the parallellogram in the complex plane formed by the four zeros of \(s(t)\) lies inside the range of integration. The algorithm also appears to produce the correct result when the crossing point is inside the range of integration but is ‘sufficiently close’ to either end of the range – the precise condition appears to be unknown. In the same paper Carlson & FitzSimmons propose a modified algorithm of the same ilk to overcome these difficulties, however it is yet to be implemented in REDUCE.

Some Examples

In this subsection some of the advantages of symmetry are illustrated; the formula for the integral of the first kind replaces the 28 formulae in Gradshteyn and Ryzhik. For the cubic case one of the factors in \(s(t)\) is simply taken to be unity whilst for cases of the form \[\int _\alpha ^\beta \frac {\mathrm {d}x}{\sqrt {(a+b x^2)(c+d x^2)}}\] one must use the substitution \(t=x^2\) to obtain \[\frac {1}{2}\int _{\alpha ^2}^{\beta ^2} \frac {\mathrm {d}t}{\sqrt {t(a+b t)(c+d t)}}\] Moreover the restriction that one limit of integration be a branch point of the integrand is eliminated without doubling the number of standard integrals in the result. Similarly the formula for the integral of the second kind replaces no less than 144 cases in Gradshteyn and Ryzhik (72 each for the quartic and cubic cases). For cases where the numerator in the integrand is unity one simply takes \(a_1=0, b_1=1\) and if there is to be no apparent term of the form \((a_4+b_4t)^{3/2}\) in its denominator one takes \(a_4=0, b_4=1\).

When \(0 \leq \phi \leq \pi /2\), the fundamental integrals of the first, second and third kinds \(\mathrm {F}(\phi ,k)\), \(\mathrm {E}(\phi ,k)\), \(\mathrm {D}(\phi ,k)\) and \(\Pi (\phi ,a,k)\) may be written as

where in all four cases use has been made of the substitution \(t=\sin ^2 \theta \) and \(s=\sin \phi \). In the last three equations one takes \(a_4=1, b_4 =0\) so that the fourth factor in \(s(t)\) is unity and in the fourth equation one also takes \(a_1=1,b_1=0\).

20.19.7 Some Numerical Utility Functions

Five utility functions are provided:

• nome2mod(q)

• nome2mod!\('\)(q)

• nome2!K(q)

• nome2!K!\('\)(q)

• nome(k)

These are only operative when the switch rounded is on and their argument is numerical. The first pair relate the nome \(q\) of the theta functions with the moduli \(k\) and \(k'=\sqrt {1-k^2}\) of the associated Jacobi elliptic functions.

The second pair return the quarter periods K and K\('\) respectively of the Jacobi elliptic functions associated with the nome \(q\).

Finally, nome(k) returns the nome \(q\) associated with the modulus \(k\) of a Jacobi elliptic function and is essentially the inverse of nome2mod.

20.19.8 Jacobi Theta Functions

These theta functions differ from those originally defined by Lisa Temme in a number of respects. Firstly four separate functions of two arguments are defined:

• elliptictheta1(u,tau)\(\qquad \vartheta _1(u, \tau )\)

• elliptictheta2(u,tau)\(\qquad \vartheta _2(u, \tau )\)

• elliptictheta3(u,tau)\(\qquad \vartheta _3(u, \tau )\)

• ellipticthetas(u,tau)\(\qquad \vartheta _4(u, \tau )\)

rather than a single function with three arguments (with the first argument taking integer values in the range 1 to 4). Secondly the periods are \(2\pi , 2\pi , \pi \) and \(\pi \) respectively (NOT 4K, 4K, 2K and 2K). Thirdly the second argument is the modulus \(\tau = a+i b\) where \(b=\Im \tau >0\) and hence the quasi-period is \(\pi \tau \).

The second parameter was previously the nome \(q\) where \(|q|<1\). As a consequence elliptictheta1 and elliptictheta2 were multi-valued owing to the appearance of \(q^{1/4}\) in their defining expansions. elliptictheta3 and elliptictheta4 were, however, single-valued functions of \(q\).

Regarded as functions of \(\tau \), elliptictheta1 and elliptictheta2 are single-valued functions. The nome is given by \(q = \exp (i\pi \tau )\) so that the condition \(\Im (\tau )>0\) ensures that \(|q| < 1\). Note also in this case \(q^{1/4}\) is interpreted as \(\exp (i\pi \tau /4)\) rather than the principal value of \(q^{1/4}\). Thus, \(\tau \), \(2+\tau \), \(4+\tau \) and \(6+\tau \) produce four different values of both elliptictheta1 and elliptictheta2 although they all correspond to the same nome \(q\).

The four theta functions are defined by their Fourier series: \begin {align*} \vartheta _1(z,\tau ) & = 2 e^{i\pi \tau /4}\sum _{n=0}^\infty (-1)^nq^{n^2+n} \sin (2n+1)z\\ \vartheta _2(z,\tau ) & = 2 e^{i\pi \tau /4}\sum _{n=0}^\infty q^{n^2+n} \cos (2n+1)z\\ \vartheta _3(z,\tau ) & = 1 +2\sum _{n=1}^\infty q^{n^2} \cos 2n z\\ \vartheta _4(z,\tau ) & = 1 +2\sum _{n=1}^\infty (-1)^n q^{n^2} \cos 2n z. \end {align*}

Utilising the periodicity and quasi-periodicity of the theta functions some generalised shift rules are implemented to shift their first argument into the base period parallelogram with vertices \[(\pi /2, \pi \tau /2),\quad (-\pi /2, \pi \tau /2),\quad (-\pi /2, -\pi \tau /2), \quad (\pi /2, -\pi \tau /2).\] Together with the relation \(\vartheta _1(0,\tau )=0\), these shift rules serve to simplify all four theta functions to zero when appropriate.

When the switches rounded and complex are on and the arguments are purely numerical and the imaginary part of \(\tau \) positive, the theta functions are evaluated numerically. Note that as \(\tau \) is necessarily complex, the switch complexmust be on.

In what follows \(a\) and \(b\) will denote the real and imaginary parts of \(\tau \) respectively and so \(|q| = \exp (-\pi b)\) and \(\arg q =\pi a\). The series for the theta functions are fairly rapidly convergent due to the quadratic growth of the exponents of the nome \(q\) – except for values of \(q\) for which \(|q|\) is near to 1 (i.e. \(b=\Im \tau \) close to zero). In such cases the direct algorithm would suffer from slow convergence and rounding errors. For such values of \(|q|\), Jacobi’s transformation \(\tau '=-1/\tau \) can be used to produce a smaller value of the nome and so increase the rate of convergence. This works very well for real values of \(q\), or equivalently for \(\tau \) purely imaginary since \(q'= q^{1/b^2}\), but for complex values the gains are somewhat smaller. The Jacobi transformation produces a nome \(q'\) for which \(|q'| = |q|^{1/(a^2+b^2)}\).

When \(\Re q < 0\), the Jacobi transformation is preceded by either the modular transformation \(\tau ' = \tau +1\) when \(\Im q < 0\), or \(\tau ' = \tau -1\) when \(\Im q > 0\), which both have the effect of multiplying \(q\) by \(-1\), so that the new nome has a non-negative real part and \(|a| \leq 1/2\). Thus the worst case occurs for values of the nome \(q\) near to \(\pm i\) where \(|q'| \approx |q|^4\).

By using a series of Jacobi transformations preceded, if necessary by \(\tau \)-shifts to ensure \(|a| <= 1/2\), \(|q|\) may be reduced to an acceptable level. Somewhat arbitrarily these Jacobi’s transformations are used until \(b > 0.6\) (i.e. \(|q| < 0.15\)). This seems to produce reasonable behaviour. In practice more than two applications of Jacobi transformations are rarely necessary.

The previous version of the numerical code returned the principal values of \(\vartheta _1\) and \(\vartheta _2\), that is the ones obtained by taking the principal value of \(q^{1/4}\) in their series expansions. The current version replaces \(q^{1/4}\) by \(\exp (i\pi \tau /4)\). If the principal value is required, it is easily obtained by multiplying by the ‘correcting’ factor \(q^{1/4}\exp (-i\pi \tau /4)\).

Derivatives of Theta Functions

Four functions are provided:

• theta1d(u,ord,tau)

• theta2d(u,ord,tau)

• theta3d(u,ord,tau)

• theta4d(u,ord,tau)

These return the \(d\)th derivatives of the respective theta functions with respect to their first argument \(u\); \(\tau \) is as usual the modulus of the theta function. These functions are only operative when the switches rounded and complex are ON and their arguments are numeric with \(d\) being a positive integer. They are provided mainly to support the implementation the Weierstrassian and Sigma functions discussed in the following subsection.

The numeric code simply sums the Fourier series for the required derivatives. Unlike the theta functions themselves the code does not use the quasi-periodicity nor modular transformations to speed up the convergence of the series by reducing the sizes of \(\Im u\) and \(|q|\). In the numerical evaluation of the Weierstrassian and Sigma functions these functions are only called after the necessary shifts of the argument \(u\) and modular transformations of \(\tau \) have been performed. These are much simpler in this context.

Nevertheless they may be used from top level and numerical experiments reveal that the rounding errors are not significant provided \(|q|\) is not near one (say \(|q|<0.9\)) and \(u\) is real or at least has a relatively small imaginary part.

20.19.9 Weierstrass Elliptic & Sigma Functions

Three main functions of three arguments are defined:

• \(\wp (u, \omega _1, \omega _3)\)  —  weierstrass(u,omega1,omega3)

• \(\zeta _w(u, \omega _1, \omega _3)\) —  weierstrassZeta(u,omega1,omega3)

• \(\sigma (u, \omega _1, \omega _3)\) —  weierstrass_sigma(u,omega1,omega3)

The notation used is broadly similar used by Lawden [Law89] which is also used in the NIST Digital Library of Mathematical Functions DLMF:NIST. However, \(\zeta _w\) is used for the Weierstrassian Zeta function to distinguish it from the Riemann Zeta function and the usual symbol \(\wp \) is used for the Weierstrassian elliptic function itself.

The two primitive periods of the Weierstrass function are \(2\omega _1\) and \(2\omega _3\) and these must satisfy \(\Im (\omega _3/\omega _1) \neq 0\). The two periods are normally numbered so that \(\tau = \omega _3/\omega _1\) has a positive imaginary part and hence the nome \(q = exp(i\pi \tau )\) satisfies \(|q| <1\).

Any linear combination \(\Omega _{m,n} = 2m\omega _1 +2n\omega _3\) where \(m\) and \(n\) are integers (not both zero) is also a period. The set of all such periods plus the origin form a lattice. In the literature \(-(\omega _1+\omega _3)\) is often denoted by \(\omega _2\) and \(2\omega _2\) is clearly also a period; this accounts for the gap in the numbering of primitive periods. The period \(\omega _2\) is not used in REDUCE the rule sets for the Weierstrassian elliptic and related functions.

The primitive periods are not unique; indeed any periods \(2\Omega _1\) and \(2\Omega _3\) defined by the unimodular integer bilinear transformation: \[\Omega _1 = a\omega _1 + b\omega _3,\qquad \Omega _3 = c\omega _1 + d\omega _3,\qquad \text { where }ad-bc = 1\] are also primitive. This fact is very useful in the numerical evaluation of the Weierstrassian and Sigma functions as a sequence of such transformations may be used to increase the size \(\Im \tau \) and so reduce the size of \(|q|\). Thus the Fourier series for the theta functions and their derivatives will converge rapidly. In theory these transformations may be used to reduce the size of \(|q|\) until \(\Im \tau \geq \sqrt 3/2\) when \(|q|<0.06\). However, in numerical evaluations in REDUCE it is sufficient to use these transformations only until \(\Im \tau > 0.7\), i.e. until \(|q| < 0.11\). In practice only two or three iterations are required and usually very much smaller values of \(|q|\) are achieved particularly when \(\tau \) is purely imaginary i.e. \(q\) is real.

In the numerical evaluations, if a result is real (or purely imaginary) it may happen that the result returned has a very small imaginary part (resp. real part). The ratio of the ‘deliquent’ part to the actual result is invariably smaller than current PRECISION and is due to rounding. Similarly if the true result is actually zero the result returned may have a very small absolute value – again smaller than the current PRECISION.

The Weierstrassian function is even and has a pole of order 2 at all lattice points. The Zeta and Sigma functions are only quasi-periodic on the lattice. Zeta is odd and has simple poles of residue 1 at all lattice points. The basic Sigma function \(\sigma (u,\omega _1,\omega _3)\) is odd and regular everywhere as is the function \(\vartheta _1(u,\tau )\) to which it is closely related. It has zeros at all lattice points. All three functions \(\wp \), \(\zeta _w\) and \(\sigma \) are homogenous of degrees -2, -1 and +1 respectively. The functions are related by \[ \wp (u) = -\zeta _w'(u),\qquad \qquad \zeta _w(u) = \sigma '(u)/\sigma (u),\] where the lattice parameters have been omitted for conciseness.

Rule sets are provided which implement all the properties such as double periodicity discussed above. For numerical evaluation the switches rounded and complex must both be ON and all three parameters must be numeric. It is not, however, necessary to ensure \(\Im (\omega _3/\omega _1) >0\) as the second and third parameters will be swapped if required.

Alternative forms of the Weierstrass Functions

Two commonly used alternative forms of the Weierstrassian functions in which they are regarded as functions of the lattice invariants \(g_2\) and \(g_3\) rather than the primitive periods \(\omega _1\) and \(\omega _3\) are provided:

• \(\wp (u \mid g_2, g_3)\)  —  weierstrass1(u,g2,g3)

• \(\zeta _w(u \mid g_2, g_3)\) —  weierstrassZeta1(u,g2,g3).

Note that for output they are distinguished from the two discussed above by separating the first and second arguments by a vertical bar rather than a comma. The rule for differentiation of the Weierstrass function is simpler when expressed in this alternative: \[ \wp '(u \mid g_2,g_3)^2 = 4\,\wp (u \mid g_2,g_3)^3 - g_2\, \wp (u \mid g_2,g_3) -g_3. \]

Other Sigma Functions

Three further Sigma functions are also provided:

• \(\sigma _1(u, \omega _1, \omega _3)\)  —  weierstrass_sigma1(u,omega1,omega3)

• \(\sigma _2(u, \omega _1, \omega _3)\) —  weierstrass_sigma2(u,omega1,omega3)

• \(\sigma _3(u, \omega _1, \omega _3)\) —  weierstrass_sigma3(u,omega1,omega3)

These are all even functions, regular everywhere, homogenous of degree zero and doubly quasi-periodic. They are closely related to the theta functions \(\vartheta _2\), \(\vartheta _3\) and \(\vartheta _4\) respectively; but note the difference in numbering. For more information on the properties these sigma functions, see Lawden [Law89]; they do not appear in the NIST Digital Library of Mathematical Functions, but are included here for completeness.

Quasi-Period Factors & Lattice Functions

Ten functions are provided:

• \(e_1(\omega _1, \omega _3)\)  —  lattice_e1(omega1, omega3);

• \(e_2(\omega _1, \omega _3)\)  —  lattice_e2(omega1, omega3);

• \(e_3(\omega _1, \omega _3)\)  —  lattice_e3(omega1, omega3);

• \(g_2(\omega _1, \omega _3)\)  —  lattice_g2(omega1, omega3);

• \(g_3(\omega _1, \omega _3)\)  —  lattice_g3(omega1, omega3);

• \(\Delta (\omega _1, \omega _3)\)  —  lattice_delta(omega1, omega3);

• \(\mathrm {G}(\omega _1, \omega _3)\)  —  lattice_g(omega1, omega3);

• \(\eta _1(\omega _1, \omega _3)\)  —  eta_1(omega1, omega3);

• \(\eta _2(\omega _1, \omega _3)\)  —  eta_2(omega1, omega3);

• \(\eta _3(\omega _1, \omega _3)\)  —  eta_3(omega1, omega3).

These are operative when the switches rounded and complex are ON and their arguments are numerical. The first three are referred to as lattice roots and are related to the invariants \(g_2, g_3\), the discriminant \(\Delta = g_2^3-27g_3^2\) and a closely related invariant \(\mathrm {G} = g_2^3/(27 g_3^2)\) of the Weierstrassian elliptic function \(\wp \). The lattice roots also appear in the numerical evaluation of the Weierstrass function. These lattice roots satisfy: \[e_1+e_2+e_3=0,\qquad g_2=2(e_1^2+e_2^2+e_3^2),\qquad g_3= 4e_1e_2e_3.\] If the discriminant \(\Delta \) vanishes or equivalently if \(\mathrm {G} = 1\), there are at most two distinct lattice roots and the elliptic function degenerates to an elementary one. The advantage of the invariant \(\mathrm {G}\) is that it is a function of \(\tau = \omega _3/\omega _1\) only.

The remaining three functions eta_1, eta_2 & eta_3 appear in the rules for the quasi-periodicity of the four sigma functions and of the Weierstrassian Zeta function. They are also used in the numerical evaluation of these functions when the switches rounded and complex are ON. The quasi-period relations are: \begin {align*} \zeta _w(u+2\omega _j) & = \zeta _w(u)+2\eta _j\\ \sigma (u+2\omega _j) & = -exp(2\eta _j(u+\omega _j))\sigma (u)\\ \sigma _k(u+2\omega _j) & = exp(2\eta _j(u+\omega _j))\sigma _k(u) \quad \text { if }j\neq k\\ \sigma _j(u+2\omega _j) & = -exp(2\eta _j(u+\omega _j))\sigma _j(u)\\ \zeta _w(\omega _j) & = \eta _j\\ \sigma _j(\omega _j) & = 0, \end {align*}

where the lattice parameters have been omitted for conciseness and \(j,k = 1\ldots 3\). The quasi-period factors satisfy \[\eta _1+\eta _2+\eta _3=0,\qquad \eta _1\omega _3-\eta _3\omega _1=\eta _2\omega _1-\eta _1\omega _2=\eta _3\omega _2-\eta _2\omega _3=i\pi /2.\] As well as the scalar-valued functions discussed above in this section, there are four functions which return a list as their value:

• lattice_roots(omega1, omega3) — returns \(\{e_1,\ e_2,\ e_3\}\);

• lattice_invariants(omega1, omega3) — returns \(\{g_2,\ g_3,\ \Delta ,\ \mathrm {G}\}\);

• quasi_period_factors(omega1, omega3) — returns \(\{\eta _1,\ \eta _2,\ \eta _3\}\);

• lattice_generators(g2, g3) — returns \(\{\omega _1,\ \omega _3\}\).

The first three are actually more efficient than calling the requisite scalar-valued functions individually and the fourth is used in the numerical evaluation of the Weierstrass functions regarded as functions of the invariants. These functions are only useful when the switches rounded and complex are ON and their arguments are all numerical. Note that the call sequence:

  lattice_generators(g2,g3);
  lattice_invariants(first ws, second ws);
  {first ws, second ws};

should reproduce the list {g2, g3}, perhaps with small rounding errors. The corresponding sequence with the function lattice_generators being called after lattice_invariants (and g2 & g3 replaced by w1 & w3), in general, will not produce the same pair of lattice generators since the generators are only defined up to a unimodular bilinear transformation.

For details of the algorithm used to calculate the lattice generators from the invariants see the DLMF:NIST chapter on Lattice Calculations.

20.19.10 Inverse Jacobi Elliptic Functions

The following inverses of the 12 Jacobi elliptic functions are available:-

• arcsn(u,k)

• arcdn(u,k)

• arccn(u,k)

• arccd(u,k)

• arcsd(u,k)

• arcnd(u,k)

• arcdc(u,k)

• arcnc(u,k)

• arcsc(u,k)

• arcns(u,k)

• arcds(u,k)

• arccs(u,k)

Thus, for example,

   jacobisn(arcsn(x, k), k)   --> x
   jacobisc(arcsc(x, k), k)   --> x

A rule list is provided to simplify these functions for special values of their arguments such \(x=0\), \(k=0\) and \(k=1\), to implement the inverse function simplification formulae illustrated immediately above and for differentiation of these functions with respect to their two arguments.

Note that \(\mathrm {arccs}\) is not defined to be an odd function of its first argument unlike \(\mathrm {cs}\). Instead it is taken to satisfy: \[ \mathrm {arccs}(-x, k) = 2\mathrm {K}(k)-\mathrm {arccs}(x, k).\] This is analogous to the situation in Reduce for \(\mathrm {acot}\) where \[ \mathrm {arctan}(-x) = -\mathrm {arctan}(x),\qquad \mathrm {arccot}(-x) = \pi -\mathrm {arccot}(x). \]

This choice means that the range of (real) principal values of \(\mathrm {arccs}\) is connected – it is the open set \((0, 2K(k))\).

When their arguments are numerical, these functions will be evaluated numerically if the rounded and complex switches are both ON. Note that in some cases the result may have an imaginary part even if both arguments are real, hence the necessity of the switch complex being ON.

Note also that for \(\mathrm {arcdn}\) and \(\mathrm {arcnd}\) a zero value of the modulus \(k\) is excluded (since \(\mathrm {dn}(x,0) = \mathrm {nd}(x,0) = 1 \quad \forall x\)).

As the Jacobi elliptic functions are doubly periodic, their inverse functions are multi-valued. The numerical value returned is the principal value \(v\) which lies in the parallelogram in the complex plane whose vertices are given in the table below. Other values of the inverse functions are indicated in the fifth column of the table below where \(m\) and \(n\) are arbitrary integers.

When both arguments are real and \(|k|<=1\) and when there are certain restrictions on the range of the first parameter \(x\) (see the table below), then the principal value of the inverse function is real. It lies in the range given in the third column of the table below. (c.f. the inverse trigonometric functions). Other real values of the inverse functions are indicated in the fourth column of the table below where \(n\) is an arbitrary integer.

The numerical values of the inverse functions are calculated by expressing them in terms of the symmetric elliptic integral: \[ R_F(x,y,z)=\int _0^\infty 1/\sqrt {(t-x)(t-y)(t-z)}\,\mathrm {d}t. \] For more details see the DLMF website: Inverse Jacobian Elliptic Functions.

20.19.11 Table of Elliptic Functions and Integrals

Hosted by