Special Functions: Elliptic Integrals, Elliptic Functions & Theta Functions

Overview

Elliptic integrals, elliptic functions, and theta functions form a closely knit family, tied together by inverse relationships and identities. Roughly, the relationships are as follows:

• An elliptic integral $F(\varphi, k)$ takes the amplitude $\varphi$ as its argument and returns an "arc length."
• The Jacobi elliptic functions $\operatorname{sn}, \operatorname{cn}, \operatorname{dn}$ are its inverses. For $u = F(\varphi, k)$ one has $\operatorname{sn}(u, k) = \sin\varphi$, and so on.
• The theta functions $\vartheta_1 \ldots \vartheta_4$ are the rapidly converging $q$-series "building blocks" from which elliptic functions are assembled. Both the Jacobi functions and the Weierstrass functions can be written as ratios of theta functions.
• The Weierstrass $\wp$ function is a doubly periodic function on the period lattice $2\omega_1 \mathbb{Z} + 2\omega_2 \mathbb{Z}$, with a double pole at each lattice point. It satisfies the differential equation $(\wp')^2 = 4\wp^3 - g_2\wp - g_3$.

The evaluation algorithms split broadly into two families. Elliptic integrals are computed via the Carlson symmetric forms + the duplication theorem, or via the iteration of the AGM (arithmetic-geometric mean). Theta functions and Weierstrass functions are computed via the $q$-series (with modular reduction as needed). In all cases the iteration or series converges quadratically or super-exponentially, reaching high accuracy in a small number of steps.

Related API: ellipticK, carlsonRF, jacobiSn and others (jacobiTheta1 …).

Carlson Symmetric Forms

Carlson expressed elliptic integrals uniformly through four basic integrals that are symmetric in their arguments:

$$R_F(x,y,z) = \frac{1}{2}\int_0^\infty \frac{dt}{\sqrt{(t+x)(t+y)(t+z)}}$$

$$R_D(x,y,z) = \frac{3}{2}\int_0^\infty \frac{dt}{(t+z)^{3/2}\sqrt{(t+x)(t+y)}}$$

$$R_J(x,y,z,p) = \frac{3}{2}\int_0^\infty \frac{dt}{(t+p)\sqrt{(t+x)(t+y)(t+z)}}, \qquad R_C(x,y) = \frac{1}{2}\int_0^\infty \frac{dt}{(t+y)\sqrt{t+x}}$$

$R_C(x,y) = R_F(x,y,y)$ is a degenerate form, used repeatedly inside $R_J$.

Iterative evaluation via the duplication theorem

The power of the Carlson forms lies in the duplication theorem. At each iteration, form $\lambda = \sqrt{x}\sqrt{y} + \sqrt{y}\sqrt{z} + \sqrt{z}\sqrt{x}$ and update the arguments as

$$x \leftarrow \frac{x+\lambda}{4}, \quad y \leftarrow \frac{y+\lambda}{4}, \quad z \leftarrow \frac{z+\lambda}{4}.$$

Repeating this drives $x, y, z$ toward one another, clustering around the common mean $A = (x+y+z)/3$. Once the relative deviations $\delta_x = (A-x)/A$, etc. are all small, stop iterating and add a single correction polynomial about $A$:

$$R_F \approx \frac{1}{\sqrt{A}}\left(1 - \frac{E_2}{10} + \frac{E_3}{14} + \frac{E_2^2}{24} - \frac{3\,E_2 E_3}{44}\right)$$

Here $E_2, E_3$ are symmetric functions of the deviations $\delta_x, \delta_y, \delta_z$. Because the deviations are small, the polynomial needs only a few terms, and overall the accuracy improves exponentially with the number of iterations.

$R_D$ and $R_J$ differ in that each iteration also accumulates a summation term $\sigma$ (the summation term of $R_J$ calls $R_C$ once). sangi's carlsonRF / carlsonRD / carlsonRJ / carlsonRC all implement this iteration (convergence judged at machine epsilon, capped at 100 iterations).

Expressing the complete and incomplete elliptic integrals in Carlson form

The Legendre standard forms can be written as combinations of the Carlson forms. This gives uniform evaluation while avoiding loss of significance near the endpoints:

$$K(k) = R_F(0,\, 1-k^2,\, 1), \qquad E(k) = R_F(0,\, 1-k^2,\, 1) - \frac{k^2}{3} R_D(0,\, 1-k^2,\, 1)$$

$$\Pi(n,k) = R_F(0,\, 1-k^2,\, 1) + \frac{n}{3} R_J(0,\, 1-k^2,\, 1,\, 1-n)$$

The incomplete forms take $\cos^2\varphi$ and the like as arguments:

$$F(\varphi, k) = \sin\varphi \cdot R_F(\cos^2\varphi,\ 1 - k^2\sin^2\varphi,\ 1)$$

$$E(\varphi, k) = \sin\varphi \cdot R_F(\cdots) - \frac{k^2}{3}\sin^3\varphi \cdot R_D(\cdots)$$

sangi provides these as ellipticK, ellipticE, ellipticPi (one argument gives the complete integral, two or three give the incomplete one). It also handles endpoint singularities, for example returning $K(\pm 1) = \infty$ when $k^2 \geq 1$. The Carlson forms extend naturally to complex arguments, and the same iteration runs for Complex<double> / Complex<Float> (with convergence judged by the scalar comparison $|\delta| < \varepsilon$).

AGM (Arithmetic-Geometric Mean)

$K(k)$ has another classical fast evaluation. The arithmetic-geometric mean $M(a, b)$ is the common limit of the iteration starting from $a_0 = a,\ b_0 = b$:

$$a_{n+1} = \frac{a_n + b_n}{2}, \qquad b_{n+1} = \sqrt{a_n b_n}.$$

Because the arithmetic mean and geometric mean approach each other at each step, $a_n - b_n$ tends to zero with quadratic convergence (roughly squared each time). The complete elliptic integral of the first kind is written in terms of the AGM as

$$K(k) = \frac{\pi}{2\, M(1,\ k')}, \qquad k' = \sqrt{1 - k^2}.$$

Owing to quadratic convergence, reaching a target precision of $B$ digits takes only $\mathcal{O}(\log B)$ iterations. Even at 1000 digits a dozen or so iterations suffice, and since each iteration is a single square root, the AGM is overwhelmingly more advantageous than series summation in high-precision computation.

Relation to the descending Landen transformation

Each AGM step is equivalent to a descending Landen transformation that shrinks the modulus $k$. The Landen transformation

$$k_1 = \frac{1 - k'}{1 + k'} \quad (k' = \sqrt{1-k^2})$$

produces a new modulus $k_1 < k$ and reduces $K(k)$ to $K(k_1)$. As the modulus approaches $0$ with each iteration, $K \to \pi/2$, so the integral value can be recovered by retracing the transformations. The AGM sequence $(a_n, b_n)$ and the Landen sequence $(k_n)$ are tied together by $b_n/a_n = k_n'$ and are essentially two views of the same iteration.

sangi's high-precision path for the theta functions (jacobiThetaZeroAGM) applies this AGM structure to the theta values at $z=0$. Using the fact that $a_n = \vartheta_3(0, q^{2^n})^2,\ b_n = \vartheta_4(0, q^{2^n})^2$ satisfy the forward AGM (with $a_n, b_n \to 1$ since $q^{2^n} \to 0$), it computes $a_N, b_N$ from the series at a sufficiently small $q_N = q^{2^N}$, then retraces the inverse AGM $d = \sqrt{a_n^2 - b_n^2},\ a_{n-1} = a_n + d,\ b_{n-1} = a_n - d$ for $N$ steps to obtain $\vartheta_3(0,q) = \sqrt{a_0}$, $\vartheta_4(0,q) = \sqrt{b_0}$, $\vartheta_2(0,q) = (a_0^2 - b_0^2)^{1/4}$.

Jacobi Elliptic Functions

The Jacobi elliptic functions are defined as the inverses of the elliptic integral. For $u = F(\varphi, k)$, writing the amplitude as $\varphi = \operatorname{am}(u, k)$,

$$\operatorname{sn}(u, k) = \sin\varphi, \qquad \operatorname{cn}(u, k) = \cos\varphi, \qquad \operatorname{dn}(u, k) = \sqrt{1 - k^2 \operatorname{sn}^2(u, k)}.$$

Evaluating the amplitude via the AGM / descending Landen

sangi computes the amplitude $\operatorname{am}(u, k)$ with the AGM sequence and an inverse transformation. The procedure is as follows:

• Forward AGM: starting from $a_0 = 1,\ b_0 = k',\ c_0 = k$, iterate $a_{n+1} = (a_n + b_n)/2,\ b_{n+1} = \sqrt{a_n b_n},\ c_{n+1} = (a_n - b_n)/2$, retaining the sequences $\{a_n\}, \{c_n\}$ up to the index $N$ at which $c_n$ becomes sufficiently small.
• Initial phase: take $\varphi_N = 2^N a_N u$.
• Inverse transformation (descent): retrace $\varphi_{n-1} = \tfrac{1}{2}\bigl(\varphi_n + \arcsin(\tfrac{c_n}{a_n}\sin\varphi_n)\bigr)$ from $n = N$ down to $1$. The terminal $\varphi_0$ is $\operatorname{am}(u, k)$.

The rest follows from $\operatorname{sn} = \sin(\operatorname{am})$, $\operatorname{cn} = \cos(\operatorname{am})$, $\operatorname{dn} = \sqrt{1 - k^2 \operatorname{sn}^2}$. The boundary moduli are handled in closed form: at $k = 0$, $\operatorname{am}(u, 0) = u$; at $k = 1$, $\operatorname{am}(u, 1) = \operatorname{gd}(u) = 2\arctan(\tanh(u/2))$ (the Gudermannian function).

Because $\frac{c_n}{a_n}\sin\varphi_n$ can exceed $|\cdot| > 1$ due to numerical error in the inverse transformation, sangi's implementation clamps the argument of $\arcsin$ to $[-1, 1]$ for safe evaluation. The provided functions are jacobiSn, jacobiCn, jacobiDn (each taking the modulus $k$ as an argument).

Theta Functions

The Jacobi theta functions are defined by power series in the nome $q$ (with $|q| < 1$). Because the exponent grows as the square of $n$, the series converges super-exponentially fast:

$$\vartheta_1(z, q) = 2\sum_{n=0}^\infty (-1)^n q^{(n+1/2)^2} \sin((2n+1)z)$$

$$\vartheta_2(z, q) = 2\sum_{n=0}^\infty q^{(n+1/2)^2} \cos((2n+1)z)$$

$$\vartheta_3(z, q) = 1 + 2\sum_{n=1}^\infty q^{n^2} \cos(2nz), \qquad \vartheta_4(z, q) = 1 + 2\sum_{n=1}^\infty (-1)^n q^{n^2} \cos(2nz)$$

The nome $q$ and the half-period ratio $\tau$

The nome is determined from the half-period ratio $\tau$ by $q = e^{i\pi\tau}$. If $\tau$ lies in the upper half-plane ($\operatorname{Im}\tau > 0$), then $|q| = e^{-\pi \operatorname{Im}\tau} < 1$ is guaranteed. sangi's qFromTau provides this conversion.

Since the series terms decay like $|q|^{n^2}$, machine precision is reached in a few terms for $|q| \leq 0.5$ or so. The implementation computes $q^{(n+1/2)^2}$ and $q^{n^2}$ directly by exponentiation and truncates once a term is sufficiently small relative to the partial sum. For a complex nome $q$, it switches between $q^k = \exp(k \log q)$ (half-integer powers) and integer powers, computing the integer powers by stable binary iteration.

Reduction of $\tau$ (modular transformation)

When $|q|$ is close to $1$ ($= \operatorname{Im}\tau$ is small), the series converges slowly. In that case the generators of $\mathrm{SL}(2,\mathbb{Z})$

• $T:\ \tau \to \tau + 1$ — $\vartheta_1, \vartheta_2$ pick up a phase $e^{i\pi/4}$, and $\vartheta_3 \leftrightarrow \vartheta_4$ are swapped
• $S:\ \tau \to -1/\tau$ — prefactors $(-i\tau)^{-1/2}$ and $\exp(i\pi z^2/\tau)$ appear, and the indices are permuted

are applied repeatedly to reduce $\tau$ to the fundamental domain ($|\tau| \geq 1$ and $|\operatorname{Re}\tau| \leq 1/2$, i.e. $\operatorname{Im}\tau \geq \sqrt{3}/2 \approx 0.866$) before evaluating the series. This always allows computation with a sufficiently small $|q|$, guaranteeing convergence. sangi's $\tau$-form API (jacobiTheta1Tau … jacobiTheta4Tau) performs this reduction internally.

For high-precision $z = 0$ values, the AGM acceleration of the previous section (jacobiThetaZeroAGM) obtains $\vartheta_2, \vartheta_3, \vartheta_4$ simultaneously in fewer iterations than increasing the number of series terms would. The provided functions are jacobiTheta1 … jacobiTheta4 (with real-argument, complex-argument, complex-nome, and arbitrary-precision overloads).

Weierstrass Elliptic Functions

The Weierstrass $\wp$ function is a doubly periodic function on the lattice spanned by the half-periods $\omega_1, \omega_2$, an even function with a double pole at each lattice point. It satisfies the differential equation

$$(\wp')^2 = 4\wp^3 - g_2 \wp - g_3$$

where $g_2, g_3$ are called the invariants. $\zeta_W$ and $\sigma_W$ are linked by $\zeta_W'(z) = -\wp(z)$, $\sigma_W'(z)/\sigma_W(z) = \zeta_W(z)$ (these are distinct from Riemann's zeta and sigma).

Evaluation via theta functions

Rather than summing the lattice directly, sangi evaluates via the theta representation (DLMF 23.6). It first checks that $\tau = \omega_2/\omega_1$ lies in the upper half-plane, then prepares the nome $q = e^{i\pi\tau}$, $\alpha = \pi/(2\omega_1)$, and $v = \alpha z$. Then $\wp$ is

$$\wp(z) = e_3 + \left(\alpha\, \vartheta_2(0,q)\, \vartheta_3(0,q)\, \frac{\vartheta_4(v,q)}{\vartheta_1(v,q)}\right)^2, \qquad e_3 = -\alpha^2\, \frac{\vartheta_2^4(0,q) + \vartheta_3^4(0,q)}{3}.$$

At the lattice points $\vartheta_1(v) \approx 0$, so those are returned as poles (infinity). $\zeta_W$ and $\sigma_W$ can likewise be written in terms of $\vartheta_1$ and its derivatives $\vartheta_1', \vartheta_1'''$:

$$\zeta_W(z) = \eta_1 \frac{z}{\omega_1} + \alpha\, \frac{\vartheta_1'(v,q)}{\vartheta_1(v,q)}, \qquad \eta_1 = -\frac{\pi^2}{12\,\omega_1}\, \frac{\vartheta_1'''(0,q)}{\vartheta_1'(0,q)}$$

$$\sigma_W(z) = \frac{2\omega_1}{\pi}\, \frac{\vartheta_1(v,q)}{\vartheta_1'(0,q)}\, \exp\!\left(\frac{\eta_1 z^2}{2\omega_1}\right)$$

Because the theta series converges quickly, this evaluates more stably and faster than truncating the lattice sum directly. sangi's weierstrassP, weierstrassZeta, weierstrassSigma take the half-periods $\omega_1, \omega_2$ as arguments and, with Complex<double>, handle arbitrary half-periods and arbitrary $z$. The real-argument version handles $\omega_1$ real and $\omega_2$ purely imaginary (a rectangular lattice) and returns real values away from the poles. To handle a general nome $q$ (= $\operatorname{Re}\tau \neq 0$), the theta series for complex $q$ is expanded directly internally.

Comparison Table

Rough guidance: for elliptic integrals, the argument-symmetric and robust Carlson forms are the default. If $K(k)$ is needed many times or at high digits, use the AGM. For the Jacobi functions, use the AGM plus the inverse descending Landen transformation. For theta functions and Weierstrass functions, use the $q$-series (with modular reduction of $\tau$ if convergence is slow).

References

• Carlson, B. C. (1995). "Numerical computation of real or complex elliptic integrals". Numerical Algorithms, 10(1), 13–26.
• Olver, F. W. J. et al. (eds.). NIST Digital Library of Mathematical Functions (DLMF), §19 (elliptic integrals), §20 (theta functions), §22 (Jacobi elliptic functions), §23 (Weierstrass elliptic functions). https://dlmf.nist.gov/
• Abramowitz, M. & Stegun, I. A. (1964). Handbook of Mathematical Functions. National Bureau of Standards. (§16 Jacobi elliptic functions, §17 elliptic integrals)
• Whittaker, E. T. & Watson, G. N. (1927). A Course of Modern Analysis. 4th ed. Cambridge University Press. (Ch. 20–22)
• Borwein, J. M. & Borwein, P. B. (1987). Pi and the AGM. Wiley.