Simon Abelard

Pierrick Gaudry

Pierre-Jean Spaenlehauer

Improved Complexity Bounds for Counting Points on Hyperelliptic Curves

We present a probabilistic Las Vegas algorithm for computing the local zeta function of a hyperelliptic curve of genus g defined over Fq. It is based on the approaches by Schoof and Pila combined with a modelling of the -torsion by structured polynomial systems. Our main result improves on previously known complexity bounds by showing that there exists a constant c > 0 such that, for any fixed g, this algorithm has expected time and space complexity O((log q) cg ) as q grows and the characteristic is large enough.

Introduction

Since the discovery of Schoof's algorithm [START_REF] Schoof | Elliptic curves over finite fields and the computation of square roots mod p[END_REF], the complexity of counting points on curves and Abelian varieties defined over finite fields has attracted a lot of attention due to its numerous applications in cryptology, number theory and algebraic geometry. In this paper, we investigate the complexity of computing the local zeta function of hyperelliptic curves of fixed large genus. We propose a probabilistic algorithm which relies on the same foundations as Schoof's [START_REF] Schoof | Elliptic curves over finite fields and the computation of square roots mod p[END_REF] and Pila's algorithms [START_REF] Pila | Frobenius maps of abelian varieties and finding roots of unity in finite fields[END_REF].

When the characteristic of the base field F q is small, Kedlaya's and Satoh's approaches [START_REF] Kedlaya | Counting points on hyperelliptic curves using Monsky-Washnitzer cohomology[END_REF][START_REF] Satoh | The canonical lift of an ordinary elliptic curve over a finite field and its point counting[END_REF] and their variants compute very efficiently the number of rational points of Jacobians of hyperelliptic curves. We can also mention Lauder-Wan's [START_REF] Lauder | Counting points on varieties over finite fields of small characteristic[END_REF] and Lauder's [START_REF] Lauder | Deformation theory and the computation of zeta functions[END_REF] methods that can handle very general varieties. The current best algorithms in this family for rather general curves are by Tuitman [START_REF] Tuitman | Counting points on curves using a map to P 1[END_REF][START_REF] Tuitman | Counting points on curves using a map to P 1 II[END_REF]. However, the complexities of these p-adic algorithms are exponential in log(p), where p is the characteristic of the base field. This dependency can be made as low as √ p, thanks to the work of Harvey [START_REF] Harvey | Computing zeta functions of arithmetic schemes[END_REF]. Another line of research aims at taking profit of extra structure of the curve, assuming that this structure is known in advance and described in a convenient way. The most popular case is the Complex Multiplication method [START_REF] Atkin | Finding suitable curves for the elliptic curve method of factorization[END_REF], and in [START_REF] Gaudry | Counting points on genus 2 curves with real multiplication[END_REF] it is shown how to exploit real multiplication for counting points on genus 2 curves. When there is no such explicitly known additional structure and the characteristic of the base field is large, Schoof-Pila's -adic algorithms are the main tools for counting points. These -adic methods were introduced for elliptic curves in [START_REF] Schoof | Elliptic curves over finite fields and the computation of square roots mod p[END_REF], and later extended to Abelian varieties in [START_REF] Pila | Frobenius maps of abelian varieties and finding roots of unity in finite fields[END_REF]. In particular, Pila showed that the local zeta function of a gdimensional Abelian variety can be computed within O(log(q) ∆ ) operations, where ∆ and the constant in the O() are functions of g (but they do not depend on q). This complexity result requires some assumptions on the presentation of the Abelian variety which are satisfied by Jacobians of hyperelliptic curves given via a Weierstrass form [START_REF] Pila | Counting points on curves over families in polynomial time[END_REF]. Complexity improvements were obtained in [START_REF] Huang | Counting points on curves over finite fields[END_REF] and [START_REF] Adleman | Counting points on curves and Abelian varieties over finite fields[END_REF]. The latter article gives a deterministic algorithm for counting points on hyperelliptic curves with complexity (log q) O(g 2 log g) . Pila's algorithm and its variants may differ from Schoof's algorithm when specialized to the case of elliptic curves, but they are nonetheless related because they all rely on computing the characteristic polynomial χ of the Frobenius endomorphism ϕ modulo a prime number for sufficiently many such primes to deduce the numerator of the local zeta function of the curve (which is in fact the reciprocal polynomial of χ).

More precisely, let C be a hyperelliptic curve of genus g and J be its Jacobian. When is a prime different from the characteristic of the base field, the -torsion group J[ ] is isomorphic to (Z / Z) 2g and the characteristic polynomial of the restriction of ϕ on J[ ] is exactly χ mod . Furthermore, χ(1) = #J(F q ). The principle of Schoof-Pila's algorithm is to pick elements D in J[ ] and to find conditions on the coefficients s 0 , . . . , s 2g-1 mod such that ϕ 2g (D)+s 2g-1 ϕ 2g-1 (D)+• • •+s 0 D is equal to 0 in J[ ]. By testing all the tuples (s 0 , . . . , s 2g-1 ) up to the symmetries coming from the functional equation of χ (and possibly many D), the number of possibilities for χ mod is reduced until only one remains. The numerator of the zeta function is then obtained by repeating this procedure for many and by using Weil's conjectures to bound the absolute value of the coefficients.

For such a strategy, it is of the utmost importance that we get a description of the -torsion for which computations are reasonably easy to perform. In the elliptic case, computations in the -torsion subgroup are achieved by computing in the ring F q [X]/ψ (X) where ψ is the -division polynomial, which has degree O( 2 ). The dominant part of the complexity is the computation of ϕ 2 (X) in this quotient ring. In the genus 2 case, the bottleneck of the algorithm is no longer the computation of the powers of ϕ but that of a convenient algebraic representation of the -torsion [START_REF] Gaudry | Genus 2 point counting over prime fields[END_REF]. This appears to be also the case for g > 2. In order to reach the desired complexity, our main task is to compute such a representation efficiently. This is the central part of the proof of the complexity bound, and it is obtained by combining a special modelling of the -torsion with the geometric resolution algorithm [START_REF] Giusti | A Gröbner free alternative for polynomial system solving[END_REF], and by using multi-homogeneous Bézout bounds. More precisely, we show how to construct a polynomial system whose solutions are the -torsion points. This system involves two sets of variables, the first containing a small number O(g) of variables, each of them occurring with a degree that is polynomial in , and the second set containing many more variables but all of them occur with a degree that can be bounded independently of . This bi-homogeneous structure is the key to obtain a complexity bound that is better than for an unstructured system with the same number of variables and the same degree.

Another important ingredient in the proof of our main result is the extension of degree bounds for the coefficients of Cantor's analogue to division polynomials [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF]. Indeed, these polynomials are involved in the modelling of the -torsion and the degrees of their coefficients have a direct impact on the complexity of solving the polynomial system representing the -torsion. We finally mention that our result is of a purely theoretical nature. In the case of genus 2 and 3, the geometric resolution algorithm is at best quadratic in the degree of the -torsion ideal, which brings no improvement over a more direct study of the polynomial systems describing the -torsion. And for curves of larger (fixed) genus, we are still far from a situation where practical experiments could be run.

Organization of the paper. Section 2 describes a general algorithm for point-counting on Abelian varieties along with its complexity, assuming that the -torsion can be efficiently computed. Section 3 establishes the complexity result for multi-homogeneous polynomial systems that is required to obtain our claimed complexity bound. Section 4 contains the modelling of the -torsion under some mild assumptions on its structure. Finally, Section 5 describes the complete modelling of the -torsion, which is faithful even if the assumptions required in Section 4 are not satisfied.

Overview of the main result

Our main result is a probabilistic algorithm and a complexity bound for solving the following problem.

Computing local zeta functions of hyperelliptic curves. Given an odd prime power q, a positive integer g and a squarefree univariate polynomial f ∈ F q [X] of degree 2 g + 1, let C be the hyperelliptic curve with Weierstrass form Y 2 = f (X). Compute the numerator P C ∈ Z[T ] of the local zeta function of C:

Z(C /F q , T ) = exp ∞ i=1 # C(F q i ) • T i i = P C (T ) (1 -T )(1 -qT ) .
The special form of the denominator of the local zeta function is a consequence of Weil's conjectures. We refer to [18, Ch. XI, Thm. 5.2] for more details. Throughout the paper, we shall assume that the characteristic of F q is sufficiently large compared to log q. This assumption is required by a variant of Bertini's theorem (Proposition 4).

Our main result is as follows.

Theorem 1. There exists an explicitly computable constant c such that for all genus g, there exists an integer q 0 (g) such that for all prime power q = p n larger than q 0 (g) with p ≥ (log q) cg and for all hyperelliptic curves C of genus g defined over F q , the numerator P C of the local zeta function of C can be computed with a probabilistic algorithm in expected time bounded by

(log q) cg .
This complexity result is summarized by the notation O g ((log q) O(g) ), keeping in mind that g is fixed and q grows to infinity. Indeed, such a complexity statement can hide any factor that depends only on g: a running time in f (g)(log q) cg can be transformed into (log q) c g by taking a value c larger than c and adjusting q 0 (g), so that |f (g)| ≤ (log q 0 (g)) (c -c)g .

A typical example used in this article is the multiplication of two polynomials of degree d = (log q) O(g) . Using FFT-based techniques, this can be done in O(d) operations, which can be rewritten as (log q) O(g) (log((log q) O(g) )) k for some constant k and is therefore again in O g (log(q) O(g) ). Here the function f (g) that has been hidden in the operation is polynomial in g, but we will have cases where it is a combinatorial factor that grows very quickly with g and we make no effort to optimize it.

A classical geometrical object associated to a genus g curve is its Jacobian variety. Over the algebraic closure of F q , it can be described as the multiset of at most g points of the curve and it is endowed with an Abelian group structure (it is isomorphic to the degree-0 subgroup of the Picard group of the curve). The Frobenius map acts in a natural way on this Jacobian and it is compatible with its Z-module structure.

Throughout this paper, C is a hyperelliptic curve defined over F q with at least one rational Weierstrass point. Hence C admits a Weierstrass model y 2 = f (x), where f is a squarefree monic polynomial of degree 2g + 1. If C does not have any rational Weierstrass point, then we can extend the base field so that there exists a rational Weierstrass point that we send to infinity. The degree of the extension does not depend on q (it is at most linear in g), so that this will not affect our complexity result.

For practical computations, we need a coordinate system to represent points on the Jacobian of C: they shall be encoded via their Mumford representation using 2g coordinates. The group law on points in the Jacobian can be performed with Cantor's algorithm [START_REF] Cantor | Computing in the Jacobian of a hyperelliptic curve[END_REF] which operates with elements in Mumford representation at a cost of O(g) base field operations.

