Fredrik Johansson

Numerical Evaluation of Elliptic Functions, Elliptic Integrals and Modular Forms

We describe algorithms to compute elliptic functions and their relatives (Jacobi theta functions, modular forms, elliptic integrals, and the arithmeticgeometric mean) numerically to arbitrary precision with rigorous error bounds for arbitrary complex variables. Implementations in ball arithmetic are available in the open source Arb library. We discuss the algorithms from a concrete implementation point of view, with focus on performance at tens to thousands of digits of precision. sage: C = ComplexBallField(100) sage: tau = (1 + C(-3).sqrt())/2 sage: z = C(2 + 2 * I) sage: z.elliptic_p(tau) [-13.7772161934928750714214345 +/-6.41e-26] + [+/-3.51e-26] * I sage: (z + 5 + 6 * tau).elliptic_p(tau) [-13.777216193492875071421435 +/-9.69e-25] + [+/-4.94e-25] * I

This text covers the algorithms used in Arb and discusses some of the implementation aspects. The algorithms are general and work well in most situations. However, we note that the code in Arb does not use the best available algorithms in all cases, and we will point out some of the possible improvements.

1 Available at http://arblib.org. The functionality for modular forms and elliptic functions can be found in the acb_modular (http://arblib.org/acb_modular.html) and acb_elliptic (http://arblib.org/acb_elliptic.html) modules. 2 Of course, for applications that do not require rigorous error bounds, all the algorithms can just as well be implemented in ordinary floating-point arithmetic.

Introduction

The elliptic functions and their relatives have many applications in mathematical physics and number theory. Among the elliptic family of special functions, we count the elliptic functions proper (i.e. doubly periodic meromorphic functions) as well as the quasiperiodic Jacobi theta functions, the closely related classical modular forms and modular functions on the upper half plane, and elliptic integrals which are the inverse functions of elliptic functions.

Our goal is to give a modern treatment of numerical evaluation of these functions, using algorithms that meet several criteria:

• Full domain. We should be able to compute the functions for arbitrary complex values of all parameters where this is reasonable, with sensible handling of branch cuts for multivalued functions. • Arbitrary precision. Precision much higher than 16-digit (or 53-bit) machine arithmetic is sometimes needed for solving numerically ill-conditioned problems. For example, extremely high precision evaluations are employed in mathematical physics to find closed-form solutions for sums and integrals using in-Fredrik Johansson, INRIA -LFANT, CNRS -IMB -UMR 5251, Université de Bordeaux, 33400 Talence, France, e-mail: fredrik.johansson@gmail.com teger relation methods [START_REF] Bailey | High-precision arithmetic in mathematical physics[END_REF]. Computations with elliptic functions and modular forms requiring hundreds or thousands of digits are commonplace in algebraic and analytic number theory, for instance in the construction of discrete data such as class polynomials from numerical approximations [START_REF] Enge | The complexity of class polynomial computation via floating point approximations[END_REF]. • Rigorous error bounds. When we compute an approximation y ≈ f (x), we should also be able to offer a bound |yf (x)| ≤ ε, accounting for all intermediate rounding and approximation errors in the algorithm. • Efficiency. The algorithms should be robust and efficient for arguments that are very small, large, or close to singularities. For arbitrary-precision implementations, a central concern is to ensure that the computational complexity as a function of the precision does not grow too quickly. At the same time, we must have in mind that the algorithms with the best theoretical asymptotic complexity are not necessarily the best in practice, and we should not sacrifice efficiency at moderate precision (tens to hundreds of digits).

It turns out that these goals can be achieved simultaneously and with reasonable implementation effort thanks to the remarkable amount of structure in the elliptic function family -in contrast to many other problems in numerical analysis!

The present author has implemented routines for elliptic functions, Jacobi theta functions, elliptic integrals and commonly-used modular forms and functions as part of the open source Arb library for arbitrary-precision ball arithmetic [START_REF] Johansson | Arb: efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF]. 1 The idea behind ball arithmetic is to represent numerical approximations with error bounds attached, as in π ∈ [3.14159265358979 ± 3.57 • 10 -15 ]. The algorithms use ball arithmetic internally for automatic propagation of error bounds, in combination with some pen-and-paper bounds mainly for truncations of infinite series. 2 The functions in Arb can be used directly in C or via the high-level wrappers in Sage [START_REF]SageMath, the Sage Mathematics Software System[END_REF] or the Julia package Nemo [START_REF] Fieker | Nemo/Hecke: computer algebra and number theory packages for the Julia programming language[END_REF]. As an example, we use the Arb interface in Sage to evaluate the Weierstrass elliptic function ℘ on the lattice (1, τ) with τ = 1 2 (1 + √ 3i). We check ℘(z) = ℘(z + 5 + 6τ) at the arbitrarily chosen point z = 2 + 2i, here using 100-bit precision:

There is a vast body of literature on elliptic functions and integrals, and we will not be able to explore the full breadth of computational approaches. We will, in particular, focus on arbitrary-precision arithmetic and omit techniques that only matter in machine precision. A good overview and a comprehensive bibliography can be found in chapters [START_REF] Molin | Numerical Integration and L-Functions Computations[END_REF][START_REF] Molin | Computing period matrices and the Abel-Jacobi map of superelliptic curves[END_REF][START_REF] Nogneng | On the evaluation of some sparse polynomials[END_REF] and 23 of the NIST Handbook of Mathematical Functions [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF] or its online counterpart, the Digital Library of Mathematical Functions3 . Cohen's book on computational number theory [START_REF] Cohen | A course in computational algebraic number theory[END_REF] is also a useful resource.

Many other packages and computer algebra systems also provide good support for evaluating elliptic and related functions, though not with rigorous error bounds; we mention Pari/GP [28] and of course Maple and Mathematica. For a nice application of Weierstrass elliptic functions in astrodynamics and a fast machine-precision implementation of these functions, we mention the work by Izzo and Biscani [START_REF] Izzo | On the astrodynamics applications of Weierstrass elliptic and related functions[END_REF].

The algorithms that we review are well known, but they are sometimes described without discussing arbitrary complex variables, variable precision, or error bounds. We attempt to provide an account that is complementary to the existing literature, and we also discuss some minor improvements to the algorithms as they are usually presented. For example, we have optimized Carlson's algorithm for symmetric elliptic integrals to reduce the asymptotic complexity at high precision (section 1.6.3), and we make several observations about the deployment of ball arithmetic.

General strategy

Algorithms for evaluating mathematical functions often have two stages: argument reduction, followed by evaluation of a series expansion [START_REF] Brent | Modern Computer Arithmetic[END_REF][START_REF] Muller | Elementary Functions[END_REF].

Minimax polynomial or rational function approximations are usually preferred for univariate functions in machine precision, but truncated Taylor series expansions are the tools of choice in arbitrary-precision arithmetic, for two reasons. First, precomputing minimax approximations is not practical, and second, we can exploit the fact that the polynomials arising from series expansions of special functions are typically not of generic type but highly structured.

Argument reduction consists of applying functional equations to move the argument to a part of the domain where the series expansion converges faster. In many cases, argument reduction is needed to ensure convergence in the first place. Argument reduction also tends to improve numerical stability, in particular by avoiding alternating series with large terms that would lead to catastrophic cancellation.

The classical elliptic and modular functions are no exception to this general pattern, as shown in Table 1.1. For elliptic integrals, the argument reduction consists of using contracting transformations to reduce the distance between the function arguments, and the series expansions are hypergeometric series (in one or several variables). For the elliptic and modular functions, the argument reduction consists of using modular transformations and periodicity to move the lattice parameter to the fundamental domain and the argument to a lattice cell near the origin, and the series expansions are the sparse q-series of Jacobi theta functions.

