Simon Abelard

Pierrick Gaudry

Pierre-Jean Spaenlehauer

COUNTING POINTS ON GENUS-3 HYPERELLIPTIC CURVES WITH EXPLICIT REAL MULTIPLICATION

We propose a Las Vegas probabilistic algorithm to compute the zeta function of a genus-3 hyperelliptic curve defined over a finite field Fq, with explicit real multiplication by an order Z[η] in a totally real cubic field. Our main result states that this algorithm requires an expected number of O((log q) 6 ) bit-operations, where the constant in the O() depends on the ring Z[η] and on the degrees of polynomials representing the endomorphism η. As a proof-of-concept, we compute the zeta function of a curve defined over a 64-bit prime field, with explicit real multiplication by Z[2 cos(2π/7)].

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 problem of computing efficiently zeta functions of curves defined over finite fields has attracted a lot of attention, as its applications range from the construction of cryptographic curves to testing conjectures in number theory. We focus on the problem of computing the zeta function of a hyperelliptic curve C of genus 3 defined over a finite field F q using -adic methods, in the spirit of Schoof's algorithm and its generalizations [START_REF] Pila | Frobenius maps of Abelian varieties and finding roots of unity in finite fields[END_REF][START_REF] Huang | Counting points on curves over finite fields[END_REF][START_REF] Adleman | Counting points on curves and Abelian varieties over finite fields[END_REF]. Although these methods are polynomial with respect to log q, the exponents in the best known complexity bounds grow quickly with the genus. Another line of research is to use p-adic methods [START_REF] Kiran | 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][START_REF] Carls | A p-adic quasi-quadratic time point counting algorithm[END_REF][START_REF] Michael | An extension of Kedlaya's algorithm for hyperelliptic curves[END_REF], which are polynomial in the genus but exponential in the size of the characteristic of the underlying finite field. Variants of these methods [START_REF] Kedlaya | Computing L-series of hyperelliptic curves[END_REF][START_REF] Harvey | Counting points on hyperelliptic curves in average polynomial time[END_REF][START_REF] Harvey | Computing Hasse-Witt matrices of hyperelliptic curves in average polynomial time, II[END_REF] allow to count the points of a curve defined over the rationals modulo many primes in average polynomial time, which is especially relevant when experimenting with the Sato-Tate conjecture.

The aim of this paper is to show -both with theoretical proofs and practical experiments -that the complexity of -adic methods for genus-3 hyperelliptic curves can be dramatically decreased as soon as an explicitly computable non-integer endomorphism η ∈ End(Jac(C)) is known. More precisely, we say that a curve C has explicit real multiplication by Z[η] if the subring Z[η] ⊂ End(Jac(C)) is isomorphic to an order in a totally real cubic number field, and if we have explicit formulas describing η(P -∞) for some fixed base point ∞ and a generic point P of C. By explicit formulas, we mean polynomials (η (u) i (x, y)) i∈{0,1,2,3} and (η (v) i (x, y)) i∈{0,1,2,3} in F q [x, y], such that, when C is given in odd-degree Weierstrass form, the Mumford coordinates of η((x, y) -∞) are

3 i=0 η (u) i (x, y)X i , 2 i=0 η (v) i (x, y)/η (v)
3 (x, y) X i , where (x, y) is the generic point of the curve. In cases where C does not have an odd-degree Weierstrass model, we can work in an extension of degree at most 8 of the base field in order to ensure the existence of a rational Weierstrass point.

The influence of real multiplication on the complexity of point counting was investigated for genus 2 curves in [START_REF] Gaudry | Counting points on genus 2 curves with real multiplication[END_REF], where the authors decrease the complexity from O((log q) 8 ) [START_REF] Gaudry | Genus 2 point counting over prime fields[END_REF] to O((log q) 5 ). For genus 2 curves, another related active line of research is to mimic the improvement of Elkies and Atkin by using modular polynomials [START_REF] Ballentine | Isogenies for point counting on genus two hyperelliptic curves with maximal real multiplication[END_REF]. However, the main difficulty of this method is to precompute the modular polynomials, which are much larger than their genus 1 counterparts.

Our main result is the following theorem.

Theorem 1. Let C be a genus-3 hyperelliptic curve defined over a finite field F q having explicit real multiplication by Z[η], where η ∈ End(Jac(C)). We assume that C is given by an odd-degree Weierstrass equation Y 2 = f (X). The characteristic polynomial of the Frobenius endomorphism on the Jacobian of C can be computed with a Las Vegas probabilistic algorithm in expected time bounded by c (log q) 6 (log log q) k , where k is an absolute constant and c depends only on the degrees of the polynomials η In this paper, we use the notation O() as a shorthand for complexity statements hiding poly-logarithmic terms: the complexity in the theorem would be abbreviated O((log q) 6 ). We insist on the fact that all the O() and the O() notation used throughout the paper should be understood up to a multiplicative constant which may depend on the ring Z[η] and on the degrees of the polynomials η (u) i and η

(v)

i . There are natural families of curves for which these degrees are bounded by an absolute constant and for which Z[η] is fixed: reductions at primes (of good reduction) of a hyperelliptic curve with explicit RM defined over a number field.

As in Schoof's algorithm and its generalizations in [START_REF] Pila | Frobenius maps of Abelian varieties and finding roots of unity in finite fields[END_REF][START_REF] Huang | Counting points on curves over finite fields[END_REF][START_REF] Adleman | Counting points on curves and Abelian varieties over finite fields[END_REF], the -adic approach consists in computing the characteristic polynomial of the Frobenius endomorphism by computing its action on the -torsion of the Jacobian of the curve for sufficiently many . In order to prove the claimed complexity bound, we consider primes ∈ Z such that Z[η] splits as a product p 1 p 2 p 3 of prime ideals. Computing the kernels of endomorphisms α i in each p i provides us with an algebraic representation of the -torsion Jac(C)[ ] ⊂ Ker α 1 + Ker α 2 + Ker α 3 . Then, we compute from this representation integers a, b, c ∈ Z/ Z such that the sum π + π ∨ of the Frobenius endomorphism and its dual equals a + bη + cη 2 mod . Once enough modular information is known, the values of a, b, c such that π + π ∨ = a + bη + cη 2 are recovered via the Chinese Remainder Theorem and the coefficients of the characteristic polynomial of the Frobenius can be directly expressed in terms of a, b and c. In fact, in practice we do not have to restrict to split primes: any partial factorization of Z[η] provides some modular information on a, b, c mod . We give an example with a ramified prime in Section 7.1; but on the theoretical side, considering non-split primes does not improve the asymptotic complexity.

The cornerstone of the complexity analysis is the cost of the computation of the kernels of the endomorphisms. This is achieved by solving a polynomial system. Using resultant-based elimination techniques and degree bounds on Cantor's polynomials, we prove that we can solve these equations in time quadratic in the number of solutions, which leads to the claimed complexity bound. For practical computations, we replace the resultants by Gröbner bases and we retrieve modular information only for small to speed up an exponential collision search which can be massively run in parallel. Although using Gröbner basis seems to be more efficient in practice, we do not see any hope of proving with rigorous arguments that it is asymptotically competitive.

As a proof-of-concept, we have implemented our algorithm and we provide experimental results. In particular, we were able to compute the zeta function of a genus 3 hyperelliptic curve with explicit RM defined over F p with p = 2 64 -59. To our knowledge the largest genus-3 computation that had been achieved previously was the computation of the zeta function of a hyperelliptic curve defined over F p with p = 2 61 -1, done by Sutherland [START_REF] Sutherland | A generic approach to searching for Jacobians[END_REF] using generic group methods.

Examples of curves with RM are given by modular curves. For instance, the genus-3 curve

y 2 = x 7 + 3x 6 + 2x 5 -x 4 -2x 3 -2x 2 -x -1 is a quotient of X 0 (284)
and therefore has real multiplication by an element of Q[x]/(x 3 -3x -1). This follows from the properties of the Hecke operators as explained in [START_REF] Shimura | Introduction to the arithmetic theory of automorphic functions[END_REF]Chapter 7]. Based on this theory, algorithms for constructing such curves are explained in [START_REF] Frey | Arithmetic of modular curves and applications[END_REF]; however the explicit expression for the real endomorphism is not given. We expect that tracking the Hecke correspondences along their construction, and using techniques like in [START_REF] Paul Van Wamelen | Proving that a genus 2 curve has complex multiplication[END_REF] to reconstruct the rational fractions describing the real endomorphism could solve this question. In any case, these are only isolated points in the moduli space. Larger families are obtained from cyclotomic covering. This line of research has produced several families of hyperelliptic genus-3 curves having explicit RM by Z[2 cos(2π/7)]. In particular, explicit such families are given in [START_REF] Mestre | Familles de courbes hyperelliptiques à multiplications réelles[END_REF] and [START_REF] Tautz | Explicit hyperelliptic curves with real multiplication and permutation polynomials[END_REF], and explicit formulas for their RM endomorphism are obtained in [START_REF] Kohel | Efficiently computable endomorphisms for hyperelliptic curves[END_REF]. We use the 1-dimensional family of curves from [28, Theorem 1 with p = 7] for our experiments. Other families of genus-3 curves (but not necessarily hyperelliptic) with RM have been made explicit in [START_REF] Boyer | Variétés abéliennes et jacobiennes de courbes hyperelliptiques, en particulier à multiplication réelle ou complexe[END_REF]Chapter 2], following [START_REF] Ellenberg | Endomorphism algebras of Jacobians[END_REF]. We would like to point out that within the moduli space of complex polarized abelian varieties of dimension 3, those with RM by a fixed order in a cubic field form a moduli space of codimension 3 [4, Sec. 9.2]. Since Jacobians of hyperelliptic curves form a codimension 1 space, we would expect the moduli space of hyperelliptic curves of genus 3 with RM by a given cubic order to have dimension 2.

We finally briefly mention how our algorithm and analysis could be extended in several directions. First, the complexity analysis leads, with small modifications, to a point-counting algorithm for general genus-3 hyperelliptic curves (i.e. without RM) with complexity in O((log q) 14 ). Second, if the curve is not hyperelliptic, the main difficulty is to define analogues of Cantor's division polynomials and get bounds on their degrees. Without them, it is still possible to use an explicit group law to derive a polynomial system for the kernel of an endomorphism, but getting a proof for its degree would require to take another path than what we did. Still, the complexities with or without RM are expected to remain the same for plane quartics as for genus-3 hyperelliptic curves. Third, if we go to higher genus hyperelliptic curves with RM, the main difficulty to extend our approach is in the complexity estimate of the polynomial system solving, because resultant-based approaches are not competitive when the number of variables grows, and a tedious analysis like in [START_REF] Abelard | Improved complexity bounds for counting points on hyperelliptic curves[END_REF] seems to be necessary.

The article is organized as follows. Section 2 gives a bird-eye view of our algorithm, along with a complexity analysis relying on the technical results detailed in Sections 3 to 6. Practical experiments are presented in Section 7.