The algorithm that allows to prove the theorem is essentially the same as the one proposed by Pila for Abelian varieties, which is itself inspired by Schoof's algorithm for counting points on elliptic curves. This algorithm relies on a few classical results for curves defined over finite fields:

• The numerator P C of the local zeta function is the reciprocal of the characteristic polynomial of the Frobenius morphism on the Jacobian variety J of C [18, Thm. 5.2];

• For prime numbers not dividing q, the -torsion J[ ] of the Jacobian variety is isomorphic (as an Abelian group) to (Z / Z) 2g [19, Sec. II.6, Prop. page 64], [START_REF] Cohen | Handbook of elliptic and hyperelliptic curve cryptography[END_REF]Thm. 4.73]; Therefore P C mod is the reciprocal of the characteristic polynomial of the Frobenius seen as an endomorphism of J[ ] ∼ = (Z / Z) 2g ;

• The Weil conjectures imply that P C has the following form over the complex numbers:

P C (T ) = 2g i=1 (1-u i T ) with |u i | = q 1/2 [18,
Ch. VIII, Thm. 6.1]. Moreover, if a 0 , . . . , a 2g denote the coefficients of P C , the functional equation implies that a 2g-i = q g-i a i . Consequently, the absolute value of the coefficients of P C ∈ Z[T ] are bounded by 2g g q g .

Data: q ∈ Z >0 a prime power, and f ∈ F q [X] a monic squarefree univariate polynomial. Result: The characteristic polynomial χ ∈ Z[T ] of the Frobenius endomorphism on the Jacobian J of the hyperelliptic curve defined over

F q with Weierstrass form Y 2 = f (X). ← 1; R ← 1; while R ≤ 2 2g g q g + 1 do ←NextPrime( ); if divides q then ←NextPrime( ); end Compute a description of J[ ];
Compute a 2g × 2g matrix F with coefficients in Z / Z representing the action of the Frobenius on J[ ] ∼ = (Z / Z) 2g ; Compute the characteristic polynomial χ mod of the matrix F ; R ← R • ; end Reconstruct χ using the Chinese Remainder Theorem. Algorithm 1: A bird's eye view of Pila's point counting algorithm for hyperelliptic curves.

Pila's algorithm reconstructs the numerator of the local zeta function of C by computing the action of the Frobenius on the -torsion for sufficiently-many prime numbers and by using the Chinese Remainder Theorem. A bird's eye view of this algorithm is given in Algorithm 1. The main difficulty resides in the step where one computes an explicit description of J[ ]. Since J[ ] is a 0-dimensional variety of degree 2g , what we will compute is a geometric resolution of the corresponding radical ideal, that is a univariate squarefree polynomial F (T ), together with 2g coordinate polynomials γ i (T ), such that the coordinates of the -torsion elements are the evaluations of the vector (γ 1 (T ), . . . , γ 2g (T )) at the roots of F .

To be more precise, the Mumford coordinates are in fact a set of g affine systems of coordinates, each corresponding to a different weight of the represented divisor (the definition is recalled in Section 4). The variety J[ ] will accordingly be represented by a set of g geometric resolutions, each encoding -torsion divisors of a given weight w ∈ [1, g]. Generically, we expect that all the elements in J[ ] have weight g, except for the neutral element which has weight 0. Most of the article is dedicated to computing efficiently this representation for J[ ]. The cornerstone of the proof of Theorem 1 relies on the following statement. Proposition 2. Let C be a hyperelliptic curve of genus g over F q with Weierstrass form Y 2 = f (X) (f monic of degree 2g + 1) and J be its Jacobian variety. Let > g be a prime not dividing q. Assuming that the characteristic of F q is sufficiently large as in Theorem 1, there is a Las Vegas probabilistic algorithm which takes as input q, , f and which computes geometric resolutions for the varieties {J w [ ]} w∈ [1,g] of -torsion points of weight w in the Jacobian variety. This algorithm can be implemented by a Turing machine with space and expected time O g ( log q) O(g) .

Assuming this complexity bound, performing a complexity analysis as done in [START_REF] Pila | Frobenius maps of abelian varieties and finding roots of unity in finite fields[END_REF] leads to a complexity bound for Algorithm 1 that corresponds to Theorem 1. We recall it here for completeness, with some simplifications due to the fact that we consider a probabilistic algorithm, so we can factor polynomials using Cantor-Zassenhaus' algorithm.

Proof of Theorem 1 assuming Proposition 2. By Weil's bounds, the absolute values of the coefficients of the characteristic polynomial χ are bounded by 2g g q g . Therefore at the end of the loop of Algorithm 1, these coefficients are completely determined by their values modulo all the primes that have been explored. It follows from [START_REF] Tenenbaum | Introduction to analytic and probabilistic number theory[END_REF]Cor. 10.1] that the largest in the loop is at most linear in g log q. From this and Proposition 2, computing the description of J[ ] as a union of geometric resolutions for all the J w [ ] can be achieved within expected complexity O g (log q) O(g) .

Factoring the univariate polynomials involved in the geometric resolutions can be done within the same time bound O g (log q) O(g) , since the sum of their degrees is 2g and factoring polynomials in finite fields can be done in time linear in log(q) and quasi-quadratic in the degree [START_REF] Zur Gathen | Modern computer algebra[END_REF]Thm. 14.14]. Therefore, it is possible to construct a Mumford representation for eachtorsion divisor within the same complexity, each of them possibly defined over a different extension of F q . In fact, due to the rationality of the group law that acts on J[ ], one of these extensions of F q contains all the others.

Using elementary linear algebra for the Frobenius endomorphism ϕ acting on J[ ] (seen as an F -vector space), we can deduce χ = χ mod . We first compute a basis of J[ ] by brute force and a dictionary of how all elements decompose on it. Then, the action of ϕ on the basis elements can be computed and the result is a matrix whose characteristic polynomial is χ . All of this fits in the O g ((log q) O(g) ) complexity bound. The loop is repeated O g (log q) times, and this additional factor does not affect the overall complexity.

Polynomial systems

This section is devoted to describing tools that we will use to estimate the complexity of computing a convenient representation of the -torsion of the Jacobian of hyperelliptic curves.

We start by fixing some notation. In the sequel, F q denotes the algebraic closure of F q . For an ideal I ⊂ F q [X 1 , . . . , X n ], we call dimension of I and note dim(I) the Krull dimension of the quotient ring F q [X 1 , . . . , X n ]/I. Moreover, by identifying a point (λ 0 , . . . , λ n ) ∈ F q n+1 with the polynomial λ

0 + λ 1 X 1 + • • • + λ n X n ∈ F q [X 1 , . . . , X n ], there is a dense Zariski open subset O ⊂ (F q n+1 ) dim(I) such that for any ( 1 , . . . , dim(I) ) ∈ O, the alge- bra F q [X 1 , . . . , X n ]/(I + 1 , . . . , dim(I)
) is a finite dimensional F q -vector space of constant dimension, which is called the degree of

I. A sequence (f 1 , . . . , f i ) ∈ F q [X 1 , . . . , X n ] i is regular if f 1 , . . . , f i = F q [X 1 , . . . , X n ] and for any j ∈ [2, i], f j does not divide zero in F q [X 1 , . . . , X n ]/ f 1 , . . . , f j-1 . The sequence (f 1 , . . . , f i ) is reduced if every intermediate ideal f 1 , . . . , f j with j ∈ [1, i] is radical.
Geometric resolutions. For describing 0-dimensional (i.e. finite) sets V ⊂ F q n where V is defined over F q , we use a data structure called a geometric resolution of V . The terminology here is borrowed from [START_REF] Cafure | Fast computation of a rational point of a variety over a finite field[END_REF], see also [START_REF] Giusti | A Gröbner free alternative for polynomial system solving[END_REF]. An F q e -geometric resolution of V is a tuple (( 1 , . . . , n ), Q, (Q 1 , . . . , Q n )) where:

• The vector ( 1 , . . . , n ) ∈ F n q e is such that the linear form

: F q n → F q (x 1 , . . . , x n ) → n i=1 i x i
takes distinct values at all points in V . The linear form is called the primitive element of the geometric resolution;

• The polynomial Q ∈ F q e [T ] equals x∈V (T -(x)); • The polynomials Q 1 , . . . , Q n ∈ F q e [T ] parametrize V by the roots of the polynomial Q, i.e. V = {(Q 1 (t), . . . , Q n (t)) | t ∈ F q , Q(t) = 0}.
We note that our definition is slightly simpler than the one in [3, Sec. 2.1] because we restrict ourselves to the 0-dimensional case in this paper (in [3, Sec. 2.1], the definition is also valid for equidimensional varieties with positive dimension).

In the following statement, if f is a polynomial in a ring

F q [X 1 , . . . , X nx , Y 1 , . . . , Y ny ], then we let deg x (f ) (resp. deg y (f )) denote the degree of f (X 1 , . . . , X nx , y 1 , . . . , y ny ) ∈ F q [X 1 , . . . , X nx ] (resp. f (x 1 , . . . , x nx , Y 1 , . . . , Y ny ) ∈ F q [Y 1 , . . . , Y ny ])
, where y 1 , . . . , y ny (resp. x 1 , . . . , x nx ) are generic values in F q .

The following proposition is a cornerstone of our complexity result for computing thetorsion of the Jacobian of a hyperelliptic curve. The statement and its proof combine three main ingredients: (1) the geometric resolution algorithm [START_REF] Giusti | A Gröbner free alternative for polynomial system solving[END_REF] and its version for finite fields [START_REF] Cafure | Fast computation of a rational point of a variety over a finite field[END_REF], which are methods for solving polynomial systems whose complexity depends mainly on geometric degrees; (2) the multi-homogeneous Bézout bound which allows us to control the geometric degrees by separating the variables in our modelling in two blocks, where the block supporting most of the degrees has small cardinality; (3) a variant of Bertini's theorem to process our polynomial system into a reduced regular sequence which is a valid input for the geometric resolution algorithm.