In the following text, we will first discuss the computation of elliptic functions starting with the special case of modular forms and functions before turning to general elliptic and Jacobi theta functions. Then, we discuss elliptic integrals, first covering the easier case of complete integrals before concluding with the treatment of incomplete integrals.

We comment briefly on error bounds. Since ball arithmetic automatically tracks the error propagation during series evaluation and through argument reduction steps, the only error analysis that needs to be done by hand is to bound the series truncation errors. If

f (x) = ∑ ∞ k=0 t k (x), we compute ∑ N k=0 t k (x)
(x)| ≤ AC k with 0 ≤ C < 1, then |R N (x)| ≤ ∑ k=N+1 AC k = AC N /(1 -C).
In some cases, further error analysis can be useful to improve the quality (tightness) of the ball enclosures.

For arbitrary-precision evaluation, we wish to minimize the computational complexity as a function of the precision p. The complexity is often measured by counting arithmetic operations. The actual time complexity must account for the fact that arithmetic operations have a bit complexity of O(p) (where the O notation ignores logarithmic factors). In some situations, it is better to use a model of complexity that distinguishes between "scalar" arithmetic operations (such as addition of two p-bit numbers or multiplication of a p-bit number by a small integer) and "nonscalar" arithmetic operations (such as multiplication of two general p-bit numbers).

The exponential function

We illustrate these principles with a commonly used algorithm to compute the exponential function e x of a real argument x to p-bit precision.

• Argument reduction. We first use e x = 2 n e t with t = xn log(2) and n = ⌊x/ log(2)⌋ which ensures that t ∈ [0, log(2)). At this point, the usual Taylor series e t = 1 + t + 1 2 t 2 + . . . does not suffer from cancellation, and we only need O(p/ log p) terms for a relative error of 2 -p independent of the initial size of |x|. As a second argument reduction step, we write e t = (e u ) 2 r with u = t/2 r , which reduces the number N of needed Taylor series terms to O(p/r).

Balancing N = O(p/r) against the number r of squarings needed to reconstruct e t from e u , it is optimal to choose r ≈ p 0.5 . This gives an algorithm for e x requiring O(p 0.5 ) arithmetic operations on p-bit numbers, which translates to a time complexity of O(p 1.5 ).

• Series evaluation. As an additional improvement, we can exploit the structure of the Taylor series of the exponential function. For example, ∑ 8 k=0 1 k! x k can be evaluated as

1 + x + 1 2 x 2 + 1 3 x 3 1 + 1 4 x + 1 5 x 2 + 1 6 x 3 1 + 1 7 x + 1 8 x 2 (1.1)
where we have extracted the power x 3 repeatedly and used the fact that the ratios between successive coefficients are small integers. As a result, we only need four nonscalar multiplications involving x (to compute x 2 , x 3 , and for the two multiplications by x 3 ), while the remaining operations are scalar divisions. With further rewriting, the scalar divisions can be replaced by even cheaper scalar multiplications.

In general, to evaluate a polynomial of degree N with scalar coefficients at a nonscalar argument x, we can compute x 2 , . . . , x m once and then use Horner's rule with respect to x m , for m ≈ N 0.5 , which reduces the total number of nonscalar multiplications to about 2N 0.5 [START_REF] Paterson | On the number of nonscalar multiplications necessary to evaluate polynomials[END_REF]. This trick is sometimes called rectangular splitting. To motivate this terminology, picture the terms of the polynomial laid out as a matrix with m columns and N/m rows.

In view of this improvement to the series evaluation, it turns out to be more efficient in practice to choose the tuning parameter r used for argument reduction e t = (e t/2 r ) 2 r slightly smaller, say about r ≈ p 0.4 for realistic p. The algorithm combining optimal argument reduction with rectangular splitting for evaluation of elementary functions such as e x is due to Smith [START_REF] Smith | Efficient multiple-precision evaluation of elementary functions[END_REF].

There are asymptotically faster algorithms that permit evaluating elementary functions using only O(log p) arithmetic operations (that is, in O(p) time), for instance based on the AGM (discussed in section 1.5 for computing elliptic integrals), but Smith's algorithm is more efficient in practice for moderate p, and in some situations still wins for p as large as 10 5 .

Modular forms and functions

A modular transformation g is a linear fractional transformation on the upper half plane H = {τ : C : Im(τ) > 0} of the form g(τ) = (aτ + b)/(cτ + d) where a, b, c, d are integers with adbc = 1. We can assume that c ≥ 0. The group of modular transformations (known as the modular group) can be identified with the projective special linear group PSL(2, Z), where g is represented by the matrix a b c d and composition corresponds to matrix multiplication.

A modular form (of weight k) is a holomorphic function on H satisfying the functional equation f (g(τ)) = (cτ + d) k f (τ) for every modular transformation g, with the additional technical requirement of being holomorphic as τ → i∞. A modular form of weight k = 0 must be a constant function, but nontrival solutions of the above functional equation are possible if we allow poles. A meromorphic function on H satisfying f (g(τ)) = f (τ) is called a modular function.

Every modular form or function is periodic with f (τ + 1) = f (τ) and has a Fourier series (or q-series)

f (τ) = ∞ ∑ n=-m a n q n , q 2πiτ (1.2)
where m = 0 in the case of a modular form. The fundamental tool in numerical evaluation of modular forms and functions is to evaluate a truncation of such a qseries. Since τ has positive imaginary part, the quantity q always satisfies |q| < 1, and provided that an explicit bound for the coefficient sequence a n is known, tails of (1.2) are easily bounded by a geometric series.

Argument reduction

The q-series (1.2) always converges, but the convergence is slow for τ close to the real line where |q| ≈ 1. However, we can always find a modular transformation g that moves τ to the fundamental domain {τ ∈ H

: |τ| ≥ 1, | Re(τ)| ≤ 1 2 } (see Fig. 1.1
). This ensures |q| ≤ e -π √ 3 ≈ 0.00433 which makes the convergence extremely rapid. Technically, the fundamental domain does not include half of the boundary (meaning that H is tiled by copies of the fundamental domain under the action of the modular group), but this does not matter for the algorithm. In fact, it is sufficient to put τ within some small distance ε of the fundamental domain, and this relaxation is especially useful in ball arithmetic since a ball may overlap with the boundary.

The well-known algorithm to construct g (see Cohen [START_REF] Cohen | A course in computational algebraic number theory[END_REF]Algorithm 7.4.2]) repeatedly applies the generators τ → τ + 1 and τ → -1/τ of the modular group:

1. Set g ← 1 0 0 1 . 2. Set τ ← τ + n, g ← 1 n 0 1 g where n = -Re(τ) + 1 2 . 3. If |τ| < 1 -ε, set τ ← -1/τ and g ← 0 -1
1 0 g and go to step 2; otherwise, stop. Exact integer operations should be used for the matrices, but we can perform all operations involving τ in this algorithm using heuristic floating-point approximations, starting with the midpoint of the initial ball. Once g has been constructed, we evaluate g(τ) and f (g(τ)) as usual in ball arithmetic.

Indeed, it is important to construct the transformation matrix g separately and then apply the functional equation for the modular form in a single step rather than applying the generating transformations iteratively in ball arithmetic. This both serves to minimize numerical instability and to optimize performance. The precision needed to construct g only depends on the size of the entries of g and not on the precision for evaluating f (τ). If τ ≈ 2 -p i, we only need about 2p bits to construct g even if the overall precision is thousands of digits. The implementation in Arb uses virtually costless machine floating-point arithmetic to construct g when 53-bit arithmetic is sufficient, switching to arbitrary-precision arithmetic only when necessary.

Standard functions

Several commonly-used modular forms and functions are implemented in Arb. The basic building blocks are the Dedekind eta function

η(τ) = e πiτ/12 ∞ ∑ n=-∞ (-1) n q (3n 2 -n)/2 , q = e 2πiτ , (1.3)
and the theta constants θ j ≡ θ j (τ),