Overview of the algorithm

Let C be a genus-3 hyperelliptic curve over a finite field F q with explicit RM, and let η be the given explicit endomorphism. We denote by µ 0 , µ 1 , µ 2 the coefficients of the minimal polynomial T 3 + µ 2 T 2 + µ 1 T + µ 0 of η over Q.

2.1. Bounds. The characteristic polynomial of the Frobenius endomorphism π is of the form χ π (T ) = T 6 -σ 1 T 5 + σ 2 T 4 -σ 3 T 3 + qσ 2 T 2 -q 2 σ 1 T + q 3 , and Weil's bounds give

|σ 1 | ≤ 6 √ q, |σ 2 | ≤ 15q, |σ 3 | ≤ 20q 3/2 .
In order to take advantage of the explicit RM, we consider the endomorphism ψ = π+π ∨ , for which we can derive the real Weil's polynomial

χ ψ (T ) = T 3 -σ 1 T 2 + (σ 2 -3q)T -(σ 3 -2qσ 1 )
, which corresponds to the characteristic polynomial of ψ viewed as an element of the real subfield of End(Jac(C)) ⊗ Q. The endomorphism ψ belongs to the ring of integers of Q(η). The ring Z[η] might be a proper suborder of the ring of integers, so let us call ∆ its index, so that ψ can be written ψ = a + bη + cη 2 , where a, b, c are rationals with a denominator that divides ∆. By computing formally the characteristic polynomial of a + bη + cη 2 in Q(η) and by equating it with the expression for the real Weil's polynomial χ ψ (T ), we obtain a direct way to compute σ 1 , σ 2 and σ 3 in terms of a, b, c:

(1) σ1 = 3 a -b µ2 -2 c µ1 + c µ 2 2 , σ2 -3q = 3 a 2 -2 a b µ2 + 2 a c (µ 2 2 -2µ1) + b 2 µ1 + 3 b c µ0 -b c µ1 µ2 - c 2 (2 µ0 µ2 + µ 2 1 ) , σ3 -2qσ1 = a 3 -a 2 b µ2 + a 2 c (µ 2 2 -2µ1) + a b 2 µ1 + a b c (3 µ0 -µ1 µ2) + a c 2 (µ 2 1 -2 µ0 µ2) -b 3 µ0 + b 2 c µ0 µ2 -b c 2 µ0 µ1 + c 3 µ 2 0 .
In Section 4, it is shown that the coefficients a, b and c can be bounded in O( √ q).

More precisely, we denote by C abc a constant that depends only on η such that their absolute values are bounded by C abc √ q. Since these bounds are much smaller than the bounds for σ 1 , σ 2 , σ 3 , it makes sense to design an algorithm that reconstruct these coefficients of ψ instead of the coefficients of χ π as in the classical Schoof algorithm, and this is what we are going to do later on.

Another important bound that we need concerns the size of small elements that can be found in ideals of Z [η]. Let be a prime that splits completely in Z[η], so that we can write = p 1 p 2 p 3 , where the p i 's are distinct prime ideals of norm . In Section 5, it is shown that each p i contains a non-zero element α i = a i + b i η + c i η 2 , where a i , b i and c i are integers and are bounded in absolute value by O( 1/3 ).

2.2.

Algorithms. The general RM point counting algorithm is Algorithm 1. We give a description of it, allowing some black-box primitives that will be detailed in dedicated sections. As mentioned above, we will work with the a, b, c coefficients of the ψ endomorphism. More precisely, we compute their values modulo sufficiently many completely split primes until we can deduce their values from the bounds of Lemma 5 by the Chinese Remainder Theorem, taking into account their potential denominator ∆. Then the coefficients of χ π are deduced by Equations [START_REF] Abelard | Improved complexity bounds for counting points on hyperelliptic curves[END_REF].

We now explain how the algorithm works for a given split . First its decomposition as a product of prime ideals Z[η] = p 1 p 2 p 3 is computed, and for each prime ideal p i , a non-zero element α i of p i is found with a small representation α i = a i + b i η + c i η 2 as in Lemma 6. In fact, p i is not necessarily principal and α i need not generate p i . The kernel of α i is denoted by J[α i ] and it contains a subgroup G i isomorphic to Z/ Z × Z/ Z, since the norm of α i is a multiple of . The two-element representation ( , η -λ i ) of the ideal p i implies that λ i is an eigenvalue of η regarded as an endomorphism of J[ ] ∼ = (Z/ Z) 6 .

On G i ⊂ J[α i ], the endomorphism η acts as the multiplication by λ i . Therefore, ψ = a + bη + cη 2 also acts as a scalar multiplication on this 2-dimensional space, and we write k i ∈ Z/ Z the corresponding eigenvalue: for any D i in G i , we have ψ(D i ) = k i D i . On the other hand, from the definition of ψ, it follows that ψπ = π 2 + q. Therefore, if such a D i is known, we can test which value of

k i ∈ Z/ Z satisfies (2) k i π(D i ) = π 2 (D i ) + qD i .
Since is a prime and D i is of order exactly , this is also the case for π(D i ). Finding k i can then be seen as a discrete logarithm problem in the subgroup of order generated by π(D i ); hence the solution is unique. Equating the two expressions for ψ, we get explicit relations between a, b, c modulo :

a + bλ i + cλ 2 i ≡ k i mod .
Therefore we have a linear system of three equations in three unknowns, the determinant of which is the Vandermonde determinant of the λ i , which are distinct by hypothesis. Hence the system can be solved and it has a unique solution modulo .

Data: q an odd prime power, and f ∈ F q [X] a monic squarefree polynomial of degree 7 such that the curve

Y 2 = f (X) has explicit RM by Z[η]. Result: The characteristic polynomial χ π ∈ Z[T ] of the Frobenius endomorphism on the Jacobian J of the curve. R ← 1; while R ≤ 2 ∆ C abc √ q + 1 do
Pick the next prime that satisfies conditions (C1) to (C4); Compute the ideal decomposition a,b,c) using the Chinese Remainder Theorem ; Deduce χ π from Equations (1).

Z[η] = p 1 p 2 p 3 , corresponding to the eigenvalues λ 1 , λ 2 , λ 3 of η in J[ ] ; for i ← 1 to 3 do Compute a small element α i of p i as in Lemma 6; Compute a non-zero element D i of order in J[α i ] ; Find the unique k i ∈ Z/ Z such that k i π(D i ) = π 2 (D i ) + qD i ; end Find the unique triple (a, b, c) in (Z/ Z) 3 such that a + bλ i + cλ 2 i = k i , for i in {1, 2, 3} ; R ← R • ; end Reconstruct (

Algorithm 1: Overview of our RM point-counting algorithm

It remains to show how to construct a divisor D i in G i , i.e. an element of order in the kernel J[α i ]. Since an explicit expression of η as an endomorphism of the Jacobian of C is known, an explicit expression can be deduced for α i , using the explicit group law. The coordinates of the elements of this kernel are solutions of a polynomial system that can be directly derived from this expression of α i . Using standard techniques, it is possible to find the solutions of this system in a finite extension of the base field (of degree bounded by the degree of the ideal generated by the system, i.e. in O( 2)), from which divisors in J[α i ] can be constructed. Multiplying by the appropriate cofactor, we can reach all the elements of G i ; but we stop as soon as we get a non-trivial one.

We summarize the conditions that must be satisfied by the primes that we work with:

(C1) must be different from the characteristic of the base field; (C2) must be coprime to the discriminant of the minimal polynomial of η; (C3) there must exist α i ∈ p i as in Lemma 6 with norm non-divisible by 3

for i ∈ {1, 2, 3}; (C4) the ideal Z[η] must split completely.
The first 3 conditions eliminate only a finite number of 's that depends only on η, while the last one eliminates a constant proportion. The condition (C3) implies that there is a unique subgroup G i of order 2 in J[α i ] (our description of the algorithm could actually be adapted to handle the cases where this is not true).

Algorithm 1 is a very natural extension of the one described in [START_REF] Gaudry | Counting points on genus 2 curves with real multiplication[END_REF] for genus 2 curves with RM. Already in [START_REF] Gaudry | Counting points on genus 2 curves with real multiplication[END_REF], the action of the real endomorphism ψ = π + π ∨ is studied on subspaces J[p i ] of the -torsion, and the corresponding eigenvalues are collected and used to reconstruct information modulo . In genus 3, we have 3 such 2-dimensional subspaces and eigenvalues to compute and recombine instead of 2 in genus 2. The main differences between the present work and [START_REF] Gaudry | Counting points on genus 2 curves with real multiplication[END_REF] are the way the -torsion elements are constructed with polynomial systems and the bounds on the coefficients of ψ. In both cases, going from dimension 2 to 3 is not immediate.

Complexity analysis.

The field Q(η) is of degree 3, so its Galois group has order at most 6 and by Chebotarev's density theorem the density of primes that split completely is at least 1/6. Therefore the main loop is done O(log q/ log log q) times, with primes that are in O(log q). All the steps that take place in the number field take a negligible time. For instance, a small generator like in Lemma 6 can be found by exhaustive search: only O( ) trials are needed since we are searching over all elements of the form a + bη + cη 2 , with |a|, |b|, |c| in O( 1/3 ).

The bottleneck of the algorithm is the computation of a non-zero element of order in the kernel J[α i ] of α i . This part will be treated in detail in Section 3, where it is shown to be feasible in O( 4 ) operations in F q . The output is a divisor D i of order in J[α i ] that is defined over an extension field F q δ , where δ is in O( 2 ).

In order to check Equation ( 2), we first need to compute π(D i ) and π 2 (D i ) which amounts to raising the coordinates to the q-th power. The cost is in O( 2 log q) operations in F q . Then, each Jacobian operation in the group generated by π(D i ) costs O( 2 ) operations in the base field, and we need O( √ ) of them to solve the discrete logarithm problem given by Equation [START_REF] Adleman | Counting points on curves and Abelian varieties over finite fields[END_REF]. The overall cost of finding k i , once D i is known is therefore O( 2 ( √ + log q)) operations in F q . Finally, the amount of work performed for each is O( 2 ( 2 + log q)) operations in the base field F q . Summing up for all the primes, and taking into account the cost of the operations in F q , we obtain a global bit-complexity of O((log q) 6 ).

Computing kernels of endomorphisms

3.1. Modelling the kernel computation by a polynomial system. Let α be an explicit endomorphism of degree O( 2 ) on the Jacobian of C, which satisfies the properties of Lemma 6. In particular, α vanishes on a subspace of J[ ]. We want to compute a triangular polynomial system that describes the kernel J[α] of α. This will provide us with a nice description of a subgroup of the -torsion on which we will be able to test the action of ψ = π + π ∨ and deduce a, b, c such that ψ = a + bη + cη 2 mod .

We first model J[α] by a system of polynomial equations that we will then put in triangular form. To do so, we consider a generic divisor D = P 1 + P 2 + P 3 -3∞, where P i is an affine point of C of coordinates (x i , y i ). We then write α(D) = 0, i.e α(P 1 -∞) + α(P 2 -∞) = -α(P 3 -∞). Generically, we expect each α(P i -∞) to be of weight 3, and we write u i , v i for its Mumford form. We derive our equations by computing the Mumford form u 12 , v 12 of α(P 1 -∞) + α(P 2 -∞) and then writing coefficient-wise the conditions u 12 = u 3 and v 12 = -v 3 . The case where the genericity conditions are not satisfied is discussed at the end of the section.

Similarly to the Schoof-Pila algorithm, we define polynomials -which are equivalent to Cantor's division polynomials -by the formulas

u 12 (X) = X 3 + 2 i=0 di (x 1 , x 2 , y 1 , y 2 ) d3 (x 1 , x 2 ) X i , v 12 (X) = 2 i=0 ẽi (x 1 , x 2 , y 1 , y 2 ) ẽ3 (x 1 , x 2 ) X i , u 3 (X) = X 3 + 2 i=0 d i (x 3 ) d 3 (x 3 ) X i , v 3 (X) = y 3 2 i=0 e i (x 3 ) e 3 (x 3 ) X i . Lemma 2.
For any i ∈ {1, 2, 3}, the degrees of di , ẽi , d i and e i are in O( 2/3 ).

Proof. Let us first remark that the di 's and ẽi 's are obtained after adding two divisors u 1 , v 1 and u 2 , v 2 such that the coefficients of the u i and v i are respectively the d j /d 3 and y i e j /e 3 evaluated at x i . Thus, since this application of the group law involves a number of operations that is bounded independently of and q, the degree stays within a constant multiplicative factor, which is captured by the O(). Therefore it is enough to prove the result for the d i 's and e i 's.

Since the endomorphism α satisfies the properties of Lemma 6, it is a linear combination of 1, η and η 2 with coefficients of size O( 1/3 ). Using the same argument about the group law, we can further reduce our proof to the case where α = nη k , with k ∈ {0, 1, 2} and n an integer in O( 1/3 ). But once again, η k does not depend on so that, provided we can prove that Cantor's n-division polynomials have degrees in O(n 2 ), we have proven that nη k (P -∞) = η k (n(P -∞)) have coefficients whose degrees are in O(n 2 ), and then so does α(P -∞). This quadratic bound on the degrees of Cantor's division polynomials is proven in Lemma 8 of Section 6 and the result follows.

3.2. Solving the system with resultants. Typical tools for solving a polynomial system are the F4 algorithm, methods based on geometric resolution, or homotopy techniques. To obtain reasonable complexity bounds, they all require some knowledge of the properties of the system, and this might be hard to prove. Since we have a system in essentially 3 variables (in fact, there are six variables x 1 , x 2 , x 3 , y 1 , y 2 , y 3 , but the y i variables can be directly eliminated by using the equation defining the curve), we prefer to stick to an approach based on resultants. It ends up having a complexity that is quasi-quadratic in the degree of the ideal, which is the best that can be hoped for anyway for all of the advanced techniques, and the complexity analysis requires only elementary tools. A complication that can occur with resultants is that Res x (f, g) is identically zero when f and g have a nonconstant GCD. This is not a problem in our case since we can divide polynomials f and g by their GCD, by factoring them at the cost of O(max(deg(f ), deg(g)) ω ) field operationswhere ω ≤ 3 is the exponent of linear algebra -using the bivariate recombination methods in [START_REF] Bostan | Complexity issues in bivariate polynomial factorization[END_REF] (the trivariate case can be reduced to the bivariate case by using the techniques in [START_REF] Zippel | Effective polynomial computation[END_REF]Sec. 21.2]). In what follows, the complexities of computing the resultants are larger than O(max(deg(f ), deg(g)) ω ), so we can forget about this complication. We also note that since the system is symmetric with respect to x 1 and x 2 , it may be possible to decrease the degrees by rewriting the system in terms of elementary symmetric polynomials in x 1 and x 2 ; however, we do not consider this symmetrization process in the analysis since it may only win a constant factor in the complexity.

