Hugo Labrande

Computing Jacobi's θ in quasi-linear time

Jacobi's θ function has numerous applications in mathematics and computer science; a naive algorithm allows the computation of θ(z, τ ), for z, τ verifying certain conditions, with precision P in O(M(P ) √ P ) bit operations, where M(P ) denotes the number of operations needed to multiply two complex P -bit numbers. We generalize an algorithm which computes specific values of the θ function (the theta-constants) in asymptotically faster time; this gives us an algorithm to compute θ(z, τ ) with precision P in O(M(P ) log P ) bit operations, for any τ ∈ F and z reduced using the quasi-periodicity of θ.

Introduction

Jacobi's θ function appears in a wide range of fields, such as non-linear differential equations (as a solution of the heat equation), the study of modular forms, and number theory, in which it is the main ingredient to convert between algebraic and analytic representations of elliptic curves. Namely, we have the embedding [13, I.4] C/(Z + τ Z) → P 3 (C) z → (θ 00 (2z, τ ), θ 01 (2z, τ ), θ 10 (2z, τ ), θ 11 (2z, τ ))

where the θ i are essentially the θ function with its z argument translated. We also have the equation ℘(z, τ ) = π 2 3 (θ 4 10 (0, τ ) -θ 4 01 (0, τ )) -π 2 θ 2 01 (0, τ )θ 2 10 (0, τ ) θ 2 00 (z, τ ) θ 2 11 (z, τ ) which allows to compute, for any point on the torus, its x-coordinate on the curve E(C) : y 2 = 4x 3 -g 2 x -g 3 .

Special values of the θ function have interesting additional properties: the theta-constants, the value of θ at points z = 0, 1 2 and τ 2 . As modular forms in τ , they are linked to other modular functions, such as the j-invariant or Dedekind's η function. Computing the value of the theta-constants allows one to compute the value of those modular functions; this has been used in record computations of class polynomials [START_REF] Enge | The complexity of class polynomial computation via floating point approximations[END_REF], which are interesting to generate safe cryptographic curves with the CM method.

The main problem we are dealing with here is to compute θ(z, τ ) with given absolute precision P , which allows us to compute the above embedding at any given precision. We will suppose throughout the article that (z, τ ) satisfy certain conditions; the general case can be deduced from this one using formulas we mention later. The θ function is defined by a rapidly convergent series; under the conditions specified on z, τ , it gives a naive algorithm that requires a running time of O(M(P ) √ P ) bit operations, where M(P ) is the number of operations needed to multiply two P -bit complex numbers. Although fast, this is a worse running time than other transcendental functions such as the exponential of a complex number, which can be computed in quasi-optimal time O(M(P ) log P ).

There is an algorithm to compute the theta-constants asymptotically faster than with the naive method, outlined in [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF]. This algorithm relies on connections between theta-constants and the arithmetico-geometric mean (AGM) of Gauss; the complex-valued AGM, when evaluated at the theta-constants, has interesting properties, and this is used along with Newton's method to, in a sense, invert the AGM and recover the values of the theta-constants. This algorithm allows computation of theta-constants for τ ∈ F with precision P in quasi-optimal time O(M(P ) log P ), independently of τ . It is faster than the naive method for precisions as low as a few thousand bits.

In this article, we provide a generalization of this algorithm which computes θ(z, τ ), for τ ∈ F and z such that Im(z) ≤ Im(τ )/2, with absolute precision P in O(M(P ) log P ) bit operations. We give two algorithms: the first one runs in quasi-optimal time in P , but its running time depends on z and τ ; we then use this algorithm as a subroutine to build a quasi-optimal algorithm with complexity independent of z and τ , provided τ ∈ F and Im(z) ≤ Im(τ )/2. An GNU MPC [START_REF] Enge | GNU MPC -A library for multiprecision complex arithmetic with exact rounding[END_REF] implementation of the algorithm was realized; it is faster than the naive method for values of P greater than a few hundred thousand digits.

Our algorithm provides the six values θ(z, τ ), θ(z + 1 2 , τ ), θ(z + τ 2 , τ ), θ(0, τ ), θ( 12 ) and θ( τ 2 ), which are sufficient to compute the projective embedding mentioned above. It can also be used to compute the Weierstrass ℘ function and its derivative in quasi-optimal time; hence, this paper provides a quasi-optimal time algorithm to compute the "Jacobi map" C/Λ → E(C) of an elliptic curve. We note that the "Abel map" E(C) → C/Λ can already be computed in quasi-optimal time using links to elliptic integrals and the Landen isogeny; see [START_REF] Bost | Moyenne arithmético-géométrique et périodes des courbes de genre 1 et 2[END_REF] and [START_REF] Cremona | The complex AGM, periods of elliptic curves over C and complex elliptic logarithms[END_REF].

We note that [START_REF] Luther | Reliable computation of elliptic functions[END_REF] gives an algorithm to compute θ with real arguments (i.e. θ(u, m) with 0 < m < 1 and 0 ≤ u ≤ K(m)), defined by its representation as an infinite product, with the same, quasi-optimal complexity. Their algorithm relies on the Landen transform for the θ function, and could perhaps be generalized to the complex setting. We had independently pursued this line of thought, but found that in the complex setting, the presence of trigonometric functions induced heavy precision losses for some inputs; however, there may be a workaround those issues, which would allow one to find an algorithm for the complex setting and with quasioptimal complexity relying on the Landen transform.

This article is organized as follows. We introduce the necessary mathematical background and the strategies we follow for argument reduction in Section 2, which justifies our choice to consider throughout the paper the case τ ∈ F and Im(z) ≤ Im(τ )/2; we then provide an analysis of the naive algorithm for θ(z, τ ) under those conditions. Section 3 introduces a sequence derived from relations between values of θ, and we prove quadratic convergence of a certain homogeneization of the sequence; this is what replaces the AGM in the general case. Section 4 gives a first algorithm for computing θ-functions, with a complexity depending on z and τ ; this is quasioptimal provided that z, τ belong to a compact set. We then get rid of the dependency in z, τ much in the same way as in the case of theta-constants, which gives a uniform algorithm with complexity O(M(P ) log P ). Section 5 shows timings for our GNU MPC implementation of this last algorithm and compare it to our implementation of the naive algorithm.

The function θ, and θ-constants

Definitions and argument reduction

We recall a few basic facts, following the presentation of [START_REF] Mumford | Tata lectures on theta[END_REF]. Proposition 2.2 (quasi-periodicity). We have θ(z+1, τ ) = θ(z, τ ) and θ(z+τ, τ ) = e -πiz-2πiτ θ(z, τ ); in fact for any integers a, b, θ(z + aτ + b, τ ) = e -iπa 2 τ -2iπaz θ(z, τ ).

(

) 1
We also define the following variants of the theta function (which are related to the definition of "theta functions with characteristics") [13, Section I.3]:

Definition 2.3. θ 00 (z, τ ) = θ (z, τ ) θ 10 (z, τ ) = exp (πiτ /4 + πiz) θ z + τ 2 , τ θ 01 (z, τ ) = θ z + 1 2 , τ θ 11 (z, τ ) = exp (πiτ /4 + πi(z + 1/2)) θ z + τ + 1 2 , τ
We define theta-constants as the values in 0 of those functions.

Those functions and their theta-constants are linked by a great number of formulas; we will give such formulas as we use them, and most of them can be found in [START_REF] Mumford | Tata lectures on theta[END_REF]Section I.5]. Note that we have:

Proposition 2.4. For any τ , the functions z → θ 00 (z, τ ), z → θ 01 (z, τ ), z → θ 10 (z, τ ) are even, while z → θ 11 (z, τ ) is odd.
The latter implies that θ 11 (0, τ ) = 0, so the only theta constants we are interested in are θ 00 (0, τ ), θ 01 (0, τ ), and θ 10 (0, τ ).

Those properties can be used for the purpose of argument reduction. For instance, we can use the parity of θ to suppose that Im(z) ≥ 0; if this is not the case, one can consider -z instead of z, which does not change the value of θ but ensures that Im(z) ≥ 0. Furthermore, Equation (1) can be used to recover θ(z, τ ) from the value of θ(z , τ ), where z is such that | Re(z )| ≤ 1 2 and | Im(z )| ≤ Im(τ ) 2 ; the added cost is the cost of computing an exponential factor. This exponential factor can become quite big; should one want to compute θ(z, τ ) with an error of at most 2 -P , they have to work with representations of at least P + C bits, with

C = log 2 (|θ(z , τ )|) + π log 2 (e)(a 2 Im(τ ) + 2a Im(z)) + 2
This is because the integral part of the result fits in C bits, while the fractional part should be coded on at least P bits to ensure a final error bounded by 2 -P . The complexity of running our algorithm and computing the exponential factor will be O(M(P + C) log(P + C)), and hence will depend on τ and z; this is inevitable. Hence, throughout the paper we suppose that z is reduced, in the sense that | Re(z)| ≤ 1 2 and 0 ≤ Im(z) ≤ Im(τ ) 2 , with the understanding that the step of argument reduction has a complexity depending on the original values of z and τ . However, as Section 2.2 shows, this hypothesis, combined with a hypothesis in τ , allows us to work with values of θ bounded by 4, which allows us to write an algorithm with complexity only depending on P for any z satisfying these conditions.

We can also reduce the second argument of θ. Define the action of SL 2 (Z) on the complex upper-half plane H by a b c d

• τ → aτ + b cτ + d .
Its fundamental domain is

F = {ω ∈ H | | Re(ω)| < 1/2, |ω| > 1} Computing τ ∈ F and M = a b c d ∈ SL 2 (Z)
such that τ = M τ can be done by finding the shortest vector in the lattice (1, τ ) using Gauss's algorithm [START_REF] Vallée | Gauss' algorithm revisited[END_REF], which (since the inertia is small) will be asymptotically negligible. The value of θ(z, τ ) can then be computed from θ(z , τ ) (for some value z ) using the following theorem:

Theorem 2.5 (extension of [START_REF] Mumford | Tata lectures on theta[END_REF]Theorem 7.1]). Let τ ∈ H and z ∈ C, and let

