Enea Milio

Damien Robert

Modular polynomials on Hilbert surfaces

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Context

Isogenies play an important role in elliptic curve cryptography. They allow to transfer the discrete logarithm problem from one curve to a possibly weaker one [START_REF] Galbraith | Extending the GHS Weil descent attack[END_REF][START_REF] Smith | Isogenies and the Discrete Logarithm Problem in Jacobians of Genus 3 Hyperelliptic Curves[END_REF]; they are used by the Schoof-Elkies-Atkin (SEA) point counting algorithm [START_REF] Schoof | Counting points on elliptic curves over finite fields[END_REF][START_REF] Morain | Calcul du nombre de points sur une courbe elliptique dans un corps fini: aspects algorithmiques[END_REF][START_REF] Elkies | Elliptic and modular curves over finite fields and related computational issues[END_REF], but also by the CRT algorithms to compute class polynomials [START_REF] Sutherland | Computing Hilbert class polynomials with the Chinese remainder theorem[END_REF][START_REF] Enge | Class invariants by the CRT method[END_REF] and modular polynomials [START_REF] Bröker | Modular polynomials via isogeny volcanoes[END_REF]. Splitting the multiplication using isogenies can improve the efficiency of the arithmetic [START_REF] Doche | Efficient scalar multiplication by isogeny decompositions[END_REF][START_REF] Gaudry | Fast genus 2 arithmetic based on Theta functions[END_REF], taking isogenies to mask the points reduces the impact of side channel attacks [START_REF] Smart | An analysis of Goubin's refined power analysis attack[END_REF], and isogenies are used to construct normal basis of a finite field [START_REF] Couveignes | Elliptic periods for finite fields[END_REF]. They have also been used to construct hash functions [START_REF] Charles | Cryptographic hash functions from expander graphs[END_REF] and to build cryptosystems [START_REF] Teske | An elliptic curve trapdoor system[END_REF][START_REF] Rostovtsev | Public-key cryptosystem based on isogenies[END_REF], in particular quantum-resistant cryptosystems [START_REF] De Feo | Towards quantum-resistant cryptosystems from supersingular elliptic curve isogenies[END_REF].

In dimension 1, the -modular polynomials φ parametrize pairs of elliptic curves E 1 and E 2 that are -isogenous over the algebraic closure. They can be computed in quasi-linear time [START_REF] Enge | Computing modular polynomials in quasi-linear time[END_REF] by the evaluation/interpolation method. More precisely the classical modular polynomials parametrize the elliptic curves from their j-invariants, so that E 1 and E 2 are -isogenous whenever φ (j(E 1 ), j(E 2 )) = 0. Other modular invariants have been proposed which yield smaller polynomials [START_REF] Enge | Comparing invariants for class fields of imaginary quadratic fields[END_REF].

Principally polarized complex abelian surfaces (which are generically Jacobians of hyperelliptic curves) are parametrized by the Siegel threefold Sp 2g (Z)\H g (with g = 2) where H g is the Siegel space of symmetric g × g complex matrices with totally positive imaginary part. The Siegel threefold is an algebraic variety birationally equivalent to the three dimensional projective space, and is parametrized by the three Igusa invariants [START_REF] Igusa | Arithmetic variety of moduli for genus 2[END_REF][START_REF] Igusa | On Siegel modular forms of genus 2[END_REF]. One can then generalize modular polynomials to this setting: the -modular polynomials classify pairs of principally polarized abelian surfaces (A, B) which admit an -isogeny A → B. More precisely the -modular polynomials evaluated on the three Igusa invariants of A describe a dimension 0 subvariety of the Siegel threefold of degree 3 + 2 + + 1 whose geometric points correspond to the triples of Igusa invariants of the -isogenous abelian surfaces B. Alternatively, these modular polynomials describe the image of X 0 ( ) inside X 0 (1) × X 0 (1) where X 0 ( ) = Γ 0 ( )\H g and Γ 0 ( ) = a b c d ∈ SL 2 (Z) : |b . These polynomials have been studied in [START_REF] Gaudry | Algorithmique des courbes hyperelliptiques et applications à la cryptologie[END_REF][START_REF] Bröker | Modular polynomials for genus 2[END_REF] and computed for = 2 in [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]. Generalizations of these modular polynomials using smaller Siegel modular invariants have more recently been computed in [START_REF] Milio | A quasi-linear time algorithm for computing modular polynomials in dimension 2[END_REF].

Unfortunately even using a quasi-linear algorithm, computing them is hard due to their size. Indeed compared to dimension 1 where modular polynomials describe a curve X 1 0 ( ) inside the plane X 1 0 (1) × X 1 0 (1), and where the degree of the projection is + 1, in dimension 2 they describe the threefold X 0 ( ) inside a dimension six space and the degree of the projection is 3 + 2 + + 1. Already these polynomials for = 7 takes 29 GB to write (even using the smaller theta invariants), so it seems hard to go much further. But having them only up to = 7 is not enough for most of the applications mentioned.

Another problem is that restricting to -isogenies does not allow one to explore the full isogeny graph of principally polarized abelian surfaces. In the CRT method to compute class polynomials, one key step of the algorithm is to take an abelian surface in the right isogeny graph, and then use isogenies to find an abelian surface with maximal complex multiplication [START_REF] Bröker | Explicit CM theory for level 2-structures on abelian surfaces[END_REF][START_REF] Lauter | Improved CRT Algorithm for Class Polynomials in Genus 2[END_REF]. But this is not always possible using only -isogenies [START_REF] Brooks | Isogeny graphs of ordinary abelian varieties[END_REF]Theorem 9.4].

We recall that an -isogeny f corresponds to a kernel V = Ker f which is maximal isotropic for the Weil pairing e on the -torsion A[ ]. The kernel of an -isogeny is then an abelian group of type ( , ). One can also consider cyclic isogenies, where the kernel is a cyclic subgroup of the -torsion. However, if A is principally polarized and V is cyclic in A[ ], then A/V is not principally polarizable in general. Indeed, the isogenous abelian surface admits a principal polarization if and only if there exists β ∈ End s,++ (A) a totally positive symmetric element of norm in the ring End s (A) of symmetric endomorphisms of A such that V ⊂ Ker β (since V is cyclic it is automatically isotropic for the β-Weil pairing). We call such an isogeny a β-isogeny, and one is naturally led to try to define β-modular polynomials parametrizing pairs of β-isogenous abelian surfaces (A, B). Generically, a complex abelian surface A has no (non trivial) symmetric endomorphisms, so to define β-modular polynomials we need to restrict to a sublocus of the moduli space of abelian surfaces with specific real multiplication.

Let O K be a maximal real quadratic order of discriminant ∆ K . The Hilbert moduli space is a surface parametrizing isomorphism classes of principally polarized abelian surfaces A together with an embedding O K → End s (A). Let β ∈ O K be a totally positive element of norm . In this article, we define β-modular polynomials on this Hilbert modular surface and we explain how to compute them by evaluation/interpolation. We use the forgetful map from the Hilbert modular surface to the Siegel space, or more precisely, to an Humbert surface, and use the tools already known there, especially those described in [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF][START_REF] Milio | A quasi-linear time algorithm for computing modular polynomials in dimension 2[END_REF] for the computation of -modular polynomials.

Results

We study several parametrizations of the Humbert surfaces. The Siegel moduli threefold is parametrized by the three Igusa functions, and in [START_REF] Milio | A quasi-linear time algorithm for computing modular polynomials in dimension 2[END_REF] a cover of the Siegel space given by level 2 theta constant is also used to give smaller modular polynomials.

Pulling back the Igusa functions to the Humbert surface gives rational coordinates which can be used to define modular polynomials. Likewise pulling back the theta functions give coordinates on a cover of the Humbert surface. Some Humbert surfaces are rational and can be parametrized by two invariants instead of the three defined above. In this paper we look in particular at the case of Humbert surfaces of discriminant 5 and 8 which can be parametrized by two Gundlach invariants.

For the modular polynomial computations, we need an algorithm for the evaluation step and one for the interpolation step. The first algorithm computes, given a period matrix τ ∈ H g 1 , the above invariants at some precision N in quasi-linear time in the precision N . The second one, given the value of the above invariants, computes a corresponding period matrix τ ∈ H g 1 in time quasi-linear in the requested precision. Instead of doing computation on Humbert surfaces, our idea is to translate back and forth between the Hilbert moduli space and the Siegel moduli space where in the latter space both algorithms have been developed by Dupont in [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]. The translations is based on the work of Igusa [START_REF] Igusa | Arithmetic variety of moduli for genus 2[END_REF][START_REF] Igusa | On Siegel modular forms of genus 2[END_REF], Gundlach [START_REF] Gundlach | Die Bestimmung der Funktionen zur Hilbertschen Modulgruppe des Zahlkörpers Q( √ 5)[END_REF][START_REF] Gundlach | Die Bestimmung der Funktionen zu einigen Hilbertschen Modulgruppen[END_REF], Resnikoff [START_REF] Resnikoff | On the Graded Ring of Hilbert Modular Forms Associated with Q( √ 5)[END_REF], Lauter-Yang [START_REF] Lauter | Computing genus 2 curves from invariants on the Hilbert moduli space[END_REF], Runge [START_REF] Runge | Endomorphism rings of abelian surfaces and projective models of their moduli spaces[END_REF], Birkenhake-Wilhelm [START_REF] Birkenhake | Humbert surfaces and the Kummer plane[END_REF]. In the case of Humbert surfaces of discriminant 5, we have inverted the formulas of [START_REF] Lauter | Computing genus 2 curves from invariants on the Hilbert moduli space[END_REF]Proposition 4.5] expressing the pullback of the Igusa invariants in function of the Gundlach invariants (see Appendix B.1). In the case of Humbert surfaces of discriminant 8, we link in Theorem A.14 and Corollary A.15 the Gundlach invariants with the pullback of the Igusa invariants, as was done by Resnikoff for discriminant 5. These algorithms are described in Section 3.2 (see also Theorem 3.4).

The main result of the paper is the computation of modular polynomials on the Hilbert (or Humbert) surface. When β ∈ O K is a totally positive prime, we define β-isogenies and β-modular polynomials in Section 4.1. There are two cases:

• If the norm of β is a prime number , then the β-isogenies correspond to isogenies with cyclic kernel

V ⊂ A[β] ⊂ A[ ].
All β-isogenies then preserve real multiplication and the β-modular polynomials parametrize all pairs of principally polarized abelian surfaces with maximal real multiplication and admitting a cyclic isogeny of degree ;

• Otherwise β is an inert prime number ∈ Z. In this case the -modular polynomials (on the Hilbert moduli space) parametrize -isogenies between abelian surfaces with maximal real multiplication. By contrast to the Siegel -modular polynomials which given A parametrize all 3 + 2 + + 1 abelian surfaces B = A/V where V ⊂ A[ ] is maximal isotropic for the Weil pairing, the Hilbert -modular polynomials parametrize all 2 + 1 abelian surfaces B = A/V where V is furthermore stable under the action of the real multiplication.

We give in Theorem 4.15 a quasi-linear algorithm for computing β-modular polynomials for a large class of invariants, like Gundlach invariants (for Q( √ 2) and Q( √ 5)), pullbacks of Igusa invariants and pullbacks of theta constants (for all real quadratic field). In the latter two cases we have three invariants for a moduli space of dimension 2 so we need to adapt the evaluation/interpolation algorithm to handle the fact that these three invariants have to satisfy a relation.

Theorem 4.15 is itself a particular case of Theorem 3.11 which gives an evaluation/interpolation algorithm to compute covers of Hilbert surfaces. Adapting this Theorem to the cover parametrizing β-isogenies then yields Theorem 4.15.

The corresponding algorithms have been implemented in Pari/GP, and we give some examples of β-modular polynomials. We mainly give examples in the case where

K = Q( √ 2) or Q( √ 5
) since this allows us to compare different kinds of invariants. Finally Martindale and Streng [START_REF] Martindale | Isogeny Graphs, Modular Polynomials, and Applications[END_REF] have also independently described an algorithm to compute modular polynomials on Hilbert moduli space. While we use evaluation/interpolation, they use linear algebra on the Fourier coefficients of the Hilbert modular form. The advantage of their method is that it works in any dimension and for any modular invariant (provided one can compute its Fourier coefficients). By contrast our evaluation/interpolation approach needs fast evaluation of modular invariants (for the complexity) and we need for the interpolation to be able to recover the period matrix from the values of the modular coefficients. We only know how to do that efficiently in dimension 2 (and 1) when the invariants are derived from theta constant (as mentioned by translating back and forth to the Siegel space and using [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]). In particular our algorithm can not be extended to higher dimension as long as the work of Dupont on the generalization of the AGM is not extended to dimension greater than 2. Work in this direction has been done by Labrande and Thomé in [START_REF] Labrande | Explicit computation of the Abel-Jacobi map and its inverse[END_REF][START_REF] Labrande | Computing theta functions in quasi-linear time in genus two and above[END_REF]. However in dimension 2, we do obtain a quasi-linear algorithm, while there is no complexity analysis in [START_REF] Martindale | Isogeny Graphs, Modular Polynomials, and Applications[END_REF], and our algorithm is practical since we have computed large polynomials. Our main limitation was the intrinsic size of the modular polynomials (see Tables 1,2,3) and not the speed of the computations.

Outline

The remainder of this article is organized as follows. In Section 2, we define the Siegel (in Section 2.1) and the Hilbert spaces (in Section 2.2) and describe the corresponding moduli data. We also give generators for the fields of modular functions on these spaces. Then in Section 2.3, we analyze the forgetful map from the Hilbert modular surface to the Siegel space. In Section 2.4, we focus on the Humbert surfaces. An Humbert surface is the image of a given Hilbert modular surface by previous map. We conclude this Section by looking at covers of an Humbert surface in Section 2.5.

Section 3 is concerned with invariants of Hilbert surfaces. In Section 3.2 we explain how to efficiently evaluate a large class of Hilbert invariants. In Section 3.3 we give an interpolation algorithm, which works even when we have relations between our invariants. Lastly we conclude the Section by giving in Section 3.4 an algorithm to compute covers of Hilbert surface.

Section 4 is concerned with modular polynomials on Hilbert surfaces. First in Section 4.1, we define the isogenies preserving real multiplication. In Section 4.3, we define the modular polynomials depending on these isogenies, explain some of their properties and give an algorithm to compute them in quasi-linear time.

Finally in Section 5, we describe some polynomials we have computed. See also Appendix B, where we apply the results of the previous sections on Gundlach invariants and pullbacks of the theta functions.

Thanks

We thank Pierre-Jean Spaenlehauer to have succesfully computed the Gröbner basis expressing the Gundlach invariants in term of the Igusa invariants for D = 2. We thank David Kohel who suggested us to look at the work of Elkies and Kumar [START_REF] Elkies | K3 surfaces and equations for Hilbert modular surfaces[END_REF] to get invariants for more Humbert or Hilbert surfaces (see Appendix A.4). We thank the anonymous reviewer for his careful reading and his many comments and suggestions.

Hilbert and Siegel modular spaces

Siegel modular space

The Siegel upper half-space in dimension 2 is the set

H 2 = {Ω ∈ M 2 (C) | Ω is symmetric and
(Ω) > 0}. It is a moduli space for principally polarized abelian surfaces with symplectic basis: such a surface is a torus C 2 /(Z 2 + ΩZ 2 ) for some Ω ∈ H 2 , and the principal polarization is induced by the Hermitian form given by (Ω) -1 (see [START_REF] Birkenhake | Complex abelian varieties[END_REF]Section 8.1]).

We define the symplectic group Sp 4 (Z) as {γ ∈ GL 4 (Z)

| t γJγ = J} where J = 0 I 2 -I 2 0
and I n is the identity matrix of size n. It acts on

H 2 by A B C D • Ω = (AΩ + B)(CΩ + D) -1 (it is a left action).
The Siegel modular threefold is the Baily-Borel compactification [START_REF] Baily | Compactification of Arithmetic Quotients of Bounded Symmetric Domains[END_REF] of the quotient space Sp 4 (Z)\H 2 and this quotient space is a moduli space for isomorphism classes of principally polarized abelian surfaces. See [START_REF] Birkenhake | Complex abelian varieties[END_REF]Section 8.2].

Let Γ be a subgroup of Sp 4 (Z) of finite index and k ∈ Z. A Siegel modular form of weight k for Γ is a holomorphic function f :

H 2 → C such that for all γ = A B C D ∈ Sp 4 (Z) and Ω ∈ H 2 , f (γΩ) = det(CΩ + D) k f (Ω).
The quotient of two Siegel modular forms for the same weight and group Γ is called a Siegel modular function. Note that if f is a Siegel modular function for Γ, and γ ∈ Sp 4 (Z), then f (γΩ) = f (ΓγΩ) so that we can consider the set of right cosets Γ\Sp 4 (Z). More generally, in the rest of the paper, we will always consider sets of right cosets for a similar reason.

Let

a, b ∈ {0, 1 2 } 2 . The classical theta constant with characteristic (a, b) is θ [ a b ] (Ω) = n∈Z g exp(iπ t (n + a)Ω(n + a) + 2iπ t (n + a)b).
To simplify the notation we define for all a = ( a 0 a 1 ) and b

= b 0 b 1 in {0, 1} 2 θ b 0 +2b 1 +4a 0 +8a 1 (Ω) := θ a/2 b/2 (Ω).
Of the 16 theta constants, 6 are identically zero and we denote P = {0, 1, 2, 3, 4, 6, 8, 9, 12, 15} the subscripts of the even theta constants (the non-zero ones). The following functions h i are Siegel modular forms of weight i for the symplectic group Sp 4 (Z)

h 4 = i∈P θ 8 i , h 6 = 60 triples (i,j,k)∈P 3 ±(θ i θ j θ k ) 4 , h 10 = i∈P θ 2 i , h 12 = 15 tuples (i,j,k,l,m,n)∈P 6 (θ i θ j θ k θ l θ m θ n ) 4
(see for example [39, Page 848], [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]Section 6.3.3] or [START_REF] Streng | Complex multiplication of abelian surfaces[END_REF]Section 7.1] for the exact definition). We define the Siegel Eisenstein series ψ k of even weight k ≥ 4 by

ψ k (Ω) = C,D det (CΩ + D) -k ,
where the sum is taken over the set of matrices A B C D in Sp 4 (Z) up to left multiplication by SL 2 (Z). Let χ 10 = -2 -12 3 -5 5 -2 7 -1 53 -1 43867(ψ 4 ψ 6 -ψ 10 ) and

χ 12 = 2 -13 3 -7 5 -3 7 -2 337 -1 131 • 593(3 2 7 2 ψ 3 4 + 2 • 5 3 ψ 2 6 - 691ψ 12
) be two Siegel modular cusps forms of weight 10 and 12 respectively. These series can be written in terms of theta constants. Indeed we have by Igusa [START_REF] Igusa | Modular Forms and Projective Invariants[END_REF]Page 848] that ψ 4 = 2 -2 h 4 , ψ 6 = 2 -2 h 6 , χ 10 = -2 -14 h 10 and χ 12 = 2 -17 3 -1 h 12 . The graded ring of holomorphic Siegel modular forms for Sp 4 (Z) is the polynomial ring of ψ 4 , ψ 6 , χ 10 and χ 12 . We define the Igusa invariants from these last modular forms:

(1) j 1 = 2 • 3 5 χ 5 .

The field of Siegel modular functions for Sp 4 (Z) is C(j 1 , j 2 , j 3 ). Generically, two principally polarized abelian surfaces are isomorphic if and only if they have the same three Igusa invariants.

(In [START_REF] Igusa | Arithmetic variety of moduli for genus 2[END_REF][START_REF] Igusa | On Siegel modular forms of genus 2[END_REF], Igusa has defined 10 absolute invariants, which allows to characterize all abelian surfaces up to isomorphisms over the algebraic closure of the base field.)

Let Γ(2) = { A B C D ∈ Sp 4 (Z) : A B C D ≡ I 4 mod 2}.
It is a normal subgroup of Sp 4 (Z) of index 720. The three following functions are Siegel modular functions for Γ(2) called the Rosenhain invariants. They are generators for the field of modular functions belonging to Γ(2) (see Mumford [55,Section 8]).

Let Γ(2, 4) = { A B C D ∈ Sp 4 (Z) : A B C D ≡ I 4 mod 2 and B 0 ≡ C 0 ≡ 0 mod 4}, where X 0 denotes the vector composed of the diagonals elements of X. It is a normal subgroup of Sp 4 (Z) of index 11520. The quotients of theta functions

(3) b i (Ω) = θ i (Ω/2)/θ 0 (Ω/2)
for i = 1, 2, 3 are Siegel modular functions for Γ [START_REF] Becker | On Gröbner bases under specialization[END_REF][START_REF] Birkenhake | Humbert surfaces and the Kummer plane[END_REF] and they are generators for the field of modular functions belonging to Γ(2, 4) (see Manni [START_REF] Manni | Modular varieties with level 2 theta structure[END_REF]Theorem 1]).

Hilbert modular space

We refer to [START_REF] Van Der | Hilbert modular surfaces[END_REF][START_REF] Bruinier | Hilbert modular forms and their applications[END_REF][START_REF] Goren | Lectures on Hilbert modular varieties and modular forms[END_REF][START_REF] Freitag | Hilbert modular forms[END_REF][START_REF] Nagaoka | On the ring of Hilbert modular forms over Z[END_REF] for more details on Hilbert modular forms and Hilbert surfaces.

Let D > 0 be a square-free integer and

K = Q( √ D) be a real quadratic field. Its discrim- inant ∆ K is D if D ≡ 1 mod 4 and 4D if D ≡ 2, 3 mod 4. Consider O K the ring of integers of K. We have that O K = Z + ωZ where ω = 1+ √ D 2 if D ≡
H + 1 × H - 1 by a b c d • (τ 1 , τ 2 ) = ( aτ 1 +b cτ 1 +d , aτ 2 +b cτ 2 +d
). The Baily-Borel compactification of the quotient space SL 2 (O K )\H + 1 × H - 1 is the Hilbert modular surface (see for example [START_REF] Van Der | Hilbert modular surfaces[END_REF]Section 7]). It parametrizes principally polarized abelian surfaces (A, θ) with real multiplication by the maximal order O K , with an explicit embedding µ : O K → End(A) (see [START_REF] Birkenhake | Complex abelian varieties[END_REF]Chapter 9]).

For computations, it can be convenient to work over another model, which we describe here. Let SL

2 (O K ⊕ ∂ K ) = { a b c d ∈ SL 2 (K) : a, d ∈ O K , b ∈ ∂ -1 K and c ∈ ∂ K }, where ∂ K is the different ideal of K. As K = Q( √ D), we have that ∂ K = √ ∆ K O K and ∂ -1 K = 1 √ ∆ K O K . We have group isomorphisms (4) φ ± : SL 2 (O K ) → SL 2 (O K ⊕ ∂ K ) a b c d → a b/ √ ∆ K c √ ∆ K d φ ± : H + 1 × H - 1 → H 2 1 (τ 1 , τ 2 ) → (τ 1 √ ∆ K , -τ 2 √ ∆ K )
and the following diagram is commutative

(5) SL 2 (O K ) × H + 1 × H - 1 / / H + 1 × H - 1 SL 2 (O K ⊕ ∂ K ) × H 2 1 / / H 2 1 (see [18, Section 3]). If τ = (τ 1 , τ 2 ) ∈ H 2 1
, the corresponding abelian surface is given by the torus

C 2 /(Φ(O K ) ⊕ τ 1 0 0 τ 2 Φ(∂ -1 K ))
where Φ : K → C 2 is given by the two real embeddings, and the polarization is induced by the symplectic form E on the lattice:

E(x 1 + x 2 τ, y 1 + y 2 τ ) = tr K/Q (x 1 y 2 -x 2 y 1 ). From the definition of ∂ -1
K we get indeed that E induces a principal polarization.

Since SL 2 (O K ) is generated by the matrices ( 1 1 0 1 ), (

1 ω 0 1 ), 0 -1 1 0 , the group SL 2 (O K ⊕ ∂ K ) is generated by the matrices 1 1/ √ ∆ K 0 1 , 1 ω/ √ ∆ K 0 1 and 0 -1/ √ ∆ K √ ∆ K 0 . For λ ∈ K and τ = (τ 1 , τ 2 ) ∈ H 2 1 , we denote λτ = (λτ 1 , λτ 2 ), N (τ ) = τ 1 τ 2 and tr(τ ) = τ 1 + τ 2 .
We define σ to be the involution σ

: (τ 1 , τ 2 ) ∈ H 2 1 → (τ 2 , τ 1 ) ∈ H 2 1 . We let σ act by conjugation on SL 2 (O K ⊕ ∂ K ) via σγσ = a b c d ∈ SL 2 (O K ⊕ ∂ K ), for γ = a b c d ∈ SL 2 (O K ⊕ ∂ K ).
= (τ 1 , τ 2 ) ∈ H 2 1 the condition f (γτ ) = N (cτ +d) k f (τ ). If moreover it satisfies f (σ(τ )) = f (τ ) for all τ ∈ H 2
1 , then we say that this form is symmetric. A Hilbert modular function is the quotient of Hilbert modular forms of the same weight and for the same group. We say it is symmetric when the forms are, for some choice of forms.

Remark 2.2.

A modular form f is then automatically holomorphic at the set of cusps SL 2 (O K )\P 1 (K) Cl(O K ) (see for example [9, Section 1.3] for more details).

For the study of Humbert surfaces in Section 2.4 we will also be interested in symmetric Hilbert modular forms and functions.

From Hilbert to Siegel

Let τ = (τ 1 , τ 2 ) ∈ H 2 1 , x ∈ K and γ = a b c d ∈ SL 2 (K). We denote τ * = τ 1 0 0 τ 2 , x * = x 0 0 x and γ * = a * b * c * d * . Fix {e 1 , e 2 } a Z-basis of O K = Z + ωZ (ω is
(K) → Sp 4 (Q) γ → Sγ * S -1 . Recall that SL 2 (O K ⊕ ∂ K ) = { a b c d ∈ SL 2 (K) : a, d ∈ O K , b ∈ 1/ √ ∆ K O K and c ∈ √ ∆ K O K }.
Proposition 2.3. The map φ e 1 ,e 2 satisfies:

• φ -1 e 1 ,e 2 (Sp 4 (Z)) = SL 2 (O K ⊕ ∂ K ); • φ e 1 ,e 2 (γ • τ ) = φ e 1 ,e 2 (γ) • φ e 1 ,e 2 (τ ) for all γ ∈ SL 2 (O K ⊕ ∂ K ) and τ ∈ H 2 1 ; • If f 1 , f 2 is another Z-basis of O K , then there exists some γ ∈ Sp 4 (Z) such that for all τ ∈ H 2 1 , φ e 1 ,e 2 (τ ) = γ • φ f 1 ,f 2 (τ );
• There exists some γ ∈ Sp 4 (Z) such that φ e 1 ,e 2 (σ(τ )) = γ • φ e 1 ,e 2 (τ ). We denote M σ this γ, and this allows us to extend φ e 1 ,e 2 to SL

2 (O K ⊕ ∂ K ) σ .
Proof. See for example [START_REF] Lauter | Computing genus 2 curves from invariants on the Hilbert moduli space[END_REF]Proposition 3.1]. For the convenience of the reader, we give more details on the key part of the proof in the reference above.

Let A = C 2 /(Φ(∂ -1 K ) ⊕ τ * Φ(O K
)) be the abelian surface corresponding to τ . Then Φ(O K ) = RZ 2 and if (f 1 , f 2 ) denotes the dual basis with respect to the trace of (e 1 , e 2 ), then

Φ(∂ -1 K ) = R ∨ Z 2 with R ∨ = f 1 f 2 f 1 f 2 . The lattice Φ(∂ -1 K ) ⊕ τ * Φ(O K ) = R ∨ Z 2 ⊕ τ * RZ 2 is isomorphic to Z 2 ⊕ R ∨,-1 τ * RZ 2 via the change of variable z → R ∨,-1 z. But an easy computation shows that t RR ∨ = 1 so R ∨,-1 = t R.
Thus, the map φ e 1 ,e 2 gives a holomorphic map from SL 2 (O K ⊕∂ K )\H 2 1 to Sp 4 (Z)\H 2 which is independent of the choice of the basis of O K . It also sends τ and σ(τ ) to the same point of Sp 4 (Z)\H 2 . Since φ e 1 ,e 2 allows us to identify SL 2 (O K ⊕ ∂ K ) and σ as subgroups of Sp 4 (Z), we will note SL

2 (O K ⊕ ∂ K ) ∪ SL 2 (O K ⊕ ∂ K )σ the group SL 2 (O K ⊕ ∂ K )
σ . We will often work with the basis e 1 = 1 and e 2 = ω. We have

φ 1,ω (τ ) = τ 1 +τ 2 τ 1 ω+τ 2 ω τ 1 ω+τ 2 ω τ 1 ω 2 +τ 2 ω 2 = Ω 1 Ω 2 Ω 2 Ω 3 ∈ H 2 and it satisfies (7) D-1 4 Ω 1 + Ω 2 -Ω 3 = 0, if D ≡ 1 mod 4; DΩ 1 -Ω 3 = 0, if D ≡ 2, 3 mod 4. Moreover, set (8)
M σ =                1 0 0 0 1 -1 0 0 0 0 1 1 0 0 0 -1 if D ≡ 1 mod 4; 1 0 0 0 0 -1 0 0 0 0 1 0 0 0 0 -1 if D ≡ 2, 3 mod 4.
We have found that the matrix M σ satisfies

(9) φ 1,ω (σ(τ )) = M σ • φ 1,ω (τ ). Consider now γ = a+a ω (b+b ω)/ √ ∆ K √ ∆ K (c+c ω) d+d ω ∈ SL 2 (O K ⊕ ∂ K ). Then (10) φ 1,ω (γ) =                       a a b b+b ( D-1 4 )a a+a b+b b+( D+3 4 )b ( D-1 4 )c -c c d ( D-1 4 )d c c d d+d    if D ≡ 1 mod 4;

Humbert surfaces

Humbert surfaces have been first studied by Humbert in [START_REF] Humbert | Sur les fonctions abéliennes singulières I[END_REF][START_REF] Humbert | Sur les fonctions abéliennes singulières II[END_REF][START_REF] Humbert | Sur les fonctions abéliennes singulières III[END_REF]. Let Ω = Ω 1 Ω 2 Ω 2 Ω 3 ∈ H 2 and a, b, c, d, e ∈ Z. We call an equation of the form: a,b,c,d,e) = 1, we say that this relation is primitive. Moreover, we define the discriminant of a singular relation to be ∆ = b 2 -4ac -4de.