As we shall see in the next sections, our polynomial system modelling the -torsion will have two blocks of variables. The first block occurs with large degree O (1) but it has very small cardinality in O(g). The second block has a larger cardinality, but the degrees of the equations with respect to this block do not depend on , but only on g. Taking this bihomogeneous structure into account is crucial to reach our claimed complexity bound. The following proposition provides a bound on the complexity of solving polynomial systems having this structure, and the sequel of this section is dedicated to its proof. Proposition 3. There exists a probabilistic Turing machine T which takes as input polynomial systems with coefficients in a finite field F q and which satisfies the following property. For any function h : Z >0 → Z >0 , for any positive number C > 0 and for any ε > 0, there exists a function ν : Z >0 → Z >0 and a positive number D > 0 such that for all positive integers g, , n x , n y , d x , d y , m > 0 such that n x < C g, n y < h(g), d x < h(g) C , d y < h(g), m < h(g), for any prime power q such that the prime number p dividing q satisfies 2 nx+ny d nx x d ny y < p, and for any polynomial system f 1 , . . . ,

f m ∈ F q [X 1 , . . . , X nx , Y 1 , . . . , Y ny ] such that • for all i ∈ [1, m], deg x (f i ) ≤ d x and deg y (f i ) ≤ d y , • the ideal I = f 1 , . . . , f m has dimension 0 and is radical, the Turing machine T with input f 1 , . . . , f m returns an F q ν(g) log -geometric resolution of the variety {x ∈ F q | f 1 (x) = • • • = f m (x) = 0}
with probability at least 5/6, using space and time bounded above by ν(g) D g (log q) 2+ε .

Proof. Postponed to the end of this section.

Since the geometric resolution requires its input to be a reduced regular sequence, we first need to ensure that we can construct such a sequence from our input system. A classical way to achieve this is to replace the input system by a generic linear combination of the polynomials. If the ideal generated by the input system is 0-dimensional and radical, then a variant of Bertini's theorem ensures that the obtained sequence is regular and reduced. Proposition 4. [START_REF] Sommese | The numerical solution of systems of polynomials arising in engineering and science[END_REF]Thm. A.8.7]

Let (f 1 , . . . , f m ) ∈ F q [X 1 , . . . , X nx , Y 1 , . . . , Y ny ] m be polyno- mials such that the ideal I = f 1 , . . . , f m has dimension 0 and is radical. Let d x , d y be two integers such that deg x (f i ) ≤ d x , deg y (f i ) ≤ d y for all i ∈ [1, m].
Let p be the characteristic of F q , and assume that 2 nx+ny d nx x d ny y < p. For M an (n x + n y ) × m matrix with entries in F q , let (f

(M ) 1 , . . . , f (M ) nx+ny ) ∈ F q [X 1 , . . . , X nx , Y 1 , . . . , Y ny ] nx+ny be defined as       f (M ) 1 f (M ) 2 . . . f (M ) nx+ny       = M •      f 1 f 2 . . . f m      .
Then there exists a nonempty open subset O ⊂ F q (nx+ny)×m of the space of (n x + n y ) × m matrices such that for any M ∈ O, for any s ∈ [1, n x + n y ], and at any point

(x, y) ∈ F q nx+ny such that f (M ) 1 (x, y) = • • • = f (M ) s (x, y) = 0, the derivatives Df (M ) 1 (x, y), . . . , Df (M ) s
(x, y) are linearly independent over F q . In particular, for any M ∈ O, the sequence (f

(M ) 1 , . . . , f (M ) nx+ny ) is reduced and regular.
Proof. This is a reformulation of [START_REF] Sommese | The numerical solution of systems of polynomials arising in engineering and science[END_REF]Thm. A.8.7] in the case of finite fields. In [START_REF] Sommese | The numerical solution of systems of polynomials arising in engineering and science[END_REF]Thm. A.8.7], this result is stated over the field C, but this statement holds true over any field k, provided that an extra separability assumption is satisfied. More precisely, set n = n x + n y and let

V s ⊂ k n × k n m be the variety of pairs ((x, y), M ) such that f (M ) 1 (x, y) = • • • = f (M ) s (x, y) = 0.
In this setting, the extra condition that is required for the proposition to hold is that the projection π of V s to k n m must be separable for all s ∈ [1, n] (this is always true in characteristic 0). We refer to [START_REF] Kleiman | Bertini and his two fundamental theorems[END_REF]Thm. 4.2] for more details on this separability argument. In our setting, the degree of a generic fiber of π is bounded by 2 n d nx x d ny y < p using the multi-homogeneous Bézout bound (see e.g. Proposition 8 below) and hence the separability condition is satisfied.

Since we are looking at polynomial systems over finite fields, we must estimate the size of the extension of the base field that is required to find with sufficiently large probability a matrix M such that f

(M ) 1 , . . . , f (M )
nx+ny is reduced and regular.

Lemma 5. Let (f 1 , . . . , f m ) ∈ F q [X 1 , . . . , X nx , Y 1 , . . . , Y
ny ] m be polynomials satisfying the assumptions of Proposition 4 and such that their total degree is bounded above by d ∈ Z ≥0 . Set n = n x + n y and e = (2n + 1) log q (d + 1) + log q [START_REF] Harvey | Computing zeta functions of arithmetic schemes[END_REF] .

If M is an n × m matrix with entries in F q e picked uniformly at random, then the probability that (f

(M ) 1 , . . . , f (M ) n
) is a reduced regular sequence is bounded below by 10/11.

Proof. Let Λ denote an n × m matrix with indeterminate entries

Λ =    λ 11 . . . λ 1m . . . . . . . . . λ n1 . . . λ nm    and let F 1 (Λ, X, Y ), . . . , F n (Λ, X, Y ) ∈ F q [X 1 , . . . , X nx , Y 1 , . . . , Y ny , λ 11 , . . . , λ nm ] be the poly- nomials defined as    F 1 (Λ, X, Y ) . . . F n (Λ, X, Y )    = Λ •    f 1 (X, Y ) . . . f m (X, Y )    .
For s ∈ [1, n], we consider the s × m matrix Λ (s) obtained by truncating Λ to its s first rows, a new set of variables {µ 1 , . . . , µ s-1 } and the following polynomial system:

F 1 (Λ (s) , X, Y ) = • • • = F s (Λ (s) , X, Y ) = 0 µ 1 • • • µ s-1 1 •        ∂F 1 ∂X 1 • • • ∂F 1 ∂X nx ∂F 1 ∂Y 1 • • • ∂F 1 ∂Y ny . . . . . . . . . . . . . . . . . . ∂F s ∂X 1 • • • ∂F s ∂X nx ∂F s ∂Y 1 • • • ∂F s ∂Y ny        = 0 • • • 0
This is a system of n + s polynomials of degree bounded above by d + 1 in n + s -1 + ms variables. By Bézout inequality (see e.g. [12, Thm. 1]), this system defines a variety V s which is either empty, or its degree is at most (d + 1) n+s . We remark that if V s is not empty, then it has dimension at least ms -1 since its vanishing ideal is generated by n + s elements. The Zariski closure of its projection W s to the space F q sm of matrices Λ (s) is either empty, the whole space or a proper sub-variety. By Proposition 4, it must be empty or a proper sub-variety. Next, we remark that the degree of the image of a variety by a linear projection cannot increase. Therefore, the sum of the degrees of the irreducible components of W s is also bounded by

(d + 1) n+s if W s = ∅.
In the sequel, we let h s (λ 11 , . . . , λ sm ) denote a polynomial vanishing on W s of degree bounded by (d + 1) n+s (we set h s (λ 11 , . . . ,

λ sm ) = 1 if W s = ∅).
The Schwarz-Zippel Lemma implies that the cardinality of the set

E =         M 11 • • • M 1m . . . . . . . . . M n1 • • • M nm    ∈ F nm q e | h 1 (M 11 , . . . , M 1m ) • • • h n (M 11 , . . . , M nm ) = 0     
is bounded above by q e /11, for the value of e given in the statement. The proof is concluded by noticing that for any M ∈ E, for any s ∈ [1, n], and for any (x, y) ∈ F q n such that f

(M ) 1 (x, y) = • • • = f (M ) s
(x, y) = 0 the derivatives Df M 1 (x, y), . . . , Df (M ) s (x, y) span the normal space at (x, y) to the variety associated with f

(M ) 1 , . . . , f (M ) s . Hence, f (M ) 1 , . . . , f (M ) n is a reduced regular sequence.
Once we have a reduced regular sequence, we can use [START_REF] Cafure | Fast computation of a rational point of a variety over a finite field[END_REF]Thm. 4.8] to solve the system. We note that in [START_REF] Cafure | Fast computation of a rational point of a variety over a finite field[END_REF] there is a general assumption that for all s ∈ [1, n] the intermediate ideals f

(M ) 1 , . . . , f (M ) s
define absolutely irreducible varieties. However, the proof of [START_REF] Cafure | Fast computation of a rational point of a variety over a finite field[END_REF]Thm. 4.8] does not require this assumption (this assumption is only required in algorithms for finding a rational point in [START_REF] Cafure | Fast computation of a rational point of a variety over a finite field[END_REF]Section 6]).

Next, we describe the data structures used in [START_REF] Cafure | Fast computation of a rational point of a variety over a finite field[END_REF] to represent polynomial systems. The algorithms take as input polynomials represented by division-free straight-line programs (DFSLP). A DFSLP defined over a field k is a sequence of polynomials

h 1 , h 2 , . . . , h ∈ k[X 1 , . . . , X n ] such that each polynomial h i is either a variable X t with t ∈ [1, n], an element in k, or h i = h j • h j ,
where j, j < i and • ∈ {+, -, ×} is an arithmetic operation. The time of a DFSLP is the total number of arithmetic operations, and its space is the minimal number of arithmetic registers required to evaluate it. A polynomial system f 1 , . . . , f m is said to be represented by a DFSLP

h 1 , . . . , h if {f 1 , . . . , f m } ⊂ {h 1 , . . . , h }. Theorem 6. [3, Thm. 4.8] Let f (M ) 1 , . . . , f (M ) n ∈ F q e [
X 1 , . . . , X n ] be a reduced regular sequence, where the polynomials are represented by a DFSLP with space S and time T . Set the following notation:

• The integer d is max i∈[1,n] (deg(f (M ) i
));

• For any real number x ≥ exp(1), U(x) = x(log x) 2 log log x;

• Let δ ∈ Z ≥0 be an integer larger than the degrees of the ideals f

(M ) 1 , f (M ) 1 , f (M ) 2 , . . . , f (M ) 1 , . . . , f (M ) n .
Assume further that q e ≥ 60 n 4 dδ 4 . There is a probabilistic Turing machine using space O((S + n + d)δ 2 log(q e δ)) and time O((nT + n 5 )U(δ)(U(dδ) + log(q e δ))U(log(q e δ))) which takes such polynomial systems as input and which outputs a F q e -geometric resolution of the algebraic set

{x ∈ F q e n | f (M ) 1 (x) = • • • = f (M ) n (x) = 0} with probability at least 11/12.
The next lemma is a first step for preparing our system in order to use Theorem 6 for bi-homogeneous systems: we need to estimate the size and space needed to represent a bihomogeneous system by a DFSLP. -1) additions.