Following our modelling, the equality of the u-coordinates gives three equations

(3) ∀i ∈ {0, 1, 2}, di (x 1 , x 2 , y 1 , y 2 )d 3 (x 3 ) = d3 (x 1 , x 2 )d i (x 3 ),
of degree O( 2/3 ) in the x i 's. By computing resultants with the equations y 2 i = f (x i ), we derive three equations E i (x 1 , x 2 , x 3 ) = 0 whose degrees are still in O( 2/3 ).

We then eliminate x 1 by computing 3 trivariate resultants R i (between the two equations E j with j = i). We get three equations R i (x 2 , x 3 ) = 0 of degrees O( 4/3 ) within a complexity in O( 10/3 ) field operations, as proven in Lemma 4 below.

Then, we compute bivariate resultants S i (between the two equations R j with j = i) to eliminate x 2 . From Lemma 3, we get three univariate equations S i (x 3 ) = 0 of degree bounded by O( 8/3 ) for a complexity in O( 4 ) field operations. And we compute the polynomial S(x 3 ) as the GCD of the S i (x 3 ), which belongs to the ideal defined by our original system.

The bound on the degree of S is much larger than 2 -1, the expected degree of the kernel. Although we can expect the actual degree to be in O( 2 ), we need to add the constraints coming from the v-coordinates to be able to prove it.

The polynomial system coming from v 12 = -v 3 has the same characteristics as the one coming from the u-coordinates. Therefore, we can proceed in a similar way and deduce, at a cost of O( 4 ) operations another univariate polynomial S(x 3 ) belonging to the ideal. Now, since all the original equations have been taken into account all common roots of S and S will correspond to a solution of the original system for which we know that there are O( 2 ) solutions. Therefore taking the squarefree part of the GCD of S and S yields a polynomial of degree O( 2).

This univariate polynomial can be factored at a cost of O( 4) operations in F q with standard algorithms [START_REF] Von | Modern computer algebra[END_REF] (there exist asymptotically faster algorithms, but we already fit in our target complexity). We then deal with each irreducible factor in turn, until one is found that leads to a genuine solution of the original system. Let δ be the degree of such an irreducible factor φ(x 3 ). In the field extension F q δ = F q [x 3 ]/φ(x 3 ), we have by construction a root x 3 of φ. We then solve again the original polynomial system where x 3 is instantiated with this root. This system is bivariate in x 1 and x 2 and there are O(1) solutions, that possibly live in another finite extension F q δ of F q δ . Since the degrees of the bivariate polynomials are in O( 2/3 ), by Lemma 3, this system solving costs O( 2 ) operations in F q δ . A solution obtained in this way must be checked, because it could come from a vanishing denominator that has been cleared when constructing the system or from non-generic situations. But given a set of candidate coordinates for a D i element of J[α i ], it is cheap to check that this is indeed an element of the Jacobian and that it is killed by α i . Also, if α i is not a generator of p i , it is necessary to check the order of D i : if this is a multiple of , then multiplying D i by the cofactor gives an order-element. But it is also possible to get an unlucky element of small order coprime to , and then we have to take another solution of the system.

Since an operation in F q δ requires a number of operations in F q that is quasilinear in δ, and since the sum of all the degrees δ of the irreducible factors of GCD(S, S) is in O( 2), the amortized cost is O( 4 ) operations in F q to deduce a divisor D i in J[α i ].

3.3.

Complexity of bi-and tri-variate resultants. In this section, the algorithms work by evaluation / interpolation, which requires to have enough elements in the base field. Were it not the case, we simply take a field extension F q δ of F q , that will add a factor O(δ) to the complexity. The complexity of the algorithms will be polynomial in the number of evaluation points, therefore, the final complexity will be logarithmic in δ, so that the cost of taking a field extension will be hidden in the O() notation. We will therefore not mention this potential complication further.

Another difficulty is that an evaluation / interpolation strategy assumes that the points of evaluation are generic enough, so that all the degrees after evaluation are generic. This is again guaranteed by taking a large enough base field. Still, the algorithm remains a Monte-Carlo one. However, the ultimate goal is to construct kernel elements, which is an easily verified property. Turning this into a Las Vegas algorithm can therefore be done with standard techniques. Lemma 3. [30, Thm. 6.22 and Cor. 11.21] Let P (x, y) and Q(x, y) be two polynomials whose degrees in x and y are bounded by d x and d y respectively. Then, R(y) = Res x (P, Q) can be computed in O(d 2

x d y ) field operations, and the degree of R is bounded by 2d x d y .

Lemma 4. Let P (x, y, z) and Q(x, y, z) be two polynomials whose degrees in each variable are bounded by d. Then, R(y, z) = Res x (P, Q) can be computed in O(d 5 ) field operations, and the degree of R in each variable is bounded by 2d 2 .

Proof. The Sylvester matrix has at most 2d columns and its entries are bivariate polynomials whose degrees in y and z are bounded by d. Thus, its determinant is a polynomial whose degrees in y and z are bounded by 2d 2 .

We first perform a Kronecker substitution by considering P (x, y) = P (x, y, y 2d 2 +1 ) and Q(x, y) = Q(x, y, y 2d 2 +1 ), which are polynomials of degrees ≤ d in x and ≤ 2d 3 + d in y. Note that the choice to replace z by y 2d 2 +1 is made to be able to invert the Kronecker substitution after the resultant computation.

Next, we compute R(y) = Res x ( P (x, y), Q(x, y)). By Lemma 3, it is a univariate polynomial of degree at most 4d 4 + 2d 2 and can be computed in O(d 5 ) operations. We can then invert the Kronecker substitution to get R(y, z), which can be done in time linear in the number of monomials, that is in O(d 4 ).

3.4.

Non-generic situations. Our analysis assumes in the first place that the -torsion elements are generic in a rather strong sense, see e.g. [START_REF] Abelard | Improved complexity bounds for counting points on hyperelliptic curves[END_REF]Def. 11] for details. This is expected to be the case with overwhelming probability, when the base field is large enough and the curve is taken at random in a large family. However, to obtain a proven complexity we must also consider the cases where there exist -torsion elements that are non-generic. We follow the strategy of [START_REF] Abelard | Improved complexity bounds for counting points on hyperelliptic curves[END_REF] where another polynomial system is designed and solved for each non-generic situation, for instance the fact that an -torsion divisor is of weight less than 3, or that some points involved in the modelling are not distinct while they generically are. We do not give all the details, but the number of polynomial systems to consider is bounded by a constant, and each of these polynomial systems describes a situation that is smaller than the generic one in the sense that it has either less variables or a lower degree, so that the complexity bound is maintained.

Bounds on the coefficients of ψ