θ 2 (τ) = e πiτ/4 ∞ ∑ n=-∞ q n(n+1) , θ 3 (τ) = ∞ ∑ n=-∞ q n 2 , θ 4 (τ) = ∞ ∑ n=-∞ (-1) n q n 2 (1.4)
in which (as a potential source of confusion) q = e πiτ . It is useful to represent other modular forms in terms of these particular functions since their q-series are extremely sparse (requiring only O(p 0.5 ) terms for p-bit accuracy, which leads to O(p 1.5 ) bit complexity) and only have coefficients ±1. We give a few examples of derived functions:

• Modular functions are precisely the rational functions of the j-invariant

j(τ) = 32 (θ 8 2 + θ 8 3 + θ 8 4 ) 3 (θ 2 θ 3 θ 4 ) 8 , j aτ + b cτ + d = j(τ). (1.5)
• The modular discriminant is a modular form of weight 12, given by

∆ (τ) = η(τ) 24 , ∆ aτ + b cτ + d = (cτ + d) 12 ∆ (τ). (1.6)
• Eisenstein series are modular forms of weight 2k for k ≥ 2, given by

G 2k (τ) = ∑ m 2 +n 2 =0 1 (m + nτ) 2k , G 2k aτ + b cτ + d = (cτ + d) 2k G 2k (τ) (1.7)
where we compute G 4 (τ), G 6 (τ) via theta constants using

G 4 (τ) = π 4 90 θ 8 2 + θ 8 3 + θ 8 4 , G 6 (τ) = π 6 945 -3θ 8 2 (θ 4 3 + θ 4 4 ) + θ 12 3 + θ 12 4
and obtain the higher-index values using recurrence relations.

The Dedekind eta function itself transforms as

η aτ + b cτ + d = ε(a, b, c, d) √ cτ + d η(τ) (1.8)
where ε(a, b, c, d) = exp(πiR/12) is a 24th root of unity. The integer R (mod 24) can be computed using Kronecker symbols [25, section 74]. The modular transformations for theta constants are a special case of the formulas for theta functions given below in section 1.4.1. However, we avoid using these transformations directly when computing the functions (1.6), (1.7) and (1.5): it is better to apply the simpler argument reductions for the top-level functions and then evaluate the series expansions (1.3) or (1.4) when τ is already reduced.

Fast evaluation of q-series

The powers of q appearing in (1.3) and (1.4) are easily generated using two multiplications per term since the exponents are successive values of quadratic polynomials.

The cost can nearly be halved using short addition sequences [START_REF] Enge | Short addition sequences for theta functions[END_REF]Algorithm 2]. The cost can be reduced even further by combining addition sequences with rectangular splitting [11, section 5]. Here, the idea is to factor out some power q m as in (1.1), but m must be chosen in a particular way -for example, in the case of the theta series ∑ ∞ n=1 q n 2 , m is chosen so that there are few distinct quadratic residues modulo m. In Arb, these optimizations save roughly a factor four over the naive algorithm.

Elliptic and theta functions

An elliptic function with respect to a lattice in C with periods ω 1 , ω 2 is a meromor- phic function satisfying f (z + mω 1 + nω 2 ) = f (z) for all z ∈ C and all m, n ∈ Z. By making a linear change of variables, we can assume that ω 1 = 1 and ω 2 = τ ∈ H. The elliptic functions with a fixed lattice parameter τ form a field, which is generated by the Weierstrass elliptic function

℘(z, τ) = 1 z 2 + ∑ n 2 +m 2 =0 1 (z + m + nτ) 2 - 1 (m + nτ) 2 (1.9)
together with its z-derivative ℘ ′ (z, τ).

The building blocks for elliptic functions are the Jacobi theta functions

θ 1 (z, τ) = ∞ ∑ n=-∞ e πi[(n+ 1 2 ) 2 τ+(2n+1)z+n- 1 2 ] = 2q 4 ∞ ∑ n=0 (-1) n q n(n+1) sin((2n + 1)πz) = -iq 4 ∞ ∑ n=0 (-1) n q n(n+1) (w 2n+1 -v 2n+1 ), θ 2 (z, τ) = ∞ ∑ n=-∞ e πi[(n+ 1 2 ) 2 τ+(2n+1)z] = 2q 4 ∞ ∑ n=0 q n(n+1) cos((2n + 1)πz) = q 4 ∞ ∑ n=0 q n(n+1) (w 2n+1 + v 2n+1 ), (1.10)
θ 3 (z, τ) = ∞ ∑ n=-∞ e πi[n 2 τ+2nz] = 1 + 2 ∞ ∑ n=1 q n 2 cos(2nπz) = 1 + ∞ ∑ n=1 q n 2 (w 2n + v 2n ), θ 4 (z, τ) = ∞ ∑ n=-∞ e πi[n 2 τ+2nz+n] = 1 + 2 ∞ ∑ n=1 (-1) n q n 2 cos(2nπz) = 1 + ∞ ∑ n=1 (-1) n q n 2 (w 2n + v 2n ),
where q = e πiτ , q 4 = e πiτ/4 , w = e πiz , v = w -1 . The theta functions are quasielliptic functions of z, having period or half-period 1 and quasiperiod τ (a shift by τ introduces an exponential prefactor). With z = 0, the theta functions θ 2 , θ 3 , θ 4 reduce to the corresponding theta constants, while θ 1 (0, τ) = 0 identically. Arb provides a complete implementation of the Jacobi theta functions themselves as well as the Weierstrass elliptic function which is computed as

℘(z, τ) = π 2 θ 2 2 (0, τ)θ 2 3 (0, τ) θ 2 4 (z, τ) θ 2 1 (z, τ) - π 2 3 θ 4 2 (0, τ) + θ 4 3 (0, τ) . (1.11)
For all these functions, Arb also allows computing an arbitrary number of zderivatives. Derivatives are handled by working with θ j (z + x, τ) and ℘(z + x, τ) as elements of C[[x]] (truncated to some length O(x D )), using power series arithmetic.

Arb also implements the quasielliptic Weierstrass zeta and sigma functions ζ (z, τ) and σ (z, τ) as well as the lattice invariants g 2 , g 3 (which are essentially Eisenstein series) and lattice roots 4z 3

-g 2 z -g 3 = 4(z -e 1 )(z -e 2 )(z -e 3 ) arising in the differential equation [℘ ′ (z, τ)] 2 = 4[℘(z, τ)] 3 -g 2 ℘(z, τ) -g 3 .
The inverse Weierstrass elliptic function is also available; see section 1.6.

The Jacobi elliptic functions sn, cn, . . . are not currently part of the library, but users can compute them via theta functions using formulas similar to (1.11).

Argument reduction

As the first step when computing theta functions or elliptic functions, we reduce τ to the fundamental domain using modular transformations. This gives us a new lattice parameter τ ′ and a new argument z ′ . As a second step, we reduce z ′ modulo τ ′ , giving an argument z ′′ with smaller imaginary part (it is not necessary to reduce z ′ modulo 1 since this is captured by the oscillatory part of exponentials). We can then compute θ j (z ′′ , τ ′ ) using the theta series (1.10).

These steps together ensure that both |q| and max(|w|, |w| -1 ) will be small. It is important to perform both transformations. Consider τ = 0.07 + 0.003i and z = 3.14 + 2.78i: without the modular transformation, the direct series evaluation would use 3710 terms for machine precision. 4 With the modular transformation alone, it would use 249 terms. With both reductions, only 6 terms are used! Depending on the arguments, the numerical stability may also be improved substantially.

Modular transformations have the effect of permuting the theta functions and introducing certain exponential prefactors. It is easy to write down the transformations for the generators τ + 1, -1/τ, but the action of a composite transformation involves a certain amount of bookkeeping. The steps have been worked out by Rademacher [25, chapter 10]. We reproduce the formulas below. 5We wish to write a theta function with lattice parameter τ in terms of a theta function with lattice parameter τ ′ = g(τ), given some g = a b c d ∈ PSL(2, Z). For j = 0, 1, 2, 3, there are R j , S j ∈ Z depending on g such that

