Improved Complexity Bounds for Counting Points on Hyperelliptic Curves

by   Simon Abelard, et al.

We present a probabilistic Las Vegas algorithm for computing the local zeta function of a hyperelliptic curve of genus g defined over F_q. It is based on the approaches by Schoof and Pila combined with a modeling 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(( q)^cg) as q grows and the characteristic is large enough.



There are no comments yet.


page 1

page 2

page 3

page 4


Counting points on hyperelliptic curves with explicit real multiplication in arbitrary genus

We present a probabilistic Las Vegas algorithm for computing the local z...

Counting points on genus-3 hyperelliptic curves with explicit real multiplication

We propose a Las Vegas probabilistic algorithm to compute the zeta funct...

Subtrajectory Clustering: Finding Set Covers for Set Systems of Subcurves

We study subtrajectory clustering under the Fréchet distance. Given one ...

Fast and Simple Methods For Computing Control Points

The purpose of this paper is to present simple and fast methods for comp...

The VC Dimension of Metric Balls under Fréchet and Hausdorff Distances

The Vapnik-Chervonenkis dimension provides a notion of complexity for sy...

Tighter Bounds for Reconstruction from ε-samples

We show that reconstructing a curve in ℝ^d for d≥ 2 from a 0.66-sample i...

Counting points on abelian surfaces over finite fields with Elkies's method

We generalize Elkies's method, an essential ingredient in the SEA algori...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Since the discovery of Schoof’s algorithm [24], 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 [24] and Pila’s algorithms [20].

When the characteristic of the base field is small, Kedlaya’s and Satoh’s approaches [14, 23] and their variants compute very efficiently the number of rational points of Jacobians of hyperelliptic curves. We can also mention Lauder-Wan’s [17] and Lauder’s [16] methods that can handle very general varieties. The current best algorithms in this family for rather general curves are by Tuitman [27, 28]. However, the complexities of these -adic algorithms are exponential in , where is the characteristic of the base field. This dependency can be made as low as , thanks to the work of Harvey [11]. 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 [2], and in [8] 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 [24], and later extended to Abelian varieties in [20]. In particular, Pila showed that the local zeta function of a -dimensional Abelian variety can be computed within operations, where and the constant in the are functions of (but they do not depend on ). 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 [21]. Complexity improvements were obtained in [13] and [1]. The latter article gives a deterministic algorithm for counting points on hyperelliptic curves with complexity . 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 be a hyperelliptic curve of genus and be its Jacobian. When is a prime different from the characteristic of the base field, the -torsion group is isomorphic to and the characteristic polynomial of the restriction of on is exactly . Furthermore, . The principle of Schoof-Pila’s algorithm is to pick elements in and to find conditions on the coefficients such that is equal to in . By testing all the tuples up to the symmetries coming from the functional equation of (and possibly many ), the number of possibilities for 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 where is the -division polynomial, which has degree . The dominant part of the complexity is the computation of 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 [9]. This appears to be also the case for . 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 [10], 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 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 [5]. 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.


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.

2 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

, a positive integer and a squarefree univariate polynomial of degree , let be the hyperelliptic curve with Weierstrass form . Compute the numerator of the local zeta function of :

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 is sufficiently large compared to . 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 such that for all genus , there exists an integer such that for all prime power larger than with and for all hyperelliptic curves of genus defined over , the numerator of the local zeta function of can be computed with a probabilistic algorithm in expected time bounded by .

This complexity result is summarized by the notation , keeping in mind that is fixed and grows to infinity. Indeed, such a complexity statement can hide any factor that depends only on : a running time in can be transformed into by taking a value larger than and adjusting , so that .

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