Lemma 7. Let d x , d y ∈ Z >0 be two positive integers. A polynomial system f 1 , . . . , f m ∈ F q [X 1 , . . . , X nx , Y 1 , . . . , Y ny ] such that for all i ∈ [1, m], deg x (f i ) ≤ d x and deg y (f i ) ≤ d
The next ingredient in order to derive Proposition 3 from Theorem 6 is an upper bound on δ. This can be obtained via the multi-homogeneous Bézout bound.

Proposition 8. Let f 1 , . . . , f m be a regular sequence in F q [X 1 , . . . , X nx , Y 1 , . . . , Y ny ] and d x , d y ∈ Z ≥0 be such that for any i ∈ [1, m], deg x (f i ) ≤ d x and deg y (f i ) ≤ d y . Then the degree of the ideal f 1 , . . . , f m is at most j1+j2=m 0≤j1≤nx 0≤j2≤ny m j 1 d j1 x d j2 y . (1)
Moreover, this degree is bounded above by

2 nx+ny d nx x d ny y .
Proof. This is a direct consequence of [22, Prop. I.1] using, with the notation of [22, Prop.

I.1], k = 1, e = 0, P = m, D i,0 = d x , D i,1 = d y , n = n x , n 1 = n y . Note that [22, Prop. I.1]
is stated when the base field is C, but the proof works without any major modification when the base field is a finite field. The last sentence of the statement follows from the fact that the regularity assumption implies that m ≤ n x + n y , and hence the sum of the binomial coefficients is bounded above by 2 m ≤ 2 nx+ny .

We now have all the ingredients needed to prove Proposition 3.

Proof of Proposition 3. Set n = n x + n y . First, we note that if f 1 , . . . , f m ∈ F q [X 1 , . . . , X nx , Y 1 , . . . , Y ny ]
is represented by a straight-line program over F q with space S and time T , then for any e ∈ Z ≥0 and any m × n matrix M with entries in F q e , the sequence f

(M ) 1 , . . . , f (M ) n ∈ F q e [X 1 , . . . , X nx , Y 1 , . . . , Y ny ]
can be represented by a straight-line program over F q e with space S and time T , where S = O(S) and T = O(T + m n). We consider the probabilistic Turing machine which performs the following steps:

1. It chooses an m × n matrix uniformly at random with entries in F q e , with e = max (2n + 1) log q (d + 1) + log q (11) , log q (60 n 4 d δ) ,

where d = d x + d y = ( C + 1) h(g), n = n x + n y = Cg + h(g), δ = 2 n d nx x d ny y = (2h(g)) Cg+h(g) C 2 g . Using the inequalities n x < C g, n y < h(g), d x < h(g) C , d y < h(g), we get that e = O g (log q );
(M ) 1 (X) = • • • = f (M ) n
(X) = 0; By Theorem 6, it returns a geometric resolution (( 1 , . . . , n ), q(T ), (q 1 (T ), . . . , q n (T ))) provided that f (M ) 1 (X), . . . , f (M ) n (X) is a reduced regular sequence; 4. It computes λ(T ) = GCD(q(T ), f 1 (q 1 (T ), . . . , q n (T )), . . . , f m (q 1 (T ), . . . , q n (T ))); 5. It computes ν 1 (T ) = q 1 (T ) mod λ(T ), . . . , ν n (T ) = q n (T ) mod λ(T ) and returns the geometric resolution (( 1 , . . . , n ), λ(T ), (ν 1 (T ), . . . , ν n (T ))).

We start by showing that the output of this algorithm is indeed a geometric resolution of the algebraic set

V = {x ∈ F q n | f 1 (x) = • • • = f m (x) = 0}
, assuming that the probabilistic algorithm in Step 3 returns the correct result. Let W be the algebraic set

{x ∈ F q n | f (M ) 1 (x) = • • • = f (M ) m (x) = 0}. Since f (M ) 1 , . . . , f (M ) n ⊂ f 1 , . . . , f m , we have V ⊂ W .
By construction, the algebraic set defined by the geometric resolution (( 1 , . . . , n ), λ(T ), (ν 1 (T ), . . . , ν n (T ))) is precisely the subset of W where all polynomials f 1 , . . . , f m simultaneously vanish.

It remains to prove that this Turing machine runs within the desired complexity. Steps 1 and 2 require negligible time. Step 3 is done within space O((S + n + d)δ 2 log(q e δ)) and time O((nT + n 5 )δ(dδ + log(q e δ)) log(q e δ)) (Theorem 6), provided that δ is an upper bound on the degrees of the intermediate ideals.

Step 4 is done within space and time bounded by O(δ e log q(T + m)) by evaluating the SLP modulo q(T ) (whose degree is bounded by δ) and then by computing m GCD using a quasi-linear algorithm. Finally, Step 5 can be done within time and space O(δ e log q).

Then, Proposition 8 shows that δ is an upper bound on the degrees of the intermediate ideals. Using the facts that nx+dx dx