θ 1+ j (z, τ) = exp(πiR j /4) • A • B • θ 1+S j (z ′ , τ ′ ) (1.12)
where if c = 0,

z ′ = z, A = 1, B = 1, (1.13)
and otherwise (if c > 0),

z ′ = -z cτ + d , A = i cτ + d , B = exp -πic z 2 cτ + d . (1.14)
We always have B = 1 when computing theta constants which have z = 0. The parameters R j , S j are computed from g as follows. If c = 0, we have θ j (z, τ) = exp(-πib/4)θ j (z, τ + b) for j = 1, 2, whereas θ 3 and θ 4 remain un- changed when b is even and swap places with each other when b is odd. For the c > 0 case, it is helpful to define the function θ m,n (z, τ) for m, n ∈ Z by θ 0,0 (z, τ) = θ 3 (z, τ), θ 0,1 (z, τ) = θ 4 (z, τ), (1.15)

θ 1,0 (z, τ) = θ 2 (z, τ), θ 1,1 (z, τ) = iθ 1 (z, τ), θ m+2,n (z, τ) = (-1) n θ m,n (z, τ) θ m,n+2 (z, τ) = θ m,n (z, τ).
With this notation, we have

θ 1 (z, τ) = ε 1 ABθ 1 (z ′ , τ ′ ), θ 2 (z, τ) = ε 2 ABθ 1-c,1+a (z ′ , τ ′ ), (1.16) θ 3 (z, τ) = ε 3 ABθ 1+d-c,1-b+a (z ′ , τ ′ ), θ 4 (z, τ) = ε 4 ABθ 1+d,1-b (z ′ , τ ′ )
where ε k is an 8th root of unity. If we denote by ε(a, b, c, d) = exp(πiR(a, b, c, d)/12) the 24th root of unity in the transformation (1.8) of the Dedekind eta function, then

ε 1 (a, b, c, d) = exp(πi[R(-d, b, c, -a) + 1]/4), ε 2 (a, b, c, d) = exp(πi[-R(a, b, c, d) + (5 + (2 -c)a)]/4), ε 3 (a, b, c, d) = exp(πi[-R(a, b, c, d) + (4 + (c -d -2)(b -a))]/4), (1.17)
ε 4 (a, b, c, d) = exp(πi[-R(a, b, c, d) + (3 -(2 + d)b)]/4).
Finally, to reduce z ′ , we compute n = ⌊Im(z ′ )/ Im(τ ′ ) + 1/2⌋ and set z ′′ = z ′nτ ′ . In this step, all theta functions pick up a prefactor exp(πi[-τn 2 -2nz]) (this data may be combined with B) while θ 1 and θ 2 pick up the additional prefactors (-1) n (this data may be combined with R j ). When computing z-derivatives of theta functions, the same formulas are applied in power series arithmetic. That is, if the initial argument consists of the formal power series z + x, then the scaling factor -1/(cτ + d) is applied coefficient by coefficient, while B = B 0 + B 1 x + . . . is obtained by squaring the power series z + x, scaling, and then evaluating a power series exponential.

As with modular forms, the transformations should be applied at the highest possible level. For example, when computing a quotient of two theta functions of the same z, τ, the prefactors A and B in (1.12) cancel out (and the leading roots of unity possibly also simplify). We should then simplify the expression symbolically and avoid computing A and B altogether, since this both saves time and improves numerical stability in ball arithmetic (in particular, e f (z) /e f (z) evaluated in ball arithmetic will not give 1 but rather a ball which can be extremely wide).

Since the description of the algorithm given above is quite terse, the reader may find it helpful to look at the code in Arb to see the concrete steps.

Theta function series evaluation

Algorithm 1 implements the expansions (1.10), with the optimization that we combine operations to save work when computing all four functions and their derivatives simultaneously (a single theta function could be computed slightly faster, but computing all four functions is barely more work than it would be to compute a pair containing either θ 1 or θ 2 and either θ 3 or θ 4 ). This is essentially the algorithm used in Arb for z = 0, while more optimized code is used for theta constants.

The main index k runs over the terms in the following order:

θ 1 , θ 2 q 0 (w 1 ± w -1 ) k = 0 θ 3 , θ 4 q 1 (w 2 ± w -2 ) k = 1 θ 1 , θ 2 q 2 (w 3 ± w -3 ) k = 2 θ 3 , θ 4 q 4 (w 4 ± w -4 ) k = 3 θ 1 , θ 2 q 6 (w 5 ± w -5 ) k = 4 θ 3 , θ 4 q 9 (w 6 ± w -6 ) k = 5 θ 1 , θ 2 q 12 (w 7 ± w -7 )
The algorithm outputs the range of scaled derivatives θ (r) j (z, τ)/r! for 0 ≤ r < D. The term of index k in the main summation picks up a factor ±(k + 2) r from r-fold differentiation of w k+2 . Another factor (πi) r /r! is needed to convert to a z-derivative and a power series coefficient, but we postpone this to a single rescaling pass at the end of the computation. In the main summation, we write the even cosine terms as w 2n + w -2n , the odd cosine terms as w(w 2n + w -2n-2 ), and the sine terms as w(w 2n -w -2n-2 ), postponing a multiplication by w for θ 1 and θ 2 until the end, so that only even powers of w and w -1 are needed.

For some integer N ≥ 1, the summation is stopped just before term k = N. Let Q = |q|, W = max(|w|, |w -1 |), E = ⌊(N + 2) 2 /4⌋ and F = ⌊(N + 1)/2⌋ + 1. The error of the zeroth derivative can be bounded as

2Q E W N+2 1 + Q F W + Q 2F W 2 + . . . = 2Q E W N+2 1 -Q F W (1.18)
provided that the denominator 1 -Q F W is positive. For the r-th derivative, including the factor (k + 2) r gives the error bound

2Q E W N+2 (N + 2) r 1 + Q F W (N + 3) r (N + 2) r + Q 2F W 2 (N + 4) r (N + 2) r + . . . (1.19)
which by the inequality (1 + m/(N + 2)) r ≤ exp(mr/(N + 2)) can be bounded as

2Q E W N+2 (N + 2) r 1 -Q F W exp(r/(N + 2)) , (1.20)
again valid when the denominator is positive.

Algorithm 1 Computation of Jacobi theta functions (using series evaluation)

Require: z, τ ∈ C with Im(τ) > 0 (can be arbitrary, but should be reduced for best performance), integer D ≥ 1 to output the D first terms in the Taylor expansions with respect to z, precision p Ensure:

θ j = [α 0 , . . ., α D-1 ] represents θ j (z + x, τ) = α 0 + α 1 x + . . . + α D-1 x D-1 , for 1 ≤ j ≤ 4 1: q 4 ← e πiτ/4 ; q ← q 4 4 ; w ← e πiz ; v ← w -1 ; Q ← |q|; W ← max(|w|, |v|) 2: Choose N with E = ⌊(N + 2) 2 /4⌋ and F = ⌊(N + 1)/2⌋ + 1 such that Q E W N+2 < 2 -p and α = Q F W exp(r/(N + 2)) < 1 3: for 0 ≤ r < D do ε[r] ← 2Q E W N+2 (N + 2) r /(1 -α) end for ⊲ Error bounds 4: w ← [1, w 2 , w 4 , . . ., w 2K-2 ]; v ← [1, v 2 , v 4 , . . ., v 2K ] for K = ⌊(N + 3)/2⌋ ⊲ Precompute powers 5: θ 1 ← [0, . . ., 0]; θ 2 ← [0, . . ., 0]; θ 3 ← [0, . . ., 0]; θ 4 ← [0, . . ., 0]; ⊲ Arrays of length D 6: for 0 ≤ k < N do 7: m ← ⌊(k + 2) 2 /4⌋; n ← ⌊k/2⌋ + 1 8:
Compute q m ⊲ Use addition sequence [11, Alg. 2] to build q m from previous powers. 9:

t ← (w[n] + v[n + (k mod 2)])q m 10: u ← (w[n] -v[n + (k mod 2)])q m ⊲ Skip when k mod 2 = 0 if D = 1. 11:
if k mod 2 = 0 then 12:

for 0 ≤ r < D do 13:

if r mod 2 = 0 then 14:

if r = 0 then t ← 4n 2 t end if 15:

θ 3 [r] ← θ 3 [r] + t; θ 4 [r] ← θ 4 [r] + (-1) ⌊(k+2)/2⌋ t 16: else 17: if r = 1 then u ← 2nu else u ← 4n 2 u end if 18: θ 3 [r] ← θ 3 [r] + u; θ 4 [r] ← θ 4 [r] + (-1) ⌊(k+2)/2⌋ u 19:
end if 20:

end for 21: else 22:

for 0 ≤ r < D do 23:

if r mod 2 = 0 then 24:

θ 1 [r] ← θ 1 [r] + (-1) ⌊(k+1)/2⌋ u; θ 2 [r] ← θ 2 [r] + t 25: else 26: θ 1 [r] ← θ 1 [r] + (-1) ⌊(k+1)/2⌋ t; θ 2 [r] ← θ 2 [r] + u 27:
end if 28:

t ← (2n + 1)t; u ← (2n + 1)u 29:

end for 30:

end if 31: end for 32: for 0 ≤ r < D do 33:

θ 1 [r] ← θ 1 [r]w + (w -(-1) r v)
⊲ Adjust power of w and add leading terms 34:

θ 2 [r] ← θ 2 [r]w + (w + (-1) r v) 35: for 1 ≤ j ≤ 4 do θ j [r] ← θ j [r] + [±ε[r]] + [±ε[r]]i end for ⊲ Add error bounds 36: C ← (πi) r /r! ⊲ Final scaling factors 37: θ 1 [r] ← -iq 4 Cθ 1 [r]; θ 2 [r] ← q 4 Cθ 2 [r]; θ 3 [r] ← Cθ 3 [r]; θ 4 [r] ← Cθ 4 [r] 38: end for 39: θ 3 [0] ← θ 3 [0] + 1; θ 4 [0] ← θ 4 [0] + 1 ⊲ Add leading terms
The time complexity of the algorithm is O(p 1.5 ) (with all inputs besides p fixed). By employing fast Fourier transforms cleverly, the complexity of evaluating theta functions from their series expansions can be reduced to O(p 1.25 ), but that method is only faster in practice for p exceeding 200 000 bits [START_REF] Nogneng | On the evaluation of some sparse polynomials[END_REF]. See also section 1.8.1 below concerning methods that are even faster asymptotically.

Complete elliptic integrals and the AGM

Complete elliptic integrals arise in period relations for elliptic functions. The complete elliptic integral of the first kind is

K(m) = 1 2 π 2 F 1 ( 1 2 , 1 2 , 1, m) and the complete elliptic integral of the second kind is E(m) = 1 2 π 2 F 1 (-1 2 , 1 2 , 1, m)
where 2 F 1 denotes the Gauss hypergeometric function, defined for |z| < 1 by

2 F 1 (a, b, c, z) = ∞ ∑ k=0 (a) k (b) k (c) k z k k! , (x) k = x(x + 1) • • •(x + k -1) (1.21)
and elsewhere by analytic continuation with the standard branch cut on [1, ∞).

The 2 F 1 function can be computed efficiently for any z ∈ C using a combination of argument transformations, analytic continuation techniques, and series expansions (where the rectangular splitting trick (1.1) and other accelerations methods are applicable). A general implementation of 2 F 1 exists in Arb [START_REF] Johansson | Computing hypergeometric functions rigorously[END_REF]. However, it is more efficient to compute the complete elliptic integrals by exploiting their connection with the arithmetic-geometric mean (AGM) described below.

A third complete elliptic integral Π (n, m) is also encountered, but this is a more complicated function that is not a special case of 2 F 1 , and we handle it later in terms of an incomplete integral without using a dedicated algorithm for the complete case.

The arithmetic-geometric mean M(x, y) of two nonnegative real numbers x, y is defined as the common limit of the sequences

a n+1 = a n + b n 2 , b n+1 = a n b n (1.22)
with initial values a 0 = x, b 0 = y. In different words, the AGM can be computed by repeatedly applying the functional equation M(x, y) = M (x + y)/2, √ xy . It is a well known fact that each step of the AGM iteration roughly doubles the number of accurate digits in the approximation a n ≈ b n ≈ M(x, y), so it only costs O(log(p)) arithmetic operations to compute the AGM with an accuracy of p bits, resulting in a bit complexity of O(p).

For complex x, y, defining the AGM becomes more difficult since there are two possible choices for the square root in each step of the iteration, and these choices lead to different limits. However, it turns out that there is an "optimal" choice which leads to a well-defined and useful extension of the AGM to complex variables. We rely on several properties of this function proved in earlier work [START_REF] Cox | The arithmetic-geometric mean of Gauss[END_REF][START_REF] Cremona | The complex AGM, periods of elliptic curves over C and complex elliptic logarithms[END_REF][START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF].

With complex variables, it is convenient to work with the univariate function M(z) = M(1, z), with a branch cut on (-∞, 0]. The general case can be recovered as M(x, y) = xM(1, y/x). The complete elliptic integrals (with the conventional branch cuts on [1, ∞)) are now given by

K(m) = π 2M( √ 1 -m) , E(m) = (1 -m)(2mK ′ (m) + K(m)). (1.23)
For implementing the function M(z), we can further assume that Re(z) ≥ 0 holds. If this is not the case, we first apply the functional equation M(z) = (z + 1)M(u)/2 where u = √ z/(z + 1). The correct square root in the AGM iteration is now always equal to √ a n √ b n , written in terms of the usual principal square root function. This can be computed as

√ a n b n , i √ -a n b n , -i √ -a n b n , √ a n √
b n respectively if both a n and b n have positive real part, nonnegative imaginary part, nonpositive imaginary part, or otherwise. When the iteration is executed in ball arithmetic, the computed balls may end up containing points with negative real part, but this just inflates the final result and does not affect correctness.

The iteration should be terminated when a n and b n are close enough. For positive real variables, we can simply take lower and upper bounds to get a correct enclosure. For complex variables, it can be shown [8, p. 87

] that |M(z) -a n | ≤ |a n -b n | if Re(z) ≥ 0,
giving a convenient error bound. However, instead of running the AGM iteration until a n and b n agree to p bits, it is slightly better to stop when they agree to about p/10 bits and end with a Taylor series. With t = (ab)/(a + b), we have which is valid at least when |t| < 1 and a, b have nonnegative real part, and where the tail (. . .) is bounded by ∑ ∞ k=10 |t| k /64. This algorithm follows the pattern of argument reduction and series evaluation. However, unlike the elementary functions and the incomplete elliptic integrals described below, there is no asymptotic benefit to using more terms of the series. The quadratic convergence of the AGM iteration is so rapid that we only get a speedup from trading O(1) of the O(log p) square roots for lower-overhead multiplications. Although there is no asymptotic improvement, the order-10 series expansion nevertheless gives a significant speedup up to a few thousand bits.