γ = a b c d ∈ SL 2 (Z)
. Suppose c > 0, or c = 0 and d > 0; if not, take -γ. Then we have:

θ i z cτ + d , aτ + b cτ + d = ζ i,γ,τ √ cτ + de iπcz 2 /(cτ +d) θ σ(i) (z, τ ) (2)
where the square root is taken with positive real part, ζ i,γ,τ is an eighth root of unity and σ is a permutation of the elements (00, 01, 10), defined by the following

e γ (τ ) = cτ + d, f γ (z, τ ) = e iπcz 2 /(cτ +d)
With a bit of care, one can prove

e γ1 (γ 2 τ )e γ2 (τ ) = e γ1γ2 (τ ), f γ1 z c 2 τ + d 2 , γ 2 τ f γ2 (z, τ ) = f γ1γ2 (z, τ )
The maps T : τ → τ + 1 and S : τ → -1 τ are generators of SL 2 (Z), and we have [13,

θ i (z, T τ ) = ζ i,T,τ e T (τ )f T (z, τ )θ σ T (i) (z, τ ), θ i z τ , Sτ = ζ i,S,τ e S (τ )f S (z, τ )θ σ S (i) (z, τ )
with square roots taken with real parts and for some ζ ∈ U 8 and σ S , σ T ∈ S 3 . Hence for all γ ∈ SL 2 (Z),

θ i z e γ (τ ) 2 , γτ = ζ i,γ,τ e γ (τ )f γ (z, τ )θ σγ (i) (z, τ )
for some root of unity and some permutation. The correspondance γ → σ γ can be determined from [13, p. 36], although one can simply notice it is independent of z and use the tables found by Gauss in the case of theta-constants [START_REF] Cox | The arithmetic-geometric mean of Gauss[END_REF]Eq. 2.15]. Finally, we could attempt to give a formula for ζ i,γ,τ , but it is more efficient to simply compute a very low precision approximation of θ(z, τ ) and compare it to the full-precision value to determine which eighth root is needed. We note that in general, because of the permutation σ γ , we need to have computed all three values θ 00,01,10 z cτ +d , aτ +b cτ +d in order to be able to use the formula to compute, say, θ 00 (z, τ ). We will occasionally talk about computing θ 11 , but this will not be the focus of the paper. Hence, the problem we consider in this paper is the following:

Compute

θ 00 (z, τ ), θ 00 (0, τ ), θ 01 (z, τ ), θ 01 (0, τ ), θ 10 (z, τ ), θ 10 (0, τ )

where |τ | > 1, | Re(τ )| ≤ 1 2 , Im(τ ) > 0, | Re(z)| ≤ 1 2 , 0 ≤ Im(z) ≤ Im(τ ) 2 (3)
in quasi-optimal time, i.e. O(M(P ) log P ).

Naive algorithm to compute θ

Partial summation of the series defining θ

We define the following partial summation for the series defining θ(z, τ ):

S B (z, τ ) = 1 + 0<n<B q n 2 (e 2iπnz + e -2iπnz )
where use the notation q = e iπτ . We have Proposition 2.6. Suppose that Im(τ ) ≥ 0.35 and 0 ≤ Im(z) ≤ Im(τ )/2; in particular, this is the case if the conditions (3) are satisfied. Then, for

B ≥ 1, |θ(z, τ ) -S B (z, τ )| ≤ 3|q| (B-1) 2 .
Proof. We look at the remainder of the series:

|θ(z, τ ) -S B (z, τ )| ≤ n≥B |q| n 2 (|e 2iπnz | + |e -2iπnz |) ≤ n≥B |q| n 2 (1 + |q| -n ) ≤ 2 n≥B |q| n 2 -n ≤ 2 n≥B |q| (n-1) 2 ≤ 2 n≥0 |q| (B-1+n) 2 ≤ 2|q| (B-1) 2 n≥0 |q| 2n(B-1)+n 2 ≤ 2 |q| (B-1) 2 1 -|q| 2B-1 (4)
A numerical calculation shows that for Im(τ ) ≥ 0.35, we have 2 1-|q| ≤ 3, which proves the proposition.

We can prove the same inequality for θ 01 , since the series that define it has the same terms, up to sign, as the series for θ. Note that, unlike the analysis of [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF] for naive theta-constant evaluation, we cannot get a bound for the relative precision: since θ( 1+τ 2 , τ ) = 0, there is no lower bound for |θ(z, τ )|. 1 If we set

B(P, τ ) = P + 2 π Im(τ ) log 2 (e) + 1.
we have 4|q| (B-1) 2 ≤ 2 -P , which means the approximation is accurate with absolute precision P . We just showed that:

Theorem 2.7. To compute θ(z, τ ) with absolute precision P bits, it is enough to sum over all k ∈ Z such that

|k| ≤ P + 2 π Im(τ ) log 2 (e) + 1
Note that this bound is larger than the one of [7, p. 5].

Naive algorithm

We then present a naive algorithm to compute not only the value of θ(z, τ ), but also the value of θ 01 (z, τ ), θ 00 (0, τ ), θ 01 (0, τ ) for only a marginal amount of extra computation; this is the algorithm we will use for comparison to the fast algorithm we propose in this article. The algorithm performs computations at a precision P, which we determine later so that the result is accurate to the desired precision P . Define the sequence (v n ) n∈N as v n = q n 2 (e 2iπnz + e -2iπnz ) so that θ(z, τ ) = 1 + n≥1 v n . This satisfies the following recurrence relation for n > 1:

v n+1 = q 2n v 1 v n -q 4n v n-1 .
We use this recursion formula to compute v n efficiently, which is similar to the trick used by [START_REF] Enge | Computing class polynomials for abelian surfaces[END_REF]Prop. 3]. This removes the need for divisions and the need to compute and store e -2iπnz , which can get quite big; indeed, computing it only to multiply it by the very small q n 2 is wasteful. The resulting algorithm is Algorithm 1.

Error analysis and complexity

We have the following theorem:

Theorem 2.8. For z, τ satisfying conditions (3), Algorithm 1 with P = P + log B + 7 computes θ 00 (z, τ ), θ 01 (z, τ ), θ 00 (0, τ ), θ 01 (0, τ ) with absolute precision P bits. This gives an algorithm which has bit complexity O M(P ) P Im(τ ) .

Algorithm 1 Compute θ 00,01 (z, τ ), θ 00,01 (0, τ ) for z, τ satisfying conditions (3).

1: prec ← P 2: B ← P +2 π Im(τ ) log 2 (e) + 1 3: θ 0,z ← 1, θ 1,z ← 1, θ 0,0 ← 1, θ 1,0 ← 1 4: q ← e iπτ , q 1 ← q, q 2 ← q 5: v 1 ← e 2iπ(z+τ /2) + e -2iπ(z-τ /2) , v ← v 1 , v ← 2 6: for n = 1..B do 7: /* q 1 = q n , q 2 = q n 2 , v = v n , v = v n-1 */ 8: θ 0,z ← θ 0,z + v, θ 1,z ← θ 1,z + (-1) n × v 9: θ 0,0 ← θ 0,0 + 2q 2 , θ 1,0 ← θ 1,0 + (-1) n × 2q 2 10:
q 2 ← q 2 × (q 1 ) 2 × q 11:

q 1 ← q 1 × q 12: temp ← v, v ← q 2 1 × v 1 × v -q 4 1 × v v ← temp 13: end for
Performing the analysis of this algorithm requires bounding the error that is incurred during the computation. We then compensate the number of inaccurate bits by increasing the precision. We use the following theorem: Theorem 2.9. For j = 1, 2, let z j = x j + iy j ∈ C and zj = xj + i ỹj its approximation. Suppose that |z j -zj | ≤ k j 2 -P and that k j ≤ 2 P/2 . Then

1. | Re(z 1 + z 2 ) -Re( z1 + z2 )| ≤ (k 1 + k 2 )2 -P 2. | Re(z 1 z 2 ) -Re( z1 z2 )| ≤ (2 + 2k 1 |z 2 | + 2k 2 |z 1 |)2 -P 3. | Re(z 2 1 ) -Re( z1 2 )| ≤ (2 + 4k 1 |z 1 |)2 -P
and the same bounds apply to imaginary parts as well; and

4. |e z1 -e z1 | ≤ |e z1 | 7k1+8.5 2 2 -P . Furthermore if |z j | ≥ 2k j 2 -P , 5. | Re z1 z2 -Re z1 z2 | ≤ 6(2+2k1|z2|+2k2|z1|) |z2| 2 + 2(4+8k2|z2|)(2|z1||z2|+1)+2 |z2| 4 2 -P
and the same bound applies to the imaginary part, and

6. | √ z 1 - √ z1 | ≤ k1 √ |z1| 2 -P .
This theorem is not very hard to prove; we refer to [START_REF] Labrande | Absolute error in complex fixed-point arithmetic[END_REF] for details.

Proof of Theorem 2.8. We first determine the size of the quantities we are manipulating; this is needed to evaluate the error incurred during the computation, as well as the number of bits needed to store fixed-point approximations of absolute precision P of the intermediate quantities.

Taking B = 1 in Proposition 2.6 gives |θ(z, τ ) -1| ≤ 3, so |θ(z, τ )| ≤ 4; actually, this also proves |S B (z, τ )| ≤ 4, which means that |θ 0,z |, |θ 1,z |, |θ 0,0 |, |θ
≥ 1 |v n | ≤ |q| n 2 +n + |q| n 2 -n ≤ (1 + |q| 2n )q n 2 -n ≤ 1.0049q n 2 -n ≤ 1.0049
Hence, storing all the complex numbers above, including our result, with absolute precision P only requires P + 2 bits, since their integral part is coded on only 2 bits. Note that, had we computed e -2iπnz before multiplying it by q n 2 , we would have needed O(Im(τ )) more bits, which worsens the asymptotic complexity.

Computing the absolute precision lost during this computation is done using Theorem 2.9. We start with the bounds |τ -τ | ≤ 1 2 2 -P and |z -z| ≤ 1 2 2 -P , coming from the hypothesis that the approximations of z and τ are correctly rounded with precision P. We then need to estimate k v1 and k q , which can be done using the formula giving the absolute error when computing an exponential from Theorem 2.9. Given that τ ∈ F, we have

