ELLIPFN: A Package for Elliptic Functions and Integrals

Additional documentation and examples can be found in the ﬁles 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

The package ELLIPFN is designed to provide algebraic and numeric manipulations of many elliptic functions, namely:

20.19.2 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 in 2021 to support Weierstrassian Elliptic functions, Sigma functions and inverse Jacobi elliptic functions.

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 and derivatives of the theta functions.

20.19.3 Jacobi Elliptic Functions

These diﬀer 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

Extensive rule lists are provided for diﬀerentiation of these functions with respect to either argument, to implement the standard addition formulae, argument shifts by multiples of the two quarter-periods K and iK′ and ﬁnally Jacobi’s transformations for a purely imaginary ﬁrst argument.

Four utility rule lists to_sn, to_cn, to_dn and no_glaisher are provided for user convenience. The ﬁrst to_sn replaces occurences of the squares of cn and dn in favour of the square of sn using the ‘Pythagorean’ identities. to_cn and to_dn apply these identities to replace squares in favour of cn and dn respectively. At most one of these rule sets should be active at any one time to avoid a potential inﬁnite recursion in the simpliﬁcation. The ﬁnal rule list no_glaisher replaces all occurences of the nine subsidiary Jacobi elliptic functions (ns, sc, sd, …by reciprocals and quotients of the ‘basic’ Jacobi elliptic functions sn, cn and dn.

When their arguments are purely numerical, these twelve functions will be evaluated numerically if the rounded switch is ON. For complex arguments it is also better if the complex switch is also ON. However, if complex is OFF, an addition rule will be applied followed by Jacobi rules for a purely imaginary ﬁrst arguments and then the result evaluated numerically. This will be considerably more expensive computationally than when switch complex is ON and rounding errors are likely to be larger.

Jacobi Amplitude Function

The amplitude of u can be evaluated using the jacobiam(u,k) command. A rule list is provided for diﬀerentiation of this functions with respect to either argument.

Arithmetic Geometric Mean (AGM)

A procedure to evaluate the AGM of initial values $a0,b0,c0$ exists as
AGM_function( $a0,b0,c0$ ) and will return
{N,AGM,{a_N,…,a₀},{b_N,…,b₀},{c_N,…,c₀}}, 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 $a0 = 1$ ; $√ ------ b0 = 1− k2$ ; $c0 = k$ .

This procedure and the following one are primarily intended for use in the numerical evaluation of the various elliptic functions and integrals rather than for direct use by users.

Descending Landen Transformation

The procedure to evaluate the Descending Landen Transformation of ϕ and α uses the following equations:

It can be called using landentrans(ϕ₀, α₀) and will return
{{ϕ₀,…,ϕ_n},{α₀,…,α_n}}.

20.19.4 Elliptic Integrals

These again diﬀer 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 ﬁrst 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 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, E(u,k) is used for the Legendre form and ℰ(u,k) for Jacobi’s form. This diﬀers from the 2019 version of this section which followed Lawden [Law89], where the notation D(ϕ,k) and 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 simpliﬁcations 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:

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:

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) and will return the value

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:

When called with one argument EllipticE(k) will return the value of the Complete Elliptic Integral of the Second Kind:

The complementary Elliptic E′ function can be used as EllipticE!′(k) and will return the value

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:

The relationship between the two forms of incomplete elliptic integrals can be expressed as

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 ℰ(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 Z(u,k) whereas the operator Zeta will invoke Riemann’s ζ function.

20.19.5 Some Numerical Utility Functions

These are only operative when the switch rounded is on and their argument is numerical. The ﬁrst pair relate the nome q of the theta functions with the moduli k and k′ = √------
1 − k2

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.6 Jacobi Theta Functions

These theta functions diﬀer from those originally deﬁned by Lisa Temme in a number of respects. Firstly four separate functions of two arguments are deﬁned:

rather than a single function with three arguments (with the ﬁrst argument taking integer values in the range 1 to 4). Secondly the periods are 2π,2π,π and π respectively (NOT 4K, 4K, 2K and 2K). Thirdly the second argument is the modulus τ = a + ib where b = ℑτ > 0 and hence the quasi-period is πτ.

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 deﬁning expansions. elliptictheta3 and elliptictheta4 were, however, single-valued functions of q.

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

Utilising the periodicity and quasi-periodicity of the theta functions some generalised shift rules are implemented to shift their ﬁrst argument into the base period parallelogram with vertices

Together with the relation 𝜗₁(0,τ) = 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 τ positive, the theta functions are evaluated numerically. Note that as τ is necessarily complex, the switch complex must be on.

In what follows a and b will denote the real and imaginary parts of τ respectively and so |q| = exp(−πb) and arg q = π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 = ℑτ close to zero). In such cases the direct algorithm would suﬀer from slow convergence and rounding errors. For such values of |q|, Jacobi’s transformation τ′ = −1∕τ 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 τ purely imaginary since q′ = q^1∕b², but for complex values the gains are somewhat smaller. The Jacobi transformation produces a nome q′ for which |q′| = |q|^{1∕(a²+b²)}.

When ℜq < 0, the Jacobi transformation is preceded by either the modular transformation τ′ = τ + 1 when ℑq < 0, or τ′ = τ − 1 when ℑq > 0, which both have the eﬀect of multiplying q by −1, so that the new nome has a non-negative real part and |a|≤ 1∕2. Thus the worst case occurs for values of the nome q near to ±i where |q′|≈|q|⁴.

By using a series of Jacobi transformations preceded, if necessary by τ-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 𝜗₁ and 𝜗₂, 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πτ∕4). If the principal value is required, it is easily obtained by multiplying by the ‘correcting’ factor q^1∕4 exp(−iπτ∕4).

Derivatives of Theta Functions

These return the dth derivatives of the respective theta functions with respect to their ﬁrst argument u; τ 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 ℑ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 τ 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 signiﬁcant 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.7 Weierstrass Elliptic & Sigma Functions

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, ζ_w is used for the Weierstrassian Zeta function to distinguish it from the Riemann Zeta function and the usual symbol ℘ is used for the Weierstrassian elliptic function itself.

The two primitive periods of the Weierstrass function are 2ω₁ and 2ω₃ and these must satisfy ℑ(ω₃∕ω₁)≠0. The two periods are normally numbered so that τ = ω₃∕ω₁ has a positive imaginary part and hence the nome q = exp(iπτ) satisﬁes |q| < 1.

Any linear combination Ω_m,n = 2mω₁ + 2nω₃ 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 −(ω₁ + ω₃) is often denoted by ω₂ and 2ω₂ is clearly also a period; this accounts for the gap in the numbering of primitive periods. The period ω₂ 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Ω₁ and 2Ω₃ deﬁned by the unimodular integer bilinear transformation:

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 ℑτ 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 ℑτ ≥ √ --
3

∕2 when |q| < 0.06. However, in numerical evaluations in REDUCE it is suﬃcient to use these transformations only until ℑτ > 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 τ 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 σ(u,ω₁,ω₃) is odd and regular everywhere as is the function 𝜗₁(u,τ) to which it is closely related. It has zeros at all lattice points. All three functions ℘, ζ_w and σ are homogenous of degrees -2, -1 and +1 respectively. The functions are related by

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 ℑ(ω₃∕ω₁) > 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₂ and g₃ rather than the primitive periods ω₁ and ω₃ are provided:

Note that for output they are distinguished from the two discussed above by separating the ﬁrst and second arguments by a vertical bar rather than a comma. The rule for diﬀerentiation of the Weierstrass function is simpler when expressed in this alternative:

Other Sigma Functions

These are all even functions, regular everywhere, homogenous of degree zero and doubly quasi-periodic. They are closely related to the theta functions 𝜗₂, 𝜗₃ and 𝜗₄ respectively; but note the diﬀerence 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

These are operative when the switches rounded and complex are ON and their arguments are numerical. The ﬁrst three are referred to as lattice roots and are related to the invariants g₂,g₃, the discriminant Δ = g₂³ − 27g₃² and a closely related invariant G = g₂³∕(27g₃²) of the Weierstrassian elliptic function ℘. The lattice roots also appear in the numerical evaluation of the Weierstrass function. These lattice roots satisfy:

If the discriminant Δ vanishes or equivalently if G = 1, there are at most two distinct lattice roots and the elliptic function degenerates to an elementary one. The advantage of the invariant G is that it is a function of τ = ω₃∕ω₁ 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:

where the lattice parameters have been omitted for conciseness and j,k = 1…3. The quasi-period factors satisfy

As well as the scalar-valued functions discussed above in this section, there are four functions which return a list as their value:

The ﬁrst three are actually more eﬃcient 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:

should reproduce the list {g2, g3}, perhaps with small rounding errors. The corresponding sequence with the calls to lattice_generators and lattice_invariants interchanged (and g2 & g3 replaced by w1 & w3), in general, will not produce the same pair of lattice generators since the generators are only deﬁned 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.8 Inverse Jacobi Elliptic Functions

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 simpliﬁcation formulae illustrated immediately above and for diﬀerentiation of these functions with respect to their two arguments.

Note that for simplicity arccs is now deﬁned to be an odd function of its ﬁrst argument like cs. This choice means that the range of (real) principal values is no longer connected as it is the open set (−K(k),K(k)) ∖ 0. It brings it into line with arcns and arcds whose principal value ranges are necessarily not connected as they have poles at zero; similarly for the range of principal values of arcdc which has a pole at K(k).

When their arguments are numerical and real, these functions will be evaluated numerically if the rounded switch is ON provided, of course that the deﬁning elliptic integral is convergent. Note that in some cases the result may be imaginary even if both arguments are real. In these cases if further computation is required involving the result, it is better if the switch complex is also ON

Note also that for arcdn and arcnd a zero value of the modulus k is excluded (since dn(x,0) = nd(x,0) = 1 ∀x).

As the Jacobi elliptic functions are doubly periodic, the inverse functions are multi-valued. When both arguments are real and |k| <= 1 and when the result exists and is real there are certain restrictions on the range of acceptable values of the ﬁrst parameter x (see the table below).

The numerical value returned is the principal value which 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.

Fn	Domain	Principal Value v	Other real values
arcsn:	\|x\| <= 1	−K(k) <= v <= K(k)	2nK(k) + (−1)ⁿv
arccn:	\|x\| <= 1	0 <= v <= 2K(k)	4nK(k) ± v
arccd:	\|x\| <= 1	0 <= v <= 2K(k)	4nK(k) ± v
arcns:	\|x\| >= 1	−K(k) <= v <= K(k) & v≠0	2nK(k) + (−1)ⁿv
arcnc:	\|x\| >= 1	0 <= v <= 2K(k) & v≠K(k)	4nK(k) ± v
arcdc:	\|x\| >= 1	0 <= v <= 2K(k) & v≠K(k)	4nK(k) ± v
arcdn:	k′ <= x <= 1	0 <= v <= K(k)	2nK(k) ± v
arcnd:	1 <= x <= 1∕k′	0 <= v <= K(k)	2nK(k) ± v
arcds:	\|x\| >= k′	−K(k) <= v <= K(k) & v≠0	2nK(k) + (−1)ⁿv
arcsd:	\|x\| <= 1∕k′	−K(k) <= v <= K(k)	2nK(k) + (−1)ⁿv
arcsc:	x ∈ ℝ	−K(k) < v < K(k)	2nK(k) + v
arccs:	x ∈ ℝ	−K(k) < v < K(k) & v≠0 2nK(k) + v

If one or both arguments are complex, the current version of the package attempts to evaluate the result numerically provided that the switch complex is also ON. This code is still under development and whilst it succeeds in many cases, it may fail with a divergent integral error message. If the function in question has a branch point or pole at the evaluation point, this is of course correct, but in other cases the value may exist, but the current choice of symmetric integral in the algorithm may diverge. This can occur when the ﬁrst parameter is a subrange of the real or imaginary axis.