M(a, b) = (a + b)π 4K(t 2 ) , π 4K(t 2 ) = 1 2 -1 8 t 2 -
For computing the second elliptic integral E(m) or the first derivative M ′ (z) of the AGM, a simple method is to use a central finite difference to compute (M(z), M ′ (z)) ≈ (M(z + h) + M(zh))/2, (M(z + h) -M(zh))/(2h). This requires two evaluations at 1.5 times increased precision, which is about three times as expensive as evaluating M once. Error bounds can be obtained using the Cauchy integral formula and the inequality |M(z)| ≤ max(1, |z|) which is an immediate consequence of the AGM iteration. This method has been implemented in Arb. A more efficient method is to compute E(m) using an auxiliary sequence related to the AGM iteration, which also generalizes to computing Π (n, m) [23, 19.8.6 and 19.8.7]. This method has not yet been implemented in Arb since it requires some additional error analysis and study for complex variables.

Higher derivatives of the arithmetic-geometric mean or the complete elliptic integrals can be computed using recurrence relations. Writing W (z) = 1/M(z) and

W (z + x) = ∑ ∞ k=0 c k x k , we have -2z(z 2 -1)c 2 = (3z 2 -1)c 1 + zc 0 , -(k + 2)(k + 3)z(z 2 -1)c k+3 = (k + 2) 2 (3z 2 -1)c k+2 + (3k(k + 3) + 7)zc k+1 + (k + 1) 2 c k when z = 1 and -(k + 2) 2 c k+2 = (3k(k + 3) + 7)c k+1 + (k + 1) 2 c k when z = 1.

Incomplete elliptic integrals

A general elliptic integral is an integral of the form b a R(t, P(t))dt where R is a bivariate rational function and P is a cubic or quartic polynomial without repeated roots. It is well known that any elliptic integral can be expressed in terms of integrals of rational functions and a finite set of standard elliptic integrals.

Such a set of standard integrals is given by the Legendre incomplete elliptic integrals of the first, second and third kind

F(φ , m) = φ 0 dt 1 -m sin 2 t , E(φ , m) = φ 0 1 -m sin 2 t dt, (1.25)
Π (n, φ , m) = φ 0 dt (1 -n sin 2 t) 1 -m sin 2 t . (1.26)
The complete elliptic integrals are the special cases

E(m) = E(π/2, m), K(m) = F(π/2, m), and Π (n, m) = Π (n, π/2, m).
The definitions for complex variables do not appear to be standardized in the literature, but following the conventions used in Mathematica [START_REF] Research | The Wolfram Functions Site -Elliptic Integrals[END_REF], we may fix an interpretation of (1.25)-(1.26) on -π/2 ≤ Re(φ ) ≤ π/2 and use the quasiperiodic extensions 6The Legendre forms of incomplete elliptic integrals are widely used by tradition, but they have some practical drawbacks. Since they have a complicated (and not standardized) complex branch structure, transforming their arguments using functional equations or using them to represent other functions often requires making complicated case distinctions. As a result, it is cumbersome both to compute the functions themselves and to apply them, outside of a restricted parameter range.

F(φ + kπ, m) = 2kK(m) + F(φ , m), E(φ + kπ, m) = 2kE(m) + E(φ , m), Π (n, φ + kπ, m) = 2kΠ (n, m) + Π (n, φ , m) for k ∈ Z.
We remark that F and E can be expressed in terms of the Appell F 1 hypergeometric function of two variables, while Π can be expressed in terms of the three-variable Lauricella hypergeometric function F

(3) D , generalizing the 2 F 1 representations for the complete integrals. Such formulas are by themselves mainly useful when the hypergeometric series converge, and provide no insight into the analytic continuations.

In the 1960s, Carlson introduced an alternative set of standard elliptic integrals in which all or some of the variables are symmetric [START_REF] Carlson | Numerical computation of real or complex elliptic integrals[END_REF]. The Carlson incomplete elliptic integrals are

R F (x, y, z) = 1 2 ∞ 0 dt (t + x)(t + y)(t + z) (1.27)
and

R J (x, y, z, p) = 3 2 ∞ 0 dt (t + p) (t + x)(t + y)(t + z) (1.28)
together with three special cases R D (x, y, z) = R J (x, y, z, z), R C (x, y) = R F (x, y, y), and

R G (x, y, z) = zR F (x, y, z) - 1 3 (x -z)(y -z)R D (x, y, z) + √ x √ y √ z . (1.29)
The Carlson forms have several advantages over the Legendre forms. Symmetry unifies and simplifies the argument transformation formulas, and the Carlson forms also have a simpler complex branch structure, induced by choosing the branch of the square root in (1.27) and (1.28) to extend continuously from +∞. We can define and compute the Legendre forms from the Carlson forms using

F(φ , m) = sR F (x, y, 1), E(φ , m) = sR F (x, y, 1) -1 3 ms 3 R D (x, y , 1), (1.30)
Π (n, φ , m) = sR F (x, y, 1) + 1 3 ns 3 R J (x, y, 1, p)
on -π/2 ≤ Re(φ ) ≤ π/2 (with the quasiperiodic extensions elsewhere) where x = c 2 , y = 1ms 2 , p = 1ns 2 and s = sin(φ ), c = cos(φ ). This is the approach used to implement the Legendre forms in Arb. The Carlson forms themselves are also exposed to users. Formulas for other elliptic integrals can be found in [START_REF] Carlson | Numerical computation of real or complex elliptic integrals[END_REF].

Elliptic integrals can also be characterized as the inverse functions of elliptic functions. For example, the inverse of the Weierstrass elliptic function, which by definition satisfies ℘(℘ -1 (z, τ), τ) = z, is given by the elliptic integral

℘ -1 (z, τ) = 1 2 ∞ z dt (t -e 1 )(t -e 2 )(t -e 3 ) = R F (z -e 1 , z -e 2 , z -e 3 ). (1.31)
The implementation in Arb simply computes the lattice roots e 1 , e 2 , e 3 using theta constants and then calls R F . The inverses of Jacobi's elliptic functions can be computed similarly, but at this time they are not implemented in Arb.

Carlson gives algorithms for computing R F and R J using argument reduction and series evaluation [START_REF] Carlson | Numerical computation of real or complex elliptic integrals[END_REF]. The algorithm for R F is correct for all complex x, y, z (Carlson restricts to the cut plane with (-∞, 0) removed, but it is clear that the algorithm also works on the branch cut by continuity). The algorithm for R J is not correct for all values of the variables, but it is always correct when computing R D (otherwise, a sufficient condition is that x, y, z have nonnegative real part while p has positive real part). Carlson also provides modifications of the algorithms for computing the Cauchy principal values of the integrals.

We will now describe Carlson's algorithm for R F and adapt it to the setting of arbitrary precision and ball arithmetic. The algorithm given in [START_REF] Carlson | Numerical computation of real or complex elliptic integrals[END_REF] for R J and R D works analogously, but we do not reproduce all the steps here since the formulas would be too lengthy (the code in Arb can be consulted for concrete details).

Argument reduction

Argument reduction for R F uses the symmetric "duplication formula"

R F (x, y, z) = R F x + λ 4 , y + λ 4 , z + λ 4 (1.32) where λ = √ x √ y + √ y √ z + √ z √
x. Each application of (1.32) reduces the distance between the arguments by roughly a factor 4. The analogous formula for R J reads

R J (x, y, z, p) = 1 4 R J x + λ 4 , y + λ 4 , z + λ 4 , p + λ 4 + 6 d R C (1, 1 + e) (1.33)
where λ is defined as above and d, e are certain auxiliary terms (see [4, (24)-( 28)]).

The formulas (1.32) and (1.33) are iterated until all the parameters are close so that a series expansion can be used, as detailed in the next subsection. It is interesting to note the similarity between (1.32) and the AGM iteration, although the convergence rate of (1.32) only is linear. When computing R C or R D , some redundant operations in the reductions for R F and R J can be avoided. R C (x, y) can also be expressed piecewise using inverse trigonometric and hyperbolic functions. The special case R C (1,

1 + t) = atan( √ t)/ √ t = 2 F 1 (1, 1 2 , 3 2 , -t
) is particularly important, as it is needed in the evaluation of R J . This function is better computed via the inverse tangent function (or a direct Taylor series for small |t|) than by invoking Carlson's general method for R F .