|q -q| ≤ 0.07 7 × 1/2 + 8.5 2 2 -P ≤ 0.42 × 2 -P |v 1 -ṽ1 | ≤ 6(|e -π(Im(τ )+2 Im(z)) | + |e π(2 Im(z)-Im(τ )) |) × 2 -P ≤ 6(|q| + 1)2 -P ≤ 6.42 × 2 -P
which means that k q ≤ 0.42 and k v1 ≤ 6.42. We then need to evaluate the loss of precision for each variable and at each step of the algorithm, which gives recurrence relations with nonconstant coefficients. Solving those is rather tedious, and we use loose upper bounds to simplify the computation; we do not detail this proof in the present article. The results obtained by this method show that the error on the computation of the theta-constants is bounded by (0.3B + 105.958)2 -P , and the one on the computation of the theta function is smaller than (5.894B + 28.062)2 -P . This proves that the number of bits lost is bounded by log 2 B + c, where c is a constant smaller than 7; hence we set P = P + log B + 7.

Finally, evaluating π and exp(z) with precision P can be done in O(M(P) log P) [START_REF] Borwein | Pi and the AGM: a study in the analytic number theory and computational complexity[END_REF], but this is negligible asymptotically. In the end, computing an approximation up to 2 -P of θ(z, τ ) can be done in O M (P + log(P/ Im(τ )) + c)

P Im(τ ) = O M(P ) P Im(τ )
bit operations.

Computing θ 10

We mentioned in section 2.1 the need to compute θ 10 (z, τ ) and θ 10 (0, τ ) as well. One could think of recovering those values using Jacobi's quartic formula and the equation of the variety:

θ 00 (0, τ ) 4 = θ 01 (0, τ ) 4 + θ 10 (0, τ ) 4
(5)

θ 2 00 (z, τ )θ 2 00 (0, τ ) = θ 2 01 (z, τ )θ 2 01 (0, τ ) + θ 2 10 (z, τ )θ 2 10 (z, τ ) (6)
that is to say, compute

θ 10 (0, τ ) = θ 00 (0, τ ) 4 -θ 01 (0, τ ) 4 1/4 θ 10 (z, τ ) = θ 2 00 (z, τ )θ 2 00 (0, τ ) -θ 2 01 (z, τ )θ 2 01 (0, τ ) θ 10 (0, τ ) .
However, this approach induces an asymptotically large loss of absolute precision for both θ 10 (0, τ ) and θ 10 (z, τ ). According to Theorem 2.9, both square root extraction and inversion induce a loss of precision proportional to |z| -1 ; since θ 10 (0, τ ) ∼ 4q 1/2 , the number of bits lost by applying those formulas is O(Im(τ )). Note that those formulas would also induce a big loss in relative precision; since θ 00 (0, τ ) and θ 01 (0, τ ) are very close when Im(τ ) goes to infinity, the subtraction induces a relative precision loss of O(Im(τ )) bits (for more details, see [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF]Section 6.3]).

Either of those analyses show that, in order to compensate precision loss, the naive algorithm should actually be run with a precision of O(P + log B + Im(τ )), which gives a running time that worsens, insteads of getting better, when Im(τ ) gets big. We do not recommend this approach.

Instead, one should compute partial summations of the series defining θ 10 , much in the same way as we did for θ(z, τ ). We outline the analysis in this case, which is very similar to the one for θ: supposing n ≥ 2, we have n 2 -2n ≥ (n -2) 2 , which can be used to prove that |θ 10 (z, τ ) -S B | ≤ 3|q| (B-2) 2 , so that the bound on B is thus just one more than for θ; the recurrence relation is the same; q 2n |v 1 | is bounded by 2 instead of 1, which in the worst case means log B more guard bits are needed. In what follows, we will refer to this algorithm as "the naive algorithm to compute θ 10 (0, τ ), θ 10 (z, τ )"; its asymptotic complexity is, just like Algorithm 1, O M(P )

P Im(τ )
bit operations, which gets better as Im(τ ) increases. We note that similar considerations apply to the problem of computing θ 11 . One can compute θ 11 (z, τ ) using the formula [13, p.22]

θ 11 (z, τ ) 2 = θ 01 (z, τ ) 2 θ 10 (0, τ ) 2 -θ 10 (z, τ ) 2 θ 01 (0, τ ) 2 θ 00 (0, τ ) 2 . ( 7
)
Using this formula loses only a few bits of precision since θ 00 (0, τ ) is bounded; however, one then needs to compute a square root, which potentially loses O(Im(τ )) bits. Hence, a summation of the series, which directly gives θ 11 , is preferable.

Fast computation of theta-constants

Recall the definition of the arithmetico-geometric mean (AGM) for two positive real numbers a, b:

a 0 = a, b 0 = b a n+1 = a n + b n 2 b n+1 = a n b n (8)
The sequences (a n ) n∈N and (b n ) n∈N both converge to the same limit, called the arithmeticogeometric mean of a and b. Furthermore, (a n ) and (b n ) are quadratically convergent, in the sense of the following definition:

Definition 2.10. A sequence (a n ) is said to be quadratically convergent (to a limit ) if there is a C > 0 such that for n large enough:

|a n+1 -a n | ≤ C|a n -a n-1 | 2
The constant C in the case of the AGM can be taken as It is possible to generalize the AGM to complex numbers, but there are two possibilities for the choice of the square root at each step. We then call an AGM sequence for a and b any sequence (a n , b n ) n∈N such that

a 0 = a, b 0 = b, 2a n+1 = a n + b n , b 2 n+1 = a n b n
Note that there are uncountably many AGM sequences for a, b. We define unambiguously the AGM of two complex numbers following [START_REF] Cox | The arithmetic-geometric mean of Gauss[END_REF]: Proposition 2.11. Let a, b ∈ C and let (a n , b n ) n∈N be an AGM sequence for a and b. We say that the choice of signs is good at the rank n if

|a n -b n | < |a n + b n | or |a n -b n | = |a n + b n | and Im b n a n > 0
We call the AGM sequence for a and b in which all the choices of signs are good the optimal AGM sequence, and define AGM(a, b) as the limit of the optimal AGM sequence for a and b.

Finally we have the following proposition: Proposition 2.12 ([4, Proposition 2.1]). Let (a n , b n ) n∈N be an AGM sequence for a, b:

• If (a n , b n ) has infinitely many bad choices of sign, lim n→∞ a n = 0 and the convergence is at least linear;

• If (a n , b n ) has only finitely many bad choices of sign (for instance if it is optimal), (a n ) and (b n ) both converge quadratically to the same non-zero limit.

The link between the complex AGM and theta-constants is well-known:

Proposition 2.13. We have the following formulas linking theta-constants: 9)

θ 2 00 (0, 2τ ) = θ 2 00 (0, τ ) + θ 2 01 (0, τ ) 2 (
θ 2 01 (0, 2τ ) = θ 00 (0, τ )θ 01 (0, τ ) (10)
This shows that (θ 2 00 (0, 2 n τ ), θ 01 (0, 2 n τ )) n∈N is an AGM sequence for θ 2 00 (0, τ ) and θ 2 01 (0, τ ), and it converges quadratically to 1. Whether or not this sequence is the optimal AGM sequence is controlled by the following result: Proposition 2.14 ([7, Theorem 2] or [4, Lemma 2.9]). Define

F k = τ ∈ H such that | Re(τ )| < 1, τ + 3 4 ≥ 1 4 , τ + 1 4 > 1 4 , τ - 1 4 ≥ 1 4 , τ - 3 4 > 1 4 ⊂ H
Let τ ∈ F k , and let (a n , b n ) n∈N the optimal AGM sequence for θ 2 00 (0, τ ) and θ 2 01 (0, τ ). Then for all n we have Algorithm 4], an algorithm relying on the AGM, and with complexity O(M(P ) log P ), is given to compute the value of the theta-constants with precision P bits. The algorithm uses the fact that k (τ ) = θ 2 01 (0,τ ) θ 2 00 (0,τ ) is a solution to the following equation:

(a n , b n ) = (θ 2 00 (0, 2 n τ ), θ 2 01 (0, 2 n τ )), which implies AGM(θ 2 00 (0, τ ), θ 2 01 (0, τ )) = 1. In [7,
i AGM(1, z) -τ AGM 1, 1 -z 2 = 0
which is a consequence of the action of SL 2 (Z) on the theta-constants, as well as of Jacobi's quartic formula (Equation ( 5)). Newton's method, when given an approximation of k (τ ) with precision P/2 as well as the knowledge of τ with precision P , computes an approximation of k (τ ) with precision P -δ, where δ is a small constant. If one carries out Newton's method while doubling the working precision with each iteration, it is asymptotically only as costly as the last iteration; this means k (τ ) can be computed with precision P in quasi-optimal running time. One can then recover the individual values of the theta-constants using the equation

AGM 1, θ 01 (0, τ ) 2 θ 00 (0, τ ) 2 = 1 θ 00 (0, τ ) 2 . ( 11
)
However, the complexity of this algorithm is not uniform -that is to say, it reaches this complexity only for τ within a compact set. A variant of the algorithm is proposed in [7, Algorithm 5] which makes the complexity uniform: if P ≤ 2 log Im(τ ), use the naive algorithm (which gives the right complexity); if not, compute the value of the theta-constants at τ 2 n for some n ≤ log Im(τ ), and use the AGM to compute the theta-constants at τ . This gives an algorithm which complexity does not depend on τ .

A sequence related to θ-functions

Definition of the F sequence

We start with the following formula:

Proposition 3.1. θ 00 (z, 2τ ) 2 = θ 00 (z, τ )θ 00 (0, τ ) + θ 01 (z, τ )θ 01 (0, τ ) 2 ( 12
)
θ 01 (z, 2τ ) 2 = θ 00 (z, τ )θ 01 (0, τ ) + θ 01 (z, τ )θ 00 (0, τ ) 2 ( 13
)
This formula is called in [3, formula 3.13, p.39] the change of basis formula from the F 2 basis to the F (n,2) 2 basis. However a direct proof can be obtained with limited effort, using the series defining θ and some manipulations and term reorganization akin to n+m≡0 (mod 2)

q n 2 +m 2 = i,j∈Z q (i+j) 2 +(i-j) 2 .
We also note that one can similarly prove the following formula, which will be used in section 4.2:

θ 10 (z, 2τ ) 2 = θ 00 (z, τ )θ 00 (0, τ ) -θ 01 (z, τ )θ 01 (0, τ ) 2 (14)
We then define the following function:

F : C 4 → C 4 (x, y, z, t) → √ x √ z + √ y √ t 2 , √ x √ t + √ y √ z 2 , z + t 2 , √ z √ t
Hence, according to Proposition 2.13 and 3.1, for some appropriate choice of roots we have F θ 2 00 (z, τ ), θ 2 01 (z, τ ), θ 2 00 (0, τ ), θ 2 01 (0, τ ) = θ 2 00 (z, 2τ ), θ 2 01 (z, 2τ ), θ 2 00 (0, 2τ ), θ 2 01 (0, 2τ ) .

Remark. One can also write rewrite F using Karatsuba-like techniques

F (x, y, z, t) = ( √ x + √ y)( √ z + √ t) + ( √ x - √ y)( √ z - √ t) 4 , ( 15
) ( √ x + √ y)( √ z + √ t) -( √ x - √ y)( √ z - √ t) 4 , ( √ z + √ t) 2 + ( √ z - √ t) 2 4 , ( √ z + √ t) 2 -( √ z - √ t) 2 4
to speed up computations. Following section 2.3, we define a good choice for square roots at the rank n as the following conditions being satisfied:

• | √ z n - √ t n | < | √ z n + √ t n |, or | √ z n - √ t n | = | √ z n + √ t n | and Im √ tn √ zn > 0.
The last condition is equivalent to |z n -t n | ≤ |z n + t n |, which means that (z n , t n ) is an AGM sequence for (z 0 , t 0 ) in which all the choices of sign are good. Note that the condition |x -y| < |x + y| is equivalent to Re y x > 0. Again, similarly to the AGM, for any x, y, z, t ∈ C we define the optimal F sequence ((x n , y n , z n , t n )) n∈N as follows:

(x 0 , y 0 , z 0 , t 0 ) = (x, y, z, t)

(x n+1 , y n+1 , z n+1 , t n+1 ) = F (x n , y n , z n , t n )
where all the choices of sign for the square roots are good. The study of this sequence and its convergence is done in Section 3.4.

Link with θ-functions

More argument reduction

We go slightly further than the conditions (3) in order to justify the forthcoming results. We wish to further reduce z, as follows:

0 ≤ Im(z) ≤ Im(τ ) 4 , | Re(z)| ≤ 1 8 (16)
The first hypothesis allows us to avoid z = τ +1 2 , which is a zero of θ(z, τ ), and hence a pole of quotients of the form θi θ00 , which we consider in our algorithm much in the same way as [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF]. We prove Lemma 3. Proof. Write

Im z τ = 1 |τ | 2 | Im(z) Re(τ ) -Re(z) Im(τ )| ≤ Im(τ ) |τ | 2 1 4 1 2 + 1 8 = 1 4 Im -1 τ .
This will be used to apply Theorem 3.4 to z τ , -1 τ in Proposition 4.1. Note that those conditions are satisfied if one takes z = z 2 with z satisfying conditions (3). One can then compute θ 2 00,01 (z, τ ) from θ 2 00,01,10 (z/2, τ ) and θ 2 00,01,10 (0, τ ) using the following z-duplication formulas [13, p. 22]:

θ 00 (z, τ )θ 3 00 (0, τ ) = θ 4 01 (z/2, τ ) + θ 4 10 (z/2, τ ) (17) θ 01 (z, τ )θ 3 01 (0, τ ) = θ 4 00 (z/2, τ ) -θ 4 10 (z/2, τ ) θ 10 (z, τ )θ 3 10 (0, τ ) = θ 4 00 (z/2, τ ) -θ 4 01 (z/2, τ ).
This requires the knowledge of θ 10 (z, τ ) and the associated theta-constant; this could be computed using Jacobi's formula (Equation ( 5)) and the equation of the variety (Equation ( 6)), but we end up using a different trick in our final algorithm.

Good choices of sign and thetas

We now prove that, for the arguments we consider, the good choices of sign correspond exactly to values of θ:

Lemma 3.3.
For any τ such that Im(τ ) ≥ 0. 345 (in particular, for τ ∈ F) and z verifying the first condition in (16) we have

|θ 00 (z, τ ) -θ 01 (z, τ )| < |θ 00 (z, τ ) + θ 01 (z, τ )|,
which also proves that Re θ01(z,τ ) θ00(z,τ ) > 0, Re θ01(0,τ ) θ00(0,τ ) > 0.

Proof. Write:

|θ 00 (z, τ ) + θ 01 (z, τ ) -2| ≤ 2 n≥2,n even |q n 2 (w 2n + w -2n )| ≤ 2 n≥2,n even |q| n 2 (1 + |q| -n/2 ) ≤ 2 n≥1 |q| 4n 2 (1 + |q| -n ) ≤ 2|q| 3 + 2|q| 4 + 2q 16 1 -q 20 + 2q 14 1 -q 19 |θ 00 (z, τ ) -θ 01 (z, τ )| ≤ 2 n≥1,n odd |q| n 2 (1 + |q| -n/2 ) ≤ 2|q| 1/2 + 2|q| + 2q 9 1 -q 16 + 2q 7.5 1 -q 19
We have 2|q| 1/2 + 2|q| + 2q 9 1-q 16 + 2q 7.5

1-q 19 2 -(2|q| 3 + 2|q| 4 + 2q 16 1-q 20 + 2q 14 1-q 19 ) ≤ 1

for Im(τ ) > 0. 345, which proves the lemma.

We are now ready to prove: Theorem 3.4. Let (x n , y n , z n , t n ) be the optimal F sequence for θ 2 00 (z, τ ), θ 2 01 (z, τ ), θ 2 00 (0, τ ), θ 2 01 (0, τ ). For any τ such that Im(τ ) ≥ 0. 345 and z satisfying the first condition of (16) we have

(x n , y n , z n , t n ) = θ 2 00 (z, 2 n τ ), θ 2 01 (z, 2 n τ ), θ 2 00 (0, 2 n τ ), θ 2 01 (0, 2 n τ )
Proof. This is true for n = 0; we prove the statement inductively. Suppose it is true for n = k.

We have for any τ :

θ 00 (0, τ ) = 1 + 2q + c, |c| ≤ 2|q| 4 1 -|q| 5
For any τ such that Im(τ ) ≥ 0. 345, 2|q| ≤ 0.676 and |c| ≤ 0.027; hence Re(θ 00 (0, 2 k τ )) > 0 for any k, which proves that √ z k = θ 00 (0, 2 k τ ). Lemma 3.3 shows that Re θ01(0,τ ) θ00(0,τ ) > 0, and we also have Re

√ t k √ z k ≥ 0 since the choice of roots is good, hence √ t k = θ 01 (0, 2 k τ ). Proposition 2.13
then proves that t k+1 = θ 2 01 (0, 2 k+1 τ ) and z k+1 = θ 2 00 (0, 2 k+1 τ ).

Similarly, given that z satisfies the first condition of (16):

|θ 00 (z, τ ) -1| ≤ |q| 1/2 + |q| + |q| 3 + |q| 4 + |q| 7/2 + |q| 9 + 2|q| 14 1 -|q| 2 (18)
For Im(τ ) ≥ 0. 345, this is strictly smaller than 1; hence Re(θ 00 (z, τ )) > 0, which proves that √ x k = θ 00 (z, 2 k τ ). Again, Lemma 3.3 proves that Re θ01(z,2 k τ ) θ00(z,2 k τ ) > 0, and since the choice of signs is good, Re

√ y k √ x k ≥ 0, necessarily √ y k = θ 01 (z, 2 k τ )
. This along with Proposition 3.1 finishes the induction.

Note that a consequence of this proposition is the following fact:

Proposition 3.5. The optimal F sequence for θ 2 00 (z, τ ), θ 2 01 (z, τ ), θ 2 00 (0, τ ), θ 2 01 (0, τ ) converges quadratically to (1, 1, 1, 1).

A function with quasi-optimal time evaluation

The strategy of [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF] is to use an homogenization of the AGM to get a function f τ : C → C, on which Newton's method can be applied. To generalize this, we homogenize the function which maps to (x, y, z, t) the limit of the optimal F sequence associated to them; it becomes a function from C 2 to C 2 . We call this function F ∞ ; this function is a major building block for the function we use to compute our two parameters z, τ using Newton's method. Proposition 3.6. Let λ, µ ∈ C. Let ((x n , y n , z n , t n )) n∈N be the optimal F sequence for (x, y, z, t), and ((x n , y n , z n , t n )) n∈N the optimal F sequence for (λx, λy, µz, µt). Put lim n→∞ z n = z ∞ and

lim n→∞ z n = z ∞ . Then we have µ = z ∞ z ∞ , λ = lim n→∞ x n z ∞ 2 n × z ∞ lim n→∞ xn z∞ 2 n × z ∞
Proof. We prove by induction that

x n = n λ 1/2 n µ 1-1/2 n x n , y n = n λ 1/2 n µ 1-1/2 n y n , z n = µz n , t n = µt n ,
where Re(λ 1/2 n ) ≥ 0, Re(µ 1-1/2 n ) ≥ 0, and n is a 2 n -th root of unity. This is enough to prove the proposition above, since then

lim n→∞ x n z ∞ 2 n = lim n→∞ λµ 2 n -1 x n z ∞ 2 n = λ µ lim n→∞ x n z ∞ 2 n .
Since this is true for n = 0, suppose this is true for n = k. We have

z k+1 = z k + t k 2 = µz k+1 .
As for t k+1 , we can write

z k = z √ µ √ z k , t k = t √ µ √ t k
where z = ±1 and t = ±1, and the square roots are taken with positive real part. But since

Re √ t k √ z k ≥ 0 and Re √ t k √ z k ≥ 0, we have z = t . Hence t k+1 = ( z √ µ √ z k ) z √ µ √ t k = µt k+1 .
As for the other coordinates, we have

x k = x √ k λ 1/2 k+1 µ 1/2-1/2 k+1 √ x k , y k = y √ k λ 1/2 k+1 µ 1/2-1/2 k+1 √ y k
where the roots are taken with positive real part, and x , y ∈ {-1, 1}. Since Re

√ y k √ x k
≥ 0 and we require Re

√ y k √ x k
≥ 0, necessarily x = y ; hence

x k+1 = x k z k + y k t k 2 = k+1 λ 1/2 k+1 µ 1-1/2 k+1 √ x k √ z k + √ y k √ t k 2 = k+1 λ 1/2 k+1 µ 1-1/2 k+1 x k+1
where k+1 = x √ z is indeed such that 2 k+1 k+1 = 1. This proves the proposition. In the case of theta-functions, however, we have:

Proposition 3.7. lim n→∞ θ(z, 2 n τ ) 2 n θ(0, 2 n τ ) 2 n = 1
Proof. It is enough to prove lim n→∞ θ(z, 2 n τ ) 2 n = 1, since it also covers the case z = 0. Write, like Equation (18):

θ(z, 2 n τ ) = 1 + q 2 n (w 2 + w -2 ) + c, |c| ≤ 2|q| 2 n +3 1 -|q|
We can write θ(z,

2 n τ ) = 1 + d n with |d n | ≤ 2|q| 2 n |w 2 + w -2 |.
We then have classically

|θ(z, 2 n τ ) 2 n -1| ∼ |2 n d n |
and since lim n→∞ 2 n |q| 2 n = 0, this proves the proposition.

Combining this with Theorems 3.6 and 3.4 proves that, for (x n , y n , z n , t n ) as in Theorem 3.4, we have

λ = lim n→∞ x n µ 2 n × µ.
In particular, if we define the following function:

F ∞ : C 4 → C 2 (x, y, z, t) → lim n→∞ x n z ∞ 2 n × z ∞ , z ∞
where z ∞ = lim n→∞ z n , then we have that, for any z, τ satisfying the hypotheses of Theorem 3.4,

F ∞ λθ 2 00 (z, τ ), λθ 2 01 (z, τ ), µθ 2 00 (0, τ ), µθ 2 01 (0, τ ) = (λ, µ). For instance, F ∞ 1, θ 2 01 (z, τ ) θ 2 00 (z, τ ) , 1, θ 2 01 (0, τ ) θ 2 00 (0, τ ) = 1 θ 2 00 (z, τ ) , 1 θ 2 00 (0, τ ) . ( 19
)
This is similar to Equation ( 11), and will play a similar role in the computation of θ(z, τ ).

Convergence

Let us start by showing that, contrary to the AGM and despite Proposition 3.5, an optimal F sequence does not always converge quadratically; for instance, the optimal F sequence for (2, 2, 1, 1) is ((2 1/2 n , 2 1/2 n , 1, 1)) n∈N , which does not converge quadratically. This is a big difference from the AGM, and this is why we are reluctant to call optimal F sequences a "generalization of the AGM". However, we now show that the sequence (λ n ) = Proof. The upper bound result follows from a trivial induction using the equations defining F . We prove in Section 4.3.5 that there is a suitable C for any z, τ we consider in our final algorithm. We now prove the existence of c. Recall that the choice of signs for the square roots are good at all steps, since we assume (x n , y n , z n , t n ) is an optimal F sequence. Thus there exists

α, β ∈ C * such that Re(x 1 /α) > 0, Re(y 1 /α) > 0, Re(z 1 /β) > 0, Re(t 1 /β) > 0.
For instance, in most cases one can take α = x 1 and β = z

| √ x n + √ y n | ≥ √ 2c | √ z n + √ t n | ≥ √ 2c
and hence

| √ x n - √ y n | ≤ |x n -y n | √ 2c , | √ z n - √ t n | ≤ |z n -t n | √ 2c
Proof. The parallelogram identity gives

| √ x n + √ y n | 2 = 2| √ x n | 2 + 2| √ y n | 2 -| √ x n - √ y n | 2 ≥ 2| √ x n | 2 + 2| √ y n | 2 -| √ x n + √ y n | 2 since
-t k | ≤ 2 -P -k-c1
, then λ k+1 is an approximation of λ with absolute precision P bits.

Proof. The point here is that once z n and t n are close enough, x n+1 and y n+1 are also close and the value of λ n does not change much after that. Let c 1 ≥ 0, and take n the first integer for which |z n -t n | ≤ η with η = 2 -P -c1-n . We then have for all k ≥ 0 [7, Theorem 1]:

|z n+k -t n+k | ≤ A 2 k -1 η 2 k with A = π 8 min(|z0|,|t0|) . Furthermore, |z n+1 -z n | = 1 2 |z n -t n |, so that |z ∞ -z n+k | ≤ 1 2 ∞ i=k A 2 i -1 η 2 i and we have |z ∞ -z n+k | ≤ 1 A (Aη) 2 k
. Finally, using Equation ( 15), one can write

|x n+k+1 -y n+k+1 | ≤ | √ x n - √ y n || √ z n - √ t n | 2 ≤ √ C √ 2 |z n+k+1 -t n+1 | since z n+k+1 -t n+k+1 = ( √ z n+k - √ t n+k ) 2 2 ≤ √ 2AC|z n+k -t n+k |.
Now, define q n = (xn/z∞) 2 xn-1/z∞ , so that λn+1 λn = q 2 n n . Note that if one makes the approximation

x n+k+1 = y n+k+1 and z n+k+1 = t n+k+1 = z ∞ , we have x n+k+2 = √ x n+k+1 z ∞ which gives q n+k+2 = 1.
We take a closer look at those approximations:

|x n+k+2 - √ x n+k+1 √ z n+k+1 | ≤ | √ y n+k+1 - √ x n+k+1 || √ z n+k+1 + √ t n+k+1 | 4 + | √ y n+k+1 + √ x n+k+1 || √ z n+k+1 - √ t n+k+1 | 4 ≤ √ C 2 (| √ y n+k+1 - √ x n+k+1 | + | √ z n+k+1 -t n+k+1 |) ≤ √ C 2 √ 2AC | √ x n+k+1 + √ y n+k+1 | |z n+k -t n+k | + √ 2A|z n+k+1 -t n+k+1 | ≤ B(Aη) 2 k so q n+k+2 -1 = (x n+k+2 /z ∞ ) 2 -x n+k+1 /z ∞ x n+k+1 /z ∞ ≤ 1 z 2 ∞ ( √ x n+k+1 √ z n+k+1 + B(Aη) 2 k ) 2 -x n+k+1 /z ∞ x n+k+1 /z ∞ ≤ x n+k+1 z n+k+1 /z 2 ∞ -x n+k+1 /z ∞ x n+k+1 /z ∞ + 2C/z 2 ∞ B(Aη) 2 k + B 2 (Aη) 2 k+1 x n+k+1 /z ∞ ≤ 1 2z ∞ ∞ i=k+1 (Aη) 2 i + CB(Aη) 2 k x n+k+1 z ∞ + B 2 (Aη) 2 k+1 x n+k+1 /z ∞ ≤ B × (Aη) 2 k
This proves that (q n ) converges quadratically to 1; using the equivalent q 2 n n -1 ∼ 2 n q n , we have that (q 2 n n ) also converges quadratically to 1, which proves the convergence of the sequence (λ n ). Finally we have

q 2 n+2 n+2 ...q 2 n+k n+k -1 ≤ k-2 i=0 1 + B (Aη) 2 i 2 n+i+2 -1 < exp 4B k-2 i=0 2 n+i (Aη) 2 i -1 ≤ 4B k-2 i=0 2 n+i (Aη) 2 i 1 -(4B k-2 i=0 2 n+i (Aη) 2 i )/2
which proves, for P large enough,

|λ -λ n+1 | ≤ 8B |λ n+1 | ∞ i=0 2 n+i (Aη) 2 i ≤ 16B |λ n+1 |Aη2 n ≤ 16AcB |λ n+1 |2 -c1 × 2 -P
This inequality proves that, at least for P large enough, |λ n+1 | ≤ 2|λ|. Hence if we suppose log 2 (32Ac|λ|) ≤ c 1 , we have that λ n+1 is an approximation of λ with P bits of absolute precision. (x, y, z, t) ← F (x, y, z, t) 6: end while 7: (x, y, z, t) ← F (x, y, z, t)

8: Return x z 2 n+1 × z, z
This gives an algorithm, Algorithm 2, to compute F ∞ (x, y, z, t). According to [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF]Theorem 12], if n = max(log | log |z 0 /t 0 ||, 1) + log(P + c 1 ), a n is an approximation of AGM 1, | z0 t0 | with relative precision P bits. This proves that at the end of the algorithm, n = O(log P ); in fact, we have more precisely n ≤ log 2 P + C with C a constant independent of P . Finally in the next subsection proves that one can take P = P + O(log P ), which means that this algorithm computes F ∞ in O(M(P ) log P ) bit operations.