aΩ 1 + bΩ 2 + cΩ 3 + d(Ω 2 2 -Ω 1 Ω 3 ) + e = 0 a singular relation. If gcd (

Theorem 2.4 (Humbert's Lemma). Let

Ω = Ω 1 Ω 2 Ω 2 Ω 3
satisfying the singular relation:

aΩ 1 + bΩ 2 + cΩ 3 + d(Ω 2 2 -Ω 1 Ω 3 ) + e = 0 of discriminant ∆ = b 2 -4ac -4de. Then there exists a matrix γ ∈ Sp 4 (Z) such that γ • Ω = Ω 1 Ω 2 Ω 2 Ω 3
satisfies a normalized singular relation of the form:

(11) kΩ 1 + Ω 2 -Ω 3 = 0
where k and are determined uniquely by ∆ = 4k + and ∈ {0, 1}.

Proof. See [START_REF] Humbert | Sur les fonctions abéliennes singulières I[END_REF][START_REF] Humbert | Sur les fonctions abéliennes singulières II[END_REF][START_REF] Humbert | Sur les fonctions abéliennes singulières III[END_REF] Let also ∆ = ∆ be non-square discriminants. Then:

• A Ω is simple if and only if Ω ∈ m>0 H m 2 ;
• If Ω ∈ H ∆ , then End(A Ω ) contains the maximal order of Q( √ ∆), i.e., equivalently, there exists a symmetric endomorphism of discriminant ∆ on A Ω ; We denote now Γ(1) = SL 2 (O K ⊕ ∂ K ). Proposition 2.3 and Equations ( 7), ( 8) and ( 9) say that the images by φ 1,ω of H 2 1 and of ( Γ(1) ∪ Γ(1)σ)\H 2 1 are in the Humbert surface of discriminant ∆ K . This is also true for any φ e 1 ,e 2 because the images of τ by φ 1,ω and by φ e 1 ,e 2 are equivalent modulo the action of Sp 4 (Z) (which means that these maps send τ to the same point of the Humbert surface). More precisely, the Hilbert surface maps onto the Humbert surface:

• if Ω ∈ H ∆ ∩ H ∆ , then either A Ω is simple and End(A Ω ) ⊗ Q is a totally indefinite quaternion algebra over Q, or A Ω is isogenous to E × E,
Proposition 2.8. The maps φ e 1 ,e 2 of Equation ( 6) give rise to the following commutative diagram :

Γ(1)\H 2 1 π ( ( H 2 1 o o φe 1 ,e 2 / / H 2 ( Γ(1) ∪ Γ(1)σ)\H 2 1 ρ / / Sp 4 (Z)\H 2
where π is a map of degree 2 and ρ is a map generically of degree 1 onto the Humbert surface H ∆ K .

Proof. See Van der Geer [70, Page 328] for the map ρ and the commutativity of the rectangular part of the diagram. The fact that π is of degree 2 is obvious. It remains to see that ρ • π is generically of degree 2. But H ∆ K is the locus of principally polarized abelian surfaces (A, θ) with real multiplication by O K , and the preimages correspond to explicit embeddings µ : O K → End(A). Generically there are only two such embeddings which differ by the real conjugation, which corresponds to the action of σ.

The analytic quotient space ( Γ(1)∪ Γ(1)σ)\H 2 1 is called a symmetric Hilbert modular surface; it is birational to the Humbert surface.

Lemma 2.9. The pullbacks by ρ of the Igusa invariants to the symmetric Hilbert modular surface ( Γ(1) ∪ Γ(1)σ)\H 2 1 generate the function field of symmetric Hilbert modular functions. (These pullbacks can also be seen as the restriction of the Igusa invariants to the Humbert surface).

Proof. This is a well-known result. We give the proof in Appendix A.1.

Symmetric and non symmetric covers of the Humbert surface

We study here the covers of the Hilbert modular surface SL

2 (O K ⊕∂ K )\H 2 1 given by a subgroup Γ of finite index in SL 2 (O K ⊕ ∂ K ).
When Γ = SL 2 (Z), the subgroups Γ(n), Γ 0 ( ) and Γ(2, 4) are standard, and of main interest for modular polynomials of elliptic curves. We want to generalize these notations to the Hilbert modular group. It is easier to define them first in the model of SL 2 (O K ) acting on H + × H - and then transport them to the model of SL

2 (O K ⊕ ∂ K ) action on H 2 via the automorphism φ ± of Equation (4). Definition 2.10. Let (12) Γ(n) = a b c d ∈ SL 2 (O K ) : a ≡ d ≡ 1 mod n, b ≡ c ≡ 0 mod n .
Define then for D ≡ 1 mod 4 and D ≡ 2, 3 mod 4

(13) Γ(2, 4) = a b c d ∈ Γ(2) : b ≡ c ≡ 0 mod 4 , (14) Γ(2, 4) = a (b 1 +b 2 ω) (c 1 +c 2 ω) d ∈ Γ(2) : b 2 ≡ c 2 ≡ 0 mod 4 respectively.
By abuse of notation, we use the same notation for their image by φ ± :

(15) Γ(n) = a b/ √ ∆ K √ ∆ K c d ∈ SL 2 (O K ⊕ ∂ K ) : a ≡ d ≡ 1 mod n, b ≡ c ≡ 0 mod n .
Define then for D ≡ 1 mod 4 and D ≡ 2, 3 mod 4

(16) Γ(2, 4) = a b/ √ ∆ K √ ∆ K c d ∈ Γ(2) : b ≡ c ≡ 0 mod 4 , (17) Γ(2, 4) = a (b 1 +b 2 ω)/ √ ∆ K √ ∆ K (c 1 +c 2 ω) d ∈ Γ(2) : b 2 ≡ c 2 ≡ 0 mod 4
respectively. Note the subtlety in the definition of Γ(2, 4) for D ≡ 2, 3 mod 4.

Remark 2.11. By Serre [START_REF] Serre | Le Problème des Groupes de Congruence Pour SL 2[END_REF] a group Γ of finite index in SL 2 (O K ⊕ ∂ K ) is necessarily a congruence subgroup, meaning that it contains the congruence subgroup Γ(n) for some positive integer n.

Lemma 2.12. Let G be a subgroup of SL

2 (O K ⊕ ∂ K ) σ of finite index. If σ ∈ G then G ⊂ SL 2 (O K ⊕ ∂ K ). Otherwise G = Γ σ for a subgroup Γ ⊂ SL 2 (O K ⊕ ∂ K )
of finite index and normalized by σ (meaning that Γ is stable under the real conjugation). In the latter case we say that G is symmetric.

Proof. Indeed as a set it is easy to see

that if σ ∈ G, then G = Γ ∪ Γσ for a subgroup Γ ⊂ SL 2 (O K ⊕ ∂ K ). It remains to check that σ normalizes Γ. But since G is a group, Γ = σ Γσ -1 ⊂ G, so Γ = Γ.
Definition 2.13. We denote by C G the field of meromorphic functions of H 2 1 invariant under the action of G. It is the function field of the Hilbert surface

H G = G\H 2 1 .
Remark 2.14. H G admits a Baily-Borel compactification [START_REF] Baily | Compactification of Arithmetic Quotients of Bounded Symmetric Domains[END_REF], which in turn admits a smooth birational model. In this article we only work with functions of the Hilbert modular function field, so only up to birational equivalence, so we do not distinguish between these models.

Consider now Γ a subgroup of Sp 4 (Z) of finite index. The projection π : 8) is not in Γ, then G = Γ, otherwise G = Γ ∪ Γσ. By Proposition 2.8 we get that the following diagram is commutative:

Γ\H 2 → Sp 4 (Z)\H 2 is a finite map. Recall that if ∆ K is the discriminant of O K , we denote by H ∆ K the Humbert surface of discriminant ∆ K . An irreducible component of H Γ ∆ K = π -1 (H ∆ K ) in Γ\H 2 is called a Humbert surface component. Let G = φ -1 1,ω (Γ) and Γ = G ∩ SL 2 (O K ⊕ ∂ K ). If the matrix M σ of Equation (
H 2 1 φ 1,ω / / H 2 G\H 2 1 ρ / / Γ\H 2
where ρ is a map generically of degree 1 onto its image, which is a In particular if M σ ∈ Γ, the pullbacks generate the symmetric Hilbert modular functions for Γ; while if M σ ∈ Γ the pullbacks generate the full function field C Γ of Hilbert modular functions for Γ.

Proof. This is identical to the proof of Lemma 2.9 (see Lemma A.1 in Appendix A.1).

We have seen that by Lemma 2.9 we can take jk = φ * 1,ω j k , for k = 1, 2, 3, as invariants on a symmetric Hilbert modular surface. These functions are algebraically dependent. Similarly, we will apply Proposition 2.15 to the functions bk = φ * 1,ω b k and rk = φ * 1,ω r k for k = 1, 2, 3. Recall Equations ( 1), ( 2) and (3) for the definitions of these invariants. Proof. By Equation [START_REF] Charles | Cryptographic hash functions from expander graphs[END_REF], we have that φ [START_REF] Birkenhake | Humbert surfaces and the Kummer plane[END_REF]. Thus, the functions bk are modular for Γ [START_REF] Becker | On Gröbner bases under specialization[END_REF][START_REF] Birkenhake | Humbert surfaces and the Kummer plane[END_REF]. Moreover, if D ≡ 2, 3 mod 4, then these functions are also modular for Γ(2, 4)σ, as the matrix M σ of Equation ( 8) belongs to Γ [START_REF] Becker | On Gröbner bases under specialization[END_REF][START_REF] Birkenhake | Humbert surfaces and the Kummer plane[END_REF]. Similarly, 2) and the rk are modular for Γ(2) and also by Γ(2)σ when D ≡ 2, 3 mod 4. We conclude using Proposition 2.15 and the fact that the b i (resp. r i ) are generators for the field of Siegel modular functions invariants by Γ(2, 4) (resp. Γ( 2)). The pullbacks are indeed well defined because the denominators of these invariants divide χ 10 , so the locus of the denominators are components above the Humbert surface H 1 .

-1 1,ω (Γ(2, 4)) ∩ SL 2 (O K ⊕ ∂ K ) = Γ(2,
φ -1 1,ω (Γ(2)) ∩ SL 2 (O K ⊕ ∂ K ) = Γ(
We refer to Appendix A.5 for the index of the subgroups Γ(2) and Γ(2, 4) of Γ( 1), the number of Humbert surface components for Γ(2) and for Γ [START_REF] Becker | On Gröbner bases under specialization[END_REF][START_REF] Birkenhake | Humbert surfaces and the Kummer plane[END_REF] and equations of the Humbert component corresponding to the image of φ 1,ω for Γ(2, 4), Γ(2) and D ∈ {2, 3, 5}.

Equations for Hilbert surfaces 3.1 Generators of the field of Hilbert modular functions

Let Γ ⊂ SL 2 (O K ⊕ ∂ K ) be a subgroup of finite index. We note H Γ = Γ\H 2
1 the corresponding Hilbert modular surface, and H Γ,σ = ( Γ ∪ Γσ)\H 2 1 the corresponding symmetric Hilbert modular surface. We let G = Γ in the first case, and G = Γ ∪ Γσ in the second one. Recall that C G is defined in Definition 2.13.

Proposition 3.1. Let H G be a Hilbert surface as above. Then

C G = C(i 1 , i 2 , i 3 ) where i 1 and i 2 are symmetric Hilbert modular functions for SL 2 (O K ⊕∂ K ) and i 3 is algebraic over C(i 1 , i 2 ). Moreover i 3 is symmetric if and only if H G is symmetric.
Proof. Since H G is a surface, the field of Hilbert modular functions C G is of transcendence degree 2. By the primitive element theorem, C G is generated by two transcendental functions i 1 , i 2 (called primary invariants) and a third one i 3 algebraic over

C(i 1 , i 2 ) (called a secondary invariant). Since C G is algebraic over C SL 2 (O K ⊕∂ K )∪SL 2 (O K ⊕∂ K )σ , we can take i 1 , i 2 ∈ C SL 2 (O K ⊕∂ K )∪SL 2 (O K ⊕∂ K )σ . They are then symmetric, so H G is symmetric if and only if i 3 is symmetric.
Usually working with symmetric Hilbert modular surface yields invariants easier to compute. For instance while H Γ(1) is not often a rational surface according to Theorem A.16, from Elkies and Kumar [START_REF] Elkies | K3 surfaces and equations for Hilbert modular surfaces[END_REF] we have that H Γ(1),σ is a rational surface for every fundamental discriminant ∆ K < 100. Hence for these surfaces we need only two birational primary invariants to define the modular polynomials. The drawback of symmetric modular surfaces is that they can not be used for all the applications of isogenies as we will see in Example 4.17.

Note that by the general theory of Shimura varieties H G has a (birational) model defined over an algebraic number field F . In fact by Van Der Geer [71, Section X.4], the Hilbert surface can be defined over Q, and its connected components over an abelian extension of Q. In particular if the invariants i 1 , i 2 , i 3 come from this model defined over F , the equation

E(i 1 , i 2 , i 3 ) = 0 can be written as E = c k (i 1 , i 2 )i k 3 , where c k ∈ F (i 1 , i 2 )
. The lemma below show that if the Hilbert invariants have their Fourier coefficients are in a number field F then they are defined over F . In practice the invariants we use for computation (pullbacks of Igusa invariants, pullbacks of theta functions, Gundlach invariants) even have Fourier coefficients in Q. Lemma 3.2. Let i 1 , . . . , i n be Hilbert modular functions generating the Hilbert modular field C G , and let E be the ideal of equations between the i k and H E the corresponding birational model of H G . Then if the Fourier coefficients of each i k are in F , then the ideal E is generated by equations with coefficients in F , so H E has a model over F .

Proof. The proof uses a similar argument of Bröker and Lauter in [START_REF] Bröker | Modular polynomials for genus 2[END_REF]Theorem 5.2]. If we fix a monomial ordering, the generators of E are uniquely determined when they form a Gröbner basis. This Gröbner basis induces a set of linear relations on the Fourier coefficients of the i k from which its coefficients (as unknown) are the unique solution. But since the Fourier coefficients lie in F , this linear system is defined over F , so the solution is defined over F . Remark 3.3. The condition on the Fourier coefficients is a sufficient condition, but far from a necessary condition. In general the field of definition of the cusps will be larger than the field of definition of the Hilbert surface, so to know if the equations between the Hilbert functions i k as in Lemma 3.2 will lie in a polynomial ring over a proper subfield of F , one needs to look at the Galois action on the Fourier coefficients (which lie in F by hypothesis).

Fast evaluation of Hilbert modular functions

We will compute modular polynomials using an evaluation/interpolation approach. To be able to compute these polynomials in time quasi linear in their size, we need two properties for the invariants used:

• For the evaluation, given τ = (τ 1 , τ 2 ) ∈ G\H 2 1 we need to be able to compute the invariants (i 1 (τ ), i 2 (τ ), i 3 (τ )) ∈ C 3 in time quasi-linear in the required precision;

• Given the value of (i 1 (τ ), i 2 (τ ), i 3 (τ )) ∈ C 3 we need to be able to recover the matrix τ ∈ G\H 2 1 in time quasi-linear in the required precision. This allows us to have well-chosen evaluation points so that we can do the interpolation of trivariate rational fractions.

Recall that the notation Õ(N ) means O(N ) ignoring logarithmic factors. By quasi-linear complexity in N , we mean a complexity in Õ(N ).

Recall that by Theorem 2.16, the functions bk , for k = 1, 2, 3, generate C Γ(2,4) when D ≡ 1 mod 4 and C Γ(2,4)∪ Γ(2,4)σ when D ≡ 2, 3 mod 4. Let Γ ⊂ SL 2 (O K ⊕ ∂ K ) be a subgroup of finite index such that Γ ⊃ Γ [START_REF] Becker | On Gröbner bases under specialization[END_REF][START_REF] Birkenhake | Humbert surfaces and the Kummer plane[END_REF]. Let G = Γ or G = Γ∪ Γσ, and i 1 , i 2 , i 3 be as in Proposition 3.1 and such that F (H G ) = F (i 1 , i 2 , i 3 ), where F is the field of definition of H G .

Then i 1 and i 2 can be expressed as rational functions in the bk :

i k = R k ( b1 , b2 , b3
). This is also true for i 3 unless G = Γ and D ≡ 2, 3 mod 4 because then the bk are symmetric while i 3 is not, and thus can not be expressed as a rational function in the bk . But in this case, i 3 + σ(i 3 ) and i 3 σ(i 3 ) are symmetric, and can be expressed as rational functions R 3,1 and R 3,2 in the bk .

Similarly, by Lemma 2.9, the pullbacks of the Igusa invariants generate the function field of symmetric Hilbert modular functions. Thus there are rational fractions in i 1 , i 2 , i 3 giving the pullbacks of the Igusa invariants. Theorem 3.4. Let Γ, G and i 1 , i 2 , i 3 be as above. Assume that we know R 1 , R 2 , and R 3 or R 3,1 and R 3,2 according to the case, and the rational fractions giving the pullbacks of the Igusa invariants from i 1 , i 2 , i 3 . Assume that in the case G = Γ and D ≡ 2, 3 mod 4, we are able to evaluate i 3 at small precision (knowing the first Fourier coefficients for example). Then, not taking rounding errors into account, both the evaluation map

G\H 2 1 → C 3 , τ → (i 1 (τ ), i 2 (τ ), i 3 (τ )

) and the inversion map (the computation of a preimage of a point of the image of the evaluation map) can be computed in time quasi-linear in the precision.

Proof. We first do the symmetric case. By Enge and Thomé in [21, Theorem 5.8, Remark 5.9], given a Siegel matrix Ω ∈ H 2 , evaluating the b k (Ω) can be done in time quasi-linear in the precision (see also [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF][START_REF] Labrande | Explicit computation of the Abel-Jacobi map and its inverse[END_REF][START_REF] Labrande | Computing theta functions in quasi-linear time in genus two and above[END_REF] on this subject), if we neglige the loss of precision. Given a period matrix τ ∈ H 2 1 , one can use the map φ 1,ω from Section 2.3 to get Ω = φ 1,ω (τ ) ∈ H 2 , then we compute the values of bk (τ ) = b k (φ 1,ω (τ )) in time quasi-linear, and finally we deduce the values of i k (τ ) = R k ( b1 (τ ), b2 (τ ), b3 (τ )).

For the inverse, the (restriction of the) Igusa invariants j1 , j2 , j3 can also be expressed as rational functions in the invariants i 1 , i 2 , i 3 . From the values (i 1 (τ ), i 2 (τ ), i 3 (τ )), one can then compute ( j1 (τ ), j2 (τ ), j3 (τ )). Using [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF] (in particular [16, Théorème 9.3] and [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]Pages 200 and 223] of Dupont), we compute in time quasi-linear in N an approximation to precision N of a matrix Ω ∈ H 2 which is equivalent modulo Sp 4 (Z) to φ 1,ω (τ ). The matrix Ω lies in the Humbert surface of discriminant ∆ K , so it satisfies a singular relation. By Section 2.4 there is a constructive algorithm to find γ ∈ Sp 4 (Z) such that γ.Ω satisfies a normalized singular relation. By Section 2.3, γ.Ω is in the image of φ 1,ω , so one can compute τ = φ -1 1,ω (γ.Ω) ∈ H 2 1 . It then only remains to compute all classes of τ under the action of the finite group G\SL 2 (O K ⊕ ∂ K ) to find a τ in the same class as τ modulo G.

For the non symmetric case, recovering τ modulo G from the values of the invariants uses the same algorithm. The only difficulty is for the evaluation in the case D ≡ 2, 3 mod 4 because in this case the bk are symmetric while i 3 is not. However, since t = i 3 + σ(i 3 ) and n = i 3 σ(i 3 ) are symmetric, one can evaluate t(τ ) and n(τ ) in time quasi-linear using the techniques above for the symmetric case. Thus i 3 (τ ) is a root of X 2 -t(τ )X + n(τ ). The two roots can be computed in quasi-linear time in the precision, and choosing the correct one only requires an evaluation with small precision of i 3 .

Remark 3.5. Computing the rational functions is just a precomputation step and this computation does not affect the asymptotic complexity. They can be computed by linear algebra on the Fourier coefficients, or by linear algebra on the evaluation of these modular functions at several period matrices τ (where the evaluation uses the slow summation series given by the Fourier coefficients).

In practice, it is important to optimize the speed of the computation of the invariants i k as rational functions of the bk to be able to do concrete computations. Rather than using linear algebra, one can use an interpolation approach. This approach requires to be able to obtain, from the values (b 1 (Ω), b 2 (Ω), b 3 (Ω)), a period matrix Ω equivalent to Ω modulo Γ(2, 4), at some given precision. We know how to do it by [START_REF] Milio | A quasi-linear time algorithm for computing modular polynomials in dimension 2[END_REF][START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]. It also requires the equation P ( b1 , b2 , b3 ) = 0 of the Humbert component described by the bk . We refer to Section 3.3 for more details about the interpolation.

Likewise, to recover τ , rather than expressing the Igusa invariants j k in terms of the Hilbert invariants i k , one could simply use Newton's method to invert the equations i k = R k ( b1 , b2 , b3 ), P ( b1 , b2 , b3 ) = 0 for k = 1, 2, 3 to recover the values of the bk hence the matrix Ω, hence the matrix τ . Remark 3.6. In the first paragraph of the second page of [START_REF] Labrande | Computing theta functions in quasi-linear time in genus two and above[END_REF] about the computation of theta functions in genus g > 1 in quasi-linear time, Labrande and Thomé explain that they deliberately omit dealing with the precision losses because they expect that the loss of precision is not significant asymptotically. Indeed, this is what the analysis done by Labrande in genus g = 1 in [START_REF] Labrande | Computing Jacobi's theta in quasi-linear time[END_REF] concludes.

Interpolation by Hilbert modular functions

Let H G be a Hilbert surface defined over F of level G = Γ or G = Γ ∪ Γσ, and c a Hilbert modular function in F (H G ). We assume that the invariants i 1 , i 2 , i 3 are such that the map τ ∈ G\H 2 1 → (i 1 (τ ), i 2 (τ ), i 3 (τ )) can be inverted in time quasi-linear (see Theorem 3.4) at the working precision. We explain how to get a fast interpolation algorithm to express c as a rational function in i 1 , i 2 , i 3 . (Without the above property, one can still do linear algebra on the Fourier coefficients or the evaluations, which gives a slow interpolation algorithm). We assume that we have a quasi-linear time algorithm in the precision of the computations to evaluate c.

This interpolation algorithm will be used to compute the modular polynomials. Indeed, we will see that to compute them, we will have to interpolate bivariate or trivariate rational functions, that we can evaluate at any point (see the proof of Theorem 3.11 for the evaluation).

Case where H G is a rational surface. In this case, F (H G ) can be written as F (J 1 , J 2 ), using only two primary invariants. Thus c is a bivariate rational function in J 1 and J 2 . The interpolation of c is explained in [51, Section 2], where a complexity analysis is given. We explain it briefly. We compute c(J 1 , J 2 ) by computing first by interpolation c(J 1 , J 1 J 2 ) and then by substituting J 2 with J 2 /J 1 . Assume that we can write

c(J 1 , J 1 J 2 ) = A(J 1 , J 2 ) B(J 1 , J 2 ) = d A J 1 m=0 d A J 2 n=0 a m,n J m 1 J n 2 d B J 1 m=0 d B J 2 n=0 b m,n J m 1 J n 2 = d A J 1 m=0 a m (J 2 )J m 1 d B J 1 m=0 b m (J 2 )J m 1 ,
such that b 0 (J 2 ) = 1 (this is always possible because we consider c(J 1 , J 1 J 2 ) and not c(J 1 , J 2 ) except when b 0 (J 2 ) = 0 in which case we interpolate c(J 1 + t 0 , J 1 J 2 + t 1 ) instead, for some integers t 0 , t 1 ). Denote by d A T the total degree of the numerator of c(J 1 , J 2 ), i.e. d A J 1 , and similarly for B. Find points z m ∈ H 2 1 for m = 1, . . . , d A T + d B T + 2 such that (J 1 (z m ), J 2 (z m )) is of the form (u m , vu m ) for a fixed v ∈ C (this is where we need the inversion map). Interpolate to find the univariate rational fraction c(J 1 , vJ 1 ) and write this fraction such that the coefficient of degree 0 of the denominator is 1. Thus we know the values a m (v) and b m (v) (this would not have been possible if we had considered c(J 1 , J 2 ) directly). Compute in this way the fractions c(J 1 , v n J 1 ) for n = 1, . . . , max (d A J 2 , d B J 2 ) + 1. Interpolate by Lagrange's or Newton's method the univariate polynomials a m and b m to obtain c(J 1 , J 1 J 2 ).

In practice for the modular polynomials the coefficients of the bivariate rational fractions will be defined over Q. So the computations are done at precision N which has to be large enough so that we can recognize the coefficients of the bivariate rational fractions as algebraic numbers in Q using a continuous fraction algorithm. We do not usually know any bounds for the precision so that in practice we double the precision until we manage to find a sufficient precision to compute the modular polynomials. The complexity of the interpolation of a bivariate rational fraction is Õ(

d T d J 2 N ), where d T = max (d A T , d B T ) and d J 2 = max (d A J 2 , d B J 2 ).
General case. Here, we have three invariants i 1 , i 2 , i 3 where i 1 and i 2 are primary, and i 3 is a secondary invariant, so there is an equation

E(i 1 , i 2 , i 3 ) = 0 describing the surface H G .
Like in the previous case and as in [51, Section 2], we would like to work with points z j ∈ H 2 1 with the property that (i 1 (z j ), i 2 (z j ), i 3 (z j )) is of the form (u m , v n u m , w r u m ), where the subscripts m, n and r vary from 1 to the maximal degree the variables i 1 , i 2 and i 3 appear. But this is not possible because of the equation E that i 1 , i 2 , i 3 have to satisfy. Indeed, for fixed values (u m , v n u m ) such that u m = i 1 (z), v n u m = i 2 (z) for some z ∈ H 2 1 , the values that i 3 (z) can take are determined (moreover, they will not be of the form w r u m and the number of values will be less than the degree in i 3 ). A solution to this problem consists in remarking that

F (i 1 , i 2 , i 3 ) = F (i 1 , i 2 )[X 3 ]/(E(i 1 , i 2 , X 3 )).
Thus the modular function c can be written as c(i

1 , i 2 , i 3 ) = d-1 i=0 c i (i 1 , i 2 )i i 3
, where d is the degree in which the variable i 3 appears in E and c i ∈ F (i 1 , i 2 ).

Assume as in the previous case that we have an algorithm to evaluate c at points in H 2 1 . The interpolation is done as follows. For sufficiently many values u m and v n , compute the d

roots w r of E(u m , v n u m , X 3 ). For r = 1, . . . , d, find z r ∈ H 2 1 such that (i 1 (z r ), i 2 (z r ), i 3 (z r )) = (u m , v n u m , w r ) and evaluate c at z r to obtain c(z r ) = d-1 i=0 c i (u m , v n u m )w i r .
Since w r = i 3 (z r ), we first interpolate c as a univariate polynomial in i 3 by interpolating on the d values w r to recover the d coefficients c i (u m , v n u m ). This allows us to evaluate the bivariate rational functions c i (i 1 , i 2 ) at a given point z r . It remains to do the interpolation of the coefficients c i to recover them as rational functions in i 1 , i 2 as was done in the previous case, where H G was a rational surface. We summarize this discussion by the following theorem.

Theorem 3.7. Let G be a subgroup of finite index in SL 2 (O K ⊕ ∂ K ) ∪ SL 2 (O K ⊕ ∂ K )σ.
Let N be a precision for the computations. Let i 1 , i 2 , i 3 generating C G be such that a preimage at precision N of the image of a point by the evaluation map τ ∈ G\H 2 1 → (i 1 (τ ), i 2 (τ ), i 3 (τ )) can be found in time quasi-linear in N . Let E(i 1 , i 2 , i 3 ) be the equation describing the Hilbert surface H G , and let d be the degree

deg i 3 (E) of i 3 in E.
Let c be a Hilbert modular function in C G . We represent it in the form c(i

1 , i 2 , i 3 ) = d-1 k=0 c k (i 1 , i 2 )i k 3 , where c k (i 1 , i 2 ) ∈ C(i 1 , i 2 )
are bivariate rational functions. We let d T be the maximal total degree of all the coefficients c k (where the degree of a rational function is the maximal of the degree of its numerator and denominator), and d i 2 the maximal degree in i 2 of the coefficients c k .

Assume that we have a quasi-linear time algorithm in the precision N taking τ ∈ H 2 1 and N as inputs and returning the evaluation of c at τ at precision N . Assume we know d T and d i 2 and N is large enough so that the interpolations give the correct coefficients c k (i 1 , i 2 ) up to the precision. Then, not taking rounding errors into account, the coefficients c k can be computed in precision N in time Õ(dd

T d i 2 N ).
Assume furthermore that the c k lie in F (i 1 , i 2 ), where F is a number field. Then the coefficients c k can be recovered exactly in time Õ(dd T d i 2 N ) if N is large enough so that the coefficients can be recovered exactly using the Lenstra-Lenstra-Lovász (LLL) algorithm.

Proof. Indeed the evaluation of c will be executed O(dd T d i 2 ) times and we will interpolate O(d) bivariate rational functions and do O(d T d i 2 ) interpolations of a univariate polynomial. The complexity is then [START_REF] Elkies | K3 surfaces and equations for Hilbert modular surfaces[END_REF] O(dd

T d i 2 ) + O(d) Õ(d T d i 2 N ) + O(d T d i 2 ) Õ(dN ) ⊂ Õ(dd T d i 2 N ). Given a coefficient c k ∈ C computed at precision N , if c k lie in F (i 1 , i 2 )
, where F ⊂ C is a number field, then one can use the LLL algorithm [START_REF] Lenstra | Factoring polynomials with rational coefficients[END_REF] to recover c k ∈ F . More precisely, if h is the maximal (logarithmic) height of the coefficients c k seen in C, then N can be bounded by a term in O(h). Using fast versions of LLL this reconstruction step can be done in time Õ(N ) (according to Novocin, Stehlé and Villard [START_REF] Novocin | An LLL-reduction algorithm with quasi-linear time complexity[END_REF]). Remark 3.8. In practice, the computation of d T , d i 2 and N is a precomputation step. We fix values for them and try to do the computations. If it fails, we double each values and try again. The exact values of d T and d i 2 can be found in this way.

More generally a similar technique could be used if we had several secondary invariants i 3 , i 4 , . . . i . There is no unique expression of c in terms of the i k due to the equations between the invariants i k . But for the interpolation to work we need to interpolate the same rational function expression. A solution is to fix a monomial ordering, since this defines a unique rational function expressing c modulo the corresponding Gröbner basis. As long as the partial evaluation of the Gröbner basis corresponds to the Gröbner basis of the partial evaluation of the equation (see [START_REF] Becker | On Gröbner bases under specialization[END_REF][START_REF] Kalkbrener | On the stability of Gröbner bases under specializations[END_REF]), the interpolation step will interpolate the correct expression of the rational function.

Equations for covers of Hilbert surfaces

Let G 2 ⊂ G 1 ⊂ SL 2 (O K ⊕ ∂ K ) ∪ SL 2 (O K ⊕ ∂ K )σ be congruence subgroups. Then H G 2 → H G 1 is a covering. Let i 1 , i 2 , i 3 be Hilbert modular functions such that C G 1 = C(i 1 , i 2 , i 3 ) and j 1 , j 2 , j 3 be Hilbert modular functions such that C G 2 = C(j 1 , j 2 , j 3 ).
To describe the cover H G 2 → H G 1 we need to give the full set of relations between i 1 , i 2 , i 3 , j 1 , j 2 , j 3 . To be more precise, as always in this text we work up to birational equivalence, and i 1 , i 2 , i 3 only give an embedding of an open subset of H G 1 , and similarly for j 1 , j 2 , j 3 . To describe the full cover we would potentially need to give the relations between more modular functions invariant by G 1 (respectively G 2 ), but these relations can be computed with the same tool as for the relations between i 1 , i 2 , i 3 , j 1 , j 2 , j 3 .

Let i 1 , i 2 , i 3 be generators of the Hilbert modular field C G 1 such that the evaluation and its inverse can be computed in time quasi-linear (see for instance Theorem 3.4).

Let j be a generator of the field extension C G 2 /C G 1 . Such a generator always exists by the primitive element theorem. The cover H G 2 → H G 1 is then (up to birationality) uniquely described by

• the minimal polynomial Φ j ∈ C G 1 [X] of j over C G 1 ; • and polynomials Q k ∈ C G 1 [X] such that j k = Q k (j).
In practice it is more convenient to use the polynomial Ψ k ∈ C G 1 [X] of minimal degree such that j k Φ j (j) = Ψ k (j) than the polynomials Q k . The polynomial Ψ k is called the Hecke representation of j k and is more convenient for computations than Q k because it has smaller coefficients [28, Section 3]. Lemma 3.9. The polynomial

Ψ k is Ψ k (X) = γ∈G 2 \G 1 j γ k Φ j (X)/(X -j γ ).
Proof. Let M/K be a finite Galois extension of Galois group G, and for f ∈ M and γ ∈ G denote by f γ the action γ.f of

γ on f . Let G 2 ⊂ G 1 ⊂ G and let K 2 = M G 2 , K 1 = M G 1 . Let j be a generator of K 2 /K 1 ; then its minimal polynomial is Φ(X) = γ∈G 2 \G 1 (X -j γ ). Let J ∈ K 2 , then there exists Q ∈ K 1 [X] of minimal degree such that J = Q(j). Since J γ = Q(j γ ),
we can use Lagrange interpolation to find Q. Indeed, evaluating δ∈G 2 \G 1 J δ δ =δ (X -j δ )/(j δ -j δ ) at j γ gives J γ . Now, this expression is equal to δ∈G 2 \G 1 J δ δ =δ (X -j δ )/Φ (j δ ) = δ∈G 2 \G 1 J δ Φ(X)/((X -j δ )Φ (j δ )) and we deduce that taking Ψ(X) = δ∈G 2 \G 1 J δ Φ(X)/(X -j δ ), we have the property J γ Φ (j γ ) = Ψ(j γ ).

We apply this to the extension

C Γ(n) /C SL 2 (O K ⊕∂ K )∪SL 2 (O K ⊕∂ K )σ where Γ(n) is a congruence subgroup included in G 2 . Indeed this is a Galois extension of Galois group Γ(n)\(SL 2 (O K ⊕ ∂ K ) ∪ SL 2 (O K ⊕ ∂ K )σ
), and we apply the result above to G 1 = Γ(n)\G 1 and G 2 = Γ(n)\G 2 with the notations of the Lemma. Theorem 3.10. Let i 1 , i 2 , i 3 , j, j 1 , j 2 , j 3 be Hilbert modular functions such that

C G 1 = C(i 1 , i 2 , i 3 ) and C G 2 = C G 1 (j) = C(j 1 , j 2 , j 3 ).
Let E be the equation E(i 1 , i 2 , i 3 ) = 0 of the surface birational to H G 1 described by i 1 , i 2 , i 3 , and let d be the degree deg i 3 (E) of i 3 in E. Assume that all these modular functions have Fourier coefficients in an algebraic number field

F ⊂ C. Let Φ(X, i 1 , i 2 , i 3 ) = γ∈G 2 \G 1 (X -j γ ) = X L + L-1 m=0 c m (i 1 , i 2 , i 3 )X m be the minimal polynomial of j over C G 1 , where L = #G 2 \G 1 . Let Ψ k ∈ C G 1 [X]
be the polynomial given in Lemma 3.9 for j k . A birational model of the cover H G 2 → H G 1 is described by the equations

(19) Φ(j) = 0, j 1 Φ (j) = Ψ 1 (j), j 2 Φ (j) = Ψ 2 (j), j 3 Φ (j) = Ψ 3 (j).
The coefficients c m of the polynomial Φ can be written as c m = d-1 n=0 c mn (i 1 , i 2 )i n 3 , and similarly for Ψ k . We have c mn ∈ F (i 1 , i 2 ).

Proof. The birational model comes from the definition of the polynomials Ψ k . As i 1 , i 2 , i 3 , j have Fourier coefficients in F , the same argument as in Lemma 3.2 or [START_REF] Bröker | Modular polynomials for genus 2[END_REF]Theorem 5.2] shows that c m ∈ F (i 1 , i 2 , i 3 ). Moreover by the same argument the equation E is defined over F , and using that

F (i 1 , i 2 , i 3 ) = F (i 1 , i 2 )[X 3 ]/(E(i 1 , i 2 , X 3 )), we can also write c m = d-1 n=0 c mn (i 1 , i 2 )i n
3 and c mn ∈ F (i 1 , i 2 ).

Theorem 3.11. With the conditions of Theorem 3.10, assume moreover that

• there is a fast evaluation algorithm in the precision to evaluate j, j 1 , j 2 , j 3 .

• a preimage of the evaluation of (i 1 , i 2 , i 3 ) at some point can be computed in time quasilinear in the precision;

We let d T be the maximal total degree of all the coefficients c mn (where the degree of a rational function is the maximal of the degree of its numerator and denominator), and d i 2 be the degree in i 2 of the coefficients c mn . Let N be the maximal height (over F ) of the coefficients of each rational function c mn ∈ F (i 1 , i 2 ). Then Φ and the Ψ k can be computed in time Õ(dd T d i 2 LN ).

Proof. We compute Φ by computing the coefficients c m ∈ F (i 1 , i 2 , i 3 ) by evaluation/interpolation. The evaluation of c m at some τ ∈ H 2 1 is done as follows.

• Compute each value j(γ.τ ) in precision N . It can be done with a complexity in L Õ(N ) time;

• Evaluate Φ(j(τ )) = γ∈G 2 \G 1 (X -j(γ.τ )) using a subproduct tree (see [START_REF] Gathen | Modern Computer Algebra[END_REF]Section 10.1]). The polynomial Φ(j(τ )) can be obtained in Õ(LN ) time.

• Separating the coefficients according to powers of X gives the values c m (i 1 (τ ), i 2 (τ ), i 3 (τ )).

By Section 3.3 and Theorem 3.7, to recover Φ, the evaluation step will be executed O(dd T d i 2 ) times and we will interpolate O(dL) bivariate rational fractions and do O(Ld T d i 2 ) interpolation of a univariate polynomial. Recall that given the coefficient c mn ∈ C computed at precision O(N ), using the LLL algorithm to recover c mn ∈ F can be done in time Õ(N ) ( [START_REF] Novocin | An LLL-reduction algorithm with quasi-linear time complexity[END_REF]). The final complexity is then [START_REF] Enge | Class invariants by the CRT method[END_REF] O(dd

T d i 2 ) Õ(LN ) + O(dL) Õ(d T d i 2 N ) + O(Ld T d i 2 ) Õ(dN ) ⊂ Õ(dd T d i 2 LN ).
The same algorithm works for the Ψ k , where at the evaluation step, Ψ k (j(τ )) is computed via a double subproduct tree on Ψ k and Φ.

Modular polynomials 4.1 Isogenies preserving real multiplication

We want now to define modular polynomials, which parametrize isogenies between principally polarized abelian surfaces with real multiplication by O K . To achieve this, we need the following subgroups of SL

2 (O K ⊕ ∂ K ) Γ0 (β) = a b/ √ ∆ K √ ∆ K c d ∈ SL 2 (O K ⊕ ∂ K ) : b ∈ βO K , Γ(β) = a b/ √ ∆ K √ ∆ K c d ∈ SL 2 (O K ⊕ ∂ K ) : a -1, d -1, b, c ∈ βO K ( 21
)
where β ∈ O K is a totally positive prime number. By abuse of notation, we use the same notation for their preimage by the isomorphism φ ± of Equation ( 4). Then the Hilbert cover Γ0 (β)\H 2 1 parametrizes pairs (A, B) of (principally polarised) abelian surfaces with real multiplication by O K together with a β-isogeny A → B, or equivalently pairs (A, K) where A has real multiplication by O K and K ⊂ A[β] is a kernel stable by O K and maximally isotropic for the β-Weil pairing.

We now explain those terms and give more details on isogenies preserving the real multiplication. Let (A, θ A ) be a principally polarized abelian surface over a field k, with real multiplication given by µ : O K → End(A). Let f : A → B be an isogeny with kernel V . Then it is easy to see that B has real multiplication by O K compatible with f if and only if V is stable under the action of µ(O K ).

It remains to see when B admits a principal polarization. If θ B is such a principal polarization, then θ = f * θ B is a polarization on A. Denote by End s (A) the subring of endomorphisms of A commuting with the Rosati involution induced by θ A . By [3, Proposition 5.2.1 and Theorem 5.2.4], the group isomorphism between the Néron-Severi group of A and End s (A) induces a bijection between the polarizations of some degree L of A and the set End s,++ (A) of totally positive symmetric endomorphisms of A with analytic norm L. When End s (A) = O K (which is the case generically for an element of the Hilbert surface), then θ comes from a totally positive element β ∈ O ++ K . Furthermore it is easy to check that V is a totally isotropic subgroup for the Weil pairing e β on A[β] = {P ∈ A(k) | βP = 0}. Looking at degrees, we also get that #V = N K/Q (β).

Conversely, let β ∈ O ++

K and note θ β the polarization on A induced from θ A by β (the polarization which arises from the composition of the polarization isogeny of θ A with β), and V ⊂ A[β] a maximal isotropic subgroup for the Weil pairing e β . Then by descent theory, θ β descends to a polarization θ B on B = A/V , and since V is maximal, θ B is principal. To emphasize the role of β, we call the isogeny f induced by V a β-isogeny.

Remark 4.1. The notation θ β comes from the fact that if θ is induced by a symmetric line bundle L and β = ∈ N, then θ is induced by the symmetric line bundle L .

For more details we refer to [START_REF] Robert | Computing cyclic isogenies using real multiplication[END_REF][START_REF] Dudeanu | Computational Aspects of Jacobians of Hyperelliptic Curves[END_REF][START_REF] Dudeanu | Cyclic Isogenies for Abelian Varieties with Real Multiplication[END_REF] and in particular in [15, Theorem 1.1] of Dudeanu, Jetchev, Robert, and Vuille. We are mainly interested with cyclic isogenies of prime degree , these are induced by β of norm . We sum up the discussion above by the following Proposition 4.2. Let (A, θ) be a principally polarized abelian surface lying on the Humbert surface H ∆ K . Then there exists cyclic isogenies of degree (possibly defined over an extension of the field of definition of (A, θ)) from A to a principally polarized abelian surface if there exists a totally positive element β ∈ O ++ K of norm . And conversely if the abelian surfaces lying on the Humbert surface admit cyclic isogenies of degree generically, then there exists such a β.

We will compute modular polynomials when β = ∈ Z is a prime number, and when β is a totally positive element of O K of norm . We first give the degree of the cover Γ0 (β 2 1 , since this will determine the degrees of the modular polynomials. When β = , the Weil pairing is the usual pairing e on A[ ], and the correspondingisogenies come from isotropic kernels of degree 2 . Over the splitting field of A[ ] over the field of definition of A, it is easy to see that there are 3 + 2 + + 1 such isogenies (this is the size of the quotient Γ 0 ( )\Sp 4 (Z)). The computation of the corresponding modular polynomials is described in [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF][START_REF] Bröker | Modular polynomials for genus 2[END_REF][START_REF] Milio | A quasi-linear time algorithm for computing modular polynomials in dimension 2[END_REF]. On the Hilbert side of things, not all such isogenies stay on the Humbert surface. Indeed this is the case if and only if the kernel is stable under the real multiplication by O K . Since the Weil pairing is compatible with endomorphisms, A[ ], as a O K -module, is given by a symplectic basis e 1 , e 2 . To such a basis one can associate the subgroup V = O K e 1 which is maximal isotropic for the Weil pairing and stable under the real multiplication by O K . All other such kernels are obtained in a similar way via the action of Γ0 ( )\SL 2 (O K ⊕ ∂ K ) on the symplectic basis (e 1 , e 2 ).

)\H 2 1 → SL 2 (O K ⊕ ∂ K )\H

Proposition 4.3. The number of -isogenies whose kernel is stable under the real multiplication is

• 2 + 1, if is inert in O K ; • ( + 1) 2 , if is split in O K ; • 2 + , if is ramified in O K .
Proof. Using φ ± , we have to find the cardinality of Γ0 ( )\SL 2 (O K ). We compute it as the quotient of # Γ( )\SL 2 (O K ) by # Γ( )\ Γ0 ( ). We use the group isomorphism ψ : Γ( )\SL 2 (O K ) SL 2 ( O K \O K ) and we compute the cardinality of Γ( )\ Γ0 ( ) looking at its image by ψ. Moreover we use that

O K \O K is isomorphic to F 2 if is split, to F 2 if is inert and to F [X]/(X 2 ) if is ramified. If is inert, #SL 2 (F 2 ) = 2 ( 4 -1) and # Γ( )\ Γ0 ( ) = 2 ( 2 -1). If is split, #SL 2 (F 2 ) = 2 ( 2 -1) 2 and # Γ( )\ Γ0 ( ) = 2 ( -1) 2 . If is ramified, #SL 2 (F [X]/(X 2 )) = 6 -4 and # Γ( )\ Γ0 ( ) = 2 ( 2 -).
Next suppose that we have β ∈ O ++ K totally positive of norm . In this case either is ramified in O K and there is only one kind of cyclic isogenies of degree , the β-isogenies, or splits as = ββ and

A[ ] = A[β] ⊕ A[β]
and there are two kind of cyclic isogenies: the β-isogenies and the β-isogenies. Proof. We have seen that β-isogenies correspond to maximally isotropic kernels of size in

A[β]. Since A[β] is of size 2 ,
such kernels are exactly the cyclic kernels of size . Since βO K \O K F , the elements of O K act by scalar multiplication on A[β] so they stabilize all the cyclic subgroups. And indeed since Γ(β)\SL 2 (O K ) SL 2 (F ) it is easy to check that Γ0 (β)\SL 2 (O K ) is of size + 1 and a set of representatives is given by the matrices ( 1 x 0 1 ) for x ∈ {0, . . . , -1} and 0 1 -1 0 . Indeed since βO K \O K Z\Z, we have that Γ0 (β)\SL 2 (O K ) Γ 0 ( )\SL 2 (Z) whose set of representatives is well known.

Furthermore it is easy to see that the composition of a β-isogeny and a β-isogeny is an -isogeny (preserving real multiplication). Conversely a counting argument shows that any -isogeny preserving real multiplication splits as a β-isogeny and a β-isogeny (which may be defined over an extension of greater degree). So when splits as ββ, β ∈ O ++ K , we only need to compute β and β Hilbert modular polynomials. Then -1 (V ) = V is the kernel of an 2 β-isogeny, and the isogenous variety A/ -1 (V ) is isomorphic to A/V (as principally polarized abelian varieties).

Proof. Let be any endomorphism of A and θ a principal polarization. Then the pullback * θ is induced by the symmetric endomorphism ˆ where • denote the Rosati involution. More generally, if β is totally positive, then * θ β = θ ˆ β .

In particular, if f : A → B is a β-isogeny, then f • is a ˆ β -isogeny. It suffices to apply this to ∈ O × K (so that ˆ = ) and f : A → B the isogeny with kernel V . If θ B is the principal polarization induced by the descent of θ β , then the descent of (A,

θ β ) induced by -1 (V ) is (B, θ -2 B ) and -1 : B = A/V → A/ -1 (V )

induces the required isomorphism of principally polarized abelian varieties.

From this Lemma we deduce that the 2 β-modular polynomial will be the same as the β-modular polynomial. Remark 4.6. For simplicity of the exposition we work with the maximal real order O K . However everything outlined above still work with a real order O that is only locally maximal at . It is also straightforward to adapt the results of Section 3 for computing invariants on the Hilbert surfaces of abelian surfaces with real multiplication by O, and also to compute modular polynomials for such an O.

Modular functions for β-isogenies

We let β ∈ O ++ K be a prime element of norm L. So L = if ∈ Z is a prime number which splits or ramifies in O K , and L = 2 if stays inert. Generalising a bit the preceding section, we let Γ ⊂ SL 2 (O K ⊕ ∂ K ) be a congruence subgroup containing Γ(n). We want to apply the results of Section 3.4 to the extension C Γ0 (β)∩ Γ/C Γ. This cover corresponds to pairs of (principally polarised) abelian surfaces (A, B) with real multiplication together with a beta isogeny A → B as in Section 4.1, which is furthermore compatible with the level structure induced by Γ on A and B.

For simplicity, we now assume that n is prime to L. We first want to give an explicit set of representatives of Γ0 (β) ∩ Γ\ Γ. Recall that there in an isomorphism 4), so that by looking at the preimage by φ ± we can assume here that Γ ⊂ SL 2 (O K ) (this is more convenient to study the quotient). Recall that in this model, Γ0