The system of equations ( 1) giving σ 1 , σ 2 and σ 3 in terms of a, b, c is homogeneous if we put weight 1/2 to a, b, c and σ 1 , weight 1 to q and σ 2 , weight 3/2 to σ 3 , and weight 0 to µ 0 , µ 1 , and µ 2 so any polynomial in a reduced Gröbner basis of the corresponding ideal will have the same property. Computing such a Gröbner basis with the lexicographical ordering a > b > c > σ 1 > σ 2 > σ 3 > µ 0 > µ 1 > µ 2 > q (we did this computation with the Magma V2.23-4 software), we get a polynomial Ψ c of degree 6 in c that does not involve a or b, and which has the following form:

Ψc(q, c, σ1, σ2, σ3, µ0, µ1, µ2) = D(µ0, µ1, µ2) 3 c 6 + 5 i=0 ψ (i) c (q, σ1, σ2, σ3, µ0, µ1, µ2) c i ,
where

D(µ 0 , µ 1 , µ 2 ) = -27 µ 2 0 + 18 µ 0 µ 1 µ 2 -4 µ 0 µ 3 2 -4 µ 3 1 + µ 2 1 µ 2 2 is the discriminant of the polynomial T 3 + µ 2 T 2 + µ 1 T + µ 0 .
By computing Gröbner bases for other lexicographical orderings (with a > c > b > σ 1 > σ 2 > σ 3 > µ 0 > µ 1 > µ 2 > q and b > c > a > σ 1 > σ 2 > σ 3 > µ 0 > µ 1 > µ 2 > q respectively), we obtain that polynomials of the following form also belong to the ideal generated by the polynomials in the system of equations (1):

Ψ b (q, b, σ1, σ2, σ3, µ0, µ1, µ2) = D(µ0, µ1, µ2) 3 b 6 + 5 i=0 ψ (i)
b (q, σ1, σ2, σ3, µ0, µ1, µ2) b i , Ψa(q, a, σ1, σ2, σ3, µ0, µ1, µ2) = D(µ0, µ1, µ2) 3 a 6 + 5 i=0 ψ (i) a (q, σ1, σ2, σ3, µ0, µ1, µ2) a i .

The polynomials ψ

(i) a , ψ (i)
b and ψ

(i)
c are homogeneous of weighted degree 3 -i/2 with respect to the grading given above.

Lemma 5. The absolute values of the coefficients a, b, c of ψ = a + bη + cη 2 are bounded above by O(q 1/2 ).

Proof. First, we consider the equation

Ψ c = 0. We write c = c q 1/2 , σ 1 = σ 1 q 1/2 , σ 2 = σ 2 q, σ 3 = σ 3 q 3/2 . Since ψ (i)
c is homogeneous and has weighted degree 3 -i/2, there is a polynomial θ

(i) c ( σ 1 , σ 2 , σ 3 , µ 0 , µ 1 , µ 2 ) such that (4) ψ (i) c (q, σ 1 , σ 2 , σ 3 , µ 0 , µ 1 , µ 2 ) • c i = q 3 c i θ (i) c ( σ 1 , σ 2 , σ 3 , µ 0 , µ 1 , µ 2 ). Weil's bounds imply that | σ i | = O(1) for i ∈ {1, 2, 3}.
Therefore, for all i ∈ {0, . . . , 5}, we obtain that |θ 3 is the minimal polynomial of a totally real algebraic number, the discriminant D(µ 0 , µ 1 , µ 2 ) must be nonzero. Equations Ψ c = 0 and (4) imply the following inequality:

(i) c ( σ 1 , σ 2 , σ 3 , µ 0 , µ 1 , µ 2 )| = O(1). For fixed µ 0 , µ 1 , µ 2 ∈ Q such that µ 0 +µ 1 T +µ 2 T 2 +T
| c| 6 - 5 i=0 |θ (i) c ( σ 1 , σ 2 , σ 3 , µ 0 , µ 1 , µ 2 )| |D(µ 0 , µ 1 , µ 2 )| 3 | c| i ≤ 0.
Then | c| must be smaller or equal to the largest root of this polynomial inequality, which can itself be bounded, for instance, with Cauchy's bound

| c| ≤ 1 + max 0≤i≤5 |θ (i) c ( σ 1 , σ 2 , σ 3 , µ 0 , µ 1 , µ 2 )| |D(µ 0 , µ 1 , µ 2 )| 3 ,
which shows that | c| = O(1), and hence |c| = O(q 1/2 ). The proof for the bounds on |a| and |b| are similar, using the equations Ψ a = 0 and Ψ b = 0.

Small elements in ideals of Z[η]

We first recall that we consider only primes that do not divide the discriminant of the minimal polynomial of η (Condition (C2)). Hence, if Z[η] is not the maximal order of Q(η), this has no consequence on the factorization properties of . Lemma 6. For any prime that splits completely in Z[η], each prime ideal p i above contains a non-zero element α i of the form

α i = a i + b i η + c i η 2 , where |a i |, |b i | and |c i | are integers in O( 1/3 ), and the norm of α i is in O( ).
Proof. The coefficients of the elements of the ideal p i represented by polynomials in η form a lattice. Applying Minkowski's bound to this lattice, we obtain the existence of a non-zero element

α i = a i + b i η + c i η 2 in p i for which the L 2 -norm of (a i , b i , c i ) is in O( 1/3
). From this bound on the L 2 -norm, we derive a bound on the L ∞ -norm, and finally on the norm of α i as an algebraic number. At each step, the constant hidden in the O() gets worse but still depends only on Z[η].

For any given η, it is not difficult to make the constants in the O() fully explicit. We do it in the particular case of Z[η 7 ], with η 7 = 2 cos(2π/7), which is the RM used in our practical experiments. Since Z[η 7 ] is a principal ring, a more direct approach leads to bounds for a generator that are tighter than what would be obtained by a naive application of the previous lemma.

Lemma 7. Every ideal p i of norm in Z[η 7 ] has a generator α i of the form a i + b i η 7 + c i η 2 7 , where a i , b i , c i ∈ Z satisfy |a i | < 2.415 • 1/3 ; |b i | < 1.850 • 1/3 ; |c i | < 1.764 • 1/3 .
Proof. By abuse of notation, we identify Q(η 7 ) with the algebraic number field

Q[X]/(X 3 + X 2 -2X -1)
and we let σ 1 , σ 2 , σ 3 be the three real embeddings of Q(η 7 ) in R. Let 1 = 1 -η 2 7 and 2 = 1 + η 7 be a pair of fundamental units, and let µ i be a generator of p i . The logarithmic embedding ϕ : x → (log|σ 1 (x)|, log|σ 2 (x)|, log|σ 3 (x)|) sends the set of generators of p i to the lattice generated by ϕ( 1 ) and ϕ( 2 ) translated by ϕ(µ i ). Solving a CVP for the projection of ϕ(µ i ) on the plane where the 3 coordinates sum-up to zero, we deduce a unit ξ i such that α i = ξ i µ i is a generator whose real embeddings are bounded by

|σ 1 (α i )| ≤ 2.247 • 1/3 , |σ 2 (α i )| ≤ 1.803 • 1/3 , |σ 3 (α i )| ≤ 2.247 • 1/3 . Writing α i = a i + b i η 7 + c i η 2
7 , the real embeddings can also be expressed as

(σ 1 (α i ), σ 2 (α i ), σ 3 (α i )) T = V • (a i , b i , c i ) T ,
where V is the Vandermonde matrix of (σ 1 (η 7 ), σ 2 (η 7 ), σ 3 (η 7 )). A numerical evaluation of its inverse allows to translate the bounds on σ 1 (α i ), σ 2 (α i ), σ 3 (α i ) into the claimed bounds on a i , b i , c i .

Bounding the degrees of Cantor's division polynomials in genus 3

The purpose of this section is to prove the following lemma on the Cantor's division polynomials, which are explicit formulas for the endomorphism corresponding to scalar multiplication [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF].

Lemma 8. In genus 3, the degrees of Cantor's -division polynomials are bounded by O( 2 ).

In [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF], there are exact formulas for the degrees of the leading and the constant coefficients d 3 and d 0 . However, there is no formula or bounds for the degrees of the other coefficients of the -division polynomials. Still, our proof strongly relies on [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF] and we do not try to make it standalone: we assume that the reader is familiar with this article and all references to expressions, propositions or definitions in this proof are taken from this paper.

For a polynomial P whose coefficients are themselves univariate polynomials, we denote by maxdeg(P ) the maximum of the degrees of its coefficients.

We first prove a bound on the degrees of the coefficients of the quantities α r and γ r defined in [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF], from which the wanted bounds will follow. The key tools are the recurrence formulas (8.31) and (8.33) that relate quantities at index r to quantities at index around r/2, in a similar fashion as for the division polynomials of elliptic curves. More precisely, the following lemma shows that when the index r is (roughly) doubled, maxdeg α r and maxdeg γ r are roughly multiplied by 4, which leads to the expected quadratic growth. Lemma 9. Let ≥ 12, and assume that for all i ≤ ( + 9)/2 the degrees maxdeg α i and maxdeg γ i are bounded by C, then maxdeg α and maxdeg γ are bounded by 4C + 36 + 108.

Proof. We first deal with the bound on maxdeg γ . Let us consider r and s around /2 such that = r + s -5: we take either r = s -3 = /2 + 1 if is even, or r = s -4 = ( + 1)/2 otherwise.

From Equations (8.30) and (8.31), the degree of γ [h]ψ s-r ψ r-2 ψ s-2 ψ r-1 ψ s-1 is that of the determinant of the matrix E rs [h] defined by:

E rs [h] =     α r-3 α s [0] α r-3 α s [1] ψ r-3 ψ s γ r-3 γ s [h] α r-2 α s-1 [0] α r-2 α s-1 [1] ψ r-2 ψ s-1 γ r-2 γ s-1 [h] α r-1 α s-2 [0] α r-1 α s-2 [1] ψ r-1 ψ s-2 γ r-1 γ s-2 [h] α r α s-3 [0] α r α s-3 [1] ψ r ψ s-3 γ r γ s-3 [h]     .
Therefore we have an expression for the degrees of the coefficients of γ in terms of objects at index around r and s:

deg γ [h] ≤ deg det E rs [h] -deg(ψ r-2 ψ s-2 ψ r-1 ψ s-1 ).
In this last formula, the factor ψ s-r has been omitted, because s -r is either 3 or 4, and by (8.17) this has non-negative degree in any case. Thus, we simply bounded it below by 0 in the previous inequality. Before entering a more detailed analysis, we use Equation (8.8) to rewrite the first column with expressions for which we have exact formulas for the degree:

E rs [h] =     ψ r-4 ψ s-1 α r-3 α s [1] ψ r-3 ψ s γ r-3 γ s [h] ψ r-3 ψ s-2 α r-2 α s-1 [1] ψ r-2 ψ s-1 γ r-2 γ s-1 [h] ψ r-2 ψ s-3 α r-1 α s-2 [1] ψ r-1 ψ s-2 γ r-1 γ s-2 [h] ψ r-1 ψ s-4 α r α s-3 [1] ψ r ψ s-3 γ r γ s-3 [h]     .
The determinant of E rs [h] is the sum of products of 4 ψ factors and 4 α or γ factors. The degrees of the former are explicitly known, while by hypothesis we have upper bounds on the latter, since all the indices are at most ( + 9)/2. We can then deduce an upper bound on the degree of this determinant. All the ψ i have indices with i in the range [r -4, s] (remember that r ≤ s), and since their degrees increases with the indices, we can upper bound the degree of the products of the four ψ factors by 4 deg ψ s . Therefore we have

deg det E rs [h] ≤ 4(deg ψ s + C).
In order to deduce an upper bound on maxdeg γ , it remains to get a lower bound on the degree of the deg(ψ r-2 ψ s-2 ψ r-1 ψ s-1 ) term, and again by monotonicity of the degree in the index, we lower bound it by 4 deg ψ r-2 . So finally, we get

maxdeg γ ≤ 4C + (deg ψ 4 s -deg ψ 4 r-2
). Using (8.16) and (8.17), we deduce that for all k, we have deg(ψ 2 k ) = 3(k 2 -9) and substituting this value and the expression of r -2 and s in term of , we obtain

deg ψ 4 s -deg ψ 4 r-2 =
30 + 90 if is even, 36 + 108 if is odd, and the result follows for maxdeg γ .

The proof for maxdeg α follows the same line. Using the matrix F rs [h] defined in (8.32) in a similar way as we used the matrix E rs [h] and with the help of the formula (8.33), we end up with the following bounds maxdeg α ≤ 4C + 30 -30 if is even, 4C + 36 -36 if is odd, which are stricter than our target. Finally, the bound ≥ 12 is necessary to ensure that the quantities r and s are at least 5, as required in [START_REF] Cantor | On the analogue of the division polynomials for hyperelliptic curves[END_REF] to apply the formulas (8.31) and (8.33).

We can now finish the proof of Lemma 8. We define two sequences ( i ) i≥0 and (C i ) i≥0 as follows: let 0 = 12 and let C 0 be a bound on the degrees of the coefficients of all the α i and γ i for i ≤ 0 . Then for all i ≥ 1, we define the sequences inductively by i+1 = 2 i -9 C i+1 = 4C i + 36 i+1 + 108. By Lemma 9, for all i and all ≤ i , the degrees maxdeg α and maxdeg γ are bounded by C i . The expression i = ( 0 -9)2 i +9 = 3•2 i +9 can be derived directly from the definition and substituted in the recurrence formula of C i+1 to get C i+1 = 4C i + 216 • 2 i + 432. This recurrence can be solved by setting Γ i = C i + 108 • 2 i + 144, so that Γ i+1 = 4 Γ i , and we obtain C i = (C 0 + 252) 4 i -108 • 2 i -144. Finally, for any , we select the smallest i such that ≤ i . This value of i is log 2 (( -9)/3) . The corresponding bound for maxdeg α and maxdeg γ is then C i , which grows like O( 2 ) (and we remark that the effect of the ceiling can make the constant hidden in the O() expression grow by a factor at most 3).

Using the expression (8.10), we have maxdeg δ ≤ maxdeg α + maxdeg γ , and therefore the bound O( 2 ) also applies to the degrees of the coefficients of δ . And using the formula (8.13), the same holds as well for the coefficients of /y. This concludes the proof of Lemma 8.

Experimental results

In order to evaluate the practicality of our algorithm, we have tested it on one of the families of genus-3 hyperelliptic curves having explicit RM given in [28, Theorem 1]. Formulas for their RM endomorphisms are described in [START_REF] Kohel | Efficiently computable endomorphisms for hyperelliptic curves[END_REF]: for t = ±2, the curve C t with equation y 2 = x 7 -7x 5 + 14x 3 -7x + t, admits an endomorphism given in Mumford representation by

η 7 (x, y) = X 2 + 11 xX/2 + x 2 -16/9, y .
The fact that this expression has degree 2 while one would generically expect a degree 3 is no accident: it comes from the construction in [START_REF] Tautz | Explicit hyperelliptic curves with real multiplication and permutation polynomials[END_REF] of the endomorphism as a sum of two automorphisms on a double cover of the curve. We have η 3 7 + η 2 7 -2η 7 -1 = 0, so that the ring Z[η 7 ] is isomorphic to the ring of integers Z[2 cos(2π/7)] of the real subfield of the cyclotomic field Q(e 2iπ/7 ). All the numerical data in this section have been obtained for the parameter t = 42, on the prime field F p with p = 2 64 -59.

In our practical computations, the main differences with the theoretical description are the following: we use Gröbner basis algorithms instead of resultants, we consider also small non-split primes and small powers, and we finish the computation with a parallel collision search. The source code for our experiments is available at https://members.loria.fr/SAbelard/RMg3.tgz.

7.1.

Computing modular information with Gröbner basis. Although the polynomial system resolution using resultants has a complexity in O( 4 ), the real cost for small values of is already pretty large. In the resolution method described in Section 3.2, each bivariate resultant is computed by evaluation / interpolation and hence requires the computation of many univariate resultants. We illustrate this by counting the number of univariate resultants to perform and their degrees for the main step of the resolution (the part that reaches the peak complexity). We also measure the cost of such resultant computations using the NTL 10.5.0 and FLINT 2.5.2 libraries, both linked against GMP 6, when the base field is F 2 64 -59 . These costs do not include the evaluation / interpolation steps which might also be problematic for large instances, because they are hard to parallelize.

#res Deg Cost (NTL) Cost (FLINT) 13 525M 16,000 1,850 days 735 days 29 12.8G 80,000 310,000 days 190,000 days We were more successful with the direct approach using Gröbner bases that we now describe. For computing the kernel of a given endomorphism, we computed a Gröbner basis of the system (3) with some small modifications. First, we observe that the only occurrences of y 1 and y 2 are within the monomial y 1 y 2 . Consequently, we can remove one variable by replacing each occurrence of y 1 y 2 by a fresh variable y. Next, we need to make the system 0-dimensional by encoding the fact that d 3 (x 3 ) and d 3 (x 1 , x 2 ) are nonzero. This is done by introducing another fresh variable t and by adding the polynomial S(x 1 , x 2 , x 3 )t -1 to the system, where S(x 1 , x 2 , x 3 ) is the squarefree part of d 3 (x 3 ) d 3 (x 1 , x 2 ). Finally, since each polynomial is symmetric with respect to the transposition of the variables x 1 and x 2 , we can rewrite the equations using the symmetric polynomials s 1 = x 1 + x 2 and s 2 = x 1 x 2 . This divides by two the degree in x 1 and x 2 of the equations. We end-up with a system in 5 variables.

The whole construction can be slightly modified to compute the pre-image of a given divisor by the endomorphism: to model α(D) = Q, we write D = P 1 + P 2 + P 3 -3∞ and solve for α(P 1 -∞) + α(P 2 -∞) = Q -α(P 3 -∞). In that case, the variable y 3 gets involved in all the equations, so that we get a system in 6 variables.

For = 2, the 2-torsion elements are easily deduced from the factorization of f , and by computing a pre-image of a 2-torsion divisor, we get a point in J [START_REF] Birkenhake | Complex Abelian varieties[END_REF] from which we could deduce a, b, c mod 4. Dividing again by 2 was too costly, due to the fact that the 4-torsion point was in an extension of degree 4. For = 3, which is an inert prime, we ran the kernel computation for the multiplication-by-3 endomorphism, without using the RM property. The norm being 27, this is the largest modular computation that we performed (and the most costly in terms of time and memory). The prime = 7 ramifies in Z[η 7 ] as the cube of the ideal generated by α 7 = -2 -η 7 + η 2 7 . The kernel of α 7 can be computed but it yields only one linear relation in a, b, c mod 7. Dividing the kernel elements by α 7 would give more information, but again, this computation did not finish due to the field extension in which the divisors are defined. The first split prime is = 13. We use the following small generators: (13) = (2-η 7 -2η 2 7 )(-2+2η 7 +η 2 7 )(3+η 7 -η 2 7 ), which seem to produce the polynomial systems with the smallest degrees. For instance, the apparently smaller element 1 + η 2 7 of norm 13 yields equations of much higher degrees 7, 71, 72, 73, 72. The next split prime is 29, which would maybe have been feasible, but was not necessary for our setting. In the following table, we summarize the data for these systems, that were obtained with Magma V2.23-4 on a Xeon E7-4850v3 at 2.20GHz, with 1.5 TB RAM. The classical square-root-complexity search in genus 3 requires O(q) group operations [START_REF] Elkies | Elliptic and modular curves over finite fields and related computational issues[END_REF]. For RM curves, this can be improved by searching for the coefficients a, b, c of ψ = π + π ∨ in Z[η]. This readily yields a complexity in O(q 3/4 ), using the equation aD + bη(D) + cη 2 (D) = (q + 1)D, that must be satisfied for any rational divisor D. While a baby-step giant-step approach is immediate to design, it needs O(q 3/4 ) space and this is the bottleneck. A low-memory, parallel version of this search can be obtained with the algorithm of [START_REF] Gaudry | A low-memory parallel version of Matsuo, Chao and Tsujii's algorithm[END_REF], where the details are given only for a 2-dimensional problem, while here this is a 3-dimensional problem. But we did not hit any surprise when adapting the parameters to our case. Also, just like in [START_REF] Gaudry | A low-memory parallel version of Matsuo, Chao and Tsujii's algorithm[END_REF], including some anterior 1 The F4 algorithm can be highly sensitive to the modelling of the problem and we refer to the source code. In particular, thanks to serendipity, we saved a factor greater than 12 in the runtime for = 7, 13 by forgetting to take the squarefree part of the saturation polynomial. We have no explanation for this phenomenon.

modular knowledge is straightforward: if a, b, c are known modulo m, the expected time is in O(q 3/4 /m 3/2 ). We wrote a dedicated C implementation with a few lines of assembly to speed-up the additions and multiplications in F p , taking advantage of the special form of p. This implementation performs 10.7M operations in the Jacobian per second using 32 (hyperthreaded) threads of a 16-core bi-Xeon E5-2650 at 2 GHz. We used the knowledge of ψ modulo 156 but not of the known relation modulo 7 for simplicity (there is no obstruction to using it and saving an additional 7 1/2 factor).

After computing about 190,000 chains of average length 32,000,000, we got a collision, from which we deduced ψ = 2551309006 + 2431319810 η 7 -847267802 η 2 7 , and the coefficients of the characteristic polynomial χ π of the Frobenius are then σ1 = 986268198, σ2 = 35389772484832465583, σ3 = 10956052862104236818770212244.

The number of group operations that were done is slightly less than 43 (p 3/4 /156 3/2 ). This factor 43 is close to the average that we observed in our numerous experiments with smaller sizes. Scaled on a single (physical) core, we can estimate the cost of this collision search to be 105 core-days.

  i and on the ring Z[η].