Loss of precision

We use Theorem 2.9 in order to evaluate the precision lost when computing F ∞ (x, y, z, t). First note that the upper and lower bounds on the terms of the sequence allow us to write 1

|z n | + 1 |t n | |z n | + |t n | ≤ b/2 1 |z n | + 1 |t n | |x n | + |y n | ≤ b 1 |x n | + 1 |y n | |z n | + |t n | ≤ b
for some b > 1; for instance, one can take b = max 1, 4 C c . We prove in Section 4.3.5 the existence of c and C, and hence of b, for any values of theta we consider as arguments.

We first evaluate a bound on the error incurred when computing F using Equation ( 15). Using those formulas allows us to get error bounds that are identical for F x and F y , and F z and F t . For simplicity, we assume that the error on z and t is the same, as well as the error on x and y. This gives:

| Re(F x ) -Re( Fx )| ≤ 1 + k z 1 |z| + 1 |t| |x| + |y| + k x 1 |x| + 1 |y| |z| + |t| 2 -P | Re(F z ) -Re( Fz )| ≤ 1 + 2k z 1 |z| + 1 |t| |z| + |t| 2 -P
We thus get the following induction relations when looking at what happens when applying F n times in a row:

k (n) x ≤ 1 + bk (n-1) z + bk (n-1) x , k (n) z ≤ 1 + bk (n-1)
z