Series expansions

Carlson's incomplete elliptic integrals are special cases of a multivariate hypergeometric function that may be written as

R -a (z 1 , . . . , z n ) = A -a ∞ ∑ N=0 (a) n ( 1 2 n) N T N (Z 1 , . . . , Z n ) (1.34)
where A = 1 n ∑ n j=1 z j , Z j = 1z j /A, and

T N (Z 1 , . . . , Z n ) = ∑ m 1 +...+m n =N m 1 ,...,m n ≥0 n ∏ j=1 ( 1 2 ) m j m j ! Z m j j . (1.35)
We have R F (x, y, z) = R -1/2 (x, y, z) and R J (x, y, z, p) = R -3/2 (x, y, z, p, p). The crucial property of this hypergeometric representation is that the expansion point is the arithmetic mean of the arguments. After sufficiently many argument reduction steps have been performed, we will have For the evaluation of (1.35), Carlson noted that it is more efficient to work with elementary symmetric polynomials E j = E j (Z 1 , . . . , Z n ) instead of the direct variables Z 1 , . . . , Z n , giving

T N (Z 1 , . . . , Z n ) = ∑ m 1 +2m 2 +...+nm n =N m 1 ,...,m n ≥0 (-1) M+N ( 1 2 ) M n ∏ j=1 E m j j m j ! . (1.37)
The key observation is that the symmetric choice of expansion variables Z j with respect to z 1 , . . . , z n implies that E 1 = 0, which eliminates most of the terms in (1.37). 7This dramatically reduces the amount of work to compute T N compared to (1.35).

For the R -1/2 series, there are (N + 1)(N + 2)/2 terms in (1.35) and roughly N/6 terms in (1.37); for example, if N = 8, there are 45 terms in the former and only two nonzero terms (with monomials E 1 E 2 2 and E 4 1 ) in the latter.

Series evaluation and balanced argument reduction

The argument reduction effectively adds 2B bits per step if a series expansion of order B is used. In other words, roughly p/(2B) argument reduction steps are needed for p-bit precision. Carlson suggests using a precomputed truncated series of order B = 6 or B = 8, which is a good default at machine precision and up to a few hundred bits. At higher precision, we make the observation that it pays off to vary B dynamically as a function of p and evaluate the series with an algorithm. Algorithm 2 gives pseudocode for the method implemented in Arb to compute R F using combined argument reduction and series evaluation. The subroutine RSum evaluates the series for R -1/2 truncated to an arbitrary order B using rectangular splitting combined with recurrence relations for the coefficients (one more optimization used in the implementation but omitted from the pseudocode is to clear denominators so that all coefficients are small integers). The exponents of E m

← (x + y + z)/3; (X,Y, Z) ← (1 -x/A, 1 -y/A, 1 -z/A); (E 2 , E 3 ) ← (XY -Z 2 , XY Z) 4: ε ← 2 ∑ ∞ k=B 9 8 max(|X|, |Y |, |Z|) k ⊲ Series error bound 5: procedure RSUM(E 2 , E 3 , B) ⊲ Compute R = ∑ B-1 N=0 ( 1 2 ) N /( 3 2 ) N T N 6: Precompute E k 2 for 2 ≤ k ≤ ⌊(B -1)/2⌋ 7: R ← 0; c 3 ← ( 1 2 ) ⌊(B-1)/3⌋ /(⌊(B -1)/3⌋)! 8:
for (m 3 ← ⌊(B -1)/3⌋; m 3 ≥ 0; m 3 ← m 3 -1) do 9:

if m 3 = ⌊(B -1)/3⌋ then 10: O(N 2 + p/N) operations, which is minimized by setting B ≈ p 1/3 ; this gives us an O(p 1.667 ) bit complexity algorithm for evaluating R F . With rectangular splitting for the series evaluation, the optimal B should be closer to B ≈ p 0.5 for moderate p. Timings in Arb show that expanding to order B = 2p 0.4 for real variables and B = 2.5p 0.4 for complex variables is nearly optimal in the range 10 ≤ p ≤ 10 6 bits. Compared to the fixed order B = 8, this results in a measured speedup of 1.5, 4.0, 11 and 31 times at a precision of 100, 1000, 10 000 and 100 000 decimal digits respectively. The algorithm for the R J series is essentially the same, except that the summation uses four nested loops instead of two to iterate over the exponents with 2m 2 + 3m 3 + 4m 4 + 5m 5 < B, with corresponding nested recurrence relations to update coefficients c 2 , c 3 , c 4 , c 5 (see the Arb source code for details). In this case, rectangular splitting is used to split the number of variables in half by precomputing a twodimensional array of powers E m 2 2 E m 3 3 and using the Horner scheme with respect to E 4 and E 5 . The speedup of combining rectangular splitting with optimal argument reduction is smaller for R J than for R F but still appreciable at very high precision.