Acknowledgements. We are grateful to Benjamin Smith for fruitful discussions and to Allan Steel for his help with memory issues with Magma. We also wish to thank anonymous referees for their comments which helped improve the paper.

Université de Lorraine, CNRS, Inria Email address: simon.abelard@loria.fr Email address: pierrick.gaudry@loria.fr Email address: pierre-jean.spaenlehauer@loria.fr


Keywords: palm-print recognition, image segmentation, feature extraction, ELM, image matching

In the contemporary period, significant attention has been focused on the prospects of innovative personal recognition methods based on palm print biometrics. However, diminished local consistency and interference from noise are only some of the obstacles that hinder the most common methods of palm-print imaging such as the grey texture and other low-level of the palm. Nevertheless, the development of the process and tackling of the obstacles faced have a potential solution in the form of high-level characteristic imaging for palm-print identification. In this study, Zernike Moments are used for acquiring superpixel features that are spiral scanned images, which is an innovative recognition method. By using the extreme learning machine, the inter-and intra-similarities of the palm-print feature maps are determined. Our experiments yield good results with an accuracy rate of 97.52 and an equal error rate of 1.47 % on the palm-print PolyU database.

Introduction

As a means of enhancing palm-print recognition precision, the development of strong algorithms has been the focus of biometric analysts over the last twenty years. Varying degrees of user tolerability, effectiveness, and precision have been observed with the use of palm-print, hand, face, voice, iris, and fingerprint biometric modalities. Owing to the recently reduced expenditure required for the sensors and high user tolerability of palm-print identification methods, the current biometric modalities are now being increasingly considered for supplementation with this original and promising biometric mechanism. Fingerprinting is often found to be more disconcerting for individuals than palm printing. Furthermore, individuals have been perceived to be concerned when an identification device requires that the eyeball of a user be scanned and iris calculations be performed. An additional advantage of palm printing is that a myriad of significant textures, ridges, creases, and lines that make for distinct characteristics exist on a palm, and it also presents a large surface for identifying such characteristics in contrast with various other biometric modalities.

The prospects for various biometric systems that make use of palm prints have been explored contemporarily; appearance-based and structural-based systems are the two primary types of biometric systems that have been identified. In the case of the structural-based systems, an array of edge identification mechanisms can be adopted to acquire a palm-print scan's distinctive lines and ridges [START_REF] Wua | Palmprint classification using principal lines[END_REF][START_REF] Han | Personal authentication using palm-print features[END_REF]. Owing to the typical characteristics of a palm of having ridges and wrinkles that intersect and the difficulty in producing an image of such features, the precision of the characteristic identification method is a highly significant aspect of performance according to existing studies. The efficiency of the use of palm-print lines for identifying a user, on having extracted its features by using an amended Radon transform, was confirmed in the contemporary research performed by (Huan et al. 2008). However, in order to resolve the issue of the deficient abstraction of the palm-print line features in the region of interest (ROI), an amended finite Radon transform with a larger training array for acquiring the orientation code was proposed by ( Jai et al. 2008).

In the case of appearance-based systems, a high-dimensional issue is regularly produced owing to the two-dimensional intensity variation features being treated and evaluated as a sole dataset in relation to the palm print. Therefore, the identification procedure should be improved in terms of precision through a crucial stage that involves transforming the high-dimensional palm-print information into a lowdimensional image. The primary statistical consistencies of the image can be determined and quantified as a means of conducting palm-print identification using such methods; both unsupervised and supervised classifications can be distinguished in this manner. The local binary pattern (LBP) [START_REF] Hanmandlu | Palm print recognition using Local Binary Pattern operator and support vector machines[END_REF]) and principal component analysis (PCA) [START_REF] Lu | Palmprint recognition using eigenpalms features[END_REF][START_REF] Zhao | Palmprint recognition with 2DPCA+PCA based on modular neural networks[END_REF][START_REF] Ribaric | A biometric identification system based on eigenpalm and eigenfinger features[END_REF] are two approaches that can be adopted to acquire the necessary data without the use of data label indexing for the unsupervised classifications. In the case of the supervised classifications, linear discriminant analysis and other common statistical means can be adopted for high-dimensionality reduction of the training phase's characteristic vectors while cataloguing on the basis of the predetermined class labels (Wue et al. 2003;[START_REF] Jing | A face and palmprint recognition approach based on discriminant DCT feature extraction[END_REF]. Maintaining the information's local structure through low-dimensional imaging has been contemporarily attempted using the LBP approach. However, the importance of the discriminant data, which is crucial to the identification phase, is reduced with the use of this method.

Individual pixel processing of representations within the spatial dimension has largely characterized appearance-orientated systems. However, if flaws exist, several disadvantages are posed by such methods. The enhancement of the discernment capabilities is possible by altering a basic characteristic selection data in the frequency domain through the use of common transform methods [START_REF] Laadjel | Palmprint recognition using fishergabor feature extraction[END_REF]Naani et al.2008). The defining characteristic data was acquired by using an enhanced form of LDA (Linear discriminant analysis) and choosing the DCT (discrete cosine transform) frequency bands using a two-dimensional separability judgment algorithm with the use of a 3040 palm-print scan database showing positive outcomes in Xiao et al.'s research (Jing et al.2004). In order to enhance the identification precision and reduce the expenditure on computation methods, an LDA was used for characteristic discernment across five representation levels of the palm print in Xiangqian et al.'s study using a Gaussian pyramid to produce the five levels and reduce the resolution (Wue et al. 2003). Furthermore, the high-dimensionality reduction of the characteristic space was also undertaken by implementing 2DPCA+PCA in Xin et al.'s study after acquiring the local characteristics by using Gabor filters, which decomposed the palm-print representations into various orientations and scales. However, the implementation of the Gabor filter bank causes the problem of exploding dimensionality, which is the major disadvantage of this approach [START_REF] Pan | Palmprint recognition using gabor feature-based (2D) 2PCA[END_REF]. Furthermore, by incorporating independent component analysis (FASTICA) with radial basis probabilistic neural networks (RBPNN), an innovative palm-print identification method was devised by Li-Shanga et al. However, the larger the database, the greater the decline observed in terms of magnitude and severity in identification precision, while the computational requirements are significant [START_REF] Shang | Palmprint recognition using fast ICA algorithm and radial basis probabilistic neural network[END_REF]. Individual identification on the basis of degraded palm-print images has not been an issue in any of the researches outlined above, where algorithms were analyzed using clean palm-print databases. Overall, the identification success is most likely to be undermined by the degradation of the palm-print images as a result of distortion or noise interference.

An array of obstacles is faced in LLF-based approaches despite various studies providing satisfactory experimental outcomes on this subject [START_REF] Wua | Palmprint classification using principal lines[END_REF]Jai W et al.2008). For example, the local consistency features of palm prints cannot be defined using low-level characteristics. Furthermore, representation discrepancies typically affect the corresponding algorithms, such as the Hamming distance algorithm, while low-level characteristic palm-print systems are vulnerable to pixel noise (Shang et al.2006).

Therefore, implementing Zernik moments in relation to the superpixel region is the novel approach presented in this research in order to overcome the aforementioned obstacles. The palm-print PolyU database showed a mere 1.47% equal error rate (EER), and thus, some excellent results were obtained in this research.

The ability of the aforementioned approach of acquiring all the palm data associated with the palm surface in order to produce a superpixel-based characteristic acquisition algorithm is the primary contribution of this study. Section 2 provides an in-depth explanation of the method adopted in this research. The experimental outcomes are outlined and evaluated in section 3, while section 4 presents the conclusions.

Methodology

Figure 1 introduces the SPC (superpixel scheme)-based original palm-print identification system. Having provided an overall summary of the devised system in figure 1, the distinct stages of the system can be identified. First, pre-processing of the representation occurs, which includes image normalization, ROI extraction, and image alignment. Second, in order to produce the superpixels, the simple linear iterative clustering (SLIC) technique [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF]) was adopted during the superpixel oversegmentation stage. Third, superpixel feature (SPF) extraction is performed using Zernike moments in order to acquire the superpixel characteristics (Kal et al.2 013). Fourth, feature normalization is performed [START_REF] Aksoy | Feature normalization and likelihood-based similarity measures for image retrieval[END_REF]. The final stage comprises the matching score calculation with template palm-print representations, and the input's standardized characteristics are used to calculate the matching score by using an extreme learning machine (ELM) [START_REF] Huang | Extreme learning machine: theory and applications[END_REF].

Extraction of the Region of Interest

As a means of straightening the palm-print representation, the key points were identified based on the technique explained for key point identification in [START_REF] Zhang | Online palmprint identification[END_REF]. The acquisition of a 128 × 128size ROI was possible after the palm-print representation had been aligned. Figure 2 presents the extracted ROI and a palm-print example.

Over-segmentation of Superpixels

Superpixels are typically acquired through superpixel over-segmentation. The production of superpixels via finger-vein representation over-segmentation was undertaken using the SLIC technique [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF] . The location data and local intensity are merged to produce small areas of the representation that were identified using SLIC. The advantages of using this technique are the reduced computational requirements and the requisite superpixel amount being produced based on a single particular parameter. There are two crucial points of departure from the k-means of the superpixel generation in the case of SLIC. Firstly, the superpixel size is equal to the constraints of the search area, thus considerably diminishing the total distance measurements of the optimization. Furthermore, the complexity is independent of the number of superpixels while also having linearity with the amount of pixels. Secondly, the regulation of superpixel density and size is possible, in addition to an amalgamation of the spatial and color closeness, by using the weighted distance measure. If the requisite amount of superpixels is represented by k, while the amount of pixels in the representation area is represented by n, then a superpixel-size equivalent area will be restricted for search in the SLIC. Thus, / s n k , and s s is the rough size of the superpixel area, while the grid interval is represented by s . Ultimately, as shown in Eq. ( 1), the weighted distance measure D is used to determine the amalgamation of the image plane and color space for the pixel grouping via SLIC. This produces superpixels that are dense and almost standardized.

2 2 c s c s d d D N N (1) 2 c i j d l l (2) 2 2 s i j i j d x x y y (3)

Region of interest extraction

With the normalized factors being C N and S N , the space proximity is represented by s d , the color proximity is represented by dc, the cluster centers are represented by k k k l x y , the pixel location is represented by x y , and colored regions of the greyscale representations are represented by T l .

Fig. 3 Over-segmentation of superpixels

Characteristics Extraction using Zernike Moments

The irregular partitions of the representation are indicated by the region adjacency graph (RAG), which has E edges and V vertices and is an undirected graph G. The feature vector F(v) is present for every participation while also possessing a vertex (v).