φ ± : SL 2 (O K ) → SL 2 (O K ⊕ ∂ K ) defined in Equation (
(β) = { a b c d ∈ SL 2 (O K ) : β|b}.
Lemma 4.7. Let N be an integer. Then the map SL

2 (O K ) → SL 2 (O K /N O K ) is surjective.
Proof. This is an application of Strong approximation theory. In this case an elementary proof is also given in Bourbaki, Algebre Commutative, VII, §2, n.4: since SL n (O K /N O K ) is a product of local rings, it is generated by elementary matrices, so it suffices to lift these matrices. Example 4.9. We describe in more details the important case Γ = SL 2 (O K ). The group Γ is generated by the three matrices M 1 = 0 -1 1 0 , M 2 = ( 1 1 0 1 ) and M 3 = ( 1 ω 0 1 ) (see Section 2.2 for the definition of ω). Note that M 2 1 0 -1 1 M 2 = -M 1 so that it will be sometimes more convenient to use the matrix 1 0 -1 1 instead of M 1 . By Lemma 4.8, the subgroup Γ0 (β) of Γ is of index L + 1 and the set of matrices

C β = M 1 , M i
2 , i ∈ {0, . . . , L -1} is a set of representatives of the classes of Γ0 (β)\ Γ. We can give a different proof using the matrices M 1 , M 2 and M 3 : the L+1 matrices of

C β are clearly in different classes of the quotient Γ0 (β)\ Γ. Remark that t M 2 = M 1 M -1 2 M -1 1 ∈ Γ0 (β) and t M 3 = M 1 M -1 3 M -1 1 ∈ Γ0 (β) and that Γ is generated by M 1 , t M 2 and t M 3 . For all i ∈ {0, . . . , L}, t M 2 M i 2 and t M 3 M i 2 are in the class of M i 2 while t M 2 M 1 and t M 3 M 1 are in the class of M 1 . Moreover, M 1 M i 2 is in the class of M 1 and M 1 M 1 = -I 2 which
shows that there can not be more than the L + 1 classes that we already know.

Example 4.10. Another important example is the case Γ = Γ [START_REF] Becker | On Gröbner bases under specialization[END_REF][START_REF] Birkenhake | Humbert surfaces and the Kummer plane[END_REF]. By the above Lemma, the subgroup Γ(2, 4) ∩ Γ0 (β) of Γ(2, 4) is of index L + 1.

If γ ∈ Γ0 (β)\SL 2 (O K ) then there exists an element γ ∈ Γ0 (β) such that γ γ ∈ Γ(2, 4). For applications it is useful to have a constructive definition of γ .

We look at γ such that γ γ ≡ 0 mod 4, namely such that γ ≡ γ -1 mod 4, and such that γ ≡ ( * 0 * * ) mod . By the Chinese Remainder Theorem, these conditions modulo 4 and gives a matrix γ which must satisfy conditions modulo 4 and by Lemma 4.7, γ can be lifted to a matrix in Γ. Now we go back to the usual model Γ ⊂ SL 2 (O K ⊕ ∂ K ). Let G be either Γ or Γ ∪ Γσ. We have G ∩ Γ0 (β) = Γ ∩ Γ0 (β). In the case that σ ∈ G, we recall that by Lemma 2.12 we can assume that the group Γ is stable under the real conjugation.

Let i 1 , i 2 , i 3 be generators of the Hilbert modular field C G . Recall that C G is defined in Definition 2.13. (Later we will assume that they are chosen such that the evaluation and its inverse can be computed in time quasi-linear, like in Theorem 3.4.)

Let j be a generator of the field extension C Γ∩ Γ0 (β) /C G . Such a generator always exists by the primitive element theorem. In fact it is easy to find such a generator: Proposition 4.11. Let i 1 , i 2 , i 3 be generators of the Hilbert modular field C Γ. Let j be a Hilbert modular function invariant by Γ ∩ Γ0 (β) but not by Γ.

Then C Γ∩ Γ0 (β) = C(i 1 , i 2 , i 3 , j).
Let G = Γ ∪ Γσ. Let i 1 , i 2 , i 3 be generators of the symmetric Hilbert modular field C G . Let j be a Hilbert modular function invariant by Γ ∩ Γ0 (β) but not by Γ. Then if j is symmetric,

C ( Γ∩ Γ0 (β)) σ = C(i 1 , i 2 , i 3 , j), otherwise C Γ∩ Γ0 (β) = C(i 1 , i 2 , i 3 , j).
Proof. Since the symmetric case is easily deduced from the non symmetric one, we only do the case G = Γ. We have seen in the proof of Lemma 3.9 that the extension

C Γ(Ln) /C Γ(1)∪ Γ(1)σ is Galois of Galois group Γ(Ln)\( Γ(1) ∪ Γ(1)σ). Let K 1 = C Γ(j) = C(i 1 , i 2 , i 3 , j) and K 2 = C Γ(Ln)\( Γ∩ Γ0 (β)) Γ(Ln)
= C Γ∩ Γ0 (β) . Then K 1 ⊂ K 2 and we want to prove the equality. By Galois theory, the subfields between K 1 and K 2 correspond to subgroups of Γ containing Γ ∩ Γ0 (β). If we show that the group Γ ∩ Γ0 (β) is maximal in Γ, then we would deduce that

K 1 = C Γ or K 1 = K 2 .
By assumption, only the last possibility can be true. Looking at the image of Γ0 (β) ∩ Γ by the isomorphism Γ(β) ∩ Γ\ Γ Γ(β)\ Γ(β) Γ = Γ(β)\ Γ(1), we can see that is is King [41,Theorem 4.1], the set of triangular matrices of SL 2 (Z/ Z) is maximal and thus π( Γ0 (β)) is maximal in SL 2 (Z/ Z) 2 . As π is surjective, we deduce that Γ0 (β) is maximal in Γ.

enough to prove that Γ0 (β) is maximal in Γ(1). Let π : Γ SL 2 (O K / O K ). If β is of norm L = prime (so that ) is split, then SL 2 (O K / O K ) SL 2 (Z/ Z) 2 and π( Γ0 (β)) = {( * 0 * * ) × ( * * * * )}. By
If β = is inert, then the image of π( Γ(β)) is given by triangular matrices of SL 2 (F 2 ) so it is also maximal.

If is ramified, then SL

2 (O K / O K ) SL 2 ((Z/ Z)[X]/(X 2
)) and π( Γ0 (β)) is the set of matrices of the form ( * xX * * ) for any x ∈ Z/ Z. Let G be a group which contains strictly π( Γ0 (β)). Then there exists some matrix

A B C D ∈ G, with B(0) = 0. If A is invertible (namely A(0) = 0) then 1 0 -AC 1 A -1 0 0 A A B C D = 1 A -1 B 0 1 ∈ G and (A -1 B)(0) = 0 so that A -1 B = x 0 + x 1 X with x 0 = 0. Finally we have 1 x 0 +x 1 X 0 1 1 -x 1 X 0 1 = 1 x 0 0 1
from which we deduce that ( 1 1 0 1 ) ∈ G. As this last matrix and the matrices 1 0 -1 1 and 1 X 0 1 are all in G and are generators for SL 2 (O K / O K ), we deduce that G is π( Γ), that π( Γ0 (β)) is maximal and thus by surjectivity that Γ0 (β) is also maximal. If A is not invertible but D is, the proof proceeds similarly. Otherwise, if both A and D are not invertible, then B and C are. Moreover,

Computing modular polynomials

Using the notations of Section 4.2, we want to compute modular polynomials classifying all β-isogenies from an abelian surface with real multiplication by O K . First consider the case Γ = Γ(1). Recall from Section 4.1 that geometrically, a point in H Γ0 (β) corresponds to a triple (A, θ, V ) with a principally polarized abelian surface (A, θ) and V the kernel of a β-isogeny (equivalently V is maximally isotropic for the e β Weil pairing on A[β]). We note π : (A, θ, V ) → (A, θ) × (A/V, θ ) where θ is the polarization induced on A/V by θ β . This defines an algebraic map (a modular correspondence) H Γ0 (β) → H Γ(1) × H Γ(1) . The β-modular polynomials describe the algebraic relations giving the image of this map.

Concretely, if i 1 , i 2 , i 3 generate C( Γ(1)), the β-modular polynomials for the invariants i k describe the locus of the modular points ((i

1 (z), i 2 (z), i 3 (z)), (i 1 (z/β), i 2 (z/β), i 3 (z/β)) for z ∈ H 2 1 .
In particular the β-modular polynomials classify the β-isogenies.

Indeed if z ∈ Γ(1)\H 2 1 , the β-isogenous varieties are 1 β γ • z for γ ∈ Γ0 (β)\ Γ(1)
. Furthermore since σ Γ0 (β)σ = Γ0 (β), the β-isogenous varieties are given by 1 β γ • z, for γ ∈ Γ0 (β)\ Γ(1). Furthermore, by Lemma 4.13, the modular function i 1 (z/β) is invariant by Γ0 (β), so by Proposition 4.11 is a generator of C(H Γ0 (β) )/C(H Γ( 1) ). This implies that the modular correspondance 1) is is an embedding (at least birationally), so that the β-modular polynomials do give equations for H Γ0 (β) .

H Γ0 (β) → H Γ(1) × H Γ(
More generally, for a group Γ containing a congruence subgroup Γ(n) with n prime to L, we would like to define β-modular polynomials describing the image of a map (a modular correspondence)

H Γ∩ Γ(β) → H Γ × H Γ. A point in H Γ∩ Γ(β) corresponds to a triple (A, θ, V )
as above together with an extra level structure G defined by Γ. For the modular correspondence to give an embedding of H Γ∩ Γ(β) we need for G to induce a unique extra level structure G on (A/V, θ ).

(z) = i(z/β). Recall that, for γ ∈ Γ(1) ∪ Γ(1)σ, i γ (z) = i(γ • z) and define i γ β (z) = i( 1 β γ • z). Then for γ ∈ Γ0 (β), i γ β (z) = i( 1 β γ • z) = i(γ β • ( 1 β z)) = i γ β ( 1 β z) i σ β (z) = i( 1 β σz) = i σ ( 1 β z)) i γσ β (z) = i( 1 β γσ • z) = i(σγ β • ( 1 β z)) = i σ γ β ( 1 β z) Corollary 4.14. Let i be a Hilbert modular function for Γ ⊂ SL 2 (O K ⊕ ∂ K ). Let Γβ = {γ ∈ SL 2 (O K ⊕ ∂ K ) | γ β ∈ Γ} ⊂ Γ0 (β). Then i β is modular for Γβ . Furthermore if i is symmetric and β = β, then i β is symmetric. Assume that for every γ ∈ Γ ∩ Γ0 (β), γ β ∈ Γ, so (22) Γ ∩ Γ0 (β) = Γ ∩ Γβ .
Then if i is a Hilbert modular function for Γ, then i β is a Hilbert modular function for Γ∩ Γ0 (β).

If Γ satisfy Equation ( 22) (such is the case when Γ = Γ(n) is a congruence subgroup), one can then define the modular correspondence as

H Γ∩ Γ0 (β) → H Γ × H Γ, z → ((i 1 (z), i 2 (z), i 3 (z)), (i 1 (z/β), i 2 (z/β), i 3 (z/β))) for z ∈ H 2 1 and i 1 , i 2 , i 3 generating C Γ.
Theorem 4.15. Non symmetric case: let Γ be a congruence subgroup such that Γ(2, 4)

⊂ Γ ⊂ SL 2 (O K ⊕ ∂ K ). Let β ∈ O ++
K be a prime of norm L, and assume that for every γ ∈ Γ ∩ Γ0 (β), γ β ∈ Γ.

Let C β be a set of representatives of Γ ∩ Γ0 (β)\ Γ.

Let i 1 , i 2 , i 3 modular functions generating C Γ and with Fourier coefficients in a number field F .

Define the modular polynomials:

(23) Φ β (X, i 1 , i 2 , i 3 ) = γ∈C β (X -i γ 1,β ), and Ψ k,β (X, i 1 , i 2 , i 3 ) = γ∈C β i γ k,β Φ β (X, i 1 , i 2 , i 3 ) X -i γ 1,β for k = 2, 3. They lie in F (i 1 , i 2 , i 3 )[X].
Then after a precomputation step described in Theorem 3.4 (which does not depend on β, only on i 1 , i 2 , i 3 ), and under the heuristics of [START_REF] Milio | A quasi-linear time algorithm for computing modular polynomials in dimension 2[END_REF]Theorem 34], the modular polynomials can be computed in quasi-linear time in their size.