≤ (n x +d x ) nx = O g ( C 2 g

Computing generic -torsion points

Let C be a hyperelliptic curve of genus g over F q with Weierstrass form Y 2 = f (X) (f monic, squarefree, and deg(f ) = 2g + 1) and J be its Jacobian. Let > g be a prime not dividing q. In this section, we define a notion of genericity for -torsion elements in J and we show that a geometric resolution for the variety they form can be computed efficiently using the tools described in Section 3 by solving a polynomial system of g 2 + g equations in g 2 + g variables. Our starting point is the modelling of the -torsion sketched by Cantor in the point (5) of Section 9 of [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF]. This section and the next one that deals with the non-generic cases rely heavily on the Mumford representation, that we recall here, and refer to [START_REF] Galbraith | Mathematics of public key cryptography[END_REF] for more details.

Definition 9 (Mumford representation). Every element in J can be uniquely represented by a pair of polynomials u(X), v(X)

, where u is monic of degree w ≤ g, the polynomial v has degree less than w, and u divides v 2 -f . We call the integer w in this definition the weight of the element. An element in Mumford representation u(X), v(X) , corresponds to a divisor (called reduced divisor) of the form 1≤i≤w (P i -∞), where the P i = (x i , y i ) are the affine points of C such that u(x i ) = 0 and y i = v(x i ), with appropriate multiplicities. This implies in particular that if two P i 's share the same x-coordinate, then they are equal.

In what follows, we often also call Mumford representation a pair of polynomials where u is not monic. In that case, unicity of the representation is no longer guaranteed, but there is no ambiguity in the element of J represented this way.

In genus 1, the -torsion points are the points whose abscissae are the roots of thedivision polynomial, which has degree O( 2 ). For higher genera, Cantor [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF] described analogous polynomials δ and ε such that, for (x, y) a generic point of the curve and > g, we have

• ((x, y) -∞) = δ x -X 4y 2 , ε x -X 4y 2 .
Lemma 10. The polynomial δ (X) has degree g and its coefficients are polynomials in F q [x] of degree bounded by

1 3 g 3 + O g ( 2 )
. The polynomial ε (X)/y has degree less than g and its coefficients are rational fractions in F q (x). The degrees of the numerators and denominators of these coefficients are bounded by 2 3

g 3 + O g ( 2 )
. Furthermore, any root of a denominator is also a root of the leading coefficient of δ (X).

Proof. An exact formula for the degree of the leading coefficient of δ (X) was given by Cantor.

For the other coefficients and the claims about the roots of the denominator, the proof is postponed to Section 6. We remark that our bounds are not tight but they are sufficient for our purpose.

We shall prove that Lemma 10 is also valid for a non-prime integer . This will be useful in Section 5 where we handle non-generic situations.

Later on, we will need explicit names for these coefficients of δ and ε , so we define the univariate polynomials d i and e i (the notation does not show the dependence on for simplicity) such that, after clearing denominators we have:

δ x -X 4y 2 = g i=0 d i (x)X i , and ε x -X 4y 2 = y g-1 i=0 e i (x) e g (x) X i .

Definition 11.

In what follows, we shall say that an element of J is -generic if it has weight g and the corresponding reduced divisor g i=1 (P i -∞) satisfies the following two properties: • For any i, the u-coordinate of the divisor • (P i -∞) in Mumford form has degree g;

• For any i = j, the u-coordinates of the divisors • (P i -∞) and • (P j -∞) are coprime. This implies that all the P i are distinct. This also implies that if an affine point P occurs in the support of a • (P i -∞) then neither P nor -P appears in the support of another

• (P j -∞).

Proposition 12. For any ε > 0, there is a constant D such that for all prime > g coprime to the base field characteristic, there is a Monte Carlo algorithm which computes an F q e -geometric resolution of the sub-variety of J[ ] consisting of -generic -torsion elements, where e = O g (log ). The time and space complexities of this algorithm are bounded by O g ( Dg (log q) 2+ε ) and it returns the correct result with probability at least 5/6.

Proof. Let D = g i=1 (P i -∞) be an -generic divisor in J. We shall consider a system equivalent to • D = 0 but let us first introduce some notation. For each point P i = (x i , y i ) in the support of D, we denote u i , v i the Mumford form of • (P i -∞) and (α ij , β ij ) 1≤j≤g the coordinates of the g points in its support counted with multiplicities, which means that for any i the g roots of u i are exactly the α ij , and that for any j, β ij = v i (α ij ). Note that using the previous notation,

u i (X) = δ xi-X 4y 2 i and v i (X) = ε xi-X 4y 2 i .
We have • D = 0 if and only if the sum of the divisors g i=1 • (P i -∞) is a principal divisor. The only pole is at infinity, so this is equivalent to the existence of a non-zero function ϕ ∈ F q (C) of the form P (X) + Y Q(X) with P and Q two polynomials such that the g 2 points (α ij , β ij ) are the zeros of ϕ, with multiplicities. Since we want ϕ to have g 2 affine points of intersection with the curve C (once again, counted with multiplicities), the polynomial

Res Y (Y 2 -f, P + Y Q) = P 2 -f Q 2 must have degree g 2 which yields 2 deg(P ) ≤ g 2 and 2 deg(Q) ≤ g 2 -2g -1.
Exactly one of those two bounds is even (it depends on the parity of g), and for this particular bound, the inequality must be an equality, otherwise the degree of the resultant would not be g 2 . Since the function ϕ is defined up to a multiplicative constant, we can normalize it so that the polynomial P 2 + f Q 2 is monic, which is equivalent to enforce that either P or Q is monic depending on the parity of g.

For a fixed i ∈ [1, g], requiring the (α ij , β ij ) to be zeros of ϕ amounts to asking for the α ij to be roots of P (X) + Q(X)v i (X), with multiplicities. Since the α ij are by definition the roots of the u i , • D = 0 is equivalent to g congruence relations P + Qv i ≡ 0 mod u i which we can rephrase using Cantor's polynomials:

P (X) + ε x i -X 4y 2 i Q(X) ≡ 0 mod δ x i -X 4y 2 i . ( 2
)
Thus, for any -generic divisor, • D = 0 is equivalent to the existence of P and Q satisfying the above g congruence relations.

The variables are the coefficients of P and Q, as well as the x i and y i . With the degree conditions and the normalization, we have g 2 -g variables coming from P and Q. Adding the 2g variables x i and y i , we get a total of g 2 + g variables. Each one of the g congruence relations (2) amounts to g equations providing a total of g 2 conditions on the coefficients of P and Q. The fact that the (x i , y i ) are points of the curve yields the g additional equations y 2 i = f (x i ). Finally, we have to enforce the -genericity of the solutions, which can be done by requiring that i d g (x i ) i<j Res(u i , u j ) = 0. Therefore, we get a polynomial system with g 2 + g equations in g 2 + g variables, together with an inequality. We remark that in principle, the denominators e g (x i ) involved in ε would generate additional conditions, but by Lemma 10 this is already covered by the condition d g (x i ) = 0.

In order to apply Proposition 3, we now estimate the degrees to which the variables occur in the equations. We start with the equations coming from [START_REF] Atkin | Finding suitable curves for the elliptic curve method of factorization[END_REF]. Each congruence relation is obtained by reducing P + Qv i , which is a polynomial of degree O(g 2 ) in X, by u i which is of degree g. We can do it by repeatedly replacing X g by -j<g (d j (x i )/d g (x i ))X j , which we will have to do at most O(g 2 ) times. Since by Lemma 10 the d j have degree in O g ( 3 ) in x i the fully reduced polynomial will have coefficients that are fractions for which the degrees of the numerators and of the denominators are at most O g ( 3 ) in the x i variables. In these equations, the degree in the y i variables and in the variables for the coefficients of P and Q is 1. The degrees in x i and y i in the curve equations are 2g + 1 and 2 respectively.

It remains to study the degree of the inequality. Each resultant is the determinant of a 2g × 2g Sylvester matrix whose coefficients are the d i , which have degrees bounded by O g ( 3 ). Since for any i there are exactly g resultants involving x i in the product, the degree of this inequality in any x i is in O g ( 3 ), and it does not involve the other variables. In order to be able to use Proposition 3, we must model this inequality by an equation, which is done classically by introducing a new variable T and by using the equation

T • i d g (x i ) i<j Res(u i , u j ) = 1.
To conclude, we have a polynomial system with two blocks of variables: the 2g variables x i and y i and the g 2 -g variables coming from the coefficients of P and Q. The degree of the equations in the first block of variables grows cubically in , while the degree in the other block of variables depends only on g. The system therefore verifies the conditions of Proposition 3 and the complexity follows, provided that we can show that the system is 0-dimensional and radical.

Let us consider the sub-variety S ⊂ J[ ] consisting of -generic -torsion elements, and I the corresponding ideal. More precisely, we see I as the ideal of a sub-scheme of the -torsion scheme, which is the kernel of a finite and étale map because is coprime to the characteristic. Therefore I is 0-dimensional and radical. Since all the elements in S have the same weight g we can use the Mumford coordinates u(X), v(X) with deg u = g and deg v < g -1 as a local system of coordinates to represent them. But the polynomial system that we have built is with the (x i , y i ) coordinates, that is, it generates the ideal I unsym obtained by adjoining to the equations defining I the 2g equations coming from u(X) = (X -x i ) and y i = v(x i ). Then we have deg I unsym = g! deg I. By the -genericity condition, all the fibers in the variety have exactly g! distinct points corresponding to permuting the (x i , y i ) which are all distinct. Therefore the radicality of I implies the radicality of I unsym and we can apply Proposition 3 to our polynomial system. We emphasize that, although the algorithm in Proposition 12 is Monte Carlo, we expect that it returns a correct and verifiable result in most of the cases. Indeed, if all the 2g -1 nonzero -torsion elements are -generic (which is the situation that we expect to happen in most of the cases) and if the algorithm returns the correct result, then we can check that these elements are indeed -torsion elements, and that we have all of them. In that favorable case, the proof of Proposition 2 is completed.

Non-generic cases

For most of the curves, we expect that for all the primes considered in Algorithm 1 the set J[ ] contains only -generic elements (apart from 0), so that the result of the previous section is sufficient. If this is not the case, then it is very likely that the orbit under the Frobenius endomorphism of the -torsion elements computed contains an F -basis of J[ ], so that we can easily recover the missing elements using the group law or the Frobenius. Still, unless we could prove otherwise, we can not exclude the case where the set of -generic -torsion elements generate a proper subgroup of J[ ] which is stable under the action of ϕ. In that unlikely case, we would maybe not be able to deduce χ . An option is then to skip this unlucky and proceed with the algorithm; this would only marginally increase the largest considered . But then, we would be left to prove that the number of unlucky 's is small enough. The number of isomorphism classes of hyperelliptic curves of genus g over F q grows like q 2g-1 . Assuming the non-genericity of the -torsion is driven by the vanishing of some polynomial of degree O (1) that we consider as random, independent events, we see that when q grows to infinity and g is fixed, the probability of having a single curve for which the -torsion is non-generic for all the small 's goes to zero. But this kind of discussion does not seem to lead to a provable result, because the polynomials involved have a highly non-generic Galois structure, and furthermore because there exist curves for which this can even become more structured, for instance those whose Jacobians are not simple.

Our only remaining option is to perform a tedious, systematic study of all the non-generic cases and to show that they can all be modelled by polynomial systems that can be solved within the target complexity. The number of these systems must also be bounded independently of , so that with our setting where g is fixed and q grows to infinity the global complexity remains the same. All this is the purpose of the present section. As a warm-up, we will first describe some simple degeneracy cases and, informally, how to deal with them. Then we will consider the most general non-generic case and construct the polynomial system to model it in order to reach our result. Note that all these systems will have more equations (O(g 4 ), see Table 2) than variables (O(g 2 ), see Table 1), which is no wonder since we expect them to have no solution in general.

Simple degeneracies

Case 1: Low weight -torsion elements. In order to compute the -torsion elements that satisfy all the conditions of -genericity except that their weight is less than g, we can proceed as in the proof of Proposition 12 with the following modifications. This time, D = w i=1 (P i -∞), and the only difference is that there are w points instead of g. Following the same method, we search ϕ of the form P (X) + Y Q(X) such that the points in the reduced divisor • (P i -∞) are exactly the zeros of ϕ. We now want ϕ to have gw points of intersection with C instead of g 2 , and we similarly deduce 2 deg(P ) ≤ gw and 2 deg(Q) ≤ gw -2g -1. By similar parity considerations we deduce that exactly one of these bounds is even, and the corresponding polynomial will be made monic to normalize the function. The number of variables from P and Q is thus gw -g, and after adding the 2w variables x i and y i , we have a total of (g +1)w +w -g variables. As for the number of equations, the number of congruence relations is now w but the relations themselves remain unchanged, and we get a total of (g + 1)w equations after adding the w equations y 2 i = f (x i ). Since we keep the degrees unchanged but reduce the number of variables, the complexity bounds are still valid in this case.

Case 2: Multiple points in the -torsion divisor. It may happen that the reduced forms of -torsion divisors contain multiple points. In that case, the u-coordinate in the Mumford representation of such a point is not squarefree. Although the modelling by the polynomial system described in Section 4 is still faithful, such multiple points will induce multiplicities since what we actually compute is the variety describing the points in the reduced divisor. Therefore, the ideal generated by the polynomial system is not radical in this case. We use the following workaround: For λ = (λ 1 , . . . , λ k ) a partition of w, we write a polynomial system generating a radical ideal whose solutions represent the reduced divisors of the form D =

λ 1 P 1 + • • • + λ k P k -w ∞.
To build this polynomial system, we do as if we were looking for elements of weight k, but instead of multiplying P i by , we multiply it by λ i , using Cantor's polynomials δ λi and ε λi . This system has the same number of variables and equations as if we were looking for elements of weight k. Since λ i is bounded above by g, the degrees of the equations are multiplied by a quantity which depends only on g but not on . Consequently, the complexity bounds are still valid in this case. To avoid multiplicity problems that could arise from subpartitions of λ, we add the inequalities x i = x j for i = j, where x i is the x-coordinate of P i . Again, this does not change our complexity estimate.

Case 3: Low weight after multiplication by .

We study here the case where thegenericity property that is not verified is that the • (P i -∞) are of weight g, all the others being satisfied. We denote by w i ≤ g the weight of •(P i -∞). Then each u i will have degree w i , so that each congruence relation (2) yields only w i equations instead of g. In Cantor's article (on top of page 141 in [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF]), it is stated that • (P i -∞) is of weight w i if and only if for any k such that w i < k ≤ g we have ψ -k+wi+1 (x i ) = 0 and ψ -g+wi (x i ) = 0, where the polynomials ψ i are efficiently computable and of degrees bounded by O g ( 2 ). Therefore the total number of equations is unchanged. Since the function ϕ will have to vanish at i w i points instead of g 2 , we also reduce the degree of P and Q accordingly. The number of variables from P and Q thus becomes i w i -g which is smaller than in the generic case, while the number of equations remains the same, and their degrees are also smaller. Thus we can still describe this non-generic situation with systems that can be handled within the same complexity bounds.

Case 4: Non semi-reduced principal divisor. We now consider the case where thegenericity property fails due to the presence of a point of abscissa ξ which appears with positive multiplicity ν i in an • (P i -∞) and with a negative multiplicity -ν j in another • (P j -∞). Let ν = min(ν i , ν j ). This event implies that (X -ξ) ν divides both P and Q so that we can write ϕ(X, Y ) = (X -ξ) ν ( P (X) + Y Q(X)), with P coprime to Q. The number of variables coming from ϕ is reduced compared to the generic case: we add one (the variable ξ), but the number of coefficients in P is reduced by ν compared to P , and the same is true for Q and Q. To write the conditions on ϕ, we write the congruences exactly like in the generic case and we add conditions to ensure that the multiplicities are respected. Namely, u i , u j and v i + v j must all be divisible by (X -ξ) ν , which adds 3ν ≤ 3g equations. The degree in ξ in these equations is bounded by g 2 . Since this does not depend on , the complexity result is maintained. The general study will cover the case where there are several ξ's at which the semi-reduction genericity assumption fails. Also, there is no reason why such a root ξ should occur in only two of the • (P i -∞)'s. Such a situation will be also taken into account in Section 5.2.

Case 5: Multiplicity in •D. The last situation that could lead to not satisfying -genericity is when the same point is shared within different • (P i -∞), which causes some trouble as the congruence relations of the generic case will not be able to handle the subsequent multiplicity. Note that if the multiplicity occurs only within a single • (P i -∞) this is already dealt within the generic case. One can view our method as using the Chinese remainder theorem on the modular conditions [START_REF] Atkin | Finding suitable curves for the elliptic curve method of factorization[END_REF] to see that multiplicities within a single congruence is handled whereas common factors within different u i -polynomials are an obstacle that needs special strategies. There is some similarities with the previous case that also implies a common factor between two different u i 's.

We devise the following workaround: instead of considering the congruences modulo the u i 's separately, we group them into a single congruence of the form P + QV ≡ 0 mod U , with U = i u i and V a polynomial whose coefficients shall be new variables such that V ≡ v i mod u i for all i. Note that if some non semi-reduced case occurs simultaneously, U must actually be divided by the aforementioned X -ξ; such situations will be dealt with later, in the general study (Section 5.2). In order for V to encode enough information and ensure that the condition P + QV ≡ 0 mod U enforces a function with exactly the correct principal divisor, we have to follow Mumford's representation and add the condition

U |V 2 -f , with deg V < deg U .
Together with the other conditions on U and V , we then have existence and unicity (up to a constant factor): they are the result of Cantor's composition algorithm.

In order to write the polynomial system modelling this situation, some care must be taken so as to stay within the scope of Proposition 3. The polynomial U is of degree g 2 and its coefficients are polynomials in the x i 's of degrees bounded by O g ( 3 ). New variables are added for the coordinates of V . For each i, the condition V ≡ v i mod u i is converted in O(g) equations, with degrees O g ( 3 ) in x i and 1 in the coordinates of V . The condition U |V 2 -f contributes to O(g 2 ) additional equations, each of them of degree 2 in the coordinates of V , and degree O g ( 3 ) in the coordinates x i . And finally, the equation P + QV ≡ 0 mod U , contributes also to O(g 2 ) equations, each of them of degree 1 in the coordinates of V , P and Q, and of degree O g ( 3 ) in the coordinates x i . Skipping the details, we can again apply Proposition 3 and get the expected complexity.

Combining all possible degeneracies

A data structure to describe each type of non-genericity. We want to describe a family of polynomial systems that covers all the possible non-generic cases, possibly mixing all kind of problems that have been listed. We begin by grouping together non-genericity situations that can be covered by the same polynomial system.

We consider an -torsion divisor D of weight w ≤ g (like in case 1). Next, a partition λ = (λ 1 , . . . , λ k ) of w is picked to represent the multiplicity pattern in the u-coordinate of the -torsion divisor, as in case 2 so that D = k i=1 λ i (P i -∞). Then, a vector t = (t 1 , . . . , t k ) is chosen, to represent the weights of the P i after multiplication by λ i as in case 3: For i in [1, k], the reduced divisor λ i • (P i -∞) is of weight t i . Then, we need to consider how many common or opposite points these divisors are in their support to take into account the cases 4 and 5. We denote by Q 1 , . . . , Q s the points in the union of the supports of all the reduced divisors λ i • (P i -∞), keeping only one point in each orbit under the hyperelliptic involution. We represent the non-genericity by a k × s matrix M such that its non-zero entries

m ij verify m ij = ord Qj (λ i • (P i -∞)) when Q j is in the support of λ i • (P i -∞) or m ij = -ord Q j (λ i •(P i -∞)) when the hyperelliptic conjugate Q j of Q j is
in the support. Note that this matrix, that we shall call the matrix of shared points, represents both multiplicities and non-semi-reduction. Since the row i represents what happens with points in the support of λ i • (P i -∞), which is of weight t i , the sum of the absolute values of the entries of the row i of M is equal to t i . Also, by construction, in each column, there is at least one non-zero entry. An additional complication arises when one of the P i is a ramification point, i.e. when its y-coordinate is zero, because this would cause multiplicities if care is not taken, leading to non-radicality of the polynomial system we build. Since this corresponds to P i -∞ being of order 2, the weight t i is equal to λ i mod 2, namely 0 or 1. If t i = 0, then the divisor D -λ i (P i -∞) is also an -torsion divisor of weight w -λ i , so that we can reconstruct D from another polynomial system. There is however no obvious way to preclude the possibility t i = 1. Therefore, we will encode the fact that P i is a ramification point by a bit i that can be set only in the cases where t i = 1 and λ i = 1.

A tuple (w, λ = (λ 1 , . . . , λ k ), t = (t 1 , . . . , t k ), = ( 1 , . . . , k ), M ) is from now on the piece of data with which we represent a non-generic situation, and a polynomial system will be associated to each tuple. Changing the order of the columns of M amounts to permuting the points Q j . Also, changing the sign of all the entries of a column j corresponds to taking the opposite of the point Q j . While it would not change the final complexity not to do so, it therefore makes sense to consider only normalized tuples, in the sense that the columns of M are sorted in lexicographical order, and the choice between a point Q j and its opposite is done so that the sum of all elements in the corresponding column is nonnegative. We remark that this is not enough to guarantee that two normalized tuples do not describe similar situations. For instance, if λ = (1, . . . , 1) and two t i values are equal, then permuting the two corresponding rows could lead to another normalized matrix that would describe the same situation. This is not a problem for the general algorithm: we might get the same -torsion elements from two different systems, but what is important to us is non-multiplicity (i.e. radicality of the ideal) in each individual system. Definition 13. A normalized non-genericity tuple is a tuple (w, λ, t, , M ), where 1 ≤ w ≤ g is an integer, λ = (λ 1 , . . . , λ k ) is a partition of w, t and are vectors t = (t 1 , . . . , t k ) and = ( 1 , . . . , k ) of the same length as λ with 1 ≤ t i ≤ g and i ∈ {0, 1}, where i can be 1 only if t i = 1 and λ i = 1, and finally M is a matrix with k rows and s columns, where 0 ≤ s ≤ g k, and its entries are integers such that:

• For all 1 ≤ i ≤ k, the sum of the absolute values of the entries on the row i is equal to t i ;

• The columns are sorted in lexicographical order;

• The sum of the rows of the matrix is a vector whose coordinates are nonnegative.

From the discussion above, any -torsion element is described by (at least) one normalized non-genericity tuple. In the following we will give a polynomial system for each normalized non-genericity tuple, so that all -torsion elements described by it are modelled by this system. Furthermore, the system will have the properties required to apply Proposition 3, so that the complexity result will follow.

Before starting this, we discuss briefly a bound on the number of normalized non-genericity tuples. Assuming everything is always of maximal size, and not sorted, we have g choices for w, then at most g g choices for λ and t, at most 2 g choices for , and finally at most (g 2g+1 ) g 2 choices for M , which gives g O(g 3 ) . As bad as it is, such a factor that depends only on g will not hinder the final complexity estimate in O g (log q) O(g) , as explained in Section 2.

Non-generic division polynomials. The expression of λ i • (P i -∞) in Mumford representation will be the same as in the generic case when its weight t i is equal to g and Lemma 10 can be applied. But when t i is strictly less than g, the weight-g coordinate system is no longer available; this is explicitly visible by the fact that the denominator e g (x i ) of the coefficients of the v-polynomial vanishes.

Therefore we need to use a weight-t coordinate system for describing a non-generic divisor λ i • (P i -∞) in Mumford representation. In this paragraph, in order to keep simple notation, we will work with • (P i -∞), keeping in mind that we do not impose any condition on , so that we can later replace by λ i .

We consider, for 1 ≤ t < g, the set V ,t of points of the curve which are mapped to a weight-t divisor after multiplication by :

V ,t = {(x, y) ∈ C | • ((x, y) -∞) is of weight t} .
This is a (possibly empty) variety of dimension 0 that can be described with the classical (generic) division polynomials of Cantor: we define ∆ ,t = GCD(ψ (x), ψ -1 (x), . . . , ψ -g+t+1 (x)), so that V ,t is precisely the set of points (x, y) for which ∆ ,t (x) = 0 and ψ -g+t (x) = 0, as stated by Cantor in [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF] on page 141. The polynomial ψ is essentially the square root of the leading coefficient of δ . It can be computed efficiently and has degree in O g ( 2 ) by Theorem 8.17 of [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF]. To avoid multiplicities, we define ∆ ,t (x) the square-free polynomial whose roots are exactly the roots of ∆ ,t (x) that are not roots of ψ -g+t (x). The degree of ∆ ,t (x) is again bounded by O g ( 2 ). Furthermore since the points of V ,t come in pairs of conjugate points sharing the same x-value, the degree of V ,t is 2 deg ∆ ,t (x). Definition 14. The non-generic division polynomials u ,t and v ,t are the polynomials in X with coefficients in F p [x, y]/( ∆ ,t (x), y 2 -f (x)) such that

• ((x, y) -∞) = u ,t (X), v ,t (X) ,
in weight-t Mumford representation: u ,t (X) is monic of degree t, v ,t (X) is of degree at most t -1 and they satisfy u ,t | v 2 -f . Just like for the classical division polynomials, the coefficients of u ,t (X) and of 1 y v ,t (X) are in F p [x]/ ∆ ,t (x) (they do not depend on y) and we can choose representatives of them that are polynomials of degree less than deg ∆ ,t (x). Hence, the bounds given in Lemma 10 are also valid for the non-generic division polynomials; and since there are no denominators in the coefficients of v ,t (X), the other part of Lemma 10 also holds trivially.

The non-generic division polynomials can be computed efficiently, once the classical division polynomials are known: the polynomial ∆ ,t (x) can be easily deduced, and then working in the quotient algebra yields the result in a time O g ( 2 ), which is negligible compared to the other parts of the algorithm.

Polynomial system derived from a normalized non-genericity tuple. We now want to write a polynomial system whose solutions are the -torsion elements following a given normalized non-genericity tuple (w, λ, t, , M ).

First, we need variables for the coordinates of the P i such that the -torsion element is D = k i=1 λ i (P i -∞), with P i = ±P j for all i = j. As a consequence, we introduce 2k variables for the coordinates (x i , y i ) of all the points P i . Since these points are on the curve, they verify y 2 i = f (x i ), however if P i is a ramification point this can be simplified into y i = 0 = f (x i ), which avoids the multiplicities. We get a first set of equations

y 2 i = f (x i ) = 0, for all i in [1, k] such that i = 0, y i = f (x i ) = 0, for all i in [1, k] such that i = 1. (Sys.1)
As we just discussed, we must model the fact that P i = ±P j for i = j. This is done via the following set of inequalities:

x i = x j , for all i, j in [1, k] such that i = j. (Sys.2)
The next step is to enforce the fact that the element λ i • (P i -∞) is of weight t i . For the indices for which t i < g, this is encoded by the equation defining V λi ,ti : ∆λi ,ti (x i ) = 0, d ti (x i ) = 0, for all i in [1, k] such that t i < g, (Sys.3) while for the indices for which t i = g, this is encoded by the non-vanishing of the leading coefficient of the Cantor polynomial in degree λ i :

d g (x i ) = 0, for all i in [1, k] such that t i = g. (Sys.4)
We now need to model the fact that the λ i • (P i -∞) satisfy the conditions given by the matrix M . We write λ i • (P i -∞) = u i (X), v i (X) in Mumford representation, where u i (X) and v i (X) are Cantor's classical division polynomials in degree λ i if t i = g or the non-generic division polynomials u λi ,ti and v λi ,ti , if t i < g. In both cases, these are polynomials in X whose coefficients are polynomials in x i and y i . Recall that the entries of M , denoted by (m ij ) i∈ [1,k],j∈ [1,s] , are such that m ij is the order of Q j in λ i • (P i -∞) if it is positive, or the opposite of the order of Q j if it is negative. To this effect, we introduce s new variables ξ j for the abscissae of the Q j , and the following equations enforce the multiplicities:

u (n) i (ξ j ) = 0, for all i, j in [1, k] × [1, s] and for all n ≤ |m ij | -1 (Sys.5) u (|mij |) i (ξ j ) = 0, for all i, j in [1, k] × [1, s] (Sys.6) v i (ξ j ) -v i (ξ j ) = 0, for all i, i , j such that m ij m i j > 0 (Sys.7)
v i (ξ j ) + v i (ξ j ) = 0, for all i, i , j such that m ij m i j < 0 (Sys.8)

ξ j = ξ j , for all j = j . (Sys.9)

In Equations Sys.5 and Sys.6, the notation u

(n) i
is for the n-th derivative of u i . This simple way of describing multiple roots is valid because the characteristic is large enough.

The next step of the construction is to consider a semi-reduced version of the divisor

• D = k i=1 λ i • (P i -∞).
This semi-reduction process can be described directly on the matrix M : if two entries in a same column have opposite signs, a semi-reduction can occur (corresponding to subtracting the principal divisor of the function (x -ξ j )), thus reducing the difference between these entries. This semi-reduction can continue until one of these two entries reaches zero. This whole process can be repeated as long as there are still columns containing entries with opposite signs. This is formalized in Algorithm 2, which takes as input a matrix M and returns a matrix M with the same dimensions such that if M describes all the multiplicities in a divisor, then M describes all the multiplicities of a semi-reduced divisor equivalent to the input divisor. More precisely, the matrix M satisfies the following properties: [START_REF] Adleman | Counting points on curves and Abelian varieties over finite fields[END_REF] In each column, all elements are nonnegative; [START_REF] Atkin | Finding suitable curves for the elliptic curve method of factorization[END_REF] The sum of the rows of M equals the sum of the rows of M ; (3) For all i, j such that M i,j is nonnegative, M ij ≤ M ij .

Data: M the k × s matrix of shared points of the system Result: M , the matrix after semi-reduction M ← k × s zero matrix for j from 1 to s do

µ j ← k i=1 M ij for i from 1 to k do if M ij > 0 then M ij ← min(M ij , µ j ) µ j ← µ j -M ij else M ij ← 0 end end end return M
Algorithm 2: Reducing the matrix of shared points

The function ϕ that we will use to model the principality of the divisor • D will have two parts: a product of "vertical lines" corresponding to semi-reductions, and a part of the form P (X) + Y Q(X), where P and Q are coprime. Modelling the existence of this second part requires to introduce new entities u i that are the u i polynomials from which we remove the linear factors coming from semi-reduction as described by M . Formally, we have the following equations, defining u i :

u i (X) = u i (X) s j=1 (X -ξ j ) |mij |-mij , for all i ∈ [1, k].
(Sys.10) Indeed, by definition of the matrix M , the factor (X -ξ j ) |mij | divides exactly u i (X), and the factor (X -ξ j ) mij divides exactly u i (X). In order to express these conditions efficiently in the polynomial system, we introduce new variables for the coefficients of the u i polynomials.

Since we are now dealing with a semi-reduced divisor, we can consider its Mumford representation, i.e. two polynomials U and V with the following properties:

U = k i=1 u i , U |V 2 -f, (Sys.11) V ≡ v i mod u i , for all i ∈ [1, k]. (Sys.12)
The expression of U is simple enough, so we do not have to introduce new variables for its coefficients. However, this will be necessary for the coefficients of the V polynomial. Finally, in order to impose that the semi-reduced part of ϕ has exactly the zeros described by this divisor, we have the equation P + QV ≡ 0 mod U, (Sys.13) which is expressed with new variables for the coefficients of P and Q.

In Table 1, we summarize all the variables used in the polynomial system and count them. A key quantity for this count is the degree of U which is the sum of the degrees of the u i 's. It can be computed directly from the tuple (w, λ, t, , M ). Then, to ensure existence and unicity of the V polynomial to represent the semi-reduced divisor, we have to impose that deg V < deg U , so that we have exactly deg U variables for the coefficients of V . For the polynomials P and Q, we need the degree of P 2 -Q 2 f to be exactly deg U . After a normalization like in Section 4 depending on the parity of deg U , we get deg U -g variables for their coefficients.

Variables

Number of variables Bound Coordinates (x i , y i ) of P i 2k 2g Abscissae ξ j of shared points s, column-size of the matrix M g 2 Coefficients of the u i polynomials deg

U = i (t i -j (|m ij | -m ij )) g 2 Coefficients of the V polynomial deg U g 2 Coefficients of the P and Q polynomials deg U -g g 2 -g Total s + 2k + 3 deg U -g 4g 2 + g
Table 1: Summary of the variables in the polynomial system corresponding to a normalized non-genericity tuple (w, λ, t, , M ).

In order to apply Proposition 3, we need to evaluate the degrees of all the equations and inequalities that we have listed, with respect to two groups of variables: The first group contains just the variables x i and y i , and we will denote deg 1 (f ) the degree of a polynomial f with respect to those variables (said otherwise, deg 1 (f ) is the degree of f if we consider only the symbols x i , y i as variables, and all the other indeterminates are considered as parameters). The second group of variables contains all the other indeterminates and the degree with respect to this group is denoted by deg 2 .

The crucial point is to ensure that each polynomial equation has a deg 1 bounded by O g ( 3 ), while deg 2 is bounded by O g [START_REF] Adleman | Counting points on curves and Abelian varieties over finite fields[END_REF]. For the inequalities, we require the same degree conditions: Indeed, an inequality f = 0 can be modeled by the equality T • f -1 = 0, where T is a fresh variable that belongs to our second group of variables. This trick requires only one more variable for each inequality and the degree of the equation T • f -1 = 0 is only one more than the degree of the inequality. Since the number of inequalities is bounded by O g [START_REF] Adleman | Counting points on curves and Abelian varieties over finite fields[END_REF], the number of extra variables required in the second group will not impact the asymptotic complexity (the second group already contains O g (1) variables). We remark that the input of the geometric resolution algorithm over fields of characteristic 0 in [START_REF] Giusti | A Gröbner free alternative for polynomial system solving[END_REF] allows inequalities. However, we use the aforementioned trick to model inequalities by equalities since the solving method that we use is the variant for positive characteristic whose complexity analysis is given in [START_REF] Cafure | Fast computation of a rational point of a variety over a finite field[END_REF]Thm. 4.8].