The last equation rewrites as k

(n) z + 1 b-1 ≤ b k (n-1) z + 1 b-1 , which gives k (n) z ≤ b n k z + 1 b-1 . The induction for x becomes k (n) x ≤ 1 + (k z + 1 b-1 )b n + bk (n-1) x
, which we solve:

k (n) x ≤ (1 + b + ... + b n )k x + nb n k z + 1 b -1 ≤ b n nk z + b + n b -1
For b > 1, we have for n large enough that k (n) ≤ 2b 2n , which ultimately means the number of bits lost when applying F n times in a row is bounded by 2n log b + 1.

Finally we need to find the number of bits lost in the computation of xn z∞

2 n
. Call E k the error made after computing k squarings in a row; we have the following recurrence relation:

E k+1 ≤ 2 + 4E k |x n /z ∞ | 2 k However, since (λ n ) converges, |λ n | ≤ ρ for some constant ρ; furthermore, for any k ≤ n, one has |x n /z ∞ | 2 k ≤ 1 + ρ z∞ . Hence the recurrence becomes E k+1 ≤ 2 + 4 1 + ρ z∞ E k , which we solve to get E n ≤ 2 C n+1 -1 C -1 ≤ 2 C -1 C n+1
with C = 4 1 + ρ z∞ . This means the number of bits lost after n successive squarings is at the most (n + 1) log C + 1 -log(C -1).

Overall, if we write that the final value of n in Algorithm 2 is bounded by log 2 P + C , we have that the number of bits lost is bounded by

(2 log 2 b + log C )(log 2 P + C ) + log C + 2 -log C -1 which is O(log P ).

Fast computation of θ

We use a similar method as [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF], that is to say finding a function F such that

F θ 2 01 (z, τ ) θ 2 00 (z, τ ) , θ 2 01 (0, τ ) θ 2 00 (0, τ ) = (z, τ ),
which can then be inverted using Newton's method. One can then compute θ(z, τ ) by, for instance, using Equation ( 19) and extracting a square root, determining the correct choice of sign by computing a low-precision (say, 10 bits) approximation of the value using the naive method; we use a different trick in our final algorithm (Algorithm 5). We build this function F using F ∞ as a building block.

Building F

Just as with the algorithm for theta-constants, we use formulas derived from the action of SL 2 (Z) on the values of θ in order to get multiplicative factors depending on our parameters; this will allow us to build a function which computes z, τ from the values θ i (z, τ ). We define the function F as the result of Algorithm 3.

Algorithm 3 Compute F (s, t)

1: b ← √ 1 -t 2
Choose the root with positive real part [4, Prop 2.9]

2: a ← 1-st b 3: (x, y) ← F ∞ (1, a, 1, b) 4: (q 1 , q 2 ) ← F ∞ (1, s, 1, t) 5:
Return log q2x q1y × q2/y -2π , i q2 y , choosing the sign of the square root so that it has positive imaginary part.

F θ 2 01 (z, τ ) θ 2 00 (z, τ ) , θ 2 01 (0, τ ) θ 2 00 (0, τ ) = (z, τ ) Proof. Equation (19) proves that (q 1 , q 2 ) = 1 θ00(z,τ ) 2 ,
1 θ00(0,τ ) 2 . Furthermore, using Jacobi's formula [START_REF] Cremona | The complex AGM, periods of elliptic curves over C and complex elliptic logarithms[END_REF] and the equation defining the variety [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF], it is easy to see that b = θ10(0,τ ) 2 θ00(0,τ ) 2 and a = θ10(z,τ ) 2 θ00(z,τ ) 2 .

20

The formulas in [START_REF] Mumford | Tata lectures on theta[END_REF]Table V,p.36] give

θ 2 00 (z, τ ), θ 2 10 (z, τ ), θ 2 00 (0, τ ), θ 2 10 (0, τ ) = λθ 2 00 z τ , -1 τ , λθ 2 01 z τ , -1 τ , µθ 2 00 0, -1 τ , µθ 2 01 0, -1 τ with λ = e -2iπz 2 /τ -iτ , µ = 1 -iτ .
From the discussion in Section 3.2.1, the conditions on z, τ allow us to apply Theorem 3.4 to z τ , -1 τ . This proves that

F ∞ θ 2 00 (z, τ ), θ 2 10 (z, τ ), θ 2 00 (0, τ ), θ 2 10 (0, τ ) = e -2iπz 2 /τ -iτ , 1 -iτ ,
and by homogeneity, (x, y)

= e -2iπz 2 /τ -iτ θ00(z,τ ) 2 , 1 -iτ θ00(0,τ ) 2 .
This means that, starting from the knowledge of z and τ with precision P and a low-precision approximation of the quotients θ01(z,τ ) θ00(z,τ ) and θ01(0,τ ) θ00(0,τ ) , one can compute those quotients with precision P using Newton's method. This is Algorithm 4.

Algorithm 4

Compute θ 2 00 (z, τ ), θ 2 01 (z, τ ), θ 2 00 (0, τ ), θ 2 01 (0, τ ) with precision P . Input: (z, τ ) with absolute precision P .

1: Compute θ 2 00,01 (z, τ ), θ 2 00,01 (0, τ ) with absolute precision P 0 using Algorithm 1.

2: s ← θ01(z,τ ) 2 θ00(z,τ ) 2 , t ← θ01(0,τ ) 2 θ00(0,τ ) 2
3: p ← P 0 4: while p ≤ P do 5:

p ← 2p 6:
Compute a 11 = ∂Fx ∂x (s, t), a 22 = ∂Fy ∂y (s, t), a 12 = ∂Fx ∂y (s, t) with precision p. We make a few remarks:

• The proof of the invertibility of the Jacobian of F is delayed until Section 4.3.2.

• Much in the same way as [START_REF] Enge | Computing class polynomials for abelian surfaces[END_REF], we find it preferable to use finite differences to compute the coefficients a 11 , a 21 , a 22 of the Jacobian, as it does not require the computation of the derivative of F, which could be tedious.

• The value of P 0 has to be large enough that Newton's method converges. We note that, in general, a lower bound on P 0 may depend on the arguments; for instance, [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF] experimentally finds 4.53 Im(τ ) to be a suitable lower bound for P 0 when computing theta-constants. However, we outline in the next section a better algorithm which only uses the present algorithm for z, τ within a compact set; hence, P 0 can be chosen to be a constant, and we use in practice P 0 = 30000.

We do not provide a full analysis for this algorithm: we outline in the next section a better algorithm, which uses this algorithm as a subroutine, and we will provide a full analysis at that time. It is enough to say that the computation of F ∞ at precision p is done in time O(M(p) log p) using Algorithm 2; however, this running time depends on z, τ , since it depends on the bounds C, c that one can write for |x n |, |y n |, |z n |, |t n |. Hence, the cost of evaluating F at precision p is O(M(p) log p) bit operations, and the fact that we double the working precision at every step means that the algorithm is as costly as the last iteration. Furthermore, one should choose P so that the final result is accurate with absolute precision P . This means compensating the loss of absolute precision incurred during the computation of F; in general, this only depends on Im(τ ) and linearly in log p. Furthermore, we have the following proposition: Proposition 4.2. Let x and y be approximations of a and b with absolute precision N , and x , y the result of Step 7 in Algorithm 4 when using finite differences to approximate the Jacobian matrix. Then x , y are approximations of a, b with precision 2P 0 -δ.

Its proof can be adapted from the proof of [START_REF] Enge | Computing class polynomials for abelian surfaces[END_REF]Theorem 12]. In practice, we found δ = 4 to be enough. Determining the number of bits lost at each step can be done in the same way as [10, p. 19]: if s (n-1) and s (n-2) agree to k bits, and s (n) and s (n-1) agree to k bits, the number of bits lost can be computed as 2k -k . In the end, working at precision P = P + c log P + d, with c, d independent of P but functions of z, τ , is enough to compensate all precision losses; this proves that the running time of this algorithm is asymptotically O(M(P ) log P ).

Computing θ(z, τ ) in uniform quasi-optimal time

We now show an algorithm with uniform (i.e. independent in z and τ ) quasi-optimal complexity that computes θ(z, τ ) for any (z, τ ) satisfying conditions [START_REF] Cosset | Applications des fonctions thêta à la cryptographie sur courbes hyperelliptiques[END_REF]. We use the same strategy as [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF]; namely, we use the naive algorithm when Im(τ ) is large; and for smaller values of Im(τ ), we put τ = τ 2 s so that τ is within a compact set, then use Algorithm 4, which complexity will be uniform since its arguments belong to a compact set. However we also need to divide z by a power of 2 so that it also belongs to a compact set, and so that (z , τ ) satisfies conditions (3) and (16). Once θ z 2 t , τ 2 s has been computed by the previous algorithm, we alternate between using Equation [START_REF] Luther | Reliable computation of elliptic functions[END_REF] to double the second argument and Equation (17) to double the first argument, until finally recovering θ(z, τ ). This is Algorithm 5.

A few notes on this algorithm:

• We note that, at several steps of the algorithm (e.g. Steps 9, 14, 16) we need to compute theta-constants from their square. The correct choice of signs is given by the proof of Theorem 3.4, which shows that Re(θ 00 (0, τ )) ≥ 0 and Re(θ 01 (0, τ )) ≥ 0; and furthermore, since Re(q 1/4 ) ≥ |q| 1/4 cos(π/8), we also have Re(θ 10 (0, τ )) ≥ 0.

• Taking τ 2 = τ 1 /2 allows us to use Equation ( 14) in step 8 instead of Equation ( 5) and Equation ( 6), which is more efficient and loses fewer bits.

• The knowledge of θ 2 10 (2 i-2 z 1 , 2 i τ ) is enough for the z-duplication formulas of step 17, and it can be computed directly from θ 00 and θ 01 using Equation ( 14).