The value currently returned may not always be the principal value (ie. the value in the fundamental period parallelogram of the elliptic function concerned). However, the value returned should be an inverse value

The numerical values of the inverse functions are calculated by expressing them in terms of the symmetric elliptic integral:

This is then evaluated using a sequence of quadratic transformations which converge rapidly to the elementary hyperbolic integral:

20.19.9 Table of Elliptic Functions and Integrals

Function	Operator

am(u,k)	jacobiam(u,k)
sn(u,k)	jacobisn(u,k)
dn(u,k)	jacobidn(u,k)
cn(u,k)	jacobicn(u,k)
cd(u,k)	jacobicd(u,k)
sd(u,k)	jacobisd(u,k)
nd(u,k)	jacobind(u,k)
dc(u,k)	jacobidc(u,k)
nc(u,k)	jacobinc(u,k)
sc(u,k)	jacobisc(u,k)
ns(u,k)	jacobins(u,k)
ds(u,k)	jacobids(u,k)
cs(u,k)	jacobics(u,k)
Inverse Functions of the above:
arcsn(u,k)	arcsn(u,k)
arccn(u,k)	arccn(u,k)
...
arccs(u,k)	arccs(u,k)
Complete Integral (1st kind) K(k)	ellipticK(k)
K^′(k)	ellipticK!’(k)
Incomplete Integral (1st kind) F(ϕ,k)	ellipticF(phi,k)
Complete Integral (2nd kind) E(k)	ellipticE(k)
E^′(k)	ellipticE!’(k)
Legendre Incomplete Int (2nd kind) E(u,k)	ellipticE(u,k)
Jacobi Incomplete Int (2nd kind) ℰ(u,k)	jacobiE(u,k)
Jacobi’s Zeta Z(u,k)	jacobiZeta(u,k)
𝜗₁(u,τ)	elliptictheta1(u,tau)
𝜗₂(u,τ)	elliptictheta2(u,tau)
𝜗₃(u,τ)	elliptictheta3(u,tau)
𝜗₄(u,τ)	elliptictheta4(u,tau)
℘(u,ω₁,ω₃)	weierstrass(u,omega1,omega3)
ζ_w(u,ω₁,ω₃)	weierstrassZeta(u,omega1,omega3)
σ(u,ω₁,ω₃)	weierstrass_sigma(u,omega1,omega3)
σ₁(u,ω₁,ω₃)	weierstrass_sigma1(u,omega1,omega3)
σ₂(u,ω₁,ω₃)	weierstrass_sigma2(u,omega1,omega3)
σ₃(u,ω₁,ω₃)	weierstrass_sigma3(u,omega1,omega3)
℘(u∣g₂,g₃)	weierstrass1(u,g2,g3)
ζ_w(u∣g₂,g₃)	weierstrassZeta1(u,g2,g3)

(1 + sinα_n+1)(1 + cosα_n) = 2	where α_n+1 < α_n,
tan(ϕ_n+1 − ϕ_n) = cosα_n tanϕ_n	where ϕ_n+1 > ϕ_n.

𝜗₁(z,τ)	= 2e^iπτ∕4 ∑ _n=0^∞(−1)ⁿq^n²+n sin(2n + 1)z
𝜗₂(z,τ)	= 2e^iπτ∕4 ∑ _n=0^∞q^n²+n cos(2n + 1)z
𝜗₃(z,τ)	= 1 + 2∑ _n=1^∞q^n² cos2nz
𝜗₄(z,τ)	= 1 + 2∑ _n=1^∞(−1)ⁿq^n² cos2nz.

ζ_w(u + 2ω_j)	= ζ_w(u) + 2η_j
σ(u + 2ω_j)	= −exp(2η_j(u + ω_j))σ(u)
σ_k(u + 2ω_j)	= exp(2η_j(u + ω_j))σ_k(u) if j≠k
σ_j(u + 2ω_j)	= −exp(2η_j(u + ω_j))σ_j(u)
ζ_w(ω_j)	= η_j
σ_j(ω_j)	= 0,

20.19 ELLIPFN: A Package for Elliptic Functions and Integrals