The number of equations and inequalities and their degrees with respect to the two groups of variables can be easily checked and are summarized in

k i=1 t i ≤ g 2 O g ( 3 ) ≤ g Eq. Sys.11 deg U ≤ g 2 0 O(g 3 ) Eq. Sys.12 k i=1 deg u i ≤ g 2 O g ( 3 ) O(g 2 ) Eq. Sys.13 deg U ≤ g 2 0 O(g 3 )
Table 2: Summary of the degrees of the equations in the polynomial system corresponding to a normalized non-genericity tuple (w, λ, t, , M ).

Finally, since we have been very careful in describing elements that are -torsion points on J, without room for parasite solutions or multiplicities, we can again appeal to the finite and étale property of multiplication by in J to deduce that the system is 0-dimensional and radical. Therefore, by Proposition 3, each system can be solved in the claimed complexity bound. To conclude the proof of Proposition 2, and hence of our main result, we need a few more observations. First, notice that the solutions of our polynomial systems can be grouped by weight of the -torsion divisor: once a geometric resolution of two 0-dimensional sets V 1 and V 2 are known, a geometric resolution of V 1 ∩ V 2 can be computed very efficiently. The strategy to do so is to change the primitive element of the geometric resolutions for a random element, so that both resolution share the same primitive element. This can done within complexity linear in the number of variables and polynomial in deg(V 1 ∪V 2 ) using Algorithm 6 in [START_REF] Giusti | A Gröbner free alternative for polynomial system solving[END_REF]. Then, computing the LCM of the univariate polynomials of the geometric resolutions and interpolating the parametrization provides a geometric resolution of V 1 ∪V 2 . Using this procedure for regrouping the solutions of all the systems derived from the non-degeneracy tuples with the same weight w provides geometric solutions of J w [ ] within the claimed complexity.