A classical geometrical object associated to a genus curve is its Jacobian variety. Over the algebraic closure of , it can be described as the multiset of at most points of the curve and it is endowed with an Abelian group structure (it is isomorphic to the degree- 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 -module structure.

Throughout this paper, is a hyperelliptic curve defined over with at least one rational Weierstrass point. Hence admits a Weierstrass model , where is a squarefree monic polynomial of degree . If 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 (it is at most linear in ), so that this will not affect our complexity result.

For practical computations, we need a coordinate system to represent points on the Jacobian of : they shall be encoded via their Mumford representation using coordinates. The group law on points in the Jacobian can be performed with Cantor’s algorithm [4] which operates with elements in Mumford representation at a cost of 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 of the local zeta function is the reciprocal of the characteristic polynomial of the Frobenius morphism on the Jacobian variety of [18, Thm. 5.2];

  • For prime numbers not dividing , the -torsion of the Jacobian variety is isomorphic (as an Abelian group) to [19, Sec. II.6, Prop. page 64], [6, Thm. 4.73]; Therefore is the reciprocal of the characteristic polynomial of the Frobenius seen as an endomorphism of ;

  • The Weil conjectures imply that has the following form over the complex numbers: with [18, Ch. VIII, Thm. 6.1]. Moreover, if denote the coefficients of , the functional equation implies that . Consequently, the absolute value of the coefficients of are bounded by .

Data: a prime power, and a monic squarefree univariate polynomial.
Result: The characteristic polynomial of the Frobenius endomorphism on the Jacobian of the hyperelliptic curve defined over with Weierstrass form .
while  do
       if  divides  then
       end if
      Compute a description of ;
       Compute a matrix with coefficients in representing the action of the Frobenius on ;
       Compute the characteristic polynomial of the matrix ;
end while
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 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 . Since is a -dimensional variety of degree , what we will compute is a geometric resolution of the corresponding radical ideal, that is a univariate squarefree polynomial , together with coordinate polynomials , such that the coordinates of the

-torsion elements are the evaluations of the vector

at the roots of .

To be more precise, the Mumford coordinates are in fact a set of affine systems of coordinates, each corresponding to a different weight of the represented divisor (the definition is recalled in Section 4). The variety will accordingly be represented by a set of geometric resolutions, each encoding -torsion divisors of a given weight . Generically, we expect that all the elements in have weight , except for the neutral element which has weight 0. Most of the article is dedicated to computing efficiently this representation for . The cornerstone of the proof of Theorem 1 relies on the following statement.

Proposition 2.

Let be a hyperelliptic curve of genus over with Weierstrass form ( monic of degree ) and be its Jacobian variety. Let be a prime not dividing . Assuming that the characteristic of is sufficiently large as in Theorem 1, there is a Las Vegas probabilistic algorithm which takes as input and which computes geometric resolutions for the varieties of -torsion points of weight

in the Jacobian variety. This algorithm can be implemented by a Turing machine with space and expected time


Assuming this complexity bound, performing a complexity analysis as done in [20] 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 . 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 [26, Cor. 10.1] that the largest in the loop is at most linear in . From this and Proposition 2, computing the description of as a union of geometric resolutions for all the can be achieved within expected complexity .

Factoring the univariate polynomials involved in the geometric resolutions can be done within the same time bound , since the sum of their degrees is and factoring polynomials in finite fields can be done in time linear in and quasi-quadratic in the degree [29, Thm. 14.14]. Therefore, it is possible to construct a Mumford representation for each -torsion divisor within the same complexity, each of them possibly defined over a different extension of . In fact, due to the rationality of the group law that acts on , one of these extensions of contains all the others.

Using elementary linear algebra for the Frobenius endomorphism acting on (seen as an -vector space), we can deduce . We first compute a basis of 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 complexity bound. The loop is repeated times, and this additional factor does not affect the overall complexity. ∎

3 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, denotes the algebraic closure of . For an ideal , we call dimension of and note the Krull dimension of the quotient ring . Moreover, by identifying a point with the polynomial , there is a dense Zariski open subset such that for any , the algebra is a finite dimensional -vector space of constant dimension, which is called the degree of . A sequence is regular if and for any , does not divide zero in . The sequence is reduced if every intermediate ideal with is radical.

Geometric resolutions.

For describing -dimensional (i.e. finite) sets where is defined over , we use a data structure called a geometric resolution of . The terminology here is borrowed from [3], see also [10]. An -geometric resolution of is a tuple where:

  • The vector is such that the linear form

    takes distinct values at all points in . The linear form is called the primitive element of the geometric resolution;

  • The polynomial equals

  • The polynomials parametrize by the roots of the polynomial , i.e.

We note that our definition is slightly simpler than the one in [3, Sec. 2.1] because we restrict ourselves to the -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 is a polynomial in a ring , then we let (resp. ) denote the degree of (resp. ), where (resp. ) are generic values in .

The following proposition is a cornerstone of our complexity result for computing the -torsion of the Jacobian of a hyperelliptic curve. The statement and its proof combine three main ingredients: (1) the geometric resolution algorithm [10] and its version for finite fields [3], 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 but it has very small cardinality in . 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 . Taking this bi-homogeneous 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 which takes as input polynomial systems with coefficients in a finite field and which satisfies the following property. For any function , for any positive number and for any , there exists a function and a positive number such that for all positive integers such that , , , , , for any prime power such that the prime number dividing satisfies , and for any polynomial system such that

  • for all , and ,

  • the ideal has dimension and is radical,

the Turing machine with input returns an -geometric resolution of the variety

with probability at least

, using space and time bounded above by .


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 -dimensional and radical, then a variant of Bertini’s theorem ensures that the obtained sequence is regular and reduced.

Proposition 4.

[25, Thm. A.8.7] Let be polynomials such that the ideal has dimension and is radical. Let be two integers such that , for all . Let be the characteristic of , and assume that . For an matrix with entries in , let be defined as

Then there exists a nonempty open subset of the space of matrices such that for any , for any , and at any point such that , the derivatives are linearly independent over . In particular, for any , the sequence is reduced and regular.


This is a reformulation of [25, Thm. A.8.7] in the case of finite fields. In [25, Thm. A.8.7], this result is stated over the field , but this statement holds true over any field , provided that an extra separability assumption is satisfied. More precisely, set and let be the variety of pairs such that . In this setting, the extra condition that is required for the proposition to hold is that the projection of to must be separable for all (this is always true in characteristic 0). We refer to [15, Thm. 4.2] for more details on this separability argument. In our setting, the degree of a generic fiber of is bounded by 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 such that is reduced and regular.

Lemma 5.

Let be polynomials satisfying the assumptions of Proposition 4 and such that their total degree is bounded above by . Set and

If is an matrix with entries in picked uniformly at random, then the probability that is a reduced regular sequence is bounded below by .


Let denote an matrix with indeterminate entries

and let be the polynomials defined as

For , we consider the matrix obtained by truncating to its first rows, a new set of variables and the following polynomial system:

This is a system of polynomials of degree bounded above by in variables. By Bézout inequality (see e.g. [12, Thm. 1]), this system defines a variety which is either empty, or its degree is at most . We remark that if is not empty, then it has dimension at least since its vanishing ideal is generated by elements. The Zariski closure of its projection to the space of matrices 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 is also bounded by if . In the sequel, we let denote a polynomial vanishing on of degree bounded by (we set if ).

The Schwarz-Zippel Lemma implies that the cardinality of the set

is bounded above by , for the value of given in the statement.

The proof is concluded by noticing that for any , for any , and for any such that the derivatives span the normal space at to the variety associated with . Hence, is a reduced regular sequence. ∎

Once we have a reduced regular sequence, we can use [3, Thm. 4.8] to solve the system. We note that in [3] there is a general assumption that for all the intermediate ideals define absolutely irreducible varieties. However, the proof of [3, Thm. 4.8] does not require this assumption (this assumption is only required in algorithms for finding a rational point in [3, Section 6]).

Next, we describe the data structures used in [3] to represent polynomial systems. The algorithms take as input polynomials represented by division-free straight-line programs (DFSLP). A DFSLP defined over a field is a sequence of polynomials such that each polynomial is either a variable with , an element in , or , where 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 is said to be represented by a DFSLP if .

Theorem 6.

[3, Thm. 4.8] Let be a reduced regular sequence, where the polynomials are represented by a DFSLP with space and time . Set the following notation:

  • The integer is ;

  • For any real number , ;

  • Let be an integer larger than the degrees of the ideals .

Assume further that . There is a probabilistic Turing machine using space and time which takes such polynomial systems as input and which outputs a -geometric resolution of the algebraic set with probability at least .

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 bi-homogeneous system by a DFSLP.

Lemma 7.

Let be two positive integers. A polynomial system such that for all , and can be represented by a DFSLP with time and space .


There are monomials in such that and . We consider the DFSLP which starts by evaluating these monomials. This costs less than multiplications, using a naive algorithm. Then we multiply each of these monomials by the corresponding coefficients, and we sum. This costs multiplications and additions. ∎

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 be a regular sequence in and be such that for any , and . Then the degree of the ideal is at most


Moreover, this degree is bounded above by .


This is a direct consequence of [22, Prop. I.1] using, with the notation of [22, Prop. I.1], , , , , , , . Note that [22, Prop. I.1] is stated when the base field is , 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 , and hence the sum of the binomial coefficients is bounded above by . ∎

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

Proof of Proposition 3.

Set . First, we note that if is represented by a straight-line program over with space and time , then for any and any matrix with entries in , the sequence can be represented by a straight-line program over with space and time , where and . We consider the probabilistic Turing machine which performs the following steps:

  1. It chooses an matrix uniformly at random with entries in , with

    where , , . Using the inequalities , we get that ;

  2. It constructs the straight-line program representing with space and time ;

  3. It applies the probabilistic Turing machine from Theorem 6 to compute a geometric resolution of the algebraic set defined by ; By Theorem 6, it returns a geometric resolution provided that is a reduced regular sequence;

  4. It computes ;

  5. It computes and returns the geometric resolution