Symmetric case:

Let G = Γ ∪ Γσ. If β = β we let C β be a set of representatives of (( Γ ∩ Γ0 (β)) ∪ ( Γ ∩ Γ0 (β))σ)\G Γ ∩ Γ0 (β)\ Γ, otherwise we let C β be a set of representatives of (G ∩ Γ0 (β))\G Γ ∩ Γ0 (β)\ Γ ∪ Γ ∩ Γ0 (β)\ Γ σ.
Then the same definition as in Equation [START_REF] Freitag | Hilbert modular forms[END_REF] applies and the corresponding modular polynomials can be computed in quasi-linear time.

Proof. This is Theorem 3.11, applied to (in the notations of the Theorem) j 1 = i 1,β , j 2 = i 2,β , j 3 = i 3,β . We only detail the non symmetric case, the adaptations to the symmetric case are obvious. Since Γ = Γ ∩ Γ0 (β), one of the i k,β is not invariant by Γ so, by Proposition 4.11, i k,β generates the field extension C Γ∩ Γ0 (β) /C Γ. Then in the notations of Theorem 3.10 we can use j = i k,β . (In Theorem 4.15 we assume k = 1). It remains to check that the i k,β can be evaluated in time quasi-linear in the precision, but this is obvious from their definition and the fact that the i k can due to Theorem 3.4. Definition 4.16. The polynomials Φ β (X, i 1 , i 2 , i 3 ) and Ψ k,β (X, i 1 , i 2 , i 3 ) for k = 2, 3 defined in Theorem 4.15 are called the β-modular polynomials for i 1 , i 2 , i 3 .

Example 4.17.

• If β = is an inert prime. Then Φ has degree 2 +1 and Ψ k, has degree 2 . If i 1 , i 2 , i 3 are symmetric, then i 1, , i 2, , i 3, also, hence they are invariant under ( Γ∩ Γ0 ( ))∪( Γ∩ Γ0 ( ))σ.

• If β has norm , so = ββ is split. Then if G = Γ is not symmetric, Φ β has degree + 1 and Ψ k,β has degree .

However if σ ∈ G, so that G = Γ σ , then since the i k,β are not symmetric, Φ β has degree 2 + 2 and Ψ k,β has degree 2 + 1. Since Γ is stable under the real conjugation, we can make explicit the action of σ as follows: if we let C β be a set of representatives of Γ ∩ Γ0 (β)\ Γ the modular polynomials are given by

Φ β (X, i 1 , i 2 , i 3 ) = γ∈C β (X -i γ 1,β )(X -i γσ 1,β ) = γ∈C β (X -i γ 1,β )(X -i γ 1,β
) and in Equation [START_REF] Bröker | Explicit CM theory for level 2-structures on abelian surfaces[END_REF]. They are defined for any D and are invariant by Γ(2, 4) (see Definition 2.10). See also Appendix B.1 and B.2 about the modular polynomials with these invariants.

Case D = 2

We have computed the β-modular polynomials with the Gundlach invariants for = 2, 7, 17, 23, 31, 41, 47 and 71. If we write, in the split case,