Finally, we need to transform the Monte Carlo algorithm from Proposition 3 in a Las Vegas algorithm. This can be easily achieved since the probability that the Monte Carlo algorithm succeeds is bounded below by a quantity which does not depend on the input size, and the output can be verified since we know that the sum of the degrees of the varieties J w [ ] for w ∈ [1, g] must equal 2g -1. Consequently, once all polynomial systems corresponding to nongeneric situations have been solved, it is easy to count the number of -torsion elements found and to check that none of them is missing by comparing their number with the theoretical value 2g -1. The Las Vegas algorithm consists in repeating the Monte Carlo algorithm until the result is verified and is correct (i.e. all elements found are -torsion elements and none of them is missing). The expected complexity of the Las Vegas variant equals the complexity of the Monte Carlo variant up to multiplication by a constant. This concludes the proof of Proposition 2.

Proof of Lemma 10

We restate Lemma 10: for a non-necessarily prime integer > g, the polynomial δ (X) of degree g in X has coefficients in F q [x] whose degrees in x are bounded by g 3 /3 + O g ( 2 ); the polynomial ε (X)/y has coefficients in F q (x) such that the degrees of the numerators and the denominators have degrees bounded by 2g 3 /3 + O g ( 2 ). Furthermore, the roots of the denominators are roots of the leading coefficient of δ (X).