A precise characteristic descriptor was produced and an array of invariants devised by means of twodimensional Zernike moments (2DZMs). Inside unit disk D, the polar coordinates (ρ, θ) are adopted to characterize the 2DZMs of the pth order possessing repetition q. , 0 , , D p q p q p p q is even (4) and * 1 , ,

pq pq x y p Z V xyf xy (5)
The definition of the Zernike polynomials can be represented as follows:

* , iq pq pq V R e (6)
The complex conjugate is denoted by the asterisk.

The vector's length from origin to pixel ( , )

x y forms the radius 2 2

x y , while the angle formed by the x-axis and the vector is represented by arctan( / ) y x . The Zernike radial polynomial's path order is pq R , which can be determined as follows:

2 , ( 1) ! 2 ! ! ! 2 2 2 p k p k pq k q p q is even p k R p k k q k q (7)

Superpixel Oversegmentation

In order to describe the representation, it is essential that the maximum order of the 2DZMs is selected. When presented in the form of geometric movements, the second-and third-order Zernike invariants were found to be equal to the Hu invariants [START_REF] Kale | Zernike moment feature extraction for handwritten Devanagari compound character recognition[END_REF]).

The Region Adjacency Graph's spiral order traversal

In this section, the results obtained when the superpixel area is subjected to spiral scanning are presented As a means of forming a spiral path of successive vertices, the central partition was selected as the initial point for spiral scanning [START_REF] Trémeau | Regions adjacency graph applied to color image segmentation[END_REF] (see figure 4).

The idea of the algorithm is simple: the regions in the image are first numbered in spiral order from the center of the image. Each node (region) of the graph (i.e., RAG) that is represented by its barycenter corresponds to a number (in spiral order) in the image. This list is then sorted in alphanumeric order starting with the index 1. Thus, all the nodes of the graph are labeled in ascending order in the spiral direction.

Feature Vector Normalization

Feature vectors are irreconcilable owing to the distribution and range variability across them. However, feature vector normalization can be implemented using a normalization system [START_REF] Aksoy | Feature normalization and likelihood-based similarity measures for image retrieval[END_REF] such as the median, z-score, or min-max system in order to remedy the aforementioned issue. The feature vectors across the range 0,1 were normalized using the min-max system. On applying the min-max normalization, the normalized feature vector '

X can be denoted as shown below when the feature vector is

1 2 3 , , ,....., n X x x x x . ( ) ' ( ) ( ) i x Min X X Max X Min X (8)

Extreme Learning Machine

Single-hidden-layer feedforward neural networks have often been learnt using ELM [START_REF] Huang | Extreme learning machine: theory and applications[END_REF] . The iterative tuning of obscured nodes is not necessary following the random initialization through ELM. Thus, learning is required to occur just for the input weight parameters. If ( , ), 1,..., j j

x y j q indicates A as the training sample, with M j

x R and M j y R , then the formula below can be provided for the ELM's output function with L obscured neurons:

1 ( ) ( ) L l i i i f x g w x x G (9)
The 1 m output nodes are associated with L hidden nodes with the output weight vector 1 ,..., L G g g . The nonlinear activation function is represented by 1 ( ) ( ),..., ( )

L x w x w x [14].
The following formula is a variation of the explicit form of the system ( ) i x :

( ) ( . ), ,

d i i i i i x x R R (10)
The obscured layer parameters ( , ) have the activation function (.) . The following least square norm was used in order to resolve the output weight and minimize the training data error during the second phase of the ELM learning:

2 * min , N M H G R (11)
The hidden-neurons layer system is identified by :

1 1 1 1 1 1 ( . ) ( . ) ( . ) ( .
)

L L N LN L x x x x (12)
While the training data matrix H can be presented as follows:

1

T T N h H h (13)
The training set ( . ., ) i e L Q values are greater than the amount of hidden neurons L, which is assumed to be the case in reality, in order for the training error in Eq. ( 11) to be reduced in the optimal solution. [START_REF] Huang | Extreme learning machine: theory and applications[END_REF] presents the optimal solution for Eq. ( 11), in the form of the Moor-Penrose generalized inverse matrix, where the inverse of is * .

Experimental results

The PPDB (PolyU Palm-print database) is an open domain, vast-scale database that was used as the basis of the comprehensive experiments undertaken in order to confirm the system's validity (). Alongside the statistical legitimacy, the recognition rate and EER are the two terms used to present the outcomes.

PolyU palm-print database

For the PolyU palm-print database, two distinct sessions were held to obtain 20 palm-print images from 352 people. The test samples were those obtained during the second stage and comprised 10 images, while the control samples were those that were obtained during the initial session with 10 images per individual [START_REF] Zhang | An online system of multispectral palmprint verification[END_REF].

Critical Evaluation in Relation to Existing Methods

In this study, the SPF extraction is conducted by using Zernike moments from the last feature vector. To obtain the best results, superpixels of various sizes were selected. The degree of accuracy versus the superpixel number is shown in figure 5. Table 1 lists the PolyU PPDB database results obtained on applying the superpixel descriptor technique, and the effectiveness of this technique can thus be determined. With an EER of 2.76%, this experimental technique was found to be superior to other feature extraction mechanism (Kong et al.2004;Kumar an Shen.2004;Sun et al.2005;[START_REF] Kong | Palmprint identification using feature-level fusion[END_REF]. The known class labels of a restricted number of palm-print images are used as training samples for the palm-print identification, while the remainder were unknown and were used as probe samples. The aim of the process is to predict the probe-image class labels by comparing images with the training samples using a one-against-many comparison. The EER (equal error rate) in the identification process is used to estimate the performance of the palm-print identification.

The probe samples are classified such that they are associated with the class of training samples to which they are most similar. In this study, the training samples were considered each palm's first Ntrain ∈ {2, 4, 6} palm-print image(s). The probe palm prints were considered to be the remainder of the palmprint images. Ntrain × Npalm matching scores were generated in a database based on the probe sample.

The proposed method was then compared with several state-of-the-art orientation-based techniques, including the ordinal code (Sun et al.2005), fusion code [START_REF] Kong | Palmprint identification using feature-level fusion[END_REF], RLOC [START_REF] Guo | Palmprint verification using binary orientation cooccurrence vector[END_REF], palmCode (Kumar an Shen.2004), BOCV [START_REF] Guo | Palmprint verification using binary orientation cooccurrence vector[END_REF], and half-orientation methods [START_REF] Fei | Half-orientation extraction of palmprint features[END_REF]) (see figure 6). Several training palm-print images (between two and six for each palm-print database) were used to evaluate the proposed method's performance. The test dataset comprised the remainder of the palm-print images. Therefore, the identification experiments were performed six times with Ntrain = 1 to 6.

Figure 6 shows the results. It is evident that the proposed method is superior to other orientation-based methods when it is applied to the PolyU databases; it achieved lower error rates in both instances. The competitive code method extracted robust orientation features from the images with a high level of identification accuracy in specific instances. However, for the majority of the conditions, the proposed method was found to be superior to the competitive method.

In order to gauge the accuracy of the technique, we used various numbers of superpixels, and the results are shown in figure 7. The EER results reveal that the most efficient number of superpixels was 50 for the PolyU databases.

Speed Effectiveness

The proposed technique was implemented in MATLAB version 10.3 using a PC with a dual-core Intel i5-6600 3.3-GHz processor and 8-GB RAM. The operating system of the PC was Windows 10.0. The method was compared to previous state-of-the-art methods based on feature extraction in terms of computational cost in order to evaluate the complexity of its computation. The average time required to perform the feature extraction and matching over 100 runs was calculated to assess the computational complexity. The results thus obtained are shown in Table 2. The proposed method was found to consume a greater amount of time during feature extraction than the competitive code, BOCV, and half-orientation methods. However, the proposed method's matching speed was slightly higher as a simple matching scheme was used.

The proposed method required 26.228 s to perform the operation. This is acceptable for practical applications, although it is slightly greater than the time consumed by several orientation-based conventional methods. Future work will be focused on reducing the feature extraction time.

Conclusions and Prospective Research

SPFs were used for palm-print identification as a means of producing an original recognition system in this study. SPFs were acquired using certain statistical methods. In the case of the standard PolyU database, a 0.0147 EER was obtained with the devised system in the experimental outcomes. In the current research, the proposed techniques were found to be superior. The devised system comprised a number of innovative features. Firstly, the high-level character depiction of the palm-print imaging was introduced through the application of SPFs. Secondly, to the best of the researchers' knowledge, SPFs have not been applied to palm-print characteristic imaging prior to this study, while significant identification effectiveness was realized through the use of this original system for palm-print identification. Two prospective research paths may now be followed: one path involves the identification of additional SPFs for palm-print identification, while the other involves the investigation of image improvement methods for improving the palm print representation quality.

Fig

  Fig. 1 Palm-print identification system summary

Fig. 2

 2 Fig. 2 Extraction of the ROI

Fig. 4

 4 Fig. 4 (a). Palm-print representation, (b). Superpixel region, and (c). Spiral order of super-pixel region

Fig. 5

 5 Fig. 5 Average accuracy for various superpixel numbers

Fig. 6

 6 Fig. 6 EERs for PolyU palm-print databases

Fig. 7

 7 Fig. 7 EERs for various superpixel numbers

Table . 1

 . . PolyU palm-print databank and the proposed technique for measuring effectiveness

	EERs	Competitive code	Palm code	Fusion code	Ordinal code	RLOC	BOCV	Half -orientation	Proposed
	PolyU	0.0261	0.0931	0.0899	0.0272	0.0360	0.0469	0.0204	0.0147

Table . 2

 . . Computational costs of various orientation-based methods in feature extraction and matching.

	Methods	Feature extraction(s)	Matching(s)
	Competitive code	18.837	0.080
	Palmcode	1.760	0.204
	Fusion code	1.607	0.254
	Ordinal code	9.716	0.991
	RLOC	65.121	3.712
	BOCV	18.896	0.963
	Half-orientation	27.194	0.150
	Proposed	26.228	0.123


Olivier Bailleux
email: olivier.bailleux@u-bourgogne.fr

Attempts to produce minimal Resolution refutations

We address the challenge of searching minimal refutations proofs of inconsistent CNF formulae using the Resolution rule. We propose two algorithms which can only afford formulae of at most 5 variables with a desktop computer. A faster but incomplete algorithm is used to produce "hard" 5 variables 3CNF formulae though a stochastic greedy search. It allowed us to find formulae that can be refuted by producing clauses of at most 3 literals, but whose all minimal refutations contain at least one clause of 4 literals.

Introduction

Almost all modern complete SAT solvers indirectly use the Resolution rule [Rob65]. It is applicable to any pair of clauses with exactly one variable occurring with opposed signs, namely x ∨ α and ¬x ∨ β such that for any other literal y of α, ¬y is not in β. This very simple inference rule is complete for refutation, meaning that applied to an inconsistent set of clauses, it always produces the empty clause after a finite number of steps. But from the same set of clauses, namely, the input CNF formula, the number of steps before the empty clause is produced can drastically vary according to the clauses which are used as premises at each application of the rule.