Φ β (X, J 1 , J 2 ) = X 2 +2 + 2 +1 i=0 c i (J 1 , J 2 )X i and Ψ β (X, J 1 , J 2 ) = 2 +1 i=0 d i (J 1 , J 2 )X i ,
then we have constated that the denominator of c i is of the form P (J 1 , J 2 ) 4 unless i = 2 + 1 where it is P (J 1 , J 2 ) 2 , and that the denominator of d i is of the form P (J 1 , J 2 ) 6 , unless i = 2 +1 where it is P (J 1 , J 2 ) 4 . We have for example for = 7

P (J 1 , J 2 ) = J 2 1 -J 1 J 2 2 + 2J 1 J 2 -81J 1 + 64J 2 2
and for = 17

P (J 1 , J 2 ) = J 7 1 -J 6 1 J 3 2 -6J 6 1 J 2 2 + J 6 1 J 2 -414J 6 1 + 428J 5 1 J 3 2 + 2387J 5 1 J 2 2 - 17760J 5 1 J 2 + 431811J 5 1 + 17728J 4 1 J 4 2 -331952J 4 1 J 3 2 -2578856J 4 1 J 2 2 + 6229197J 4 1 J 2 -80515134J 4 1 -6145536J 3 1 J 4 2 + 52974272J 3 1 J 3 2 + 535037040J 3 1 J 2 2 + 6116816412J 3 1 J 2 + 37822859361J 3 1 -91648000J 2 1 J 5 2 - 6502153216J 2 1 J 4 2 -75793205760J 2 1 J 3 2 -197144611776J 2 1 J 2 2 - 17565696000J 1 J 5 2 -7812042752J 1 J 4 2 + 110592000000J 6 2 .
Table 1 contains some informations about these polynomials. The first column is the prime number, the second the size of the β-modular polynomials, then we have put the total degree and the degree in J 1 and in J 2 of the denominator P (J 1 , J 2 ), and then similarly for the maximal degrees appearing in the numerators. The last column is the number of decimal digits of the largest coefficient appearing in the polynomials.

β (X, b2 , b3 ) = X 2 +2 + 2 +1 i=0 c i ( b2 , b3 )X i and Ψ β (X, b2 , b3 ) = 2 +1 i=0 d i ( b2 , b3 )X i .
when is split, then we have constated that the denominators of c i and of d i are of the form P (J 1 , J 2 ) 4 except for i = 2 + 1 where it is P (J 1 , J 2 ) 2 . We have for example for = 11

P (J 1 , J 2 ) = 4J 7 1 + (-12J 2 2 -19236J 2 + 119497519)J 6 1 + (12J 4 2 + 56972J 3 2 -387805052J 2 2 - 278163835056J 2 + 35953243171744)J 5 1 + (-4J 6 2 -55980J 5 2 + 449730698J 4 2 + 943837290960J 3 2 -133230692691392J 2 2 + 6651010132099840J 2 + 13001634695104256)J 4 1 + (18500J 7 2 -215193500J 6 2 -1170430882000J 5 2 + 388324233980000J 4 2 -32395226716512000J 3 2 )J 3 1 + (32609375J 8 2 + 635091750000J 7 2 -718632513000000J 6 2 + 34620677424000000J 5 2 )J 2 1 + (-124875000000J 9 2 + 601911000000000J 8 2 )J 1 -182250000000000J 10 2 .
We have computed the β-modular polynomials with theta invariants for = 5, 11, 19, 29, 31, 41 and 59. These polynomials are

Φ β (X, b1 , b2 , b3 ) = X +1 + i=0 ( 4 j=0 c i,j ( b1 , b2 ) bj 3 )X i and Ψ k,β (X, b1 , b2 , b3 ) = X +1 + i=0 ( 4 j=0 d k,i,j ( b1 , b2 ) bj 3 )X i ,
by Equation ( 28) and Appendix B.2. Table 3 contains the same informations as Table 1, but the first part concern the polynomials with the Gundlach invariants and the second the polynomials with the bi invariants. We can see that there is a gain in terms of memory space, except for = 5, which corresponds to the ramified case. The degrees are larger with the bi but there also are relations modulo 4 between the exponents.

Examples of isogenous curves

The modular polynomials allow one to compute hyperelliptic curves with isogenous Jacobians. In particular over finite field as the β-polynomials found can be reduced modulo a prime number p = ββ without loosing their meaning ([6, Section 6, page 511]).

We begin with examples of curves found when working on Q( √ 2) and taking the Gundlach invariants. The Jacobians of the following curves are (

Y 2 = 13477X 5 + 6136X 4 + 35146X 3 + 28148X 2 + 7150X + 19730, Y 2 = 2953X 5 + 26725X 4 + 14100X 3 + 6565X 2 + 22149X + 19740
and second example for (5 + 2(1 + √ 5)/2)-isogenies over F 6728947 :

Y 2 = 3739712X 6 + 4881762X 5 + 6611129X 4 + 5775262X 3 + 521647X 2 + 2066678X + 350732, Y 2 = 2707309X 6 + 1535264X 5 + 311501X 4 + 2965267X 3 + 3507011X 2 + 101110X + 5795310.
Finally, we give pairs of curves, whose Jacobians are (7 + 2 √ 2)-isogenous over F 562789 , computed using the β-modular polynomials with the bi for Q( √ 2):

Y 2 = 540913X 5 + 353915X 4 + 118050X 3 + 355166X 2 + 424096X + 379433, Y 2 = 231396X 5 + 474300X 4 + 200176X 3 + 335056X 2 + 345222X + 464702
and a pair for (5 -(1 + √ 5)/2)-isogenies over F 5362789 , computed using the polynomials with the bi for Q( √ 5):

Y 2 = 2531476X 5 + 900554X 4 + 1248025X 3 + 440959X 2 + 912166X+ 4367293, Y 2 = 1772175X 5 + 3557482X 4 + 848889X 3 + 4562893X 2 + 146681X+
475016.

The motivated reader can check that the curves are indeed isogenous in verifying that the curves have the same zeta function (by a result of Tate [START_REF] Tate | Endomorphisms of Abelian Varieties over Finite Fields[END_REF]Theorem 1]).

A Exemple of invariants on Hilbert and Humbert surfaces

A.1 Pullback of Siegel invariants

We prove here Lemma 2.9 which states that we can always use the pullback of Igusa's invariants to get invariants for the Humbert surface.

Lemma A.1. Let X be a subvariety of Y , with both X and Y irreducible and defined over a field F . Then the restriction map (which is not defined everywhere) on the functions fields

F (Y ) -→ • F (X) is surjective.
Proof. This result is well-known. We give a proof for convenience. Since X is a subvariety of Y , it is a closed variety of an open locus U of Y . The inclusion ι : X → U then yields an epimorphism of sheaves ι * : O U → O X . Looking at the stalks of the generic points we deduce that the map F (Y ) → F (X) (defined for functions f ∈ F (Y ) which are defined on the generic point of X) is surjective.

Proof of Lemma 2.9. By the theory of Shimura varieties, both ( Γ(1)∪ Γ(1)σ)\H 2 1 and Sp 4 (Z)\H 2 are algebraic, and so is ρ, the application of Proposition 2.8. This last proposition says that the map from the Symmetric Hilbert modular surface ( Γ(1) ∪ Γ(1)σ)\H 2 1 to the Siegel space is birational to its image, the Humbert surface H ∆ K . Its field of functions are the symmetric Hilbert modular functions. So, by Lemma A.1, any symmetric Hilbert modular function (seen by birationality as a rational function on the Humbert surface) can be lifted to a Siegel modular function. Since the Igusa invariants generate the field of the Siegel modular functions, it suffices to check that the restriction of these invariants to H ∆ K is well defined (on an open set). But the denominators of these functions is (up to a scalar multiple) χ 10 (see Equation ( 1)) whose locus is exactly H 1 , the set of abelian surfaces isomorphic to a product of elliptic curves. By Proposition 2.7 the intersection of H 1 and H ∆ K is a (union of) curves, so the Igusa invariants are well defined on

H ∆ K \ H 1 .
As seen in Section 2.5, we can also always take pullback of theta functions (to get invariants in some level higher than one). These pullback will be studied further in Section A.5. By Theorem 2.16 these give non symmetric invariants when D ≡ 1 mod 4. More generally to get a non symmetric invariant of level 1 we could take traces of pullback of theta functions of higher levels.

In practice we used pullback of theta functions or Gundlach invariants to compute the (symmetric) Hilbert modular polynomials. The reader interested in computing them with pullback of Igusa invariants is advised to use a variant defined by Streng in [START_REF] Streng | Complex multiplication of abelian surfaces[END_REF] to get smaller coefficients.

A.2 Gundlach invariants for

Q( √ 2) and Q( √ 5
)

When K = Q( √ 2) or K = Q( √ 5
) there is a fundamental unit of norm -1 and such that > 0 as a real number. Let α = diag(1,

√ ∆ K ). Then (25) φ 0 : H 2 1 → H 2 1 τ → √ ∆ K τ and φ 0 : SL 2 (O K ) → SL 2 (O K ⊕ ∂ K ) γ → αγα -1
are group isomorphisms which induce a commutative diagram similar as the one in Equation [START_REF] Bröker | Modular polynomials for genus 2[END_REF]. Note that when > 0 has norm -1, then < 0 so that √ ∆ K is totally positive and φ 0 (τ ) ∈ H 2 1 . Let {e 1 , e 2 } be a Z-basis of O K and q j = e 2iπ( e j τ 1 -e j τ 2 )/ √ ∆ K for j = 1, 2 and τ = (τ 1 , τ 2 ). Proposition A.2. Let a ⊂ K be a fractional ideal and SL

2 (O K ⊕ a) = { a b c d ∈ SL 2 (K) : a, d ∈ O K , b ∈ a -1 and c ∈ a}. Let a * = ∂ -1 K a -
g(τ ) = a g (0) + t=ae 1 +be 2 ∈O ++ K a g (t)q a 1 q b 2 .
Proof. See [START_REF] Lauter | Computing genus 2 curves from invariants on the Hilbert moduli space[END_REF]Proposition 3.2]. By Proposition A.2 the Fourier coefficients are indexed by

∂ -1,++ K , but if v ∈ ∂ -1,++ K , then v = √ D t with t = ae 1 + be 2 ∈ O ++ K .
We denote by A Z (SL 2 (O K )) k the Z-module of symmetric Hilbert modular forms of even weight k with rational integral Fourier coefficients and put

A Z (SL 2 (O K )) = A Z (SL 2 (O K )) k .
Define the Hilbert Eisenstein series of even weight k ≥ 2:

G k (τ ) = 1 + t=ae 1 +be 2 ∈O ++ K b k (t)q a 1 q b 2 , where b k (t) = κ k tO K ⊂µO K | O K /µO K | k-1 and κ k = ζ K (k) -1 (2π) 2k ((k -1)!) -2 ∆ 1/2-k K
(by [56, Equation (1.5)]).

Lemma A.4 (Gundlach).

• If K = Q( √ 2), let = 1+ √ 2. Then κ 2 = 2 4 •3, κ 4 = 2 5 •3•5•11 -1 and κ 6 = 2 4 •3 2 •7•19 -2 ; • If K = Q( √ 5), let = 1+ √ 5 2 . Then κ 2 = 2 3 • 3 • 5, κ 4 = 2 4 • 3 • 5, κ 6 = 2 3 • 3 2 • 5 • 7 • 67 -1 and κ 10 = 2 3 • 3 • 5 2 • 11 • 412751 -1 .
Proof. See [START_REF] Gundlach | Die Bestimmung der Funktionen zur Hilbertschen Modulgruppe des Zahlkörpers Q( √ 5)[END_REF][START_REF] Gundlach | Die Bestimmung der Funktionen zu einigen Hilbertschen Modulgruppen[END_REF]. See also [START_REF] Nagaoka | On the ring of Hilbert modular forms over Z[END_REF]Lemma 1.1].

The Hilbert Eisenstein series are symmetric Hilbert modular forms for SL 2 (O K ) with coefficients in Q. We focus now on the cases D = 2, 5 and we fix the basis {1, }, which gives a nice expression of q 1 and q 2 . We have Theorem A.5 (Nagaoka). In the case K = Q( √ 2), we put

F 4 = 2 -6 • 3 -2 • 11(G 2 2 -G 4 ) and F 6 = -5 • 7 2 2 8 3 3 13 G 3 2 + 11 • 59 2 8 3 2 5 • 13 G 2 G 4 - 19 2 2 7 3 3 5 • 13 G 6 .
Then Theorem A.8 (Nagaoka). In the case K = Q( √ 5), we put

F 6 = 67 2 5 3 3 5 2 (G 3 2 -G 6 ), F 10 = 2 -10 3 -5 5 -5 7 -1 (412751G 10 -5 • 67 • 2293G 2 2 G 6 + 2 2 3 • 7 • 4231G 5 2 ), and
F 12 = 2 -2 (F 2 6 -G 2 F 10
). The four modular forms G 2 , F 6 , F 10 and F 12 are in A Z (SL 2 (O K )) k for k = 2, 6, 10 and 12 respectively. Furthermore, they form a minimal set of generators of A Z (SL 2 (O K )) over Z.

Proof. See [START_REF] Nagaoka | On the ring of Hilbert modular forms over Z[END_REF]Theorem 2]. See also [START_REF] Gundlach | Die Bestimmung der Funktionen zur Hilbertschen Modulgruppe des Zahlkörpers Q( √ 5)[END_REF]Satz 5].

Theorem A.9 (Gundlach). In the case K = Q( √ 5), the field of symmetric meromorphic Hilbert modular functions for SL 2 (O K ) are rational functions of

J 1 = G 5 2 F 10 and J 2 = F 6 G 2 2 F 10 .
Proof. See [START_REF] Gundlach | Die Bestimmung der Funktionen zur Hilbertschen Modulgruppe des Zahlkörpers Q( √ 5)[END_REF]Satz 6]. See also [START_REF] Müller | Hilbertsche Modulformen und Modulfunktionen zu Q( √ 5)[END_REF]Page 249] or the proof in Appendix A.3. Note that it is usual to take the invariants

G 5 2 F 10 and F 6 G 3 2
. We have substituted the last one by the product of the two. As explained in Section 4.3 these invariants will give smaller modular polynomials. Indeed we will see that the denominators of the invariants determine the denominators of the modular polynomials so that it is better to have fewer factors.

Definition A.10. The functions J 1 and J 2 of Theorem A.9 are called the Gundlach invariants for Q( √ 5).

A.3 Gundlach and pullbacks of the Igusa invariants

For D = 2, 5, a fundamental unit has norm -1 and it will be more convenient to work with the basis {1, }, which was used to define the Fourier coefficients of the symmetric Hilbert modular forms in Appendix A.

(Ω) = a f (0) + T ∈Sym 2 (Z) ∨,++
a f (T )q T be a holomorphic Siegel modular form for Sp 4 (Z) of weight k. Then its pullback g = φ * e 1 ,e 2 f is a symmetric Hilbert modular form with the following Fourier expansion:

g(τ ) = f (φ e 1 ,e 2 (τ )) = a g (0) + t=ae 1 +be 2 ∈O ++ K a g (t)q a 1 q b 2 , with a g (0) = a f (0) and a g (t) = T ∈Sym 2 (Z) ∨,++ Q T (e 1 ,e 2 )=t a f (T ).
Here,

Q T (x 1 , x 2 ) = (x 1 , x 2 )T ( x 1
x 2 ) is the positive definite quadratic form associated to T and

Sym 2 (Z) ∨ = T = m 1 1 2 m 1 2 m m 2 : m i , m ∈ Z
is the dual of Sym 2 (Z). Finally, q T = e 2iπtr(T Ω) .

Proof. See [START_REF] Lauter | Computing genus 2 curves from invariants on the Hilbert moduli space[END_REF]Proposition 3.2].

We are interested in the pullbacks of the Igusa invariants (defined in Equation ( 1)). They are already known in the case D = 5.

Theorem A.12 (Resnikoff). Theorem A.16. The Hilbert modular surface is rational for D = 2, 3, 5, 13, 17.

Proof. See [START_REF] Hirzebruch | Classification of Hilbert Modular Surfaces[END_REF]Theorem 3].

The paper [START_REF] Elkies | K3 surfaces and equations for Hilbert modular surfaces[END_REF] of Elkies and Kumar allows to get the corresponding birational isomorphisms. For instance, for Q( √ 5), they give an explicit birational isomorphism between the Humbert surface H 5 and P 2 (different from the one induced by the Gundlach invariants), and show that a birational model over Q of the non symmetric Hilbert modular surface is given by the double cover of P 2 z 2 = 2(6250h 2 -4500g 2 h -1350gh -108h -972g 5 -324g 4 -27g 3 ).

As this surface is also rational, a parametrization is obtained, given by the modular functions m and n. We have by [ from which we can compute m, n explicitly. The functions g and h are easy to evaluate from the Gundlach invariants. For z, we use the expression of F i in terms of theta constant or the equation of the double cover given above and the first coefficients of the Fourier series of z for the choice of the square root. More generally, equations for Humbert and Hilbert surfaces are given in [START_REF] Elkies | K3 surfaces and equations for Hilbert modular surfaces[END_REF] for every quadratic field Q( √ ∆) for all thirty fundamental discriminants ∆ with 1 < ∆ < 100. We can then use their results to get non symmetric invariants for all these fields. Furthermore, the equation of the Humbert surface for these fields is always rational, so we only need two symmetric invariants (like the Gundlach case for Q( √ 5) and Q( √ 2)) and don't need to take the pullback of the Igusa invariants for these fields. The difficulty resides in the optimization of these invariants: for instance for computing modular polynomials it is better if they have the same denominator.

A.5 Components of Humbert surfaces of level 2 and (2, 4)

This is a complement to Section 2.5.

Proposition A.17 Let L be the locus of the principally polarized abelian surfaces with real multiplication by O K which are β-or β-isogenous to a product of elliptic curves (and which are not isomorphic to a product of elliptic curves because when this happens, the Gundlach invariants are not always defined). Proposition B.6. In the case where D = 5, the denominators of the modular polynomials Φ β and Ψ β are divisible by a polynomial L in Q[J 1 , J 2 ] describing L .

Proof. We adapt the proof of [START_REF] Bröker | Modular polynomials for genus 2[END_REF]Lemma 6.3]. Let z ∈ H 2 1 which is β-or β-isogenous to a product of elliptic curves and let c i be a coefficient of Φ β . The cusp form χ 10 vanishes at products of elliptic curves and by Theorem A.12, we have F 10 = -4φ * χ 10 so that F 10 also vanishes at product of elliptic curves. Thus J 1 and J 2 have poles at these values and there exists some γ ∈ Γ0 (β)\ Γ(1) such that J γ 1,β (z) or J γ 1,β (z) is infinite. The evaluation of c i at z is a symmetric expression in the J γ 1,β (z) and in the J γ 1,β (z). Generically, there is no algebraic relation between these values and the evaluation of c i at z is therefore infinite. Since J 1 (z) and J 2 (z) are finite, the numerator of c i is finite. The denominator of c i must vanish at z which means that c i is divisible by L . The proof for Ψ β is similar.