• Computing θ 11 (z, τ ) is also possible; one should use a partial summation if P ≤ 25 Im(τ ).

In the other case, since all the z-duplication formulas for θ 11 (z, τ ) involve a division by θ 10 (0, τ ) [13, p.22], it is just as efficient to simply use Equation [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF] after Step 20, then extract the square root. The square root extraction loses O(Im(τ )) = O(P ) bits, and this also gives a quasi-optimal algorithm.

Proving the correctness of the algorithm

This section is devoted to proving the following theorem:

Theorem 4.3. For any τ, z satisfying conditions (3), Algorithm 5 with P = 2P computes θ 00,01,10 (z, τ ), θ 00,01,10 (0, τ ) with absolute precision P in O(M(P ) log P ) bit operations.

As we discussed in Section 2, this also gives an algorithm that computes θ(z, τ ) for any (z, τ ) ∈ C × H; one simply needs to reduce τ in τ ∈ F, then reduce z in z , and deduce θ(z, ) from θ(z , τ ) using Equations ( 1) and ( 2). This causes a loss of absolute precision which depends on z and τ , and this algorithm is no longer uniform.

We need to perform an analysis of the number of bits lost by the algorithm; once again, we use Theorem 2.9. For each step, we proceed as follows: assuming the error on all the quantities is bounded by k, determine a factor x such that the error on the quantities we get after the computation is bounded by xk, then declare the number of bits lost in this step to be log x; this gives a very loose upper bound, but simplifies the process.

Finally, we also need to prove that the hypotheses made in Sections 3.4 and 3.5 are verified in Step 7 of the algorithm. This is necessary to prove that the sequence (λ n ) we consider is quadratically convergent, and that the number of bits lost is only O(log P ). We prove this in Section 4.3.5, which then completes the proof that the running time is indeed uniform and quasi-optimal.

Naive algorithm

As we showed in Theorem 2.8, the number of bits lost when using the naive algorithm is log B +7, although this constant could be made even smaller when taking into account that P ≤ 25 Im(τ ). Furthermore, P Im(τ ) ≤ 25, which means the running time of this step is asymptotically dominated by the cost of the computation of π, q and w with precision P = P + log B + 7, which takes O(M(P ) log P ) bit operations.

1 = ∂(F • u) 1 ∂z (z, τ ) = a × ∂u 1 ∂z (z, τ ) + ∂F ∂u 2 θ 2 1 (z, τ ) θ 2 0 (z, τ ) , θ 2 1 (0, τ ) θ 2 0 (0, τ ) × ∂u 2 ∂z (z, τ )
This already proves that a = 0, but we can actually give a more explicit form:

a = ∂u 1 ∂z (z, τ ) -1 = 2 θ 1 (z, τ ) θ 0 (z, τ ) θ 1 (z, τ ) θ 0 (z, τ ) = 2 θ 1 (z, τ ) θ 3 0 (z, τ ) (θ 1 (z, τ )θ 0 (z, τ ) -θ 0 (z, τ )θ 1 (z, τ )) = 2 θ 1 (z, τ ) θ 3 0 (z, τ ) πθ 2 2 (0, τ )θ 2 (z, τ )θ 3 (z, τ ) = 0
the last equality deriving from Formula 10 in [15, Section 23]. This finishes the proof.

Square root extraction

Steps 9, 16 and 14 require extracting square roots, which multiply the error by

Duplication formulas and finishing the proof of correctness

The algorithm uses both τ -duplication formulas and z-duplication formulas, and we need to analyse how many bits are lost for each application of those formulas. The τ -duplication formulas are nothing more than applying F to θ 2 00,01 (z, τ ) and θ 2 00,01 (0, τ ). However, the analysis here is simpler than in section 3.5, because we do not need to compute the square roots of θ 00,01 (z, τ ), since they are directly given by step 17. Hence we just need to account for the error of the additions, subtractions and multiplications in Equation ( 15); since all the quantities are bounded, this means each step loses a constant number g of bits (our analysis shows that g ≤ 10.48). In the end, the τ -duplication formulas account for the loss of g × s ≤ g log P bits of precision.

As for the z-duplication formulas, we need to perform several analyses. Looking at Equation (17), one needs to evaluate the fourth power of theta functions, then add them; then evaluate the third power of theta constants, then perform a division. Computing the error using the formulas from 2.9 is rather straightforward when one has bounds on those quantities, which are given by the following theorem: We also have: which means losing at most 16 more bits of precision.

• Suppose that 0 ≤ Im(z) ≤ Im(τ ) 8 ,
Step 19 causes the loss of a greater number of bits. We use Equation ( 6) instead of the third z-duplication formula, because dividing by θ 10 (0, τ ) 2 loses less bits than dividing by θ 10 (0, τ ) 3 , and we only need the knowledge of θ 2 10 (2 s-1 z, 2 s τ ) for the next step anyway. This amounts to computing: θ 2 10 (z, τ ) = θ 2 00 (z, τ )θ 2 00 (0, τ ) -θ 2 01 (z, τ )θ 2 01 (0, τ ) θ 2 10 (0, τ ) Computing the numerator multiplies the error by a factor at most 60, and the norm of this numerator is bounded by 4.557; we then get from Theorem 2.9 that the error is bounded by m |θ10(0,τ )| 8 ∼ m|q| -2 , with m ≤ 1600. In the end, we lose at most 2π log 2 e Im(τ ) + 11 bits.

Finally, we also lose a great number of bits during the last application of the z-duplication formulas in step 20, since the formula for θ 10 (z, τ ) requires dividing by θ 10 (0, τ ) 3 . The error is thus multiplied by |q| -3 up to a constant factor; this means a loss of 3π log 2 e Im(τ ) bits, plus a constant.

In the end, we see that the number of lost bits is bounded by (2π + π/8 + 3π) Im(τ ) log 2 e + c log P + d; given that P ≥ 25 Im(τ ), and that 5.125π log 2 e ≤ 23.3, the number of bits lost is thus less than P . This means that P = 1.93P + c log P + d ≤ 2P is enough to get a result which is accurate to absolute precision P ; this also means that we indeed never have an error k bigger than 2 (2P )/2 , which is necessary to apply Theorem 2.9.

Proof of quadratic convergence and quasi-optimal running time

It remains to prove that the complexity is the right one. If P ≥ 25 Im(τ ), log 2 P > log 2 Im(τ ) + 4.7, which means s ≤ log P and the cost of Steps 11 to 18 is O(M(P) log P ). We verify that conditions (3) and ( 16) hold:

| Re(τ 2 )| ≤ 1/2 s+2 ≤ 1/4, 0. 345 ≤ √ 3 4 ≤ Im(τ 2 ) ≤ 1, | Re(z 2 )| ≤ 1/2 s+3 ≤ 1/8, 0 ≤ Im(z 2 ) ≤ Im(τ 2 ) 4 .
This means the choices of signs are always good, and hence our result is indeed (squares of) theta-functions and theta-constants. We also need prove that there is a C > 1 such that, for all z 2 , τ 2 that we consider, . Furthermore, with a careful analysis, one can prove that c 1 = 55 is enough in Theorem 3.10.

θ
In any case, this proves that (λ n ) is quadratically convergent. We note that the fact that z 2 , τ 2 are within a compact shows that the constants b 1 , b 2 , b 3 exist and are independent of z, τ . This makes the running time of Step 7 only dependent in P , which was the point of the uniform algorithm. In particular, the number of bits lost during the computation of F ∞ or in F can be written as c 1 log P + c 2 , with c 1 , c 2 constants independent in z, τ . Hence, the number of bits that are lost in the whole of Step 7 is This means the computations in step 7 should be carried out at precision P = P + G log P + H, so that the result is accurate with P bits. This gives a running time of O(M(P ) log P ), independently of z and τ . All the other steps cost no more than O(M(P)) bit operations. Given the formula for P in subsections 4.3.1 to 4.3.4, this indeed gives us a running time of O(M(P ) log P ).

Implementation

An implementation using the GNU MPC library [START_REF] Enge | GNU MPC -A library for multiprecision complex arithmetic with exact rounding[END_REF] for arithmetic on multiprecision complex numbers was developed; we compared our algorithm to our own implementation of Algorithm 1 using MPC2 . The code is distributed under the GNU General public license, version 3 or any later version (GPLv3+); it is available at the address http://www.hlabrande.fr/pubs/fastthetas.tar.gz

We compared those implementations to MAGMA's implementation of the computation of θ(z, τ ) (function Theta). Each of those implementations computed θ(z, τ ) for z = 0.123456789 + 0.123456789i and τ = 0.23456789 + 1.23456789i at different precisions; the computations took place on a computer with an Intel Core i5-4570 processor. The results are presented in Figure 1 and Table 1.

Our figures show that our algorithm outperforms Magma even for computations at relatively low (i.e. 1000 digits) precision3 , and the naive algorithm for more than 325000 digits of precision. Hence, a combined algorithm which uses the naive method for precisions smaller than 325000 digits, and our method for larger precision, will yield the best algorithm, and outperform Magma in all cases, as shown in Table 1.

Definition 2 . 1 .

 21 Define, for z ∈ C and τ ∈ H (i.e. Im τ > 0) θ(z, τ ) = n∈Z exp πiτ n 2 + 2πinz .

π 8 min

 8 (|a|,|b|) [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF] Thm. 1]. Quadratic convergence implies that the number of exact digits approximately doubles with each iteration, so that one only needs O(log P ) iterations to compute AGM(a, b) with precision P ; hence the total cost of computing AGM(a, b) is O(M(P ) log P ) bit operations.

Lemma 3 . 2 .- 1 τ

 321 3 and Theorem 3.4 under this assumption. The second condition complements the first one as follows: Let z, τ such that | Re(τ )| ≤ 1 2 and the conditions of (16) are verified. Then z τ , verify the first condition of (16).

Algorithm 2 1 :

 21 Compute F ∞ (x, y, z, y) Work at precision P. 2: n ← 0 3: while |z -t| ≤ 2 -P -n-c1 do