This section strongly relies on [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF]. From here, all the notations, references to propositions, expressions or definitions are taken from this paper, except for the notation of the base field that we choose to be F q , and the torsion level is denoted instead of r.

Proof. Technicalities arise from the normalizations required to manipulate entities that are polynomials in x (and not rational fractions), without odd power of y involved. In Cantor's article, this normalization often depends on the parity of -g. We will concentrate on the case where g is even; for the other case some formulae must be adapted, multiplying or dividing by 2y at various places.

We recall that ν = ( 2 --g 2 + g)/2 as defined in (8.5), so that ν = ν -1 + -1. By combining (8.7) and (2.3), we obtain δ (z) = (2y) 2ν (4y 2 z) A (4y 2 z) 2 -B (4y 2 z) 2 E(4y 2 z) ,

where A and B are unnormalized versions of α and β defined in (8.7) and E(z) is defined by E(z) = f (x -z). For our purpose, it is easier to deal with non-truncated versions of α and β . Let us then introduce the following quantities, inspired from (8.7): ᾱ (z) = 2(2y) ν -1 -1 A (4y 2 z), and β (z) = (2y) ν -1 B (4y 2 z), so that δ can be rewritten as δ (z) = 1 4z ᾱ (z) 2 -1 y 2 β (z) 2 E(4y 2 z) .

By Theorem (8.15), the coefficients of α (z) and β (z) are polynomials in F q [x], and the proof is also valid for the non-truncated versions ᾱ (z) and β (z). Note that here we use the fact that g is even, so that the potential adjusting factor (2y) g is an even power of y that can be rewritten in terms of f (x). The polynomial E(4y 2 z) has coefficients which are polynomials in x of degree bounded by (2g +1) 2 . Therefore, in order to obtain a degree bound for the coefficients of δ (z), it is sufficient to bound the coefficients of ᾱ (z) and β (z). For convenience, for any polynomial P in z, with coefficients in F q [x], we write degmax(P ) the maximum degree in x of its coefficients. We are interested in a bound at fixed genus g and when grows to infinity and we use the O g notation. For k in [1, ], we will use an induction to bound degmax(ᾱ k (z)) and degmax( βk (z)). For k ≤ g + 1, none of these quantities depends on , so that all the degrees can be bounded by an expression in g only, i.e. in O g [START_REF] Adleman | Counting points on curves and Abelian varieties over finite fields[END_REF]. For k ≥ g + 1, we start from the equation (3.14) where we substitute k for r, we evaluate it at 4y 2 z, and we multiply by (2y) 2ν k -1 , so that we obtain: (2y) ν k f k-1 ᾱk+1 (z) = (2y) ν k +k-1 f k ᾱk (z) -(2y) ν k +2k-1 f k+1 z ᾱk-1 (z), where all the polynomials have coefficients in F q [x]. The expression for βk is exactly the same, but we have to multiply the expression in (3.14) by (2y) 2ν k in that case. By (8.7), (8.16) and Theorems (8.15) and (8.17), for any k, the quantity (2y) ν k f k is a polynomial in x of degree g(k 2 -g 2 )/2. Therefore the right-hand-side of the recurrence relation have coefficients with degrees bounded by an expression of the form max degmax(ᾱ k (z)), degmax(ᾱ k-1 (z)) +gk 2 /2 up to a term linear in k and cubic in g. And we finally get: degmax(ᾱ k+1 (z)) ≤ max degmax(ᾱ k (z)), degmax(ᾱ k-1 (z)) + gk 2 /2 + Err g (k), where Err g (k) is a polynomial linear in k and cubic in g. Again, this inequality is also valid for βk . By induction, we then get the following bounds: degmax(ᾱ (z)) ≤ g 3 6 + O g ( 2 ), and degmax( β (z)) ≤ g 3 6 + O g ( 2 ).

We can then propagate these bounds in the expression of δ and we get degmax(δ (z)) ≤ max(2 degmax(ᾱ (z)), 2 degmax( β (z)) + degmax(E(4y 2 z)), so that we get the claimed result concerning δ .

The fact that ε (z)/y has coefficients in F q (x) follows directly from Equation (8.13) that we recall here:

ε (z) = y z ψ 2 -1 δ +1 (z) -ψ 2 +1 δ -1 (z) ψ -1 ψ 2 ψ +1 mod δ (z).
By (8.11), the leading coefficient of δ (z) is -(4y 2 ) g ψ 2 , so that the property on the denominator of ε can not be easily deduced from this equation, due to the presence of ψ -1 and ψ +1 before the reduction modulo δ (z) occurs. We will prove it below, with a direct geometric argument, but we first give bounds on the degrees of the coefficients of the numerator and the denominator. The polynomial δ (z) is of degree g in z, so that at most two steps of reduction are required to reduce the degree of ε to strictly less than g. In fact, it can be checked that the leading coefficient of ψ 2 -1 δ +1 (z)-ψ 2 +1 δ -1 (z) is zero, so that there is at most only one reduction step. This reduction accounts for an increase of the coefficients' degrees in x by at most degmax(δ ) in the numerator and an increase of the degree of the leading coefficient of δ in the denominator. Since deg ψ = g 2 /2 + O g ( ), the degrees of the coefficients of the numerator of ε (z) are bounded by 2 3 g 3 + O g ( 2 ), and the degree of the denominator is bounded by 3g 2 + O g ( ). It remains to prove the claim on the roots of the denominator of the coefficients of ε (z)/y. For this, we consider the map from the affine part of the curve C aff to J seen as a projective Abelian variety, that sends a point (x, y) to [ ]((x, y) -∞). One of the main points of Cantor's article is that if ψ (x) = 0, then the image by this map is in J \ Θ, where Θ ⊂ J is the subvariety of elements of weight less than g. On this open subset, Mumford coordinates with a monic u of degree g and v of degree at most g -1 give a local set of coordinates that we use to describe the map. The i-th coefficient of v is y times a rational fraction c i in x that gives a finite value at any x for which ψ (x) = 0. Therefore, any root of the denominator of c i is a root of ψ . By Theorem (8.35), the Mumford v-polynomial that we are considering is ε up to a renormalization that will only introduce additional powers of 4y 2 in the denominator. Therefore, any root of the denominator of the coefficients of ε is a root of ψ or of 4y 2 , and both divide the leading coefficient of δ by (8.11).

Remark. The bounds that we obtain are not tight: from [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF], we know that the leading and constant coefficients are in O g ( 2 ) instead of O g ( 3 ). We ran experiments that allow us to conjecture precise degrees for the other coefficients. In these experiments, instead of developing δ x-X 4y 2 and ε x-X 4y 2 to compute the d i 's and e i 's, we computed • ((x, y) -∞) over the function field of the curve. This does not exactly yield the d i 's and e i 's because we actually get d i /d g and e i /e g , thus possibly missing a common factor in all the d i 's and e i 's. We denote di and ẽi the numerators and denominators of the aforementioned fractions, and we compute their degrees for each pair (g, ) with g ≤ 8 and g < ≤ g +20 (which includes non prime values of ). We found that the degrees of the di are consecutive from deg( dg ) up to deg( d0 ) = deg( dg ) + g, with the following values for deg( d0 ).

g 2 -g 3 + g if g -is even g 2 -g 3 + 2g 2 -1 if g -is odd

2 .

 2 It constructs the straight-line program representing f (M ) 1 , . . . , f (M ) n with space S = O(S) and time T = O(T + m n); 3. It applies the probabilistic Turing machine from Theorem 6 to compute a geometric resolution of the algebraic set defined by f

  ) and ny+dy dy ≤ (n y +d y ) ny = O g (1), Lemma 7 provide bounds on S and T . Summing the complexities leads to the claimed complexity estimate. Finally, the probability of success is bounded below by the probability that the sequence f n is reduced and regular (Lemma 5) multiplied by the probability of success of the probabilistic Turing machine in Theorem 6, namely 10/11 • 11/12 = 5/6.

  y can be represented by a DFSLP with time and space O (d x + d y + m) nx+dx in F q [X 1 , . . . , X nx , Y 1 , . . . , Y ny ] such that deg x (µ) ≤ d x and deg y (µ) ≤ d y .We consider the DFSLP which starts by evaluating these monomials. Then we multiply each of these monomials by the corresponding coefficients, and we sum. This costs

				nx	ny+dy ny	.
	Proof. There are nx+dx nx monomials µ This costs less than nx+dx ny+dy ny ny+dy nx ny
	m nx+dx nx	ny+dy ny	multiplications and m( nx+dx nx	ny+dy ny

(d x + d y -1) multiplications, using a naive algorithm.

Table 2 .

 2

	Equations reference Number of equations (and bound) deg 1	deg 2
	Eq. and Ineq. Sys.1 2k ≤ 2g		2g + 1 0
	InEq. Sys.2	k(k -1)/2 ≤ g(g -1)/2	1	0
	Eq. and Ineq. Sys.3 ≤ 2g		O g ( 3 ) 0
	InEq. Sys.4	≤ g		O g ( 3 ) 0
	Eq. Sys.5 InEq. Sys.6	k i=1 ks ≤ g 3	s j=1 |m ij | ≤ g 4	O g ( 3 ) ≤ g O g ( 3 ) ≤ g
	Eq. Sys.7 and Sys.8 ≤ k 2 s ≤ g 4	O g ( 3 ) ≤ g
	InEq. Sys.9	≤ s 2 ≤ g 4	0	1
	Eq. Sys.10

Acknowledgements. We are grateful to Éric Schost and Guillermo Matera for fruitful discussions and for pointing out important references. We also wish to thank anonymous referees for their comments which helped improve the paper.

Concerning the ẽi , the degrees are consecutive from deg(ẽ g-1 ) up to deg(ẽ 0 ) = deg(ẽ g ), the latter being equal to

Cantor [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF] gave simple expressions for the leading term and constant term of δ (respectively -(4y 2 ) g ψ 2 and (-1) g+1 ψ -1 ψ +1 ), from which we can deduce the degrees of d 0 and d g by evaluating δ at (x -X)/4y 2 . Assuming that there is no common factor to all the d i 's when g -is even, while the GCD of all the d i 's is f g-1 when g -is odd, these theoretical degrees are consistent with our experiments.