If D = 2, the Gundlach invariants J 1 and J 2 have poles when F 4 (z) = 0. Since by Theorem A.14, we have that φ * χ 10 = -1 4 F 4 F 6 , the set of poles is a subset of the products of elliptic curves. We have thus to consider the subset L of L of the surfaces z such that We have proved that we have in the denominators of the modular polynomials a subset of the set H β of abelian surfaces which are β-isogenous to a product of elliptic curves (and which are not isomorphic to a product of elliptic curves).

B.2 Modular polynomials in theta invariants

In this section, we define modular polynomials for any D square-free by using theta constants. They illustrate nicely the different possibilities of Theorem 4.15. We use the action of Γ(2, 4)\(SL 2 (O K ⊕ ∂ K ) ∪ SL 2 (O K ⊕ ∂ K )σ) to prove symmetries of these polynomials and accelerate their computations.

The invariants we use are the pullbacks bi = φ * 1,ω b i for i = 1, 2, 3, of the generators b 1 , b 2 , b 3 (defined in Equation ( 3)) for the group Γ [START_REF] Becker | On Gröbner bases under specialization[END_REF][START_REF] Birkenhake | Humbert surfaces and the Kummer plane[END_REF]. They are modular functions for the group Γ(2, 4), which is defined in Equations ( 16) and [START_REF] Elkies | Elliptic and modular curves over finite fields and related computational issues[END_REF] and are generators for a corresponding field of modular functions by Theorem 2.16. Recall that we define Γ(1) = SL 2 (O K ⊕ ∂ K ).

Evaluation and inversion:

We now outline efficient procedures for the computation of the values bi (τ ) of Hilbert modular functions at any τ ∈ H 2 1 and for finding some τ ∈ H 2 1 from the bi (τ ). The first one is similar to Algorithm B.3, the third step being trivial as bi = φ * 1,ω b i , and has the same complexity. For the second procedure, we also proceed as in Algorithm B.1, the first step being also trivial. For the second, it is possible to find Ω modulo Γ(2, 4) in Õ(N ) time (see [START_REF] Milio | A quasi-linear time algorithm for computing modular polynomials in dimension 2[END_REF]). The difficulty is in the third step. Indeed, we are able to find γ such that φ 1,ω (τ ) = γΩ, but γ is not necessarily in Γ(2, 4) so that we only find τ modulo Γ(1) ∪ Γ(1)σ instead of τ modulo Γ(2, 4), if D ≡ 1 mod 4, or modulo Γ(2, 4) ∪ Γ(2, 4)σ, if D ≡ 2, 3 mod 4. A solution consists in computing beforehand all the classes of the quotient Γ(2, 4)\ Γ(1) and of Γ(2, 4)σ\ Γ(1) and see how they are sent to the classes of Γ(2, 4)\Sp 4 (Z). It suffices to find in which class of Γ(2, 4)\Sp 4 (Z) the matrix γ belongs in to find a corresponding matrix γ in Γ(2, 4)\ Γ(1) or in Γ(2, 4)σ\ Γ(1). Then we have φ 1,ω (γ -1 τ ) = φ 1,ω (γ -1 )φ 1,ω (τ ) = γ -1 γΩ = Ω.

Corollary B.8. We can evaluate the three bi (τ ) at precision N for τ ∈ H 2 1 in Õ(N ) time and we can find τ at precision N modulo Γ(2, 4), or modulo Γ(2, 4) ∪ Γ(2, 4)σ according to the cases, from the values bi (τ ) with this same complexity.

Note that when we use the function bi to define modular polynomials, for the interpolation step we need the equations of the Humbert component defined by the bi , as explained in Section 3.3. We refer to Equation [START_REF] Gaudry | The 2-adic CM method for genus 2 curves with application to cryptography[END_REF] for the equations for D = 2, 3, 5 and to [START_REF] Gruenewald | Explicit algorithms for Humbert surfaces[END_REF] for larger discriminants.

12 χ 6 10 , j 2 = 2 -3 3 3 ψ 4 χ 3

 1210223

Theorem 2 . 16 .

 216 The functions rk and bk for k = 1, 2, 3 are generators for the field of Hilbert modular functions invariants by Γ(2) and Γ(2, 4), if D ≡ 1 mod 4, and by Γ(2) ∪ Γ(2)σ and Γ(2, 4) ∪ Γ(2, 4)σ, if D ≡ 2, 3 mod 4, respectively.

Proposition 4 . 4 .

 44 Let β be a totally positive element of norm . The number of β-isogenies is + 1. They correspond to cyclic kernels of size in A[β], which are stable by O K .

Lemma 4 . 5 .

 45 Let = ββ be a splitting of into totally positive elements ofO ++ K . Let V ⊂ A[β]be the kernel of a β-isogeny. Let ∈ O × K be a unit, so that 2 is totally positive and we have another splitting of as = ( 2 β)( 2 β).

Lemma 4 . 8 .

 48 The quotient Γ ∩ Γ0 (β)\ Γ is of cardinality L + 1.Proof. The sets Γ ∩ Γ0 (β)\ Γ and Γ0 (β)\ Γ0 (β) Γ are in bijection so by Propositions 4.3 and 4.4 it suffices to prove that Γ0 (β) Γ = Γ[START_REF] Baily | Compactification of Arithmetic Quotients of Bounded Symmetric Domains[END_REF]. So it suffices to prove that Γ(L) Γ(n) = Γ(1), which is obvious by the Chinese Reminder Theorem.Indeed by Lemma 4.7 it suffices to check that π : SL 2 (O K ) → SL 2 (O K /LnO K ) is surjective on Γ(L) Γ(n) (since this group contains the kernel). But since n is prime to L, SL 2 (O K /LnO K ) SL 2 (O K /LO K ) × SL 2 (O K /nO K ) and π( Γ(n)) contains the left factor while π( Γ(L)) contains the right factor.

  A+B B C+D D and (A + B)(0) = 0, which ends the proof.

Definition 4 . 12 .Lemma 4 . 13 .

 412413 Let γ ∈ Γ0 (β) = a b c d . We denote γ β = a b/β cβ d ∈ Γ(1). Let i be a meromorphic function H 2 1 → C, and define i β

2

 2

G 3 2 -2 5 3 3 F 6 ; 6 -

 2366 -4φ * χ 10 = F 10 ; 12φ * χ 12 = 3F 2 2G 2 F 10 .Proof. See[START_REF] Resnikoff | On the Graded Ring of Hilbert Modular Forms Associated with Q( √ 5)[END_REF] Theorem 1].

18 ,J 5 2 and z = -F 3 5 F

 1825 Section 6] m = -(5g 2 + 3g/2 -125h/9 + 3/25)/(g 2 + 13g/30 + 1/25), n = z/(18(g 2 + 13g/30 + 1/25))andg = (m 2 -5n 2 -9)/30, k = 3m(10g + 3)(15g + 2)/6250, h = k + 9(250g 2 + 75g + 6)/6250, z = 3n(10g + 3)(15g + 2)/25.Using these equations,[START_REF] Elkies | K3 surfaces and equations for Hilbert modular surfaces[END_REF] Corollary 15] (linking the Igusa-Clebsh invariants with g and h) and Theorem A.12, we have found the relations(27) g = -J 1 /(6J2 2 ), h = J 2 1 /15 /(2F 5 6 )

.

  The subgroups Γ(2) and Γ(2, 4) of Γ(1) are of index      36 and 576, if D ≡ 1 mod 8; 60 and 960, if D ≡ 5 mod 8; 48 and 192, if D ≡ 2, 3 mod 4.

F 4 Proposition B. 7 .

 47 ( 1 β γ • z) = 0 or F 4 ( 1 β γ • z) = 0 for some γ ∈ C β .In the case where D = 2, the denominators of the modular polynomials Φ β and Ψ β are divisible by a polynomial L in Q[J 1 , J 2 ] describing L .

. 9 .Proposition B. 10 .

 910 Modular polynomials: Recall that we denote for i = 1, 2, 3, β ∈ O ++ K and γ ∈ Γ(1)∪ Γ(1)σ: bi,β :H 2 1 → C τ → bi ( 1 β τ ) and bγ i,β : H 2 1 → C τ → bi ( 1 β γ • τ ).For a matrix γ ∈ Γ(2, 4) ∩ Γ0 (β), we would like to write bγi,β (τ ) = bi ( 1 β γ • τ ) = bi (γ β • ( 1 β τ )) = bi ( 1 β τ ) = bi,β(τ )so that the functions bi,β for i = 1, 2, 3 would be modular for the group Γ(2, 4)∩ Γ0 (β). However the third equality is true only if the matrix γ β is in Γ(2, 4) (see Corollary 4.14). A simple calculation shows that this is always the case when D ≡ 1 mod 4. When D ≡ 2, 3 mod 4, this happens only when β is of the form a + bω with b even. If D ≡ 2 mod 4, this is equivalent to ask that ≡ 1 mod 4 and else if D ≡ 3 mod 4, must necessarily satisfy ≡ 1 mod 4. In particular, in the last case, 0, 1 or 2 modular polynomials with Γ(2, 4) structure can exist for a given prime which splits in totally positive factors, according to the fundamental unit . Thus Proposition BThe functions bi,β for i = 1, 2, 3 are modular functions for Γ(2, 4) ∩ Γ0 (β) when • D ≡ 1 mod 4; • D ≡ 2 mod 4 and β = a + bω with b even, or, equivalently, ≡ 1 mod 4; • D ≡ 3 mod 4 and β = a + bω with b even; this implies that ≡ 1 mod 4. Let be a prime number. Write = β if is inert and = ββ if is split with β ∈ O ++ K . Let C β be a set of representatives of ( Γ(2, 4) ∩ Γ0 (β))\ Γ(2, 4). If D ≡ 1 mod 4, then the polynomials Φ β (X, b1 , b2 , b3 ) = γ∈C β (Xbγ 1,β ), and Ψ k,β (X, b1 , b2 , b3 ) = γ∈C β bγ k,β Φ β (X, b1 , b2 , b3 ) X -bγ 1,β

  Let A Ω be the principally polarized abelian surface associated to Ω ∈ H 2 .

	or for example Birkenhake-Wilhelm [4, Proposition 4.5] or Runge [61,
	Theorem 2].
	Remark 2.5. Let Ω ∈ H 2 . It can satisfy many singular relations of different discriminants.
	If it satisfies a singular relation of discriminant ∆, a constructive algorithm to find γ as in
	Humbert's Lemma can be found in [4, 61]. Conversely, let Ω ∈ H 2 be a matrix equivalent to
	a matrix satisfying (11). Then Ω satisfy necessarily a singular relation of discriminant ∆. By
	(7)
	Proposition 2.6. For any ∆ ≡ 0 or 1 mod 4, ∆ > 0, the set H ∆ := {{Ω} ∈ Sp 4 (Z)\H 2 : Ω
	satisfies a primitive singular relation of discriminant ∆} is a surface which we call the Humbert
	surface of discriminant ∆.
	Proof. See [4, Corollary 4.6 and Proposition 4.7] or [30, Proposition 2.11].
	Proposition 2.7.

  where E is an elliptic curve.

	Proof. See [4, Proposition 4.9] or [30, Corollary 2.10, Proposition 2.15].

Table 1 :

 1 Informations about the modular polynomials for D = 2 (Gundlach invariants)

	8.5 KB	3	0	3	4	4	2	8
	7 172 KB 3	2	2	25	23	13	66
	17 5.8 MB 9	7	6	65	61	36	196
	23 21 MB 12 11 8	87	85	48	280
	31 70 MB 17 14 10 117 111 61	401
	41 225 MB 23 21 14 157 153 84	560
	47 400 MB 26 25 16 179 177 96	665
	71 2.2 GB 42 37 24 275 265 144 1078
	We have computed the β-modular polynomials with theta invariants for = 17, 41, 73,
	89 and 97 (which are 1 modulo 4, see Proposition B.9). By Remark B.11, the β-modular
	polynomials are
	Φ

Table 3 :

 3 Informations about the modular polynomials for D = 5

	5	22 KB	5	3	5	10	10	10	53
	11 3.5 MB 10	7	10	40	40	40	252
	19 33 MB	16	12	16	64	64	64	513
	29 188 MB 25	20	25 100 100 100 830
	31 248 MB 26	21	26 104 104 104 885
	41 785 MB 35	29	35 140 140 140 1191
	59 3.6 GB	50	43	50 200 200 200 1820
	5	26 KB	16	8	8	31	19	22	5
	11 308 KB 72	40	40	84	52	52	11
	19 3.6 MB 128 96	96 132 103 108	25
	29 21 MB 200 152 152 212 160 168	44
	31 28 MB 216 160 160 224 173 172	47
	41 115 MB 288 240 240 324 272 272	69
	59 470 MB 424 352 352 440 373 370 109

  3 + √2)-isogenous over F 2333 :Y 2 = 356X 6 + 116X 5 + 1589X 4 + 986X 3 + 178X 2 + 1094X + 1229, Y 2 = 144X 6 + 2096X 5 + 387X 4 + 1562X 3 + 478X 2 + 486X + 1718

	while the Jacobians of the following ones are (5 + 2	√ 2)-isogenous over F 345267203 :
	Y 2 = 288618938X 5 + 208826828X 4 + 73681500X 3 + 329580565X 2 +
	193693317X + 328425210,
	Y 2 = 229859713X 5 + 180037958X 4 + 95105703X 3 + 68631100X 2 +
	32660205X + 107566399
	and the Jacobians of the curves hereafter are (7 +	√	2)-isogenous over F 3526982779 :
	Y 2 = 3476666651X 5 + 2997006123X 4 + 2343918968X 3 + 1313289865X 2 +
	1251164949X + 1521154595,
	Y 2 = 2390845907X 6 + 2649299485X 5 + 3307186776X 4 + 2143442296X 3 +
	1448110737X 2 + 918458873X + 1476608496.
	We also give two examples of pairs of curves computed with the β-modular polynomials with the Gundlach invariants for Q( √ 5). First example of curves for (4 -(1 + √ 5)/2)-isogenies
	over F 56311 :

  1 be the dual of a with respect to the trace, and a * ,++ the set of totally positive elements of a * .Let g be a holomorphic Hilbert modular form for SL 2 (O K ⊕ a) of weight k.

		Then it has a
	Fourier expansion
	g(τ ) = a g (0) +	a g (v)e 2πitr(vτ ) .
	v∈a * ,++
	Proof. See [9, Example 1.6 and Corollary 1.21].

Corollary A.3. Let g be a holomorphic Hilbert modular form for SL 2 (O K ) of weight k. Then it has a Fourier expansion

  G 2 , F 4 and F 6 are in A Z (SL 2 (O K )) k for k = 2, 4, 6 respectively. Furthermore, they form a minimal set of generators of A Z (SL 2 (O K )) over Z. The functions J 1 and J 2 of Theorem A.6 are called the Gundlach invariants

	Proof. See [56, Theorem 1]. See also [32, Page 119].
	Theorem A.6 (Gundlach). In the case K = Q( √	2), the field of symmetric meromorphic
	Hilbert modular functions for SL 2 (O K ) are rational functions of
	J 1 =	G 2 2 F 4	and J 2 =	G 2 F 6 4 F 2	.
	Proof. See [32, Page 119]. See also [53, Lemma 6] or Appendix A.3.
	Definition A.7. for Q( √ 2).

  where φ 0 denote the isomorphisms of Equation[START_REF] Gathen | Modern Computer Algebra[END_REF]. The map φ satisfies similar equalities as in Proposition 2.3 between the action of SL 2 (O K ) on H2 1 and the action of Sp 4 (Z) on H 2 . It also maps to the Humbert surface and Proposition 2.8 apply with φ instead of φ e 1 ,e 2 .For a basis {e 1 , e 2 }, we give in Proposition A.11 the relation between the Fourier coefficients of a Siegel modular form f and the coefficients of its pullback φ * e 1 ,e 2 f , which is a symmetric Hilbert modular form. Recall Equation (6) for the definition of φ e 1 ,e 2 .

	2. Let
	(26)	φ 1 := φ 1,	and	φ := φ 1 • φ 0 ,
	Proposition A.11. Let
	f

Deduce j 1 (Ω), j

(Ω), j

(Ω) from J 1 (τ ), J 2 (τ ), where Ω ∈ H 2 such that Ω = φ (τ ); 2 Deduce a period matrix Ω at precision N equivalent to Ω modulo Sp

(Z), up to the precision, from the three Igusa invariants; 3 Find some γ ∈ Sp 4 (Z) such that γΩ satisfies Equation[START_REF] Couveignes | Elliptic periods for finite fields[END_REF] and deduce τ ;

.

In this case the β-modular polynomials parametrize both β and β-isogenies (so they are equal to the β-modular polynomials). This is the drawback for the applications of Appendix ??, hence the interest to also have non symmetric invariants, even if they are harder to compute. , then there are also β-isogenies. (We only consider totally positive units to guarantee the fact that z ∈ H 2 1 ). If there exists n ∈ Z such that = 2n , then the matrix γ = n 0 0 -n is in Γ and γ • z = z. Thus, in this case, i k ( z) = i k (z), and, in particular, a β-isogeny is also a β-isogeny. (For a more intrinsic proof see Lemma 4.5.) When D = 2 or 5, a fundamental unit has norm -1 while ∈ O ×,+ K has norm 1, so that the latter can always been written as an even power of . Thus, the choice of the splitting of does not matter.

Remark 4. [START_REF] Enge | Computing modular polynomials in quasi-linear time[END_REF] (General modular polynomials). Let Γ ⊂ SL 2 (O K ⊕ ∂ K ) that does not satisfy Equation [START_REF] Enge | Comparing invariants for class fields of imaginary quadratic fields[END_REF]. Then from a level structure G associated to a triple (A, θ, V ) corresponds several level structures G on (A/V, θ ). From Corollary 4.14 the modular functions i k,β are modular for the group Γβ = {γ ∈ SL 2 (O K ⊕ ∂ K ) | γ β ∈ Γ} ⊂ Γ0 (β). So we can define modular polynomials in a similar way as in Theorem 4.15 except that we act by Γ ∩ Γβ \ Γ. The fibers correspond to β-isogenies together with an extra structure determined by the action of Γ ∩ Γβ \ Γ ∩ Γ0 (β). So we loose the corresponding factor in the degree of the modular polynomials. A possible solution would be to replace i 1,β by its trace under the action of Γ ∩ Γβ \ Γ ∩ Γ0 (β) to get a modular function invariant by Γ ∩ Γ0 (β).

Also, if Γ does not contain a congruence subgroup Γ(n) of level n prime to , then Γ ∩ Γ(β)\ Γ may not be isomorphic to Γ(β)\ Γ [START_REF] Baily | Compactification of Arithmetic Quotients of Bounded Symmetric Domains[END_REF], but only isomorphic to a subgroup. We can still compute modular polynomials, but they will not parametrize all β-isogenies, only those who are compatible with the structure induced by Γ.

Finally if β ∈ O ++ K is totally positive but not prime, it is easy to adapt Theorem 4.15 (if we suppose that i 1 is not invariant by Γ ∩ Γ0 (I) for strict divisors ideal I of (β)). The only difference is on the degree of the polynomials, Φ β will not be of degree the norm of β. Rather the degree depends on the factorization of (β) into prime ideals.