Proposition 4 . 1 .

 41 Let τ be such that | Re(τ )| ≤ 0.5, Im(τ ) ≥ 0. 345 and Im -1 τ ≥ 0. 345, and let z be such that the conditions (16) are satisfied. Then

7 :a 12 0 a 22 - 1 8: end while 9 :

 719 (s, t) ← (s, t) -(F(s, t) -(z, τ )) a 11 (a, b) ← F ∞ (1, s, 1, t) 10: (a, b) ← (1/a, 1/b),(s, t) = (sa, tb) 11: Return (a, s, b, t).

Theorem 4 . 5 .

 45 Assume Im(τ ) > √ 3/2. Then 0.859 ≤ |θ 00,01 (0, τ )| ≤ 1.141, |θ 10 (0, τ )| ≤ 1.018

δ

  + h + log(p/2 i ) ≤ G log P + Hsince the number n of steps in Newton's method is O(log P ).

Figure 1

 1 Figure 1: Timing results

table :

 :

		a	b	c	d	σ(00, 01, 10)
		odd	even	even	odd	(00, 01, 10)
		odd	odd	even	odd	(01, 00, 10)
		odd	even	odd	odd	(10, 01, 00)
		even	odd	odd	even	(00, 10, 01)
		odd	odd	odd	even	(10, 00, 01)
		even	odd	odd	odd	(01, 10, 00)
	Proof. Define for any γ =	a b c d	∈ SL 2 (Z):

Table V

 V

	, p.
	36]

  1,0 | are bounded by 4. We also have |q| ≤ 0.07, and |q 2 | ≤ |q| n 2 ≤ |q| n = |q 1 | ≤ |q| ≤ 0.07. As for the v i , we have v 0 = 2, and for n

  Let (x 0 , y 0 , z 0 , t 0 ) ∈ C 4 . Put (x n+1 , y n+1 , z n+1 , t n+1 ) = F (x n , y n , z n , t n ) for any integer n ∈ N,and suppose this is an optimal F sequence. Then there exists positive real constants c, C such ∀n ≥ 1, c ≤ |x n |, |y n |, |z n |, |t n | < C.

	xn z∞	2 n	× z ∞	converges
	quadratically, whence F ∞ can be computed in O(M(P ) log P ) bit operations.	n∈N
	Lemma 3.8.

  The sequence (λ n ) converges, to a limit λ. Furthermore, for P large enough, there exists a constant c 1 > 0, depending on C, c and |λ|, such that, if k is the first integer such that |z k

	Theorem 3.10.
						choice of signs are good
	and hence | √ x n +	√ y n | 2 ≥ 2c. The proof is the same for | √ z n +	√	t n |.
	We now prove that (λ n ) =	xn z∞	2 n	× z ∞	converges quadratically, by proving the
	following theorem:				n∈N

  This means the number of bits lost during this step is bounded by log |q|

	8	= π 8 log 2 e Im(τ ).

1

√ |z| . We prove in the next subsection that |θ 00,01 (0, 2 i τ 1 )| ≥ 0.859 for i = [1 . . . s]. Hence, each extraction of square root loses at most 4 bits: step 9 loses 4 bits, and step 16 loses 4s ≤ 4 log P bits.

Step 14 loses more bits since θ 10 (0, τ ) is smaller; indeed, |θ 10 (0, τ )| ∼ |q| 1/4 .

  as in Steps 10 and 17. Then |w| -2n ≤ e nπ Im(τ )/4 and 0.8038 ≤ |θ 00,01 (z, τ )| ≤ 1.1962 |θ 10 (z, τ )| ≤ 1.228 • Suppose that 0 ≤ Im(z) ≤ Im(τ ) 4 , as in Steps 20. Then |w| -2n ≤ e nπ Im(τ )/2 and 0.6772 ≤ |θ 00,01 (z, τ )| ≤ 1.3228 |θ 10 (z, τ )| ≤ 1.543 Proof. The bounds on the theta-constants come from [7, p. 5], which proves |θ 00,01 (0, τ ) -1| ≤ 2|q| 1-|q| . The techniques are the same as the proof of Lemma 3.3 or Theorem 2.6. This gives in the first case |θ 00,01 (z, τ ) -1| ≤ |q| 3/4 + |q| + |q| 7/2 + |q| 4 + |q| 8.25 + |q| 9 + -q 3/8 ≤ 0.009 so |θ 10 (z, τ )| ≤ |q| 1/4 (|w| + |w| -1 ) + 0.009 ≤ 1.228. In the second case: |θ 00,01 (z, τ ) -1| ≤ |q| 1/2 + |q| +

						|q| 15
	≤ 0.1962	since Im(τ ) ≥	√	3/2	1 -|q|
		|q| 3 1 -|q|	≤ 0.3228
	|θ q 5/4 1 -q	≤ 0.0357
	Combining these bounds with formulas from Theorem 2.9 gives the following bounds:
	error(θ 00 (2 i z 1 , 2 i+1 τ 1 )) ≤ (20050.518 + 1818.032k θ z 01 + 1966.823k θ z 10 + 33516k θ 0 00 )2 -P
	error(θ 01 (2 i z 1 , 2 i+1 τ 1 )) ≤ (20050.518 + 1818.032k θ z 00 + 1966.823k θ z 10 + 33516k θ 0 01 )2 -P

|θ 10 (z, τ ) -q 1/4 (w + w -1 )| ≤ q 15/8 1 10 (z, τ ) -q 1/4 (w + w -1 )| ≤ |q| 5/4 + |q| 9/4 + ... ≤

  This is a direct consequence of the fact that z 2 , τ 2 are within a compact set that does not contain any zero of θ(z, τ ); hence one can write (non-zero) lower and upper bounds for any of the values of theta. If one wants to be a bit more precise, using the same reasoning as in Theorem 4.5, we have for √ 3/4 ≤ Im(τ ) ≤ 1:|θ 00,01 (z, τ ) -1| ≤ |q| 1/2 + |q| +

	2 01 (0, τ 2 ) θ 2 00 (0, τ 2 )	≤ C,	θ 2 10 (0, τ 2 ) θ 2 00 (0, τ 2 )	≤ C,	θ 2 01 (z 2 , τ 2 ) θ 2 00 (z 2 , τ 2 )	≤ C,	θ 2 10 (z 2 , τ 2 ) θ 2 00 (z 2 , τ 2 )	≤ C.
			|q| 3 1 -|q|	≤ 0.7859,	|θ 10 (z, τ )| ≤ 1 + |q| 1/4 +	|q| 5/4 1 -|q|	≤ 1.958.
	This gives C ≤ 83.64 and c ≥ 0.042 2 1.7859 2	1 1808

Table 1 :

 1 Times (in s) of different methods

		10 7
			Our algo (low prec=9000)
		10 5		Naive Magma
	Time (s)	10 1 10 3
			Base 10 precision
		Prec This work Naive Magma
		4000	0.092	0.032	0.1740
		8000	0.298	0.112	0.8719
		16000	0.868	0.384	4.358
		32000	2.399	1.347	22.70
		64000	6.778	4.598	116.4
		128000	18.32	15.29	606.1
		256000	45.54	41.56
		325000	62.74	63.90
		512000	111.78	129.8
		1024000	263.7	390.3
		2048000	625.4	1275
		4096000	1468	3921

Incidentally, this is why we consider only absolute precision in this paper.

√ yn√ xn > 0;

(0,τ ) . Then (F•u)(z, τ ) = (z, τ ) and we thus have

The naive algorithm which only computes θ(z, τ ) is only 5% faster; furthermore since Algorithm 5 computes all 4 values, it is fair to compare it to Algorithm 1.

This is even though Magma only returns θ(z, τ ), when our algorithm returns

values.

Compute θ 00,01,10 (z, τ ), θ 00,01,10 (0, τ ) with precision P using the naive method (Algorithm 1 + Section 2.2.4).

3: else

4:

Take s ∈ N such that 1 ≤ |τ |/2 s < 2

5:

Put τ 1 = τ 2 s and z 1 = z 2 s , so that Im(z 1 ) ≤ Im(τ 1 )/2.

6:

Put z 2 = z 1 /4 and τ 2 = τ 1 /2.

7:

Compute approximations of absolute precision P of θ 2 00 (z 2 , τ 2 ), θ 2 01 (z 2 , τ 2 ), θ 2 00 (0, τ 2 ), and θ 2 01 (0, τ 2 ) using Algorithm 4.

8:

Compute θ 2 00 (z 2 , τ 1 ), θ 2 01 (z 2 , τ 1 ), θ 2 00 (0, τ 1 ), θ 2 01 (0, τ 1 ) using Equation ( 12), and θ 2 10 (z 2 , τ 1 ) using Equation ( 14) and θ 2 10 (0, τ 1 ) using its equivalent in z = 0.

9:

Compute θ 00,01,10 (0, τ 1 ).

10:

Compute θ 00,01 (z 1 /2, τ 1 ) using Equation (17).

11:

for i = 1 .. s do

12:

Compute θ 2 00 (0, 2 i τ 1 ), θ 2 01 (0, 2 i τ 1 ) using the AGM.

13:

Compute θ 2 00 (2 i-2 z 1 , 2 i τ 1 ), θ 2 01 (2 i-2 z 1 , 2 τ 1 ) using Equation [START_REF] Luther | Reliable computation of elliptic functions[END_REF].

14:

If i = s, compute also θ 2 10 (0, 2 i τ 1 ) using the equivalent of Equation ( 14) in z = 0, then θ 10 (0, 2 i τ 1 ) by taking the square root.

15:

Compute θ 2 10 (2 i-2 z 1 , 2 i τ 1 ) using Equation [START_REF] Vallée | Gauss' algorithm revisited[END_REF].

16:

Compute θ 00,01 (0, 2 i τ 1 ).

17:

Compute θ 00 (2 i-1 z 1 , 2 i τ 1 ), θ 01 (2 i-1 z 1 , 2 i τ 1 ) using Equation (17).

18:

end for 19:

Compute θ 2 10 (2 s-1 z 1 , 2 s τ 1 ) using Equation [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF].

20:

Compute θ 00,01,10 (z, τ ) using Equation (17). 21: end if