In this report, the following question is addressed in a practical perspective : given an inconsistent CNF formula Σ, what is the smallest number of deductions needed to produce the empty clause ? Any refutation of Σ with such a number of deductions will be called a minimal refutation. Except in special cases, for example when all the clauses have a length of 2 [BOM07], producing a minimal refutation is a very hard problem [START_REF] Alekhnovich | Minimum propositional proof length is np-hard to linearly approximate[END_REF].

Producing minimal refutations

We tried two ways for producing minimal refutations. The first one consists in producing all the possible Resolution inferences in a single graph. The second one consists in producing all the possible refutations one by one.

Parallel search

This first approach consists in producing the graph of all possible derivations, one layer after another. The layer 0 contains the initial clauses. Each layer i > 0 contains the clauses occurrences derived from one premise of the preceding layer and another premise belonging to any layer j < i. This concept of clause occurrences is important because some clauses can be derived many times in different ways. In the following, any clause occurrence will be called a deduction.

The key point of the efficiency of such an approach is to avoid producing (or to remove) as many unnecessary deductions as possible. To this end, we compute two values related to each deduction q : the weight, defined as the number of inferences used to produce q, and the cover, defined as the number of initial clauses used to produce q. The weights are given in the gray disks in figure 1. Let min denote the weight of the refutation of the lowest weight produced so far. Clearly, any deduction q of size k and weight at least min -k can be rejected because any path from q to an empty clause would have length at least k. Figure 2 shows the clauses that can be removed thanks to this criterion. Another rejection criterion can be used if the input formula is critically inconsistent, i.e., each initial clause is used in any refutation. In such a case, any deduction q with weight w(q) and cover c(q) can be rejected if w(q) ≥ min -m + c(q), where m is the number of initial clauses. Subsumption and clause duplication can be used to reject some deductions in a restricted way, namely, if a deduction q is an ancestor of a deduction p and q subsumes p or is equals to p, then p can be removed. The reason is given in [START_REF] Jean | Proof Theory and Automated Deduction[END_REF] (chapter 7, corollary 7.16).

Property 2.1. [JI97] Let S be a set of clauses, and S be a subset of S. Assume that every clause C in S is in S or is subsumed by a clause in S . If the empty clause is derivable from S in k resolution steps, then it is derivable from S in at most k resolution steps.

Unfortunately, we found no other obvious way to use subsumption or clause duplication as a rejection criterion, because the deductions belonging to a minimal refutation are not necessarily themselves minimal. Consider the following example, where

Σ i (x, y) is the formula (s i ∨ t i ∨ x ∨ y) ∧ (¬s i ∨ t i ∨ x ∨ y) ∧ (s i ∨ ¬t i ∨ x ∨ y) ∧ (¬s i ∨ ¬t i ∨ x ∨ y)
, which produces the clause x ∨ y with 3 resolution steps, and σ(z) is a formula with variables not in Σ i , from which the clause (z) can be minimally produced with 14 deductions. Figure 3 shows a minimal refutation producing the deduction () 1 with 24 resolution steps, including the deduction (a) 1 of weight 15. The clause (a) can also be produced with only 14 inferences thanks to the deduction (a) 2 , but using this deduction instead of (a) 1 would not allow to produce a minimal refutation. The reason is that some deductions used to produce (a) 1 are also used to produce (¬a).

Incomplete parallel search

As we will see in the experimental part, beyond 4 variable, the time required to produce a minimal refutation becomes prohibitive. In order to produce shorts refutation, but without any guarantee of minimality, we propose to modify the previous algorithm in the following way: if any deduction q with weight w(q) is subsumed by or is equal to an existing deduction p with w(p) ≤ w(q), then q is rejected. The example of figure 3 shows that this new algorithm can fail to produce a minimal refutation. We failed to find such examples with critically inconsistent formulae, so we leave open the question of the minimality of the refutations produced in such a case.

Sequential search

As with the parallel search algorithm, any refutation graph is organized in layers, were at least one premises of each deduction of the layer i is in the the layer i -1.

The algorithm we used enumerates all the possible refutation graphs with a given number of resolution steps. Let the depth of an inference graph denote its number of layers excepts the one containing the initial clauses. For any possible depth of the derivation graph, we try to populate the first layer with any admissible sequence of deductions at this level. Theses sequences of deductions will be called clusters. For each cluster of the first layer, we try to populate the second layer, and so on. Some simple properties allow us to restrict the possible numbers of deductions in any layer i, given the sizes of the preceding ones and the depth of the derivation graph. Clearly, the layer depth must contain exactly one deduction corresponding to the empty clause. The the layer depth -1 cannot contain more that 2 deductions, and any layer depth -k cannot contain more than 2 k deductions. In addition, each layer must contain at least one deduction. It follows that if r denotes the number of deductions that remains to be done after the level i -1, and k denotes the number of layers after the layer i, then the size of the layer i is at least max(1, r -2 k + 1) and at most min(r -k, 2 k ). The search space can also be reduced by tacking into account the number s of sinks at each stage of the construction of the current derivation graph C, as well as the number r of deductions that remains to be done. This current graph C must be connected to the graph R consisting of the r deductions to be produced. Then R must have one sink, s input arcs, and r nodes, including the sink. On the example of figure 4, we can see that if s = 5, the graph can be completed only if r ≥ 4. More generally, each new node in R remove one input arc and add two ones, then we must have r ≥ s -1.

Of course, in addition to the preceding criteria, any new occurrence of an already produced (or initial) clause can be rejected, as well as any deduction of a clause that is subsumed by an already produced one in the current, partially built inference graph.

The algorithm used for enumerating all the inference graphs is summarized figure 5.

• cur represents the current level, i.e., the level of the layer in construction.

• layer[cur] represents the layer of depth cur.

• init layer [cur].n and init layer [cur].max set the current (minimal) and the maximal size of the current layer, according to the previously seen criteria.

• init layer[cur].cluster puts the first admissible cluster into the current cluster, and layer[cur].cluster++ produces the next admissible cluster for the current layer, according to the enumeration order of the clusters.

All these procedures can either succeed (blue arrows) or fail (red arrows).

Figure 5: Enumerating all the derivation graphs.

3 Experimental investigations

Comparing sequential and parallel search

Our first experiment compares the three algorithms applied to a minimalist pebbling formula [START_REF] Ben-Sasson | Near optimal separation of tree-like and general resolution[END_REF] with 4 variables, namely 13), ( 14), ( 34), ( 23), ( 24), ( 12 The third experiment is about the following small 3CNF formula with 5 variables. This formula is intractable using the parallel search algorithm with our desktop computer, and requires Billions of Resolution tests using the sequential algorithms. The incomplete parallel search algorithm quickly produces an optimal refutation.

Σ 1 = {(
Σ 3 = {(123), ( 12

Producing difficult formulae

As proved in [BSW01], the largest Refutation proofs also are the widest ones, in the sense that if the initial clauses have a fixed width (i.e. number of literals), the size of any minimal refutation is exponentially related to the width of the widest clauses produced in such a refutation, namely,

S(F ) = exp Ω w(F 0) -w(F ) n
were S(F ) is the size of any minimal refutation of a formula F , w(F ) is the width of the widest clause of of F , w(F 0) is the width of the widest clause of any minimal refutation of F , and n is the number of variables in F . Then, asymptotically, the kCNF formulae whose minimal refutations produce "wide" clauses are "hard" for Resolution. On the other hand, because there are only a polynomial amount of clauses of width k on n variables, "large" refutations of kCNF formulae do contain "wide" clauses.

But aside from these asymptotic facts, what are the smallest 3CNF formulae that cannot be refuted without producing at least a clause with more that 3 literals ? In [Urq87], such a formula with 6 variable is given. One of our aims was to try to found such a 3CNF formulae with 5 variables. To this aim, we implemented a randomized greedy algorithm with objective function provided by the incomplete parallel search of minimal refutations (IPS for short). First, a random 3CNF formula with 5 variable is produced. Then this formula is iteratively perturbed by replacing one randomly picked clause by another randomly generated one. The score of any consistent formula is the opposite of the number of deductions produced by IPS. The score of an inconsistent formula is the size of the refutation produced by IPS.

Let us denote w-refutation any refutation producing clauses of width at most w. A formula will be said w-refutable if and only if a w-refutation exists for this formula. We spent hundreds of billions Resolution attempts to try to produce 5 variables formulae with "large" minimal refutation proof, and checking them for a 3-refutation. We did not find any inconsistent but not 3-refutable formula. The most "hard" formula we found consists in 18 clauses and can be refuted with 17 deductions. The size of the minimal refutation proof for this formula was confirmed using the (complete) sequential search algorithm at nearly no cost because the formula is critically inconsistent and contains 18 clauses, then it cannot be refuted with less than 17 refutations according to the criterion illustrated figure 4.

Σ 4 ={( 345 ), ( 3
Interestingly, the formula Σ 4 is 3-refutable, but doe not admit any 3-refutation of size 17, as verified by the sequential search algorithm (at the cost of 15 billions resolution attempts). The smallest 3-refutation we found with the incomplete parallel search have size 21. This means that, sometimes, producing wider clauses can allow to produce shorter refutations, which is somewhat counter-intuitive.

Below is another example of formula with minimal 3-refutations longer than minimal 4-refutations (11 and 12 deductions, respectively).

Perspectives

The complexity of finding short resolution proofs has already been addressed [START_REF] Iwama | Complexity of finding short resolution proofs[END_REF] [FSW04]. Nevertheless, to the best of our known, if algorithms for finding minimal refutations have been implemented and experimentally evaluated, no such facts have been yet published. So, even if our results are somewhat disappointing, they are a first open contribution to this area. We hope this contribution will be followed by many more, and we wonder how far can be pushed the practical tractability of this challenging problem.

Figure 1 :

 1 Figure 1: Producing all the deductions, one level after another.

Figure 2 :

 2 Figure 2: Removing unnecessary deductions.

Figure 3 :

 3 Figure 3: Some duplicated clauses cannot be rejected.

Figure 4 :

 4 Figure 4: Dealing with the current number of sinks and the number of remaining nodes.

  3), (1 23 ), ( 12 3), (345), ( 345), (3 45), (34 5), (14), (1 4), ( 25), (2

  45), (2 35 ), ( 12 3), (145), (134), (1 23), (2 45 ), (1 25), ( 134 ), ( 145 ), ( 135 ), ( 234), ( 135), (2 34 ), (1 24), ( 12 5), ( 245)}

Σ 5 =

 5 Figure 6: A minimal 3-refutation (left, size 12) and a minimal 4-refutation (right, size 11) of Σ 5 . An example of 3-refutation and an example of 4-refutation of Σ 5 are given figure 6. The numbers of the initial clauses correspond to their rank in the formula. The numbers listed in each node are the ones of the initial clauses used for the corresponding deduction.