Elliptic integrals, elliptic functions and modular forms

Elliptic integrals and elliptic functions reduce to a few primitives that converge quickly: the arithmetic-geometric mean, which gives the complete elliptic integrals; Carlson's duplication algorithm, which gives the incomplete integrals; and theta functions, which give the Weierstrass elliptic functions and the modular forms after a modular transformation has moved \tau to the fundamental domain. The sections below describe these methods as implemented in acb (the AGM), acb_elliptic and acb_modular, following [Joh2019]; the summation of theta series is described in Theta functions.

The arithmetic-geometric mean#

For a, b > 0, the arithmetic-geometric mean \operatorname{agm}(a, b) is the common limit of a_{n+1} = (a_n + b_n)/2, b_{n+1} = \sqrt{a_n b_n}, with a_0 = a, b_0 = b. The convergence is quadratic: |a_{n+1} - b_{n+1}| \le |a_n - b_n|^2/(8 b_{n+1}), so O(\log p) steps give p bits once a_n and b_n agree to a few bits. With complex variables, it is convenient to use the univariate function M(z) = \operatorname{agm}(1, z); then \operatorname{agm}(a,b) = a M(b/a) (acb_agm), and M is analytic on \mathbb{C} \setminus (-\infty, 0].

Functional equation. If the real part of z is not nonnegative, acb_agm1 first applies M(z) = \frac{z+1}{2} M(u) with u = 2\sqrt{z}/(z+1), which has nonnegative real part (rounding errors that make the ball for u cross the imaginary axis only inflate the result). For the derivative this becomes M'(z) = [M(u) - (z-1) M'(u) / ((1+z) \sqrt{z})] / 2.

Square roots. In the iteration the correct square root is the continuous branch \sqrt{a_n}\sqrt{b_n}, which is computed as \sqrt{a_n b_n}, i\sqrt{-a_n b_n} or -i\sqrt{-a_n b_n} when a_n and b_n both have positive real part, both have nonnegative imaginary part, or both have nonpositive imaginary part, and as \sqrt{a_n}\sqrt{b_n} otherwise.

This gives a rigorous enclosure when the iteration is stopped, and for real positive variables the lower and upper bounds b_n \le M \le a_n can be used directly (arb_agm). It is slightly faster to stop when a_n and b_n agree to about p/10 bits and finish with a series: with w = (a - b)/(a + b),

valid for |w| < 1 and a, b with nonnegative real part, where the tail is bounded by \sum_{k \ge 10} |w|^k/64.