c 3 ← c 3 • (2m 3 + 2)/(2m
= √ 2 + √ 3i, y = √ 3 + √ 5i, z = √ 5 + √ 7i,t = √ 7 + i/ √ 11.

Arb implementation benchmarks

Table 1.2 compares the performance of different functions implemented in Arb.

The complete elliptic integrals of the first and second kind are about as fast as the elementary functions at high precision due to the O(p) AGM algorithm. 8 The modular forms and functions which use O(p 1.5 ) algorithms with very low overhead are nearly as fast as complete elliptic integrals in practice.

Elliptic functions and Jacobi theta functions, also implemented with O(p 1.5 ) algorithms, are some 5-10 times slower than the special case of theta constants or modular forms. The incomplete elliptic integrals based on the R F function implemented with O(p 1.667 ) complexity have similar performance to the elliptic functions at moderate precision with a slight divergence becoming visible only at several thousand digits. Indeed, ℘(x,t) and ℘ -1 (x,t) have virtually identical performance although the algorithms are completely independent.

The incomplete elliptic integrals based on the R J function stand out as being noticeably slower than the other functions, as a result of the more complicated argument reduction and high-dimensional series expansion.

Other methods

Many numerical techniques apart from those covered in this text are useful in connection with elliptic functions and modular forms. Without going into detail, we sketch a few important ideas.

Quadratically convergent methods and Newton iteration

The algorithms described above for complete elliptic integrals have quasioptimal O(p) bit complexity owing to the quadratically convergent AGM iteration, while the algorithms for all other functions have O(p 1.5 ) or worse bit complexity. In fact, it is possible to compute general elliptic functions, modular forms and incomplete elliptic integrals with O(p) bit complexity using generalizations of the AGM iteration together with Newton's method for inverse functions. We have omitted these methods in the present work since they are more complicated, especially for complex variables, and not necessarily faster for p encountered in practice.

The asymptotically fast computation of modular forms and modular functions is discussed by Dupont [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF], and Labrande [START_REF] Labrande | Computing Jacobi's theta in quasi-linear time[END_REF] has given algorithms for general theta functions and elliptic functions. An important special case is the inverse Weierstrass elliptic function in the form of the elliptic logarithm, which can be computed using a simple AGM-type algorithm [START_REF] Cremona | The complex AGM, periods of elliptic curves over C and complex elliptic logarithms[END_REF]. For the Legendre incomplete elliptic integrals, algorithms based on the quadratic Landen transformations are classical and have been described in several other works; they have the disadvantage of involving trigonometric functions, not having a straightforward extension to complex variables, and in some regions suffering from precision loss.

Numerical integration

Direct numerical integration is a viable way to compute elliptic integrals. Numerical integration is generally slower than the more specialized algorithms already presented, but with a robust general-purpose integration algorithm, we can just plug in the formula for any particular elliptic integral. Specifying an explicit contour of integration also provides full control over branch cuts.

The double exponential or tanh-sinh quadrature method [START_REF] Bailey | High-precision numerical integration: Progress and challenges[END_REF][START_REF] Takahasi | Double exponential formulas for numerical integration[END_REF] is ideal for elliptic integrals since it is extremely simple and converges rapidly even if the integrand has algebraic singularities of unknown type at one or both endpoints. The quadrature error in the double exponential method can be estimated quite reliably using heuristics, and effective rigorous error bounds are also known [START_REF] Molin | Numerical Integration and L-Functions Computations[END_REF]. Alternatively, Gauss-Jacobi quadrature can be used for integrals with known algebraic singularities. Recently, Molin and Neurohr have studied use of both double exponential and Gauss-Jacobi quadrature with rigorous error bounds for integration of algebraic functions in the context of computing period matrices for hyperelliptic curves [START_REF] Molin | Computing period matrices and the Abel-Jacobi map of superelliptic curves[END_REF]. Rigorous numerical integration code also exists in Arb [START_REF] Johansson | Numerical integration in arbitrary-precision ball arithmetic[END_REF], but endpoint singularities require manual processing.

For integrals of smooth periodic functions, including integrals of analytic functions on closed circular contours, direct application of the trapezoidal rule is often the best choice. We conclude with the anecdote that Poisson already in the 1820s demonstrated the use of the trapezoidal rule to approximate the elliptic integral 1 2π 2π 0 1 -0.36 sin 2 (θ )dθ which is equal to 2 π E(0.36) in the Legendre notation. Poisson derived an error bound for the N-point trapezoidal approximation and showed that N = 16 gives an error less than 4.84 • 10 -6 for this integral (in fact, nine digits are correct). Due to symmetry, just three nontrivial evaluations of the integrand are required for this level of accuracy! Trefethen and Weideman [START_REF] Trefethen | The exponentially convergent trapezoidal rule[END_REF] discuss this example and provide a general error analysis for the trapezoidal rule applied to periodic functions.

  and then add the ball [±ε] or [±ε] + [±ε]i where ε is an upper bound for |R N (x)| = | ∑ ∞ k=N+1 t k (x)|. Such a bound is often readily obtained by comparison with a geometric series, i.e. if |t k

Fig. 1 . 1

 11 Fig. 1.1 The shaded region shows the canonical fundamental domain for the action of the modular group on the upper half plane.

z 1 ≈

 1 z 2 ≈ . . . ≈ z n which ensures |Z 1 |, . . . , |Z n | ≪ 1 and rapid convergence of the series. A trivial bound for the terms is |T N (Z 1 , . . . , Z n )| ≤ p(N) max(|Z 1 |, . . . , |Z n |) N , where p(N) denotes the number of partitions of N which is bounded by O(c N ) for any c > 1. An explicit calculation shows, for example, that the error when computing either R F or R J is bounded by 2A |Z 1 |, . . . , |Z n |) N (1.36) if the summation in (1.34) includes the terms of order N < B.

  ing in the series (Fig.1.2) are the lattice points m 2 , m 3 ∈ Z ≥0 with 2m 2 + 3m 3 < B: we compute powers of E 2 and then use Horner's rule with respect to E 3 .We now consider the choice of B. Since the series evaluation costs O(B 2 ) operations and the argument reduction costs O(p/B) operations, the overall cost is Algorithm 2 Computation of R F (x, y, z) 1: Choose series truncation order B optimally depending on the precision p 2: Apply argument reduction (1.32) until x, y, z are close (until ε ≈ 2 -p below) 3: A

Table 1 . 1

 11 Methods for computation of elliptic functions and integrals. This table illustrates the analogies between the elliptic function and elliptic integral cases, and the simplifications between the general (arbitrary parameters) and special (some parameters fixed at special values) cases.

		Elliptic functions	Elliptic integrals
	General case	Elliptic functions,	Incomplete elliptic integrals
		Jacobi theta functions
	Argument reduction Reduction to standard domain	Contraction of parameters
		(modular transformations,	(linear symmetric transformations)
		periodicity)
	Series expansions	Theta function q-series	Multivariate hypergeometric series
	Special case	Modular forms & functions,	Complete elliptic integrals,
		theta constants	arithmetic-geometric mean
	Argument reduction Reduction to standard domain	Contraction of parameters
		(modular transformations)	(quadratic transformations)
	Series expansions	Theta constant and eta function	Classical 2 F 1 hypergeometric series
		q-series

  5 128 t 4 -11 512 t 6 -469 32768 t 8 + . . . (1.24)

Table 1 .

 1

	Function	d = 10	d = 10 2	d = 10 3 d = 10 4 d = 10 5
	Elementary exp(x) functions log(x) Modular η(t) forms & j(t) functions (θ i (0,t)) 4 i=1 Elliptic (θ i (x,t)) 4 i=1 and theta ℘(x,t) functions (℘(x,t),℘ ′ (x,t)) 5.6 • 10 -5 7.7 • 10 -7 2.9 • 10 -6 0.00011 0.0062 8.1 • 10 -7 2.8 • 10 -6 0.00011 0.0077 6.2 • 10 -6 1.99 • 10 -5 0.00037 0.015 6.3 • 10 -6 2.3 • 10 -5 0.00046 0.022 7.6 • 10 -6 2.7 • 10 -5 0.00044 0.022 0.0016 0.089 2.8 • 10 -5 8.1 • 10 -5 0.00012 0.0021 0.11 3.9 • 10 -5 0.00017 0.0026 0.13 ζ (x,t) 0.00022 0.0028 0.14 7.5 • 10 -5 σ (x,t) 0.00022 0.0030 0.14 7.6 • 10 -5 Complete K(x) 5.4 • 10 -6 2.0 • 10 -5 0.00018 0.0068 elliptic E(y) 1.7 • 10 -5 6.1 • 10 -5 0.00072 0.025 integrals Π (x, y) 0.00046 0.014 3.6 7.0 • 10 -5 Incomplete -1 (x,t) 0.00014 0.0025 0.20 3.1 • 10 -5 elliptic F(x, y) 0.00011 0.0022 0.19 2.4 • 10 -5 integrals E(x, y) 0.00030 0.0070 0.76 5.6 • 10 -5 Π (x, y, z) 0.00017 0.00098 0.030 5.6	0.24 0.27 0.69 1.1 1.1 5.4 6.6 7.3 7.8 8.1 0.23 0.71 563 20 19 97 895
	R F (x, y, z) R G (x, y, z) R D (x, y, z) R J (x, y, z,t)	1.6 • 10 -5 9.5 • 10 -5 0.00027 4.7 • 10 -5 0.00016 2.1 • 10 -5 3.4 • 10 -5 0.00031	0.0020 0.0067 0.0046 0.012	0.18 0.75 0.57 2.6	18 95 78 428

2 Time in seconds to evaluate the function (or tuple of function values simultaneously) at d decimal digits of precision (p = ⌈d log 2 10⌉ bits) for d between 10 and 100 000. The arguments are set to generic complex numbers x

https://dlmf.nist.gov/

The number is somewhat smaller if the series is truncated optimally using a relative rather than an absolute tolerance.

We give the inverse form of the transformation.

For Π , Mathematica restricts this quasiperiodicity relation to hold only for -1 ≤ n ≤ 1.

This is an algebraic simplification, so we can take E 1 = 0 even if the input argument are represented by inexact balls.

At precision up to about 1000 digits, the elementary functions in Arb are significantly faster than the AGM due to using precomputed lookup tables and many low-level optimizations[START_REF] Johansson | Efficient implementation of elementary functions in the medium-precision range[END_REF].

Acknowledgment. The author thanks the organizers of the KMPB Conference on Elliptic Integrals, Elliptic Functions and Modular Forms in Quantum Field Theory for the invitation to present this work at DESY in October 2017 and for the opportunity to publish this extended review in the post-conference proceedings.