Of course this whole discussion is easily extended to the symmetric case.

Results

The aim of this section is to present some polynomials we have computed and to compare the polynomials with the different invariants when this comparison makes sense. All the polynomials computed are accessible at https://members.loria.fr/EMilio/ . In particular, we do not present here the polynomials computed when is inert or when D = 3.

The first invariants we used are the Gundlach invariants J 1 , J 2 , which are defined in Appendix A.2. They are defined for D = 2, 5 and are invariant by SL 2 (O K ). And we also used the pullbacks bi = φ * 1,ω b i of the b i functions defined in Equation ( 3) by the map φ 1,ω defined

We have constated that the denominators of c i and d i are of the form P ( b2 , b3 )

where it is P ( b2 , b3 ). For example, we have for = 17 and β = 5 + 2 √ 2 P ( b2 , b3 ) = b6 For = 17 and 41, the degrees of the coefficients c i and d i in the variables b2 and b3 are close to the degrees in the variables J 1 and J 2 . But with the bi , some relations between the exponents occur. The numerator of c i can be written as m n c i,m,n bm 2 bn 3 (and similarly for d i ). We have then for = 17 and

for c i and d i respectively. In the case = 41 and β = 7 + 2 √ 2, these equations are the same except the last which is m + n ≡ -i mod 4 for c i and d i . Comparing Tables 1 and2, we can see that taking the invariants based on the theta functions give better results. But, here, this is the case only when ≡ 1 mod 4. For = 7 ( ≡ 3 mod 4), we have done as explained at the end of Appendix B.2. On the one hand, we have computed the polynomials using the subgroup of index 4( + 1) and on the other hand, we have computed the polynomials using the Rosenhain invariants. The first solution give better results in terms of degree, sparsity and the whole polynomials fill 930 KB in the first case while 70 MB in the second. In both cases, the polynomials are bigger than those using the Gundlach invariants. This is also true for = 23, where using the first method, the polynomials fill 110 MB.

Case D = 5

We have computed the β-modular polynomials with the Gundlach invariants for = 5, 11, 19, 29, 31, 41 and 59. If we write

Corollary A.13 (Lauter-Yang). One has

We did not find in the literature a similar result as Theorem A.12 for Q( √ 2). Using Proposition A.11 and comparing the different Fourier series (as done in [START_REF] Resnikoff | On the Graded Ring of Hilbert Modular Forms Associated with Q( √ 5)[END_REF] in the case D = 5) we have found Theorem A.14.

. Proof of Theorems A.6 and A.9. By Lemma 2.9, any symmetric Hilbert modular function is a rational fraction with complex coefficients in the pullbacks of the Igusa invariants. By Corollaries A.15 and A.13, the pullbacks of the Igusa invariants can be expressed in terms of the Gundlach invariants for Q( √ 2) and Q( √ 5) respectively. Thus each symmetric Hilbert modular function can be expressed in terms of the Gundlach invariants.

A.4 Other invariants

By Müller [START_REF] Müller | Hilbertsche Modulformen und Modulfunktionen zu Q( √ 5)[END_REF], non-symmetric Gundlach invariants for Q( √ 5) can be obtained considering the Hilbert modular forms

F 2 5 = F 10 and by defining the modular function J 3 = F 15 /F 3 5 . To use interpolation to compute nonsymmetric Hilbert modular polynomials for J 1 , J 2 and J 3 , we need the equation of the Hilbert modular surface, which is given by

2 )/(16J 2 1 ). Müller gives an expression of F 5 and of F 15 in terms of theta constants so J 3 can be efficiently evaluated at τ ∈ H 2 1 . The polynomials obtained are smaller than the symmetric ones with J 1 and J 2 only, as they are not symmetric and thus parameterize only β-isogenies. We refer to Martindale [START_REF] Martindale | Isogeny Graphs, Modular Polynomials, and Applications[END_REF] for more details on the polynomials coming from these invariants.

Similarly for Q( √ 2), see the work of Müller [START_REF] Müller | Hilbertsche modulformen und modulfunktionen zu Q( √ 8)[END_REF]. In fact in these two cases the Hilbert surface is rational (that is birationnal to P 2 ):

Proof. We do the proof for Γ [START_REF] Becker | On Gröbner bases under specialization[END_REF][START_REF] Birkenhake | Humbert surfaces and the Kummer plane[END_REF] as the other one is similar. Note that Γ(1)/ Γ(4) SL 2 (O K /4O K ). We have then that O K /4O K is isomorphic to

The cardinality of SL 2 (O K /4O K ) is then 48 2 , 3840 and 3072 respectively. Moreover, the index of the subgroup Γ(4) of Γ(2, 4) is 4 when D ≡ 1 mod 4 and 16 when D ≡ 2, 3 mod 4. As these two sets are normal subgroups of Γ(1), the third isomorphism theorem of groups gives us the desired results.

Proposition A. 18 (Besser,Runge). The number of Humbert surface components for Γ(2) and for Γ(2, 4) is respectively

Proof. The numbers for Γ(2) are due to Besser and the other ones to Runge. See [START_REF] Runge | Endomorphism rings of abelian surfaces and projective models of their moduli spaces[END_REF]Page 293]. An heuristic argument for Γ(2, 4) is that given This is easier to see via the modular interpretation. Let Γ = Γ(2) (respectively Γ(2, 4)). Then an element of Γ\H 2 corresponds to a principally polarized abelian surface with a symplectic basis of the 2-torsion (resp. a symmetric theta structure of level 2). The cover Γ\H 2 → Sp 4 (Z)\H 2 corresponds to forgetting this extra structure, and the fibers form a torsor under the isomorphisms of this extra structure, which are equal to Γ(2)\Sp 4 (Z) (resp. Γ(2, 4)\Sp 4 (Z)).

The same is true for the map

1 and the action of G\ Γ(1)∪ Γ(1)σ on the fibers, where G is Γ(2) (resp. Γ(2, 4)) when D ≡ 1 mod 4 and Γ(2)∪ Γ(2)σ (resp. Γ(2, 4) ∪ Γ(2, 4)σ) when D ≡ 2, 3 mod 4. Except that here the extra structure has to be compatible with the action of O K . (For instance a symmetric theta structure of level 2 is induced by a symplectic basis of the 2-torsion and a compatible symplectic decomposition of the 4-torsion into maximal isotropic subgroups. For this symmetric theta structure to be compatible with the action of O K , these maximal isotropic subgroups have to be stable under the action of O K ).

In particular on the Humbert component H G ∆ K , then the action of G\ Γ(1) ∪ Γ(1)σ permute the fibers. Since this quotient is isomorphic to Γ

We give the equations of the Humbert component corresponding to the image of φ 1,ω for Γ(2, 4) and D = 2, 3, 5

and similarly for Γ(2) and D = 2 only, as for D = 3 the equations are too big to be put in the paper,

These equations have been computed by Gruenewald in [START_REF] Gruenewald | Explicit algorithms for Humbert surfaces[END_REF], where he gives equations for components of Humbert surfaces for many discriminants and many model.

B Examples of Hilbert modular polynomials

We apply now the results of the Sections 3 and 4 on Gundlach invariants for D = 2, 5 and on pullbacks of theta functions. We also propose non symmetric invariants.

B.1 Modular polynomials in Gundlach invariants

Evaluation and inversion: We first illustrate Theorem 3.4 for the Gundlach invariants J 1 , J 2 defined for Q( √ 2) and Q( √ 5) in Definition A.7 and A.10. The only small difference is that we use the map φ defined in Equation ( 26) rather than the map φ 1,ω (Equation ( 6)) to map Hilbert matrices τ ∈ H 2 1 to Siegel matrices Ω ∈ H 2 . In this case we have already seen how to express the pullbacks of the Igusa invariants in terms of the Gundlach invariants in Section 2.3 (see Corollaries A.13 and A.15). The expression is easier than the method outlined in Theorem 3.4 because the Gundlach invariants are expressed in terms of symmetric Hilbert modular forms whose relation to the pullbacks of the Siegel modular form defining the Igusa invariants are very simple (see Theorems A.12 and A.14).

We outline the algorithm (Algorithm B.1) to find τ ∈ (SL

The first step can be done using Corollary A.13 or A.15. The second is explained in [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF][START_REF] Milio | A quasi-linear time algorithm for computing modular polynomials in dimension 2[END_REF] and can be done in Õ(N ). For the third step, remark that for

The second step provides Ω ∈ H 2 which is more precisely in the Humbert surface H 5 . Thus by Humbert's Lemma we know there exists a matrix γ ∈ Sp 4 (Z) such that Ω = γΩ =

satisfies Ω 1 +Ω 2 -Ω 3 = 0 (see Remark 2.5 for the computation of γ). We have then τ * = ((

3 for the notations). For D = 2, φ (τ ) satisfies Ω 1 + 2Ω 2 -Ω 3 = 0 and we can adapt the algorithm to find the matrix γ. Thus Corollary B.2. Given J 1 (τ ) and J 2 (τ ), where J 1 and J 2 are the Gundlach invariants for D = 2 or 5 evaluated at some τ ∈ H 2 1 , then we can find an approximation of τ ∈ (SL

For the evaluation of the Gundlach invariants, using their definition as Fourier series would not give a good enough complexity. Instead Theorem 3.4 suggests to express J 1 and J 2 in term of the bk . Here, since the Gundlach invariants are invariants for the full modular group SL 2 (O K ⊕ ∂ K ), we can also express them directly in terms of the (pullbacks of the) Igusa invariants j 1 , j 2 , j 3 . Rather than doing an interpolation using Section 3.3, Equations (A.13) and (A.15) expressing the Igusa invariants in term of the Gundlach invariants are sufficiently simple to be inverted using a Gröbner basis.

In the case D = 5 we have found:

2 -42/5φ * j 7 2 φ * j 3 -7776/5φ * j 7 2 + 117/5φ * j 6 2 φ * j 2 3 -108/5φ * j 5 2 φ * j 3 3 ); In the case D = 2, the equations are too large to be included in the paper. But see https://members.loria.fr/EMilio/ .

We then have the following algorithm:

Data: τ ∈ H 2 1 and a working precision N Result: J 1 (τ ) and J 2 (τ ) at precision N 1 Compute Ω = φ (τ ); 2 Compute j 1 (Ω), j 2 (Ω) and j 3 (Ω) at precision N ; 3 Deduce J 1 (τ ) and J 2 (τ ) from the Igusa invariants;

For the first step we only have to use the definition of φ . For the second, we refer to [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]. The evaluation of the Igusa invariants can be done in Õ(N ) by [START_REF] Enge | Computing class polynomials for abelian surfaces[END_REF]. For the third, we use the equations above.

Corollary B.4. We can evaluate the Gundlach invariants J 1 (τ ) and J 2 (τ ) at precision N for D = 2 or 5 at any point τ ∈ H 2 1 with a complexity in Õ(N ) time.

Modular polynomials: Since we only have two invariants, this simplifies the definition of the modular polynomials:

Proposition B.5. Let D = 2 or 5 and be a prime number. If

. These polynomials depend only on and can be computed in time quasilinear in their size.

Proof. This is a corollary of Theorem 4.15. These polynomials depend only on as Q( √ D) for D = 2 or 5 has a fundamental unit of norm -1 (see the discussion in Remark 4.18). By construction, for any z ∈ H 2 1 , the modular polynomials satisfy Φ β (X, J 1 (z), J 2 (z)) = 0 when X is the evaluation of J 1 in one of the β-or β-isogenous points z . Then J

, where Φ β is the derivative of Φ β with respect to the variable X. Thus, given J 1 (z) and J 2 (z), the β-modular polynomials allow one to compute all the Gundlach invariants at the isogenous points of z.

. They can be computed in time quasi-linear in their size.

Proof. This is a corollary of Theorem 4.15. The difference between the cases D ≡ 1 mod 4 and D ≡ 2, 3 mod 4 comes from Equations ( 8) and ( 9): in the first case, by Proposition 2.8, the map Γ(2, 4)\H 2 1 → Sp 4 (Z)\H 2 is injective while in the second it is the map ( Γ(2, 4)∪ Γ(2, 4)σ)\H 2 1 → Sp 4 (Z)\H 2 which is injective. The coefficients of the Fourier series of the bi are in Q because it is the case of the Hilbert theta series (see [START_REF] Lauter | Hilbert theta series and invariants of genus 2 curves[END_REF]).

Note that there is three polynomials so that given b1 , b2 and b3 , one can obtain the values bγ 1,β , bγ 2,β and bγ 3,β for any γ ∈ C β . If D ≡ 1 mod 4 we are in the non symmetric case, so we compute non symmetric modular polynomials.

Remark B.11. When D = 2, Equation ( 28) says that we have to consider only two modular functions as b1 is determined by b2 and b3 . In particular the corresponding Humbert component is a rational surface.

β-modular polynomials: As Φ β is a minimal polynomial, it is the unique irreducible and monic polynomial which satisfies, for any τ ∈ H 2 1 , Φ β ( b1,β (τ ), b1 (τ ), b2 (τ ), b3 (τ )) = 0. We can look at what happen on σ(τ ). The matrix M σ of Equation ( 8 However ( bσ 1 , bσ 2 , bσ 3 ) = ( b3 , b2 , b1 ) if D ≡ 1 mod 4. The irreducible and monic polynomial Φ β ( bσ 1,β , bσ 1 , bσ 2 , bσ 3 ) has the same roots as Φ β ( b1,β , b1 , b2 , b3 ) and thus by unicity, these polynomials have to be equals. Thus, if D ≡ 1 mod 4, Φ β ( b3,β , b3 , b2 , b1 ) = Φ β ( b1,β , b1 , b2 , b3 ) and it is possible to obtain the value b3,β (τ ) for any τ ∈ H 2 1 using the β-modular polynomials. We have then, still acting by σ, b2,β (τ ) = Ψ 2,β ( b3,β (τ ), b3 (τ ), b2 (τ ), b1 (τ ))/Φ β ( b3,β (τ ), b3 (τ ), b2 (τ ), b1 (τ )) and b1,β (τ ) = Ψ 3,β ( b3,β (τ ), b3 (τ ), b2 (τ ), b1 (τ ))/Φ β ( b3,β (τ ), b3 (τ ), b2 (τ ), b1 (τ )).

We conclude that once we have the β-modular polynomials, we get the β-modular polynomials for free.

Changing β by a unit: Note that in the case where two pairs (β, β) and (β , β ) of totally positive elements, whose product is , differ by an even factor of (this always happens when has norm -1), we have that β = 2n β = n 0 0 -n β. Thus for any τ ∈ H 2 1 , if we compute bi,β (τ ), for i = 1, 2, 3, from bi (τ ) and using the β-modular polynomials, then we have bi,β (τ ) = bi -n 0 0 n 1 β τ and knowing how the matrix -n 0 0 n acts on the bi,β , we can compute the bi,β from the bi,β . In this case, it is useless to compute the β -modular polynomials.

Example B.12. When D = 2, 5 or 13, a fundamental unit > 1 has norm -1.

• If D = 2, we have that ( b1, 2 , b2, 2 , b3, 2 ) = ( b1 , b3 , b2 );

• If D = 5, we have that ( b1, 2 , b2, 2 , b3, 2 ) = ( b3 , b1 , b2 );

• If D = 13, we have that ( b1, 2 , b2, 2 , b3, 2 ) = ( b2 , b3 , b1 ).

When the norm of > 0 is 1, then if = ββ, we also have = β β , where β = β. The multiplication by does not come from the action of a matrix and the previous argument does not work.

Symmetries:

We can proceed in the same way with matrices γ ∈ Γ(2, 4)\ Γ(1) having special properties. If γ permutes the bi and the bi,β , this says that there are symmetries in the modular polynomials. In particular, if γ satisfies ( bγ

and consequently that Ψ 2,β (X, b1 , b3 , b2 ) = Ψ 3,β (X, b1 , b2 , b3 ) so that we only need to compute the two first β-modular polynomials, as the third one is deduced from the second one. For example, this happens for D = 6, = 73, β = 13 -4 √ 6 and for D = 10, = 41, β = 9 -2 √ 10. Moreover, if γ satisfies bγ k = i α k bk and bγ k,β = i β k bk,β , for k = 1, 2, 3 and α k , β k ∈ {0, 1, 2, 3} (i is the imaginary unit), then the exponents of the bk at each coefficient of the modular polynomials satisfy some relations modulo 4. As we compute the modular polynomials by evaluation/interpolation, this can be used to decrease the number of evaluations.

The existence of these matrices depend on D and β. They can be searched before the computation of the polynomials. We give some examples of relations between the exponents in Section 5 (see Equation [START_REF] Galbraith | Extending the GHS Weil descent attack[END_REF]). Similar arguments have already been used in [51, Sections 5.2 and 5.3] for the computation of -modular polynomials.

Denominator: Let L β be the locus of the principally polarized abelian surfaces z modulo Γ(2, 4) with real multiplication by O K for which z, or σ(z) in the case D ≡ 2, 3 mod 4, is β-isogenous to z such that φ 1,ω (z ) is isogenous to a product of elliptic curves by the 2-isogeny φ 1,ω (z ) → φ 1,ω (z )/2 and such that θ 0 (φ 1,ω (z )/2) = 0. Proof. Let z ∈ L β and let c i be a coefficient of Φ β . Then there is some γ ∈ ( Γ(2, 4) ∩ Γ0 (β))\ Γ(2, 4) such that bγ 1,β , or bγσ 1,β if D ≡ 2, 3 mod 4, is infinite. Indeed, recall that b i = θ i θ 0 (Ω/2) and that by [16, Proposition 6.5 and Corollary 6.1], exactly one theta constant vanishes at Ω if and only if Ω is isomorphic to a product of elliptic curves. We conclude using the same arguments as in the proof of Theorem B.6.

The reason for which we have introduced modular polynomials with the bi invariants was to obtain smaller polynomials compared to the ones with the Gundlach invariants or with the pullbacks of the Igusa invariants. But by Theorem B.9, the β-modular polynomials are not defined for all splitting in totally positives factors. We have two ways to deal with this problem, as explained in Remark 4.19. The first one consists in finding a subset of Γ(2, 4) for which bi,β is invariant (we are in the case D ≡ 2, 3 mod 4). For example, the group Γ defined as Γ [START_REF] Becker | On Gröbner bases under specialization[END_REF][START_REF] Birkenhake | Humbert surfaces and the Kummer plane[END_REF] in the case D ≡ 1 mod 4 (see Equations ( 16) and ( 17)). This subgroup is of index 4 in Γ(2, 4) and we consider the quotient ( Γ ∩ Γ0 (β))\ Γ(2, 4), containing 4( + 1) classes, to define our polynomials. The second one consists in taking other invariants, in particular the Rosenhain invariants ri = φ * 1,ω r i . We have already seen that they are generators for the field of Hilbert modular functions invariants by Γ(2) (see Theorem 2.16) and ri,β for i = 1, 2, 3 is always invariant by Γ(2) ∩ Γ0 (β). All the results of this section can be adapted to these invariants.