First derivative. acb_agm1_cpx computes (M(z), M'(z)) for exact z with |\arg z| \le 3\pi/4 by a central finite difference. Since |M(z)| \le \max(1, |z|), Cauchy's integral formula gives |M^{(k)}(z)/k!| \le C D^k with C = \max(1, |z| + r) and D = 1/r for 0 < r < |z| (FLINT takes r \approx |z|/4). For hD < 1,

so two evaluations at about 1.5 times the precision give both values, against twice the precision for a one-sided difference. For inexact z, the midpoint is used and the propagated error is bounded with |M'(z+\varepsilon)| \le \max(1, |z|+|\varepsilon|+r)/r and |M''(z+\varepsilon)| \le 2\max(1, |z|+|\varepsilon|+r)/r^2, valid if the disk of radius |\varepsilon| + r around z avoids the negative real axis (with r \approx |z|/2).

Higher derivatives. W(z) = 1/M(z) = \frac{2}{\pi} K(1 - z^2) is D-finite. The coefficients of W(z+x) = \sum_{k} c_k x^k satisfy -2 z (z^2-1) c_2 = (3z^2-1) c_1 + z c_0 and

which reduces at z = 1 to -(k+2)^2 c_{k+2} = (3k(k+3)+7) c_{k+1} + (k+1)^2 c_{k}. Given (M(z), M'(z)), the remaining coefficients follow by the recurrence and one power series inversion.

Complete elliptic integrals#

With M as above,

where K' is the derivative; both cost O(\log p) square roots (acb_elliptic_k, acb_elliptic_e). The Taylor coefficients of K(m + x) (acb_elliptic_k_jet) satisfy the recurrence of the hypergeometric equation of K = \frac{\pi}{2}\,{}_2F_1(\frac12, \frac12; 1; m):

The complete integral of the third kind, \Pi(n, m), is computed with Carlson's integrals, \Pi(n, m) = R_F(0, 1-m, 1) + \frac{n}{3} R_J(0, 1-m, 1, 1-n).

Carlson's symmetric integrals#

The symmetric integrals R_F, R_J (with R_C, R_D as special cases) and R_G are defined on Carlson symmetric integrals. All three Legendre integrals reduce to them on the strip |\operatorname{Re}\phi| \le \pi/2 [NIST2012, §19.25], and are extended outside it by quasi-periodicity.

Duplication. For \lambda = \sqrt{x}\sqrt{y} + \sqrt{y}\sqrt{z} + \sqrt{z}\sqrt{x},

and the differences between the arguments are divided by about 4 at each step, while their mean A = (x+y+z)/3 is preserved up to the factor introduced by the step. After j steps, with X = 1 - x/A, Y = 1 - y/A, Z = -X - Y, the elementary symmetric polynomials E_2 = XY - Z^2 and E_3 = XYZ of the deviations are of order 4^{-2j} and 4^{-3j}, and

Carlson [Car1995] uses the terms of total degree 2a + 3b \le 7 and enough duplication steps to make the remainder negligible, which costs O(p) square roots at precision p. FLINT instead sums the series up to degree N \approx 2.1\,p^{0.4} (real case) or 2.5\,p^{0.4} (complex case), at most 500, so that only O(p/N) duplication steps are needed [Joh2019]: the iteration stops as soon as (\max(|x - y|, |x - z|)/|x|)^N < 2^{-p}. The coefficients are computed exactly with a common denominator, the polynomial in E_2 is evaluated with precomputed powers and the polynomial in E_3 by Horner's rule. The tail is bounded using that the coefficient of X^{n_1} Y^{n_2} Z^{n_3} with n_1 + n_2 + n_3 = n in the expansion has absolute value at most 2 (9/8)^n, so that the truncation error is at most 2 \sum_{n \ge N} (\tfrac98 \max(|X|, |Y|, |Z|))^n.

R_J. The duplication theorem for R_J adds a term 6 \cdot 4^{-j} R_C(1, 1 + \delta_j) (up to scaling) at each step, and the final Taylor series is in the symmetric polynomials E_2, \ldots, E_5 of four deviations [Car1995]; R_C(1, 1+x) = \operatorname{atan}(\sqrt{x})/\sqrt{x} is computed by acb_elliptic_rc1. The principal square roots in the duplication can cross branch cuts of R_J for complex arguments, so the algorithm is only used where it is known to be valid, and otherwise R_J is shifted into the right half-plane by numerical integration; see Correctness of the duplication algorithm for R_J. For R_D the duplication algorithm is valid everywhere.

R_G is computed from 2R_G(x,y,z) = z R_F(x,y,z) - \frac13 (x-z)(y-z) R_D(x,y,z) + \sqrt{xy/z}.

Weierstrass elliptic functions#

For the lattice \mathbb{Z} + \tau\mathbb{Z} with half-periods \frac12 and \frac\tau2, the Weierstrass functions are quotients of theta functions [NIST2012, §23.6]. With \theta_j = \theta_j(0, \tau) and \theta_j(z) = \theta_j(z, \tau) in the notation of acb_modular,

where \eta(\tau) is the Dedekind eta function. The lattice roots are e_1 = \frac{\pi^2}{3}(\theta_3^4 + \theta_4^4), e_2 = \frac{\pi^2}{3}(\theta_2^4 - \theta_4^4) and e_3 = -\frac{\pi^2}{3}(\theta_2^4 + \theta_3^4), the invariants are g_2 = 60 G_4 and g_3 = 140 G_6, and the inverse function is the elliptic integral \wp^{-1}(z) = R_F(z - e_1, z - e_2, z - e_3).

All of these need theta values at one \tau, and the theta functions are evaluated after reducing \tau to the fundamental domain, so the cost is O(\sqrt{p}) multiplications for every \tau.

Modular forms#

A modular form of weight k satisfies f(g\tau) = (c\tau + d)^k f(\tau) for g = \begin{pmatrix} a & b \\ c & d\end{pmatrix} \in \mathrm{SL}_2(\mathbb{Z}), so it can be evaluated at the image \tau' of \tau in the fundamental domain \mathcal{F} = \{|\operatorname{Re}\tau| \le \frac12, |\tau| \ge 1\}, where \operatorname{Im}(\tau') \ge \sqrt{3}/2 and the q-series converge geometrically with ratio |q| \le e^{-\pi\sqrt 3} \approx 0.0043 (for q = e^{2\pi i \tau}).

Each inversion multiplies \operatorname{Im}(\tau) by 1/|\tau|^2 > 1, which guarantees termination. acb_modular_fundamental_domain_approx runs the algorithm on the midpoint in double arithmetic when possible and in arf arithmetic at the working precision otherwise, and then applies g to the ball \tau once, which does not need to be exact: reducing to |\tau'|^2 \ge 63/64 instead of 1 is enough for fast convergence.

Transformation formulas. The eta function transforms as \eta(g\tau) = \varepsilon(a,b,c,d)\sqrt{c\tau + d}\,\eta(\tau) with a 24th root of unity \varepsilon given by Dedekind sums (acb_modular_epsilon_arg), and the four theta functions are permuted with 8th roots of unity and the factors of acb_modular_theta_transform [Rad1973].

q-series. On \mathcal{F}, the theta constants and the eta function are lacunary series,

(with q = e^{2\pi i\tau} and q_1 = e^{\pi i\tau}), which need O(\sqrt{p}) terms. The powers of q with exponents in the sequence of generalised pentagonal numbers (or squares and pronic numbers) are computed with an addition sequence in which most terms cost one multiplication; Enge, Hart and Johansson [EHJ2016] construct short addition sequences for these exponents (acb_modular_addseq_eta, acb_modular_addseq_theta). For many terms, rectangular splitting is used instead: the series is written as a polynomial in q^m for a well-chosen m such that the exponents modulo m take few values, which replaces most full multiplications by cheaper additions and scalar operations. Each term is computed at a precision reduced according to its magnitude (_acb_modular_mul).

Forms computed from theta constants. With \theta_j = \theta_j(0,\tau) at the reduced point:

and the higher Eisenstein series follow from the recurrence (2k+1)(k-3) c_k = 3 \sum_{j=2}^{k-2} c_j c_{k-j} for c_k = (2k-1) G_{2k} [Apostol1997, Theorem 1.13]. j is invariant under \mathrm{SL}_2(\mathbb{Z}) and \lambda under \Gamma(2) (for \lambda, the transformation of the theta functions is applied).

Class polynomials. The Hilbert class polynomial H_D is the product of x - j(\tau_i) over the CM points \tau_i = (-b + \sqrt{D})/(2a) of the reduced forms of discriminant D. Its coefficients have about \pi\sqrt{|D|}\sum_i 1/a_i bits; computing the roots to that precision and expanding the product with a binary product tree gives the polynomial with a certificate when every coefficient ball contains a unique integer [Enge2009]. See Hilbert class polynomials.

Use in FLINT#

• acb_agm1, acb_agm1_cpx, acb_agm, arb_agm (acb, arb): the AGM.
• acb_elliptic_k, acb_elliptic_k_jet, acb_elliptic_e: complete integrals by the AGM; acb_elliptic_pi: via R_F and R_J.
• acb_elliptic_rf, acb_elliptic_rj_carlson, acb_elliptic_rj, acb_elliptic_rj_integration, acb_elliptic_rg, acb_elliptic_rc1: Carlson's integrals; acb_elliptic_f, acb_elliptic_e_inc, acb_elliptic_pi_inc: Legendre forms.
• acb_elliptic_p, acb_elliptic_p_prime, acb_elliptic_p_jet, acb_elliptic_zeta, acb_elliptic_sigma, acb_elliptic_roots, acb_elliptic_invariants, acb_elliptic_inv_p: Weierstrass functions.
• acb_modular_fundamental_domain_approx, acb_modular_theta, acb_modular_eta, acb_modular_j, acb_modular_delta, acb_modular_lambda, acb_modular_eisenstein, acb_modular_hilbert_class_poly: modular forms.