Bruno Grenet
email: bruno.grenet@lirmm.fr

Joris Van Der Hoeven

Grégoire Lecerf
email: lecerf@lix.polytechnique.fr

Deterministic root finding over finite fields using Graeffe transforms

We design new deterministic algorithms, based on Graeffe transforms, to compute all the roots of a polynomial which splits over a finite field F q . Our algorithms were designed to be particularly efficient in the case when the cardinality q -1 of the multiplicative group of F q is smooth. Such fields are often used in practice because they support fast discrete Fourier transforms. We also present a new nearly optimal algorithm for computing characteristic polynomials of multiplication endomorphisms in finite field extensions. This algorithm allows for the efficient computation of Graeffe transforms of arbitrary orders.

Motivation

In a previous paper [20], we presented efficient randomized root finding algorithms. Such algorithms were needed for the efficient interpolation, into the standard monomial basis, of polynomials that are given through evaluation functions. This latter problem is also known

Introduction

Let F q represent the finite field with q = p k elements, where p is a prime number, and k 1. Throughout this article, such a field is supposed to be described as a quotient of F p [x] by a monic irreducible polynomial. Let f ∈ F q [x] represent a separable monic polynomial of degree d 1 which splits over F q , which means that all its irreducible factors have degree one and multiplicity one. In this article we are interested in computing all the roots of f . under the name sparse interpolation, and root finding often turns out to be a bottleneck, as indicated in [START_REF] Van Der Hoeven | Sparse polynomial interpolation in practice[END_REF]. In fact, in this case, the ground field can be chosen to be F p with p = M 2 m + 1, and where 2 m is taken to be much larger than the number of terms to be discovered. In order to minimize the size of p, making it fit into a machine register, we usually take M = O(log p) as small as possible. A typical example is p = 7 • 2 26 + 1. We informally refer to such primes as FFT primes.

While working on [START_REF] Grenet | Randomized root finding over finite fields using tangent Graeffe transforms[END_REF], it turned out that some of the new methods could also be used for the design of fast deterministic methods for computing all the roots of f . Even though randomized algorithms are more useful in practice, the existence of efficient deterministic algorithms remains interesting from a theoretical point of view. The goal of this paper is to report on our progress on this issue.

Notations and prerequisites

We will use the same notations as in [START_REF] Grenet | Randomized root finding over finite fields using tangent Graeffe transforms[END_REF], which we briefly recall now. The multiplicative group F q \ {0} of F q is written F q * . Complexity bounds will frequently be expressed using the soft-Oh notation: 1) g(n). Given x ∈ R, we write ⌊x⌋ = max {k ∈ Z: k x} and ⌈x⌉ = min {k ∈ Z: k x}. The remainder of a division of g by f is denoted by g rem f . We write M(d) for the complexity of multiplication of polynomials of degree <d over an arbitrary ring. We also write I(n) and M q (d) for the bit complexities (in the Turing machine model) of n-bit integer multiplication and multiplication of polynomials of degrees <d over F q . We make the customary assumption that M(d) / d, I(n) / n and M q (d) / d are increasing functions. Currently best known bounds [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF][START_REF] Harvey | Even faster integer multiplication[END_REF][START_REF] Harvey | Faster polynomial multiplication over finite fields[END_REF] are M(d) = O(d log d log log d), I(n) = O(n log n 8 log * n ) and M q (d) = O d log q log(d log q) 8 log * (dlog q) , where log * stands for the iterated logarithm function log * n = min k ∈ N: log ... k× log x 1 . We freely use the following classical facts: ring operations in F p cost O(I(log p)) and one division or inversion in F p costs O(I(log p) log log p) [START_REF] Zur Gathen | Modern computer algebra[END_REF]Chapter 11]. The ring operations in F q cost at most O(M p (k)), and inversions take at most O(M p (k) log k + I(log p) log log p). For convenience, m q and d q will respectively denote cost functions for the product and the inverse in F q . The gcd of two polynomials of degrees at most d over F q can be computed in time O(M q (d) log d) [START_REF] Zur Gathen | Modern computer algebra[END_REF]Algorithm 11.4]. Given polynomials f and g 1 , ..., g l monic with deg f = d and deg g 1 + ••• + deg g l = O(d), all the remainders f rem g i can be computed in time O(M q (d) log l) [START_REF] Zur Gathen | Modern computer algebra[END_REF]Chapter 10]. The inverse problem, called Chinese remaindering, can be solved within a similar cost O(M q (d) log l + d d q ).

f (n) ∈ O ˜(g(n)) means that f (n) = g(n) log O(
For an integer n, the largest prime dividing n is written S 1 (n), and the second one S 2 (n). Pollard and Strassen have given a deterministic algorithm for factoring p -1, of bit complexity O M S 2 (p -1) 1/2 I(log p) log p (see for instance [START_REF] Zur Gathen | Modern computer algebra[END_REF]Corollary 19.4]). Ignoring smoothness considerations, the latter result has been improved into O I p 1/4 log p in [START_REF] Bostan | Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator[END_REF], and further refined to O I p 1/4 (log p)/ log log p √ in [START_REF] Costa | Faster deterministic integer factorization[END_REF]. In this article, when needed, we consider that the factorization of q -1 and a primitive element of F q * have been precomputed once, and we discard the necessary underlying costs. In practice, if the factorization of q -1 is given, then it is straightforward to verify whether a given element is primitive. For known complexity bounds and historical details on these tasks, we refer to [START_REF] Bach | Comments on search procedures for primitive roots[END_REF][START_REF] Mullen | Handbook of Finite Fields[END_REF][START_REF] Shoup | Searching for primitive roots in finite fields[END_REF].

Related work

It is now classical that polynomial factorization over F q reduces to finding roots over F p in deterministic polynomial time, uniformly in d and log q. But it is still unknown whether root finding can be solved in deterministic polynomial time or not, even under classical conjectures in number theory. This problem can nevertheless be solved efficiently by randomized algorithms in average polynomial time. Concerning related results and historical surveys on this topic, the reader might consult [START_REF] Zur Gathen | Modern computer algebra[END_REF][START_REF] Gathen | Factoring polynomials over finite fields: A survey[END_REF][START_REF] Kaltofen | Polynomial factorization: a success story[END_REF][START_REF] Mullen | Handbook of Finite Fields[END_REF].

Seminal algorithms for polynomial factorization over finite fields are classically attributed to Berlekamp [START_REF] Berlekamp | Factoring polynomials over finite fields[END_REF][START_REF] Berlekamp | Factoring polynomials over large finite fields[END_REF], and Cantor and Zassenhaus [START_REF] Cantor | A new algorithm for factoring polynomials over finite fields[END_REF], but central earlier ideas can be found in works of Gauss, Galois, Arwins, Faddeev and Skopin. Cantor-Zassenhaus' algorithm is randomized and well suited to compute roots of polynomials of degree d that split over F q in average time O(M q (d) (log q + log d) log d + d d q ). Of course, if q = O(d) then an exhaustive search can be naively performed in time O(M q (d) log d) (the factor log d can be discarded if a primitive element of F q * is given, by means of [8, Proposition 3]), so that the cost of root finding simplifies to O(M q (d) log q log d + d d q ). This classical approach is for instance implemented in the NTL library written by Shoup [START_REF] Shoup | NTL: A Library for doing Number Theory[END_REF]. However neither Berlekamp's nor Cantor-Zassenhaus' algorithm seems to benefit of particular prime numbers such as FFT primes. Instead, alternative approaches have been proposed by Moenck [START_REF] Moenck | On the efficiency of algorithms for polynomial factoring[END_REF], von zur Gathen [START_REF] Gathen | Factoring polynomials and primitive elements for special primes[END_REF], Mignotte and Schnorr [START_REF] Mignotte | Calcul déterministe des racines d'un polynôme dans un corps fini[END_REF], and then by Rónyai [START_REF] Rónyai | Factoring polynomials modulo special primes[END_REF]. Recent slight improvements of Moenck's and Mignotte-Schnorr's algorithms may also be found in [START_REF] Grenet | Randomized root finding over finite fields using tangent Graeffe transforms[END_REF].

Deterministic polynomial time root finding has been tackled by several authors. Camion showed that, for a fixed finite field, one can pre-compute a suitable subset which allows to factor any polynomial in deterministic polynomial time [START_REF] Camion | A deterministic algorithm for factorizing polynomials of F q [X][END_REF]. Nevertheless the construction of such a subset is purely theoretical, and it is not clear that there exists an efficient algorithm to compute it, even for FFT prime fields.

Schoof, Rónyai, Huang, and Źrałek designed different methods for particular types of input polynomials according to their syntax or to properties of the Galois group of the lifted input polynomial over Q [START_REF] Ming-Deh | Generalized Riemann hypothesis and factoring polynomials over finite fields[END_REF][START_REF] Rónyai | Factoring polynomials modulo special primes[END_REF][START_REF] Rónyai | Galois groups and factoring polynomials over finite fields[END_REF][START_REF] Schoof | Elliptic curves over finite fields and the computation of square roots mod p[END_REF][START_REF] Źrałek | Using partial smoothness of p -1 for factoring polynomials modulo p[END_REF]. Some of their algorithms require the general Riemann hypothesis (GRH) or the extended Riemann hypothesis (ERH). Sub-exponential algorithms are also known from the work of Evdokimov [START_REF] Evdokimov | Factorization of polynomials over finite fields in subexponential time under GRH[END_REF], which has recently led to new conjectures in terms of graphs [START_REF] Arora | Deterministic polynomial factoring and association schemes[END_REF]. Other conjectural algorithms can be found in [START_REF] Gao | On the deterministic complexity of factoring polynomials[END_REF][START_REF] Saha | Factoring polynomials over finite fields using balance test[END_REF].

Another series of articles concern complexity bounds which are uniformly polynomial in the degree and in S 1 (q -1). Such algorithms are practical when q -1 = p 1 m 1 ••• p r m r is sufficiently smooth. Assuming a primitive element of F q * is given, Mignotte and Schnorr [START_REF] Mignotte | Calcul déterministe des racines d'un polynôme dans un corps fini[END_REF] proposed a method based on a cascade of gcds, that needs O(M(d) i=1 r m i (log q + p i log d)) operations in F q . The computation of a primitive element might of course be expensive but it can be seen as a pre-computation. In fact, von zur Gathen proved the deterministic polynomial time equivalence (in terms of S 1 (p -1) and input size) between the computation of primitive elements and polynomial factorization [START_REF] Gathen | Factoring polynomials and primitive elements for special primes[END_REF].

Rónyai [START_REF] Rónyai | Factoring polynomials modulo special primes[END_REF] obtained a polynomial complexity bound in S 1 (p -1), d and log p by means of linear algebra techniques, using primitive roots of orders lower than in [START_REF] Gathen | Factoring polynomials and primitive elements for special primes[END_REF], but he did not explicit the exponents in the bound. For F p [x], Shoup [START_REF] Shoup | On the deterministic complexity of factoring polynomials over finite fields[END_REF] reached the bound p 1/2 log 2 p O ˜(d 2 ) in terms of the number of operations in F p , and then refined it to

p 1/2 log p O ˜(d k min (d, k)) + O ˜(d k 2 ) log p for F q [x]
, still in terms of operations in F p [START_REF] Shoup | A fast deterministic algorithm for factoring polynomials over finite fields of small characteristic[END_REF]. Finally Shoup proved the bound S 1 (p -1) 1/2 (d log p) O (1) in F p [x] under ERH [START_REF] Shoup | Smoothness and factoring polynomials over finite fields[END_REF].

Our contributions

The main contribution of this article is an efficient deterministic root finding method. The new algorithm is based on generalized Graeffe transforms and it is particularly efficient for finite fields F p such that p is an FFT prime. Roughly speaking, the generalized Graeffe transform replaces modular exponention as used in Cantor-Zassenhaus' algorithm, and gcd computations are changed to multi-point evaluations.

Thanks to a slight improvement of Shoup's adaptation of Pollard-Strassen's algorithm to find roots deterministically [START_REF] Shoup | Smoothness and factoring polynomials over finite fields[END_REF], we obtain a new deterministic complexity bound in terms of S 1 (p -1) 1/2 . In addition, in the smooth case S 1 (q -1) = O(log q), the cost of our deterministic algorithm is softly equivalent to the one of Cantor-Zassenhaus.

The complexity of generalized Graeffe transforms is studied in Section 2. A first ingredient is the transposed modular composition algorithm of Kedlaya and Umans [START_REF] Kedlaya | Fast polynomial factorization and modular composition[END_REF], also called modular power projections algorithm. The second ingredient is a special variant for finite fields of Newton-Girard's formulas for recovering a polynomial from the power sums of its roots. Such a variant had been previously designed by Bostan, Gonzalez-Vega, Perdry and Schost [START_REF] Bostan | Complexity issues on Newton sums of polynomials[END_REF] in the case of prime finite fields. We extend their results to finite fields from scratch, using very simple arguments, and independently of the framework designed in [START_REF] Caruso | Tracking p-adic precision[END_REF].

Generalized Graeffe transforms over finite fields

Classically, the Graeffe transform of a polynomial g ∈

F q [x] of degree d is the unique polynomial h ∈ F q [x] satisfying h(x 2 ) = g(x) g(-x). If g(x) = i=1 d (α i -x), then h(x) = i=1 d (α i 2 -x)
. This construction can be extended to higher orders as follows: the generalized Graeffe transform of g of order π, written G π (g), is defined as the resultant

G π (g)(x) = (-1) πd Res z (g(z), z π -x). If g = i=1 d (α i -x), then G π (g)(x) = i=1 d (α i π -x).
Equivalently, G π (g) is the characteristic polynomial of multiplication by x π in F q [x]/(g) (up to the sign). In this section we show how to compute Graeffe transforms efficiently. We start with a particularly simple and efficient algorithm for the smooth case, repeated from [START_REF] Grenet | Randomized root finding over finite fields using tangent Graeffe transforms[END_REF].

Smooth case

For our root finding algorithm, the most important case is when q -1 is smooth and the order π of the generalized Graeffe transform divides q -1.

Proposition 1. Let π 1 , ..., π m be integers 2, such that χ = π 1 ••• π m divides q -1, and let ξ i be given primitive roots of unity of order π i , for all i ∈ {1, ..., m}. If g is a monic polynomial in F q [x] of degree d, then the generalized Graeffe transforms of orders

π 1 , π 1 π 2 , π 1 π 2 π 3 , ..., χ of g can be computed in time O(m M q (µ d)) or σ O ˜(d log q log µ), where µ = max (π 1 , ..., π m ) and σ = π 1 + ••• + π m . Proof. Writing g(x) = c j=1 d (α j -x) in an algebraic closure of F q , the Graeffe transform of g of order π i is h i (x) = c π i j=1 d (α j π i -x)
. Consequently this leads to

h i (x π i ) = g(x) g(ξ i x) g(ξ i 2 x) ••• g ξ i π i -1
x . Using the latter formula, by Lemma 2 below, the transform can be obtained in time O(M q (π i d)). Taking the sum over i concludes the proof.

Lemma 2. Let g be a polynomial of degree d 1 in F q [x], let α ∈ F q , and let l be an integer. Then the product

P l (x) = g(x) g(α x) g(α 2 x) ••• g(α l-1 x) can be computed in time O(M q (l d)).
Proof. Let h := ⌊l / 2⌋. If l is even, then P l (x) = P h (x) P h (α h x), otherwise we have P l (x) = P h (x) P h (α h x) g(α l-1 x). These formulas lead to an algorithm with the claimed cost.

Notice that this lemma generalizes to arbitrary fields.

General case

In the general case, let g ∈ F q [x] be monic of degree d. Given a polynomial e of degree <d, we are interested in computing the characteristic polynomial h of the multiplication by e in R := F q [x]/(g). Our strategy essentially follows Le Verrier's method: we first compute the traces Tr(e i ) of the multiplication by e i in R for all i ∈ {1, ..., d}, which are also called the power sums of the roots of the characteristic polynomial h. The generating series

τ (x) = i 0 Tr(e i ) x i ∈ F q [[x]
] of these power sums relates to the logarithmic derivative of the reverse polynomial η = x d h(1/x) of h via a first order differential equation

η ′ (x) η(x) = τ (x) + O(x d ). (1)
This equation can be seen as a compact form of the classical Newton-Girard identities which relate power sums and symmetric polynomials. After computing τ it thus suffices to solve this equation in order to deduce h. If p d, then this method is classical. If p < d, then h cannot directly be retrieved in this way [START_REF] Pan | New techniques for the computation of linear recurrence coefficients[END_REF]Fact 2.1]. In the rest of this section, we will overcome this problem by computing with p-adic integers at a suitable finite precision κ.

Power sums

Let Z / p κ Z be the ring of truncated p-adic integers at precision κ. It is convenient to write GR(p κ , k) for the Galois ring of characteristic p κ and degree k, which is classically defined as GR(p κ , k) = ((Z/ p κ Z)[z]) / (µ), where µ ∈ (Z/ p κ Z)[z] is a monic polynomial of degree k and is irreducible modulo p.

Let g represent a given monic polynomial of degree d, let A = GR(p κ , k)[x] / (g). We write Tr: A → GR(p κ , k) the classical trace function on A. In order to compute the power sums Tr(e), Tr(e 2 ), ..., Tr(e d ), we will use the following proposition, which is essentially a consequence of a result due to Kedlaya and Umans in [START_REF] Kedlaya | Fast modular composition in any characteristic[END_REF][START_REF] Kedlaya | Fast polynomial factorization and modular composition[END_REF]:

Proposition 3. Let δ > 0 be a fixed constant. If d 1+δ q -1, then the traces Tr(e),
Tr(e 2 ), ..., Tr(e d ) for a given e ∈ A can be computed in time (d log(q κ )) 1+δ .

Proof. Let γ(x) = x d g(1 / x) denote the reverse polynomial of g. Since g is monic we have γ(0) = 1. In the canonical basis made of the classes of 1, x, ..., x d-1 in A, the trace function can be computed as the first d terms of the generating series

γ ′ (x) γ(x) = Tr(x) + Tr(x 2 ) x + Tr(x 3 ) x 2 + ••• + Tr(x d ) x d-1 + O(x d ),
in softly linear time, namely O ˜(κ d log q). The cardinality of A is q κ , and A contains at least q -1 invertible elements. Then, for every constant δ > 0 such that d 1+δ q -1, Theorem 7.7 of [27, p. 1792] provides us with an algorithm to compute Tr(e), Tr(e 2 ), ..., Tr(e d ) in time d 1+δ log 1+o(1) (q κ ).

Remark 4. For d close to q the proposition does not apply. A simple solution is to lift the computations in F q 2. Nevertheless, for finding the roots of g, this case is in fact favorable: it suffices to evaluate g at all the elements of F q . Therefore we will not need to work in F q 2.

Newton-Girard's identities

We carry on with the notations of the previous subsection. Let e be a given element in A, and let τ (x) = i 0 Tr(e i ) x i ∈ GR(p κ , k) [[x]] represent the generating series of the power sums of the roots of the characteristic polynomial h of e. If p d, then h can be deduced from Tr(e), Tr(e 2 ), ..., Tr(e d ) in softly linear time, using Newton-Girard's identities (1), as for instance in [5, Corollary 1]. If k = 1, then Theorem 1 of [START_REF] Bostan | Complexity issues on Newton sums of polynomials[END_REF] solves the equation modulo p whenever κ ⌊log d/ log p⌋ + 1. This subsection is dedicated to the remaining case when p < d and k is general. Our new alternative proofs are also shorter than in [START_REF] Bostan | Complexity issues on Newton sums of polynomials[END_REF].

Lemma 5. Let θ ∈ GR(p κ , k)[x] be of degree d satisfying θ(0) = 1, and θ ′ (x) / θ(x) = τ (x) + O(x d ). If κ > ⌊log d/log p⌋, then θ coincides with η modulo p at precision O(x d+1 ).
Proof. Let ψ = θ / η. Writing ψ as i 0 ψ i x i , we deduce from ψ ′ = O(x d ) that i ψ i = 0 for all i ∈ {1, ..., d}, whence p κ-val p i divides ψ i , where val p i denotes the p-adic valuation of i. Since κ > ⌊log d/ log p⌋ is equivalent to d < p κ , we have κ -val p i > 0, which concludes the proof. Proof. We will show how to compute a series θ ∈ GR(p κ , k) [[x]] satisfying θ(0) = 1 and θ ′ (x)/θ(x) = τ (x) + O(x d ), using Newton's method. In view of the above lemma, this will yield η followed by h modulo p.

Denote cost of multiplication in GR

(p κ , k)[x] of polynomials of degrees <d by M p κ ,k (d). We have M p κ ,k (d) = O k κ d log p log(k κ d log p) 8 log * (kκdlog p) = O ˜(κ d log q) according to [22].
Suppose that we are given θ such that θ(0) = 1 and θ ′ (x) / θ(x) = τ (x) + O(x n-1 ) for some integer n 1. Then we claim that there exists θ ˜∈ GR(p κ , k)[[x]] of valuation in x at least n and such that

θ ˜′(x) = τ (x) - θ ′ (x) θ(x) + O(x 2n-1 ). (2)
Indeed, using that τθ ′ /θ = O(x n-1 ) and θ ˜= O(x n ), this equation is equivalent to

θ ˜′ 1 + θ ˜= τ - θ ′ θ 1 -θ ˜+ θ ˜2 + ••• + O(x 2n-1 ) = τ - θ ′ θ + O(x 2n-1 ),
and therefore equivalent to

1 + θ ˜ θ ′ 1 + θ ˜ θ = θ ˜′ 1 + θ ˜+ θ ′ θ = τ (x) + O(x 2n-1 ).
This latter equation admits η / θ -1 as a solution. Having proved our claim, we can compute θ ˜as being any integral of τθ ′ / θ modulo O(x 2n ), and observe that θ ˆ= 1 + θ ˜ θ satisfies θ ˆ(0) = 1 and θ ˆ′/θ ˆ= τ + O(x 2n-1 ).

The cost of one iteration (2) is bounded by O(M p κ ,k (n)). Starting with n = 1 and θ = 1, we recursively apply the iteration for n = 1, 2, 4, ..., 2 ⌈log κ/log 2⌉-1 , the end result being

a θ ∈ GR(p κ , k)[[x]] with θ ′ /θ = τ + O(x d ).
It is classical that the cost of the last iteration dominates the total cost, which is therefore bounded by O(M p κ ,k (d)).

Characteristic polynomials

Now we come back to the original computational problem over F q .

Theorem 7. Let g ∈ F q [x] be a monic polynomial of degree d, and let e ∈ F q [x] be of degree <d. For every constant δ > 0 such that d 1+δ q -1, the characteristic polynomial h of the multiplication by e in F q [x]/(g) can be computed in time (d log q) 1+δ .

Proof. We embed

F q [x]/(g) into A = GR(p κ , k)[x]/(g), where GR(p κ , k) is constructed from (Z/ p κ Z)[z]/ (µ), and µ ∈ (Z/ p κ Z)[z] is the monic defining polynomial of degree k of F q over F p .
Let δ > 0 be such that d 1+δ q -1. Proposition 3 allows us to compute the traces of the successive powers of e seen in A in time (d log(q κ )) 1+δ . Then Proposition 6 allows to deduce h modulo p in softly linear time.

Remark 8. If the characteristic p is larger than the multiplicities of the roots of h, then Pan's algorithm [START_REF] Pan | New techniques for the computation of linear recurrence coefficients[END_REF] does not involve p-adic arithmetic, but may require higher order traces.

Corollary 9. Let g be a monic polynomial in F q [x] of degree d. For every constant δ > 0 such that d 1+δ q -1, the generalized Graeffe transform G π (g) can be computed in time (d log q) 1+δ + O ˜(d log q log π).

Proof. We apply the latter theorem to e(x) = x π rem g(x), that can be computed in time O(M q (d) log π).

Tangent Graeffe transforms

Whenever a Graeffe transform preserves the number of distinct simple roots, the so called tangent Graeffe transform can be used to directly recover the original simple roots from the simple roots of the transformed polynomial. Our randomized root finding algorithms from [START_REF] Grenet | Randomized root finding over finite fields using tangent Graeffe transforms[END_REF] heavily rely on this fact. Unfortunately, we have not found a way to exploit tangent Graeffe transforms in the deterministic setting. For completeness, we will nevertheless show that the results of this section extend to the tangent case.

Introducing a formal parameter ε with ε 2 = 0, we define the generalized tangent Graeffe transform of g of order π as being G π (g

(x + ε)) ∈ (F q [ε] / (ε 2 ))[x].
For any ring R, computations with "tangent numbers" in R[ε] / (ε 2 ) can be done with constant overhead with respect to computations in R (in the FFT model, the overhead is asymptotically limited to a factor of two).

In the context of the Graeffe method, the tangent transform is classical (for instance, see [START_REF] Malajovich | Tangent Graeffe iteration[END_REF][START_REF] Pan | Solving a polynomial equation: Some history and recent progress[END_REF] for history, references, and use in numerical algorithms). The generalized tangent Graeffe transform can also be seen as the tangent characteristic polynomial of x π modulo g(x + ε), and this construction is attributed to Kronecker in algebraic geometry [START_REF] Giusti | A Gröbner free alternative for polynomial system solving[END_REF][START_REF] Kronecker | Grundzüge einer arithmetischen Theorie der algebraischen Grössen[END_REF].

Proposition 1 from Section 2.1 admits a straightforward generalization to tangent Graeffe transforms with a constant overhead. In order to generalize Proposition 3, we introduce Proof. The value Tr ε ((e + ε e ¯)i ) equals the trace of (e(xε)

A ε = GR(p κ , k)[[ε]][x] / (g(x + ε)), where GR(p κ , k) = ((Z / p κ Z)[z]) / (µ), where µ ∈ (Z / p κ Z)[z] is
+ ε e ¯(x -ε)) i computed in GR(p κ , k)[[ε]][x]/(g).
Therefore we obtain:

Tr ε ((e(x) + ε e ¯(x)) i ) = Tr((e(x -ε) + ε e ¯(x -ε)) i ) = Tr(e i (x) + i ε (e ¯(x) -e ′ (x)) e i-1 (x)).
We therefore use Theorem 7.7 of [27, p. 1792] twice over A = GR(p κ , k)[x] / (g): once with Tr, and once for Tr • ϕ, where ϕ is the product by e ¯e ′ .

Theorem 11. Let g ∈ F q [x] be a monic polynomial of degree d, and let e, e ¯∈ F q [x] be of degrees <d. For every constant δ > 0 such that d 1+δ q -1, the characteristic polynomial h(x) + ε h ¯(x) of the multiplication by e(x)

+ ε e ¯(x) in F q [[ε]][x]/(g(x + ε)) can be computed in time (d log q) 1+δ . Proof. We embed F q [[ε]][x]/(g) into A ε = GR(p κ , k)[[ε]][x]/(g(x + ε)), where GR(p κ , k) is constructed from (Z / p κ Z)[z] / (µ), and µ ∈ (Z / p κ Z)[z] is the monic defining polynomial of degree k of F q over F p .
Let δ > 0 be such that d 1+δ q -1. Proposition 10 allows us to compute the traces of the successive powers of e + ε e ¯seen in A ε in time (d log(q κ )) 1+δ . Then Proposition 6 allows us to deduce h + ε h ¯modulo p in softly linear time, mutatis mutandis.

Corollary 12. Let g be a monic polynomial in F q [x] of degree d. For every constant δ > 0 such that d 1+δ q -1, the generalized tangent Graeffe transform G π (g(x + ε)) can be computed in time (d log q) 1+δ + O ˜(d log q log π).

Proof. We apply the latter theorem to e(x) + ε e ¯(x) = x π rem g(x + ε), that can be computed in time O(M q (d) log π).

Root finding based on Graeffe transforms

Our new root finding algorithm is technically reminiscent of numerical solvers: we make use of the Graeffe transform instead of the modular exponentiation, and replace gcds by multi-point evaluations. We begin this section with a first simple version of the algorithm, and then introduce new ingredients in order to improve it. Throughout this section, the quantity δ > 0 represents a constant that can be fixed arbitrarily small.

Simple version of the deterministic algorithm

We assume given integers

π i 2 such that q -1 = ρ π 1 ••• π m . We let χ = π 1 ••• π m , and σ = π 1 + ••• + π m as previously.
The core algorithm of this section computes a sequence of Graeffe transforms, and then deduces inductively their roots starting from the last transform and finishing to the input polynomial.

Algorithm 13.

Input. f ∈ F q [x] of degree d, monic, separable, which splits over F q , and such that f (0) = / 0; a primitive root ξ of unity of order χ. Output. e 1 , ..., e s in {0, ..., χ -1}, such that f divides i=1 s (x ρξ e i ).

1. Compute the Graeffe transform h 0 of order ρ of f , and recursively compute h i as the Graeffe transform of order π i of h i-1 , for all 1 i m -1. b. If E has cardinality more than d then remove the elements e from E such that h i-1 (ξ e ) = / 0, using a multi-point evaluation to compute the h i-1 (ξ e ). 4. Return E.

Initialize

Proposition 14. Algorithm 13 is correct. If d 1+δ q -1 then it costs σ O(M q (d) log d) + (d log 2 q) 1+δ + O(m d I(log χ) + m m q d log χ) = σ O ˜(d log q) + (d log 2 q) 1+δ + O ˜(d log 3 q).
Whenever a primitive root of unity of order ρ is given, and µ := max (ρ, π 1 , ..., π m ) = O(log q), this bound further reduces to O ˜(d log 3 q). Proof. We prove the correctness by descending induction on i from m down to 1. At the end of iteration i of the loop of step 3, we claim that f divides e∈E (x ρπ 1 •••π i-1ξ e ), and that the integers e in E are in {0, ..., χ -1} and are divisible by ρ π 1 ••• π i-1 . These properties all hold for i = m + 1 when entering step 3. By induction, we may assume that these properties hold with i + 1 when entering the loop at level i. Let e ∈ E, and let α be a root of f such that α ρπ 1 •••π i = ξ e . Since e is divisible by π i , there exists j ∈ {0, ..., π i -1} such that α ρπ 1 •••π i-1 = ξ (e+ jχ)/π i . This proves the correctness.

From now on we assume that d 1+δ q -1. As to the complexity, we can compute the χ / π i for i ∈ {1, ..., m} in time O(m I(log χ)) as follows: we first compute

π 1 , π 1 π 2 , π 1 π 2 π 3 , ..., π 1 π 2 ••• π m-1 , and π m , π m-1 π m , π m-2 π m-1 π m , ..., π 2 ••• π m-1 π m , and then deduce each χ/π i by multiplying π 1 ••• π i-1 and π i+1 ••• π m .
Step 1 requires time (d log 2 q) 1+δ by Corollary 9. In step 3.a, we compute all the e/π i in time O(d I(log χ)), and then all the ξ e/π i and ξ χ/π i by means of O(d log χ) operations in F q . Deducing all the ξ (e+jχ)/π i takes O(π i d) additional products in F q . The multi-point evaluations can be done in time O(π i M q (d) log d).

Finally, when µ = O(log q), the term (d log 2 q) 1+δ can be replaced by

(ρ + σ) O ˜(d log q) = O ˜(σ d log q + d log 2 q) = O ˜(d log 3 q) in view of Proposition 1.
The cost in the latter proposition is asymptotically higher than the one of Mignotte-Schnorr's algorithm [20, Section 2], because of the term (d log 2 q) 1+δ . In the smooth case when µ = O(log q), the asymptotic complexity bounds become similar.

The rest of this section is devoted to improve Algorithm 13. In order to diminish the dependence on log q, we stop representing roots by their logarithms as in Mignotte-Schnorr's algorithm. Instead we shall compute discrete logarithms only when necessary. Furthermore, we will show how to decrease the dependence on the π i to π i √ , by using the "baby-step giant-step" technique.

Reverse splitting of Graeffe transforms

Let g be a monic polynomial which splits over F q . Once the roots a 1 , ..., a l of the Graeffe transform of order π of g have been recovered, we may decompose g into l factors g 1 , ..., g l such that the π-th powers of the roots of g i all coincide to a i , for all i ∈ {1, ..., l}. In short, we have Gr π (g i )(x) ∽ (xa i ) deg g i . This decomposition may be achieved efficiently via the following divide and conquer algorithm.

Algorithm 15.

Input. g ∈ F q [x] of degree d, monic, separable, which splits over F q , and such that g(0) = / 0; the pairwise distinct roots a 1 , ..., a l of the Graeffe transform of a given order π 2 of g. Output. The sequence of the monic gcd(x πa i , g(x)) for all i ∈ {1, ..., l}.

1. If l = 1 then return g. 2. Let h := ⌊l /2⌋.

Compute e(x) :=

i=1 h (x -a i ).
4. Compute g 1 (x) := gcd(e(x π ), g(x)) and g 2 (x) := g(x)/ g 1 (x). 5. Call recursively the algorithm with input g 1 , a 1 , ..., a h and g 2 , a h+1 , ..., a l , and return the concatenation of the polynomial sequences returned so.

Proposition 16. Algorithm 15 is correct and costs π O ˜(d log q log π) or (d log q) 1+δ + O ˜(d log π log q), for any constant δ > 0.

Proof. If l = 1, then the π-th powers of the roots of g are all equal to a 1 , whence the correctness. If l > 1, then the roots of g 1 are exactly those of g whose π-th powers are in the set {a 1 , ..., a h }. Consequently, the π-th powers of the roots of g 2 are all in the set {a h+1 , ..., a l }. This completes the proof of the correctness by induction.

Step 3 can be performed in softly linear time by the subproduct tree technique (for the sake of efficiency the subproduct tree should be shared between all the recursive calls).

On the one hand, computing g 1 can be done naively in time O(M(π d) + M(d) log d) = π O ˜(d log q log π). On the other hand, one may proceed as follows: compute r(x) := x π rem g(x) in time O ˜(d log π log q), and then use [27, Corollary 7.2, p. 1789] to obtain e(r(x)) rem g(x) in time (d log q) 1+δ .

The sum of the degrees of the input polynomials g at the same depth of the recursive calls is at most d, and the maximal depth is O(log l) = O(log d). This yields the complexity bounds of the algorithm.

Finding roots of several polynomials in a given set

Let g 1 , ..., g l be polynomials of respective degrees d 1 , ..., d l , and let

d := d 1 + ••• + d l .
For given subset A of F q of cardinality d, we are interested in finding the roots of g i in A for each i ∈ {1, ..., l} simultaneously. We do not make additional assumptions, and in particular the g i are allowed to share common roots. A natural divide and conquer approach applies, as explained in the following algorithm.

Algorithm 17.

Input. g 1 , ..., g l ∈ F q [x] of degrees 1; a subset A of F q of cardinality at most d := Σ i=1 l deg g i . Output. The sequence of the sets of the roots of g 1 , ..., g l in A.

1. If l = 1, then evaluate g 1 at all the points of A, and return those which are roots. 2. Let h := ⌊l /2⌋. 3. Compute e 1 := i=1 h g i and e 2 := i=h+1 l g i . 4. Compute the subset A 1 of points of A which are roots of e 1 . 5. Compute the subset A 2 of points of A which are roots of e 2 . 6. Call recursively the algorithm with input g 1 , ..., g h , A 1 and g h+1 , ..., g l , A 2 , and return the concatenation of the returned sequences of sets of roots.

Proposition 18. Algorithm 17 is correct and costs

O(M q (d) log d log l) = O ˜(d log q).
Proof. For the correctness, it is clear that the roots of g 1 , ..., g h in A are to be found in A 1 -idem for A 2 . The recursive calls are legitimate since the cardinality of A 1 (resp. of A 2 ) cannot exceed deg g 1 (resp. deg g 2 ).

Steps 1 to 5 take O(M q (d) log d). The sum of the degrees of the input polynomials at the same depth of the recursive calls is at most d, and the maximal depth is O(log l). This yields the complexity bound, thanks to the super-additivity M q .

Notice that this lemma generalizes to arbitrary fields.

Simultaneous baby-step giant-step root finding

Instead of using multi-point evaluations in step 3.b of Algorithm 13, we may split the separable part of h i-1 thanks to Algorithm 15. We will detail later how this situation further reduces to computing all the roots of several polynomials of degree sum at most d, and having roots of order dividing π. To perform efficient root finding in this case, we adapt Pollard-Strassen's aforementioned strategy for integer factoring as follows.

Algorithm 19.

Input. g 1 , ..., g l ∈ F q [x] of degrees 1, monic, separable, which split over F q , and such that g j (0) = / 0 holds for all j ∈ {1, ..., l}; a multiple π of the orders of the roots of all the g j ; a primitive root ξ of unity of order π. Output. The sets L j of the ξ-logarithms of the roots of g j for all j ∈ {1, ..., l}.

1. Let d := j =1 l deg g j , let s := π d √
, and t := ⌈π /s⌉.

Compute h(x) := j =0

s-1 (ξ j x -1).

3. For all j ∈ {1, ..., l} and all i ∈ {0, ..., t -1}, compute r j,i (x) := h(x) rem g j (ξ si x).

4. For all j ∈ {1, ..., l}, compute the set I j of indices i ∈ {0, ..., t -1} such that r j ,i (x) has a proper gcd e j,i (x) with g j (ξ si x). 5. For all j ∈ {1, ..., l} and all i ∈ I j , compute all the ξ-logarithms E j ,i of the roots of e j ,i in A = {1, ξ, ξ 2 , ..., ξ s-1 } by calling Algorithm 17 ⌈s/d⌉ times on subsets of A of cardinalities at most d. 6. Return the sets L j := i∈I j (s i + E j,i ), for all j ∈ {1, ..., l}.

Proposition 20. Algorithm 19 is correct and costs

O ˜ π d √ log q + O ˜(d log q).
Proof. Of course, the returned values in L j are ξ-logarithms of roots of g j . Conversely, let ξ u be a root of g j with 0 u π -1. There exists a unique i in {0, ..., t -1} such that i s u < (i + 1) s, since s t π. Then ξ u-is is a common root of g j (ξ si x) and h(x). This proves the correctness. By Lemma 2, step 2 executes in time

O(M q (s)) = O ˜ π d √ log q . Since d t = s + O(d), step 3 costs O(M q (s) log(s) + M q (d) log d) = O ˜ π d √ log q + O ˜(d log q).
Using the super-additivity of M q , step 4 needs

O(t M q (d) log d) = O ˜ π d √ log q + O ˜(d log q).
Since the sum of the degrees of all the polynomials e j ,i is at most d, Proposition 18 implies a cost O ˜(s log q + d log q) for step 5.

If l = 1, then Algorithm 19 slightly improves Shoup's variant of Pollard-Strassen's algorithm described in [START_REF] Shoup | Smoothness and factoring polynomials over finite fields[END_REF], which takes s = π √ independently of d; taking s = π d √ turns out to be a better trade-off between the baby steps and the giant steps. We also notice that this trade-off might be further improved by introducing logarithm factors in s.

Fast discrete logarithm in smooth finite fields

Let us recall that the term O ˜(d log 3 q) in the complexity of Algoritm 13 is due to computing roots from their logarithms. In the next subsection we will proceed differently and will compute logarithms only when needed. This is why we now address the discrete logarithm problem. Let us recall that the seminal discrete logarithm method in the smooth case of F p * goes back to [START_REF] Pohlig | An improved algorithm for computing logarithms over GF(p) and its cryptographic significance (corresp.)[END_REF] and performs S 1 (p -1) O(log p) operations in F p , making use internally of Shanks' baby-step giant-step paradigm. We recall this algorithm for completeness, adapting it to our context of Turing machines.

Algorithm 21.

Input. A primitive root ξ ∈ F q * of order π 2 of unity; a in the multiplicative subgroup generated by ξ. Output. The ξ-logarithm of a in {0, ..., π -1}. 1. Let s := ⌊ π √ ⌋ and t := ⌈π /s⌉.

2. Compute 1, ξ, ξ 2 , ..., ξ s-1 .

3. Compute a, a ξ -s , a ξ -2s , ..., a ξ -(t-1)s . 4. Sort the array made of the pairs (i, ξ i ) for i ∈ {0, ..., s -1} and (j , a ξ -js ) for j ∈ {0, ..., t -1}, accordingly to the lexicographic order of the second entries of the pairs seen as k-tuples in {0, ..., p -1} k . 5. Find the first two consecutive pairs (e 1 , ξ e 1 ) and (e 2 , a ξ -e 2 s ) whose second entries coincide. [START_REF] Bostan | Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator[END_REF]. Return e 1 + s e 2 rem π. Proposition 22. Algorithm 21 is correct and requires a running time π √ O ˜(log q) + π √ log π O(log q). Proof. Since s t π, the discrete logarithm e = log ξ a can be uniquely written as e = e 2 s + e 1 with 0 e 1 < s and 0 e 2 < t. This proves the correctness.

Steps 2 and 3 perform O( π √ ) products in F q and one inversion. Using the merge-sort algorithm, step 4 costs O( π √ log q log π).

The following divide and conquer approach allows us to compute logarithms for composite orders χ = π 1 ••• π m . Independently from the previous algorithm, it leads to a good complexity bound in terms of σ = π 1 + ••• + π m . This approach might already be known, but we could not find a reference to it in the literature.

Algorithm 23.

Input. ξ ∈ F q * of order χ = π 1 ••• π m , with given π i 2; a in the multiplicative subgroup generated by ξ. Output. The ξ-logarithm of a in {0, ..., χ -1}.

1. If m = 1, then find and return e ∈ {0, ..., π 1 -1} such that a = ξ e .

Let m

1 := ⌊m/2⌋, m 2 := m -m 1 , χ 1 := π 1 ••• π m 1 and χ 2 := χ/ χ 1 . 3.
Recursively compute e 1 := log ξ χ 1 a χ 1 . 4. Recursively compute e 2 := log ξ χ 2(a / ξ e 1 ). 5. Return e 1 + χ 2 e 2 .

Proposition 24. Algorithm 23 is correct and executes in time O(σ m q + d q log χ + m q logχ log log χ + I(log χ) log log χ) = σ O ˜(log q) + O ˜(log q log χ). By using Algorithm 21 in step 1, the cost becomes µ √ O ˜(log 2 q), where µ := max (π 1 , ..., π m ).

Proof. We prove the correctness by induction. The case m = 1 is clear. Assume m 2. In step 3, we have e 1 ∈ {0, ..., χ 2 -1}. By definition of e 1 , we have (a/ ξ e 1 ) χ 1 = a χ 1 / ξ χ 1 e 1 = 1 in step 4, whence a / ξ e 1 indeed lies in the multiplicative subgroup of order χ 1 generated by ξ χ 2 . By definition of e 2 , we have a/ ξ e 1 = (ξ χ 2 ) e 2 , whence a = ξ e 1 +χ 2 e 2 , which completes the induction.

As to the complexity bound, let T (χ) denote the cost of the algorithm. We share the subproduct tree built from π 1 , ..., π m , in time O(I(log χ) log m), between the recursive calls. Using binary powering in steps 3 and 4, we obtain the recursion relation

T (χ)
O(m q log χ + d q ) + T (χ 1 ) + T (χ 2 ).

which leads to the bound T (χ) Σ i=1 m T (π i ) + O(m d q + m q logχ log m). Using a naive exhaustive search in step 1 yields T (π i ) = π i m q , whence T (χ) = O(σ m q + m d q + m q logχ log m + I(log χ) log m). The first bound thus follows using m = O(log χ). The second bound follows from Proposition 22 that gives T (π i ) = π i √ log π i O ˜(log q). Remark 25. As a byproduct of the complexity analysis, we remark that, given a primitive root of unity ξ of order χ = π 1 ••• π m , where the π i are given integers 2, one may compute primitive roots of unity ξ i for all orders π i in time O ˜(log q log χ). In the next subsection we will need the following bound, taking the number of simple and multiple roots into account.

Squarefree part

Proposition 26. Let g ∈ F q [x] be of degree d. In an algebraic closure of F q , let δ 1 be the number of simple roots of g, let δ 2 be the number of its multiple roots, and let ∆ 2 = dδ 1 be the number of multiple roots counted with multiplicities. Then the squarefree factorization and the squarefree part of f may be computed in time O ˜(d log q) + ∆ 2 O ˜(log 2 q).

Proof. Let (g 1 , q 1 , ν 1 ), ..., (g l , q l , ν l ) in (F q [x] \ F q ) × {1, p, p 2 , ..., } × N represent the separable factorization of g [START_REF] Lecerf | Fast separable factorization and applications[END_REF]Introduction], so that we have g

(x) = g 1 (x q 1 ) ν 1 ••• g l (x q l ) ν l .
Since F q is perfect, the factors g i (x q i ) with q i = ν i = 1 contain all the simple roots of g. Letting h 1 = q i =ν i =1 g i and h 2 = g / h 1 , we have deg h 1 = δ 1 and deg h 2 = ∆ 2 . Following [29, Section 4], the actual number of p-th root extractions is O(∆ 2 ). Each such extraction amounts to O(log q) products in F q . In this way we obtain the squarefree factorization of g.

Improved root finding algorithm

We are now ready to combine all the previous algorithms of this section in order to improve Algorithm 13. Recall that we are given integers π i 2 such that q -1 = ρ π 1 ••• π m , and that we let χ = π 1 ••• π m . The following algorithm is parametrized by subroutines used internally. We distinguish three cases of use, leading to complexity bounds in terms of σ, of µ, and in the smooth case of S 1 (q -1) = O(1).

Algorithm 27.

Input. f ∈ F q [x], monic, separable, which splits over F q , of degree d < q, such that f (0) = / 0; a primitive root ζ of F q * . Output. The roots of Gr ρ (f ).

1. Compute h 0 as the Graeffe transform of f of order ρ. Compute h ˜0 as the squarefree part of h 0 . Recursively compute h i as the Graeffe transform of order π i of h ˜i-1 , and h ˜i as the squarefree part of h i , for all i from 1 to m -1. 2. Initialize Z with {1}, and compute ξ = ζ ρ . 3. For i from m down to 1 do a. Use Algorithm 15 with input h ˜i-1 and Z = {a 1 , ..., a l }, and let g 1 , ..., g l be the polynomials obtained in return.

b. Initialize W 0 with the roots of all the g j of degree 1, for j ∈ {1, ..., l}. c. For all j ∈ {1, ..., l} such that deg g j 2, compute the ξ-logarithm e j of a j with Algorithm 23 (resp. with Algorithms 23 and 21). d. Use Algorithm 17 (resp. Algorithm 19) to compute the sets of roots W j of g j ξ e j /π i x for all j ∈ {1, ..., l} such that deg g j 2. e. Z := W 0 ∪ j ξ e j /π i W j . 4. Return Z.

Proposition 28. Algorithm 27 is correct and costs (σ + ρ) O ˜(d log q log (µ ρ)) + O ˜(d log 2 q). In addition, for any fixed δ > 0, the cost is also bounded by σ O ˜(d log q) + (d log 2 q) 1+δ or O ˜( µ √ d log 2 q) + (d log 2 q) 1+δ , whenever d 1+δ q -1 and where µ = max (π 1 , ..., π m ).

Proof. Let introduce h m := Gr π m h ˜m-1 , which does not need to be computed. When entering the loop in step 3, the set Z contains the single root of h m . We prove by descending induction on i from m down to 1 that Z contains the roots of h i when entering step i of the loop. Assuming this hypothesis holds for i m, conditions of Algorithm 15 are met, and finding the roots of h ˜i-1 reduces to finding the roots of g 1 , ..., g l . The g j of degree 1 contribute to roots of h ˜i-1 in a straightforward manner. For the other g j , of degree at least 2, it turns out that the roots of g j ξ e j /π i x have order dividing π i , and thus Algorithms 17 or 19 produce the remaining roots of h ˜i-1 . This shows the correctness.

First all, the primitive roots of unity of orders π i and ρ needed during the execution of the algorithm can be obtained in time O ˜(log 2 q) by Remark 25. The Graeffe transforms in step 1 may execute in time (σ + ρ) O ˜(d log q log (µ ρ)) by Proposition 1.

Let d i := deg h i , let δ 1,i be the number of simple roots of h i , let δ 2,i be the number of multiple roots of h i , and let ∆ 2,i = d iδ 1,i be the number of multiple roots of h i counted with multiplicities. By Proposition 26, computing the squarefree part h ˜i of h i takes time O ˜(d log q) + ∆ 2,i O ˜(log 2 q). From deg h ˜i = δ 1,i + δ 2,i and deg h ˜i-1 = d i , we obtain deg h ˜i-1 -deg h ˜i = d iδ 1,iδ 2,i = ∆ 2,iδ 2,i . Summing these equalities over i leads to deg h ˜0 -deg h ˜m = i=1 m ∆ 2,i -i=1 m δ 2,i . Using ∆ 2,i 2 δ 2,i , deg h ˜0 = d, and deg h ˜m = 1, we deduce that i=1 m δ 2,i d -1, and then that i=1 m ∆ 2,i 2 (d -1). Consequently, the total cost for computing all the squarefree parts drops to O ˜(d log 2 q).

Step 3.a may take π i O ˜(d log q log π i ) by Proposition 16, which yields a total cost σ O ˜(d log q log µ). In step 3.c, Algorithm 23 is called O(δ 2,i ) times. Consequently, the total cost of this step is σ d O ˜(log q) + d O ˜(log 2 q) by Proposition 24. If we use Algorithm 17 in step 3.d, then the cost is ⌈π i / d⌉ O ˜(d log q) = π i O ˜ log O(1) d log q + O ˜(d log q) by Proposition 18. The sum of these costs is bounded by σ O ˜(d log q). So far, this establishes the first complexity bound.

From now on assume that d 1+δ q -1. The cost for the Graeffe transforms in step 1 may drop into (d log 2 q) 1+δ by Corollary 9. Step 3.a may run in time (d log q) 1+δ + O ˜(d log π i log q) by Proposition 16. The total cost of this step thus amounts to (d log 2 q) 1+δ . Putting these bounds together yields the second complexity.

For the third complexity bound, we use Algorithm 23 combined with Algorithm 21 in step 3.c, so that Proposition 24 ensures a time µ √ d O ˜(log 2 q). In addition we use Algorithm 19 in step 3.d, and Proposition 20 gives a total cost O ˜( µ d √ log 2 q) + O ˜(d log 2 q). Recall that S 1 (q -1) represents the largest prime number in the irreducible factorization of q -1. If S 1 (q -1) = O(log q) and if σ + ρ = O(log q), then the first complexity bound of Proposition 28 rewrites into O ˜(d log 2 q), which is thus softly equivalent to the average cost of Cantor-Zassenhaus' algorithm. We are now able to state the main result of this section.

Theorem 29. Given a constant δ > 0, the irreducible factorization of q -1, a primitive element of F q * , and f ∈ F q [x] of degree d, monic, separable, which splits over F q , the roots of f can be computed in time O ˜ S 1 (q -1) d log 2 q + (d log 2 q) 1+δ .

Proof. Without loss of generality we can assume that f (0) = / 0 and that d q -1. If d 1+δ > q -1, then it suffices to evaluate f at all the elements of F q * in time O ˜(d + q) = O(d 1+δ ). Consequently, for when d 1+δ q -1, we appeal to the preceding proposition.

Let us finally mention that, in general, the present method improves [16, Theorem 4.1] (see complexity details at the end of the proof). Nevertheless, if k = 1, and if S 1 (p -1) is close to p, then Theorem 29 does not improve [START_REF] Shoup | A fast deterministic algorithm for factoring polynomials over finite fields of small characteristic[END_REF]Theorem 4.1(2)].

Bibliography

Proposition 6 .

 6 Let e ∈ A, and let κ = ⌊log d/log p⌋ + 1. From Tr(e 1 ), ..., Tr(e d ), we can compute the characteristic polynomial h of e modulo p in time O(M p κ ,k (d)) = O ˜(d log q).

  monic of degree k and irreducible modulo p. Let Tr ε represent the trace function in A ε . Proposition 10. Let δ > 0 be a fixed constant. If d 1+δ q -1, then the traces Tr ε (e + ε e ¯), Tr ε ((e + ε e ¯)2 ), ..., Tr ε ((e + ε e ¯)d ) modulo ε 2 , for a given e + ε e ¯∈ A ε , can be computed in time (d log(q κ )) 1+δ .

  E with the list [0]. 3. For i from m down to 1 do a. Replace E by the concatenation of [(e + j χ)/ π i | e ∈ E] for j from 0 to π i -1.

  Over any field K the separable factorization of a univariate polynomial g of degree d can computed by means of O(M(d) log d) operations in K[START_REF] Lecerf | Fast separable factorization and applications[END_REF]. If K is the finite field F q then it is possible to deduce the squarefree factorization of g with O(d) additional extractions of p-th roots[START_REF] Lecerf | Fast separable factorization and applications[END_REF] Section 4]. This result is also left as an exercise in [17, Chapter 14, Exercise 14.30]. In the special case when k = 1 and p > d, no root extraction is necessary and the squarefree part of g is obtained as g / gcd(g, g ′ ).

Bruno Grenet, Joris van der Hoeven, Grégoire Lecerf

A. Partially supported by a LIX-Qualcomm -Carnot postdoctoral fellowship.


Bruno Grenet
email: bruno.grenet@lirmm.fr

Joris Van Der Hoeven

Grégoire Lecerf
email: lecerf@lix.polytechnique.fr

Randomized root finding over finite fields using tangent Graeffe transforms

Keywords: F.2.1 [ANALYSIS OF ALGORITHMS AND PROBLEM COMPLEXITY]: Numerical Algorithms and Problems-Computations in finite fields; G.4 [MATH-EMATICAL SOFTWARE]: Algorithm design and analysis Algorithms, Theory Finite fields, polynomial root finding, algorithm, Mathemagix

Consider a finite field F q whose multiplicative group has smooth cardinality. We study the problem of computing all roots of a polynomial that splits over F q , which was one of the bottlenecks for fast sparse interpolation in practice. We revisit and slightly improve existing algorithms and then present new randomized ones based on the Graeffe transform. We report on our implementation in the Mathemagix computer algebra system, confirming that our ideas gain by an order of magnitude in practice.

INTRODUCTION

Let F q represent the finite field with q = p k elements, where p is a prime number, and k 1. Throughout this article, such a field is supposed to be described as a quotient of F p [x] by a monic irreducible polynomial. Let f ∈ F q [x] represent a separable monic polynomial of degree d 1 which splits over F q , which means that all its irreducible factors have degree one and multiplicity one. In this article we are interested in computing all the roots of f .

Motivation

One of our interests in root finding came from the recent design of efficient algorithms to interpolate, into the standard monomial basis, polynomials that are given through evaluation functions. This task is briefly called sparse interpolation , and root finding often turns out to be a bottleneck, as reported in [START_REF] Van Der Hoeven | Sparse polynomial interpolation in practice[END_REF]. In fact, in this case, the ground field can be chosen to be F p with p = M 2 m + 1, and where 2 m is taken to be much larger than the number of terms to be discovered. In practice, to minimize the size of p, so that it fits a machine register, we take M = O(log p) as small as possible. A typical example is p = 7 • 2 26 + 1. We informally refer to such primes as FFT primes.

Root finding over prime finite fields critically occurs during the computation of integer and rational roots of polynomials in Q[x], both for dense and lacunary representations. Yet other applications concern cryptography and error correcting codes. Nevertheless, practical root finding has received only moderate attention so far, existing algorithms with good average complexity bounds often being sufficient [START_REF] Kaltofen | Subquadratic-time factoring of polynomials over finite fields[END_REF][START_REF] Kedlaya | Fast modular composition in any characteristic[END_REF].

In this article, we focus on fast probabilistic root finding algorithms, targeting primarily FFT prime fields and, more generally, finite fields whose multiplicative group has smooth cardinality. At a second stage, we report on practical efficiency of our new algorithms within the Mathemagix computer algebra system [START_REF] Van Der Hoeven | [END_REF].

Notations and prerequisites

The multiplicative group F q \ {0} of F q is written F q * . In order to simplify the presentation of complexity bounds, we use the soft-Oh notation: (1) g(n) (we refer the reader to [START_REF] Zur Gathen | Modern computer algebra[END_REF]Chapter 25,Section 7] for technical details). The least integer larger or equal to x is written ⌈x⌉. The largest integer smaller or equal to x is written ⌊x⌋. The remainder of g in the division by f is denoted by g rem f .

f (n) ∈ O ˜(g(n)) means that f (n) = g(n) log O
We write M: N → Z for a function that bounds the total cost of a polynomial product algorithm in terms of the number of ring operations performed independently of the coefficient ring, assuming a unity is available. In other words, two polynomials of degrees at most d over such a ring A can be multiplied with at most M(d) arithmetic operations in A. The schoolbook algorithm allows us to take M(d) = O(d 2 ). On the other hand the fastest known algorithm, due to Cantor and Kaltofen [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF], provides us with M(d) = O(d log d log log d). In order to simplify the cost analysis of our algorithms we make the customary assumption that M(d 1 ) / d 1 M(d 2 ) / d 2 for all 0 < d 1 d 2 . Notice that this assumption implies the super-additivity of M, namely

M(d 1 ) + M(d 2 ) M(d 1 + d 2 ) for all d 1 0 and d 2 0.
For operations in Z and in finite fields, we are interested in the Turing machine model , with a sufficiently large number of tapes. In short, we use the terms bit-cost and bit-complexity to refer to this model whenever the context is not clear. For randomized algorithms, we endow Turing machines with an additional instruction which generates random bits with a uniform distribution [START_REF] Papadimitriou | Computational Complexity[END_REF].

We write I(n) for a function that bounds the bit-cost of an algorithm which multiplies two integers of bit-sizes at most n, viewed in classical binary representation. Recently, the best bound for I(n) has been improved to I(n) = O(n log n 8 log * n ), where log * represents the iterated logarithm function [START_REF] Harvey | Even faster integer multiplication[END_REF]. Again, we make the customary assumption that I(n 1 ) / n 1 I(n 2 ) / n 2 for all 0 < n 1 n 2 . We freely use the following classical facts: ring operations in F p cost O(I(log p)) and one division or inversion in F p costs O(I(log p) log log p) [START_REF] Zur Gathen | Modern computer algebra[END_REF]Chapter 11].

For polynomial operations over F q , we let M q (d) represent a function that bounds the bit-cost of an algorithm that multiplies two polynomials of degrees at most d, with the same kind of assumptions as for M and I. According to [START_REF] Harvey | Faster polynomial multiplication over finite fields[END_REF], we may take M q (d) = O(d log q log(d log q) 8 log * (dlog q) ). The ring operations in F q cost at most O(M p (k)), and inversions take at most O(M p (k) log k + I(log p) log log p). For convenience, m q and d q will respectively denote cost functions for the product and the inverse in F q .

Let us recall that the gcd of two polynomials of degrees at most d over F q can be computed in time O(M q (d) log d): One can for instance use pseudo-remainders in [START_REF] Zur Gathen | Modern computer algebra[END_REF]Algorithm 11.4], mutatis mutandis. Given monic polynomials f and g 1 , ..., g l with deg f = d and deg

g 1 + ••• + deg g l = O(d)
, all the remainders f rem g i can be computed in time O(M q (d) log l) using a subproduct tree [START_REF] Zur Gathen | Modern computer algebra[END_REF]Chapter 10]. The inverse problem, called Chinese remaindering , can be solved within a similar cost O(M q (d) log l + d d q ). In this article, when needed, we consider that the factorization of q -1 and a primitive element of F q * have been precomputed once, and we discard the necessary underlying costs. In practice, if the factorization of q -1 is given, then it is straightforward to verify whether a given element is primitive. For known complexity bounds and historical details on these tasks, we refer the reader to [START_REF] Bach | Comments on search procedures for primitive roots[END_REF][START_REF] Mullen | Handbook of Finite Fields[END_REF][START_REF] Shoup | Searching for primitive roots in finite fields[END_REF].

Related work and our contributions

Seminal algorithms for polynomial factorization over finite fields are classically attributed to Berlekamp [START_REF] Berlekamp | Factoring polynomials over finite fields[END_REF][START_REF] Berlekamp | Factoring polynomials over large finite fields[END_REF], and Cantor and Zassenhaus [START_REF] Cantor | A new algorithm for factoring polynomials over finite fields[END_REF], but central earlier ideas can be found in works of Gauss, Galois, Arwins, Faddeev and Skopin. Cantor-Zassenhaus' algorithm is randomized and well suited to compute roots of polynomials of degree d that split over F q in average time O(M q (d) (log q + log d) log d + d d q ). Of course, if q = O(d) then an exhaustive search can be naively performed in time O(M q (d) log d) (the factor log d can be discarded if a primitive element of F q * is given, by means of [6, Proposition 3]), so that the cost of root finding simplifies to O(M q (d) log q log d + d d q ). This classical approach is for instance implemented in the NTL library written by Shoup [START_REF] Shoup | NTL: A Library for doing Number Theory[END_REF]. However neither Berlekamp's nor Cantor-Zassenhaus' algorithm seems to benefit of particular prime numbers such as FFT primes. Instead, alternative approaches have been proposed by Moenck [START_REF] Moenck | On the efficiency of algorithms for polynomial factoring[END_REF], von zur Gathen [START_REF] Gathen | Factoring polynomials and primitive elements for special primes[END_REF], Mignotte and Schnorr [START_REF] Mignotte | Calcul déterministe des racines d'un polynôme dans un corps fini[END_REF], and then by Rónyai [START_REF] Rónyai | Factoring polynomials modulo special primes[END_REF].

In Section 2, for the sake of comparison, we first revisit Cantor-Zassenhaus' approach and propose a practical trick to slightly speed it up. We also briefly recall the complexity bound of Mignotte-Schnorr's algorithm [START_REF] Mignotte | Calcul déterministe des racines d'un polynôme dans un corps fini[END_REF].

We then design and analyze fast randomized algorithms based on tangent Graeffe transforms, leading to an important speed-up for FFT primes. The practical efficiency of the new algorithms is discussed on the basis of implementations in the Mathemagix computer algebra system [START_REF] Van Der Hoeven | [END_REF], thus revealing that our fastest variant gains an order of magnitude over other existing software.

Let us finally mention that our present randomized algorithms admit deterministic counterparts which are studied in [START_REF] Grenet | Deterministic root finding over finite fields using Graeffe transforms[END_REF].

KNOWN ALGORITHMS REVISTED

In this section we revisit efficient root finding algorithms previously described in the literature, and slightly improve their complexity analysis. The first stream, followed by Berlekamp [START_REF] Berlekamp | Factoring polynomials over large finite fields[END_REF], Cantor and Zassenhaus [START_REF] Cantor | A new algorithm for factoring polynomials over finite fields[END_REF], consists in splitting the input polynomial by computing suitable modular exponentiations and gcds. The second stream, followed by Moenck [START_REF] Moenck | On the efficiency of algorithms for polynomial factoring[END_REF], and Mignotte and Schnorr [START_REF] Mignotte | Calcul déterministe des racines d'un polynôme dans un corps fini[END_REF], takes advantage of the smoothness of q -1.

Cantor-Zassenhaus' algorithm

In this subsection we suppose that q has the form q = χ ρ + 1, with χ 2. The following randomized algorithm extends Cantor-Zassenhaus' one, which corresponds to χ = 2, when q is odd. We need a primitive root ξ of unity of order χ.

Randomized algorithm 1.

Input. f ∈ F q [x] of degree d 1, monic, separable, which splits over F q , and such that f (0) = / 0; a primitive root ξ of unity of order χ. Output. The roots of f . 1. If d = O(χ) then compute the roots of f by calling a fallback root finding algorithm. 2. Pick up g ∈ F q [x] at random of degree at most d -1. 3. Compute h := g ρ rem f , and set f 0 := f . 4. For all i from 1 to χ -1, compute f i := gcd(hξ i , f 0 ), make it monic, and replace f 0 by f 0 / f i . 5. Recursively call the algorithm with those of f 0 , ..., f χ-1 which are not constant, and return the union of the sets of their roots.

In the classical case when χ = 2, then we recall that step 3 costs O(M q (d) log q) and step 4 takes O(M q (d) log d) plus one inversion in F q . Since the average depth of the recursion is in O(log d) [START_REF] Ph | A branching process arising in dynamic hashing, trie searching and polynomial factorization[END_REF], the total average cost amounts to O(M q (d) (log q + log d) log d + d d q ). For the case of arbitrary χ, we propose an informal discussion, so as to justify the interest in the algorithm. If α is a root of f , then h(α) = g(α) ρ has order dividing χ, hence is a power of ξ, or zero with a very small probability. Therefore we have

f = f 0 f 1 ••• f χ-1 .
Since g is taken at random, we might expect that the values of h(α) are distributed uniformly among the powers of ξ (we discard cases when h(α) = 0). In other words, the depth of the recursive calls is expected to be in O(log d/log χ).

If we further assume that log χ = o(log q), then step 3 takes approximately C 1 M q (d) log q for a certain constant C 1 , and step 4 costs C 2 χ M q (d) log d. Discarding the cost of the inversions, and compared to the case χ = 2, we achieve an approximate speed-up of

(C1 log q + 2 C2 log d) Mq(d) (log d)/log 2 (C 1 log q + χ C 2 log d) M q (d) (log d)/log χ = log χ log 2 C1 log q /log d + 2 C2 C 1 log q /log d + χ C 2 . Whenever C 1 log q / log d ≫ χ C 2 , this speed-up is of order log χ / log 2. In general, the speed-up is maximal if χ (log χ -1) = C1 log q C2 log d .

Mignotte-Schnorr's algorithm

Assuming given a primitive element of F q * , Mignotte and Schnorr proposed a general deterministic root finding algorithm in [START_REF] Mignotte | Calcul déterministe des racines d'un polynôme dans un corps fini[END_REF], that is efficient when q -1 is smooth.

a fair comparison with our new algorithm presented in the next section (and since their article is written in French) we recall their method in a different and concise manner.

Let π 1 , ..., π m be integers 2, such that q

-1 = ρ π 1 ••• π m , and let χ = π 1 ••• π m and σ = π 1 + ••• + π m . For instance, if the irreducible factorization of q -1 = p 1 m1 ••• p r mr is known, then we may take ρ = 1 and π 1 = π 2 = ••• = π m1 = p 1 , π m1+1 = π m1+2 = ••• = π m1+m2 = p 2 , ..., and set m = m 1 + ••• + m r .
In order to split f into factors of degrees at most ρ, we may use the following algorithm.

Algorithm 2. Input. f ∈ F q [x
] of degree d 1, monic, separable, which splits over F q , and such that f (0) = / 0; a primitive root ξ of unity of order χ.

Output. (f 1 , e 1 ), ..., (f s , e s ) in (F q [x] \ F q ) × {0, ..., χ -1} such that f = f 1 ••• f s ,
and the f i are monic, separable, and divide x ρξ ei , for all i.

1. Let h 0 (x) := x ρ rem f (x), and

compute h i (x) := h i-1 (x) πi rem f (x), for all 1 i m -1.
2. Initialize F with the list [(f , 0)].

3. For i from m down to 1 do a. Compute h i-1 rem g for all pairs (g, e) in F using a subproduct tree.

b. Initialize G with the empty list, and for all j from 0 to π i -1 do For each pair (g, e) in F do Compute g j := gcd((h i-1 rem g)ξ (e+ jχ)/πi , g). If g j is not constant, then make it monic and append (g j , (e + j χ)/π i ) to G.

Let F := G. 4. Return F . Lemma 3. Algorithm 2 is correct and executes in time O(σ M q (d) log d + M q (d) log q + m d I(log χ) + m m q d log χ + d d q ) = σ O ˜(d log q) + O ˜(d log 3 q).
Proof. We prove the correctness by descending induction on i from m down to 1. In fact, at the end of iteration i of the loop of step 3, we claim that f = (g,e)∈F g, that g divides x ρπ1•••πi-1ξ e for all (g, e) ∈ F , and that the integers e in F are divisible by π 1 ••• π i-1 and bounded by χ -1. By construction, these properties are all true with i = m + 1 when entering step 3, so that by induction we can assume that these properties hold for i + 1 when entering the loop at level i. Let (g, e) ∈ F , and let α be a root of g. From α ρπ1•••πi = ξ e , since e is divisible by π i , there exists j ∈ {0, ..., π i -1} such that α ρπ1•••πi-1 = ξ e+ jχ π i . We are done with the correctness.

As to the complexity, we can compute the χ / π i for i ∈ {1, ..., m} in time O(m I(log χ)) as follows: we first compute

π 1 , π 1 π 2 , π 1 π 2 π 3 , ..., π 1 π 2 ••• π m-1 , and π m , π m-1 π m , π m-2 π m-1 π m , ..., π 2 ••• π m-1 π m , and then deduce each χ/π i by multiplying π 1 ••• π i-1 and π i+1 ••• π m .
Step 1 requires time O(M q (d) log q). In step 3.a, since the sum of the degrees of the polynomials in F is at most d, and since these polynomials are monic, this step can be done in time O(M q (d) log d). The cost of the gcds in step 3.b is O(π i M q (d) log d) by the super-additivity of M q . We compute all the e/ π i in time O(d I(log χ)), and then all the ξ e/πi and ξ χ/πi by means of O(d log χ) operations in F q . Deducing all the ξ (e+ jχ)/πi takes O(π i d) additional products in F q .

Since the cardinality of F cannot exceed d, the total number of proper factors g j in step 3.b is at most d. Therefore, we need O(d) inversions in order to make all the polynomials in F monic.

Notice that the dependence on σ is good when q -1 is smooth. For example, if ρ = 1 and σ = O(log q), then the cost of root finding via the latter proposition drops to O ˜(d log 3 q). This is higher than the average O ˜(d log 2 q) bitcost of Cantor-Zassenhaus' algorithm. Nevertheless since the term O ˜(d log 3 q) corresponds to O(d log 2 q) products in F q , with a small constant hidden in the O, Mignotte-Schnorr's algorithm is competitive for small values of log q.

the original work of Mignotte and Schnorr [START_REF] Mignotte | Calcul déterministe des racines d'un polynôme dans un corps fini[END_REF], the cost of the algorithm was estimated to

O(M(d) i=1 r m i (log q + p i log d)) operations in F q , if q -1 = p 1 m1 ••• p r mr .
Our presentation slightly differs by the use of a subproduct tree in step 3.a.

Let us briefly mention that a better bound for splitting f was achieved in [START_REF] Gathen | Factoring polynomials and primitive elements for special primes[END_REF]. Nevertheless, for finding all the roots, the method of [START_REF] Gathen | Factoring polynomials and primitive elements for special primes[END_REF] does not seem competitive in general.

Moenck's algorithm

Moenck's algorithm [START_REF] Moenck | On the efficiency of algorithms for polynomial factoring[END_REF] deals with the special case when q -1 = ρ 2 m . Let ζ be a primitive root of F q * . Let 1 i m, let f be a polynomial whose roots have orders dividing ρ 2 i , and let α ∈ F q * be one of these roots. Then either this order divides ρ 2 i-1 , or we have α ρ2 i-1 + 1 = 0, since x ρ2 i -1 = (x ρ2 i-1 + 1) (x ρ2 i -1 -1). In the latter case, we obtain

(α/ ζ 2 m-i ) ρ2 i-1 = 1. The polynomials f 1 = gcd(x ρ2 i-1 -1, f ) and f 2 (x) = (f / f 1 )(ζ 2 m -i
x) have all their roots of order dividing ρ 2 i-1 . In this way, the roots of f can be computed inductively starting from i = m down to i = 1. At the end, we are reduced to finding roots of several polynomials whose orders divide ρ. If ρ is small, then an exhaustive search easily completes the computations. Otherwise we may use a fallback algorithm. Moenck's algorithm summarizes as follows.

Algorithm 4. Input. f ∈ F q [x] of degree d 1, monic, separable, which splits over F q , and such that f (0) = / 0; a primitive element

ζ of F q * . Output. (f 1 , γ 1 ), ..., (f s , γ s ) in (F q [x] \ F q ) × F q such thatf = f 1 ••• f s ,
the f i are monic, separable, and the roots of f i (γ i x) have orders dividing ρ.

1. Compute h i (x) := x ρ2 i rem f (x) for all i ∈ {0, ..., m -1}. 2. Initialize F with the list [(f , 1)].
3. For i from m down to 1 do a. Compute the remainders of h i-1 rem g for all triples (g, γ) in F using a subproduct tree. b. Initialize G as an empty list. c. For all (g, γ) in F do i. Compute g 1 (x) := gcd((h i-1 rem g)(x)γ ρ2 i-1 , g(x)). If g 1 is not constant then make it monic and append (g 1 , γ) to G. ii. Compute g 2 := g / g 1 . If g 2 is not constant then make it monic, and append (g

2 , γ ζ 2 m-i ) to G. d. Replace F by G. 4. Return F .
Our presentation differs from [START_REF] Moenck | On the efficiency of algorithms for polynomial factoring[END_REF]. We also introduced step 3.a, which yields a sensible theoretical speed-up, similarly to Mignotte-Schnorr's algorithm. In fact Moenck's and Mignotte-Schnorr's algorithms are quite similar from the point of view of the successive splittings of f , and the intermediate operations. We thus leave out the details on correctness and cost analysis. As an advantage, the logarithms of the successive roots are not needed. Nevertheless computing all the powers of γ in steps 3.c amounts to a bit-cost O ˜(d log 3 q).

USING GRAEFFE TRANSFORMS

One advantage of Cantor-Zassenhaus' algorithm is its average depth in the recursive calls in O(log d). This is to be compared to the O(log q) iterations of Algorithms 2 and 4.

In this section we propose a new kind of root finder based on the tangent Graeffe transform, which takes advantage of a FFT prime field, with an average cost not exceeding the one of Cantor-Zassenhaus.

Generalized Graeffe transforms

Classically, the Graeffe transform of a polynomial g ∈

F q [x] of degree d is the unique polynomial h ∈ F q [x] sat- isfying h(x 2 ) = g(x) g(-x). If g(x) = i=1 d (α i -x), then h(x) = i=1 d (α i 2 -x)
. This construction can be extended to higher orders as follows: the generalized Graeffe transform of g of order π, written G π (g), is defined as the resultant

G π (g)(x) = (-1) πd Res z (g(z), z π -x). If g = i=1 d (α i -x), then G π (g)(x) = i=1 d (α i π -x).
Equivalently, G π (g) is the characteristic polynomial of multiplication by x π in F q [x]/(g) (up to the sign). For our root finding algorithms, the most important case is when q -1 is smooth and the order π of the generalized Graeffe transform divides q -1.

Proposition 5. Let π 1 , ..., π m be integers 2, such that χ = π 1 ••• π m divides q -1, and let ξ i be given primitive roots of unity of orders π i , for all i ∈ {1, ..., m}. If g is a monic polynomial in F q [x] of degree d, then the generalized Graeffe transforms of orders π 1 , π 1 π 2 , π 1 π 2 π 3 , ..., χ of g can be computed in time O(m M q (µ d)) or O(M q (σ d)), where µ = max (π 1 , ..., π m ) and

σ = π 1 + ••• + π m . Proof. Writing g(x) = c j =1 d (α j -x) in an algeb- raic closure of F q , the Graeffe transform of g of order π i is h i (x) = c πi j =1 d (α j πi -x)
. Consequently this leads to

h i (x πi ) = g(x) g(ξ i x) g(ξ i 2 x) ••• g(ξ i πi-1 x).
Using the latter formula, by Lemma 6 below, the transform can be obtained in time O(M q (π i d)). Taking the sum over i concludes the proof. Lemma 6. Let g be a polynomial of degree d 1 in F q [x], let α ∈ F q , and let l be an integer. Then the product P l (x) = g(x) g(α x) g(α 2 x) ••• g(α l-1 x) can be computed in time O(M q (l d)).

Proof. Let h := ⌊l / 2⌋. If l is even, then we have P l (x) = P h (x) P h (α h x), otherwise we have P l (x) = P h (x) P h (α h x) g(α l-1 x). These formulas lead to an algorithm with the claimed cost.

Let us mention that good complexity bounds for Graeffe transforms for general finite fields can be found in [START_REF] Grenet | Deterministic root finding over finite fields using Graeffe transforms[END_REF].

Tangent Graeffe transforms

Introducing a formal parameter ε with ε 2 = 0, we define the generalized tangent Graeffe transform of g of order π as being G π (g

(x + ε)) ∈ (F q [ε] / (ε 2 ))[x].
For any ring R, computations with "tangent numbers" in R[ε] / (ε 2 ) can be done with constant overhead with respect to computations in R (in the FFT model, the overhead is asymptotically limited to a factor of two). Whenever a Graeffe transform preserves the number of distinct simple roots, the tangent Graeffe transform can be used to directly recover the original simple roots from the transformed polynomial, as follows:

Lemma 7. Let g ∈ F q [x
] be separable of degree d, let π be coprime to p, and let h(x)

+ h ¯(x) ε + O(ε 2 ) = G π (g(x + ε)).
A nonzero root β of h is simple if, and only if, h ¯(β) = / 0. For such a root, π β h ′ (β) / h ¯(β) is the unique root α of g such that α π = β.

Proof. Let α be a root of g in an algebraic closure of F q . The lemma follows from the formula h ¯(α π ) = π α π -1 h ′ (α π ), obtained by direct calculation.

In the context of the Graeffe method, the tangent transform is classical (for instance, see [START_REF] Malajovich | Tangent Graeffe iteration[END_REF][START_REF] Pan | Solving a polynomial equation: Some history and recent progress[END_REF] for history, references, and use in numerical algorithms). The generalized tangent Graeffe transform can also be seen as the tangent characteristic polynomial of x π modulo g(x + ε), and this construction is often attributed to Kronecker in algebraic geometry [START_REF] Kronecker | Grundzüge einer arithmetischen Theorie der algebraischen Grössen[END_REF].

Overview of our methods

Let us write q -1 = ρ π 1 ••• π m . Let h 0 := f , and let h 1 , ..., h m be the Graeffe transforms of orders π 1 , π 1 π 2 , ..., π 1 ••• π m of f . The roots of h m are to be found among the elements of order ρ. One may then find the roots of h m-1 among the π m -th roots of the roots of h m , and, by induction, the roots of h i-1 are to be found among the π i -th roots of the roots of h i . This is the starting point of the efficient deterministic algorithm designed in [START_REF] Grenet | Deterministic root finding over finite fields using Graeffe transforms[END_REF]. But in order to make it much more efficient than Cantor-Zassenhaus's algorithm we introduce and analyze two types of randomizations in the next sections. First we avoid computing π i -th roots by combining random variable shifts and tangent transforms. Second we analyze the behaviour for random polynomials, which leads us to a very efficient heuristic algorithm.

A randomized algorithm

Our randomized algorithms can be competitive to Cantor-Zassenhaus' algorithm only for very smooth values of q -1. For simplicity we focus on the important case of an FFT field where q = M 2 m + 1, with M = O(log q).

We introduce the parameter τ ∈ F q . Let χ = M 2 m-l be a divisor of q -1, and let ρ = (q -1)/ χ = 2 l . Let α 1 , ..., α d denote the roots of f . The Graeffe transform g τ of f (xτ ) of order ρ equals i=1 d

((α i + τ ) ρ -x).
Given a pair of distinct roots α i and α j of f , we have (τ + α i ) ρ = / (τ + α j ) ρ for all but at most ρ values of τ . Therefore for all but at most

d (d -1) 2
ρ values of τ , the polynomial g τ has no multiple root. Considering that τ is taken uniformly at random in F q , the probability that g τ has multiple roots is at most

d (d -1) 2 1
χ . This yields the following randomized algorithm.

Randomized algorithm 8.

Input. f ∈ F q [x] of degree d 1, monic, separable, which splits over F q , and such that f (0) = / 0; a primitive root ζ of F q * . Output. The roots of f . 1. If d (d -1) q then set l to 0. Otherwise set l to the largest integer in {0, ..., m} such that d (d -1) M 2 m-l . Compute χ := M 2 m-l and ρ := (q -1)/ χ = 2 l . 2. Pick τ ∈ F q uniformly at random, and compute h 0 (x) + ε h ¯0(x) := f (xτ + ε), where ε 2 = 0.

Recursively compute the Graeffe transform h

i + ε h ¯i of order 2 of h i-1 + ε h ¯i-1 , for all 1 i m. 4. Compute the list E := [j 2 m | j ∈ {0, ..., M -1}, h m (ζ j 2 m ) = 0] of ζ-
:= [e| e ∈ E , h ¯l(ζ e ) = / 0]. 8. Compute Z 1 := [2 l ζ e h l ′ (ζ e )/h ¯l(ζ e ) -τ | e ∈ E 1 ]. Add τ to Z 1 if f (τ ) = 0. 9. Compute f 2 (x) := f (x)/ a∈Z1 (x -a). 10. Compute recursively the roots Z 2 of f 2 . 11. Return Z 1 ∪ Z 2 .
Proposition 9. Algorithm 8 is correct, and takes an average time O ˜(d log 2 q), if d < q and q = M 2 m + 1 with M = O(log q).

Proof. The correctness follows from Lemma 7, which asserts that Z 1 is a subset of the roots of f .

Step 3 takes O ˜(d log 2 q) by Proposition 5. Steps 2, 5.a, 5.b, 6, 7, 8, 9 execute in time O ˜(d log 2 q). Step 4 costs

O ˜((d + M ) log q) = O ˜(d log q + log 2 q). If d (d -1)
q then the number of iterations in step 5 is m = O(log q) = O(log d). Otherwise the number of iterations is ml = O(log d). Consequently, the total cost of all steps but step 10 is O ˜(d log 2 q). From the choice of l, we have already seen that the degree of f 2 equals 0 with probability at least 1 / 2. Writing T (d) for the average execution time, we thus have

T (d) O ˜(d log 2 q) + 1 2 T (d), which implies T (d) = O ˜(d log 2 q).
Since the base field supports FFTs it is important to perform all intermediate computations by taking advantage of sharing transforms. This technique was popularized within the NTL library for many algorithms. In addition we notice that the first few iterates of the loop of step 5 can be reduced to direct FFT computations instead of generic multi-point evaluations.

Random polynomials which split over = q

Taking a monic random polynomial of degree d which splits over F q , for a uniform distribution, is equivalent to taking its d roots at random in F q . Let α 1 , ..., α d represent the roots of f . We examine the behavior of tangent Graeffe transforms of such random polynomials.

Let P d,i be the number of monic polynomials of degree d which split over F q and have i distinct nonzero roots of order dividing χ. Such a polynomial uniquely corresponds to the choice of i distinct values among χ, namely the roots, and then of the multiplicities of these roots with sum d. Therefore we have

P d,i = χ i d -1 i -1 . One can check that i=0 d P d,i = χ + d -1 χ -1
, which corresponds to choosing a multiplicity in {0, ..., d} for each element of order dividing χ, under the constraint that the sum of the multiplicities is exactly d. The average number of roots of such polynomials is

i=0 d i Pd,i i=0 d P d,i = χ χ + d -2 χ -1 χ + d -1 χ -1 = d 1 + (d -1)/ χ .
The latter formulas are direct consequences of the classical Chu-Vandermonde identity. Assuming that χ d, the average number of distinct roots is at least d / 2. When χ is much greater than d then this lower bound is getting close to d.

Proposition 10. Let ρ be a fixed divisor of q -1, and let χ = (q -1) / ρ. Let f ∈ F q [x] be a monic separable polynomial of degree d q -1, which splits over F q , and such that f (0) = / 0. Let g be the Graeffe transform of f of order ρ. The average number of simple roots of g over all such polynomials f endowed with a uniform distribution is at least d

1 -d / χ 1 + d / χ .
Moreover, assuming that χ > 3 d the probability that g has at least 2 -

1 + d / χ 1 -d / χ d simple roots is at least 1/2.
Proof. The polynomial g is uniformly distributed in the set of monic polynomials with roots of order χ. Let s and m be the respective numbers of simple and multiple roots of g, so that d s + 2 m. We have just seen that the average number s + m of distinct roots of g is at least

d 1 + d / χ . Hence, s + (d -s)/2 d 1 + d / χ and s 2 d 1 + d / χ -d = d 1 -d / χ 1 + d / χ . Setting K = 1 -d / χ
1 + d / χ and assuming that χ > 3 d, so that 2 > 1 / K > 1, let P be the probability that g has less than (2 -K) d simple roots. Then the average number of roots K d is at most

P (2 -1 / K) d + (1 -P ) d, whence (K -1) d P (1 -1 / K) d. We conclude that P (K -1)/(1 -1/K) = K > 1/2.

A heuristic randomized algorithm

Let χ be a divisor of q -1, ρ = (q -1) / χ and consider a polynomial f ∈ F q [x] of degree d. In Section 3.4, we have that for all but at most

d (d -1) 2
ρ values of τ , the Graeffe transform g τ of order ρ of f (xτ ) has no multiple root. Taking χ > d 2 , this implies that g τ has no multiple roots with probability at least 1 / 2, when picking τ at random. Now Proposition 10 implies that at least one third of the roots remain simple with probability at least 1 / 2 under the weaker condition that χ > 4 d and for random f . It is an interesting question whether this fact actually holds for any f if τ is randomly chosen:

Heuristic. Let q -1 = ρ χ with χ > 4 d and let f ∈ F q [x]
be a monic, separable polynomial of degree d which splits over F q . Then there exist at least q / 2 elements τ ∈ F q [x] such that the Graeffe transform of order ρ of f (xτ ) has at least d/3 simple roots.

From now on, assume that f is picked at random or that the above heuristic holds and τ is chosen at random. Taking χ > 4 d, we may then find the roots of g τ by computing a few discrete Fourier transforms of length O(d), instead of more expensive multi-point evaluations. Using tangent Graeffe transforms, the simple roots can be lifted back for almost no cost. These considerations lead to the following efficient randomized algorithm.

Randomized algorithm 11.

Input. f ∈ F q [x] of degree d 1, monic, separable, which splits over F q , and such that f (0) = / 0; a primitive element ζ of F q * .

Output. The roots of f .

1. If d M 2 m-3 then evaluate f on all the elements of F q * , return the roots.

2. Otherwise set l to be the largest integer in {1, ..., m -2} such that d < M 2 m-l-2 . Compute χ := M 2 m-l and ρ := (q -1)/ χ = 2 l .

3. Pick a random τ in F q , and compute the Graeffe transform h l (x) + ε h ¯l(x) of order ρ of f (xτ + ε) with ε 2 = 0.

4. Compute [h l (ζ j 2 l )| j ∈ {0, ..., χ -1}], [h l ′ (ζ j 2 l )| j ∈ {0, ..., χ -1}], [h ¯l(ζ j 2 l )| j ∈ {0, ..., χ -1}]. 5. Compute Z 1 := 2 l ζ j 2 l h l ′ ζ j 2 l / h ¯l ζ j 2 l -τ | j ∈ {0, ..., χ -1}, h l (ζ j 2 l ) = 0, h ¯l(ζ j 2 l ) = / 0 . Add τ to Z 1 if f (τ ) = 0.
6. Compute f 2 (x) := f (x)/ a∈Z1 (xa).

Compute recursively the roots

Z 2 of f 2 . 8. Return Z 1 ∪ Z 2 .
Proposition 12. Assume the heuristic or that f is chosen at random among the monic, separable polynomials with f (0) = / 0 and which split over F q . Then Algorithm 11 takes an average time O(M q (d) log q) = O ˜(d log 2 q), whenever d < q and q = M 2 m + 1 with M = O(log q). Proof. Of course if d > M 2 m-3 then the cost is O(M q (d) log d). Now suppose we arrive in step 2. The roots of h l have order dividing χ, and we have χ / 2 4 d < χ. If the algorithm terminates then it is correct by the same arguments as for Algorithm 8.

The cost of step 3 is bounded by O(M q (d) log d) (and even by O(M q (d)) if p > d according to [4, Chapter 1, Section 2]). When computing h l using l tangent Graeffe transforms of order 2, the cost of step 3 is bounded by O(M q (d) l) = O(M q (d) log q). Steps 4 and 5 can be done in time O(M q (d) log d). The average execution time T (d) thus satisfies

T (d) O(M q (d) log q) + P T (d) + (1 -P ) T (2 d/3),
where P 1 / 2 is the probability that Z 2 contains more than 2 d / 3 elements in step 7. It follows that T

(d) = O(M q (d) log q) + T (2 d / 3), whence T (d) = O(M q (d) log q).
The main bottleneck for Algorithm 11 is step 3. It would be possible to improve our complexity bound if a better algorithm were known to compute Graeffe transforms of order 2 l . In the FFT model, the tangent Graeffe transform of a polynomial of degree d can be computed in time (4 / 3 + o(1)) M q (d). In practice, this means that the asymptotic complexity of Algorithm 11 is equivalent to 4 / 3 M q (d) log 2 q.

In order to reduce the cost of the Graeffe transforms, it is relevant to choose l such that M 2 m -l -η -1 < d M 2 m-l-η for some small η 2, which is to be determined in practice as a function of d and q. In this way more time is spent in multi-point evaluations. For the multi-point evaluations in step 4, one may use M FFTs of size 2 m-l for small values of M , or appeal to Bluestein's chirp transform [START_REF] Bluestein | A linear filtering approach to the computation of discrete Fourier transform[END_REF].

TIMINGS

In this section we report on timings for the algorithms presented above. We use a platform equipped with an Intel Core T M i7-4770 CPU at 3.40 GHz, and 8 GB of 1600 MHz DDR3 . It runs the Jessie GNU Debian operating system with a Linux kernel version 3.14 in 64 bit mode. We compile with GCC [12] version 4.9.1.

We use revision 9738 of Mathemagix. For the sake of comparison to other software, we installed the NTL library version 8.0.0 [START_REF] Shoup | NTL: A Library for doing Number Theory[END_REF] (configured to benefit from GMP [START_REF] Granlund | the GNU multiple precision arithmetic library[END_REF] version 6.0.0), and Flint version 2.4.4 [START_REF] Hart | FLINT: Fast Library for Number Theory[END_REF]. For our benchmark family we simply build polynomials f of degree d from random pairwise distinct roots.

Table 1 shows timings for the FFT prime field F q = F p with p = 7 • 2 26 + 1. The row NTL corresponds to the function FindRoots, which contains an optimization with respect to Algorithm 1 with χ = 2: in fact, the polynomials g in step 2 are taken of degree 1 (see [START_REF] Shoup | A fast deterministic algorithm for factoring polynomials over finite fields of small characteristic[END_REF] for details). We implemented the same optimization in our versions of Cantor-Zassenhaus in Mathemagix. The row Flint cor-responds to the function nmod_poly_factor_equal_deg, that does not take advantage of this optimization, which mainly explains a loss of efficiency. The rows "Alg. 1" to "Alg. 11" correspond to our own implementations in Mathemagix of the above algorithms (see files algebramix/polynomial_modular_int.hpp and algebramix/bench/polynomial_modular_int_bench.cpp). Our modified version of Cantor-Zassenhaus' algorithm, with χ > 2, does not reveal to be of practical interest for this prime size.

Concerning the deterministic algorithms, let us precise that we take π

1 = ••• = π m-1 = 2 and π m = M in Algorithm 2.
In addition, in our implementation of Algorithm 4, when g and i are small, for efficiency reasons, we compute h i-1 rem g using binary powering instead of simultaneous remainders. These deterministic methods turn out to be quite competitive.

Our Table 2 is similar to Table 1 for the FFT prime field F q = F p with p = 5 • 2 55 + 1. Nevertheless NTL does not support this larger prime within its type zz_p. This larger prime starts to reveal the interest in our modified version of Cantor-Zassenhaus with χ > 2. The speed-up measured with the modified version even more increases with larger primes: With p = 7 • 2 120 + 1, it reaches a factor 1.6 in size 2 14 . Graeffe transforms as in Algorithm 8 is not very competitive. However, in final, it is satisfactory to observe that our Algorithm 11, exploiting the tangent transform, gains by an order of magnitude over other existing methods in both tables. Coming back to our initial motivation for sparse polynomial interpolation, let us mention that plugging Algorithm 11 into our implementations reported in [START_REF] Van Der Hoeven | Sparse polynomial interpolation in practice[END_REF] leads to quite good speed-ups. The total time for Example 1 of [START_REF] Van Der Hoeven | Sparse polynomial interpolation in practice[END_REF], concerning the product of random sparse polynomials, decreases from 12 s to 7.6 s with the "coefficient ratios" technique. In fact, the time spent in root finding decreases from 50 % to 8 % during the sole determination of the support.

d
Example 2 of [START_REF] Van Der Hoeven | Sparse polynomial interpolation in practice[END_REF] concerns the sparse interpolation of the determinant of the generic n × n matrix. For n = 8, the time for the "Kronecker substitution" technique with a prime of 68 bits decreases from 95 s to 43 s. The time elapsed in root finding during the sole determination of the support drops from 70 % to 20 %. Thanks to additional optimizations, the "coefficient ratios" technique now runs this example in 46 s.

  logarithms of the roots of h m . 5. For i from m down to l + 1 do a. Replace E by the concatenation of [(e + j χ) / 2| e ∈ E] for j ∈ {0, 1}. b. If E has cardinality more than d then remove the elements e from E such that h i-1 (ζ e ) = / 0. 6. If l = 0 then return {ζ eτ | e ∈ E }. 7. Compute E 1

Table 1 .

 1 Randomized root finding in degree d over Fp, with p = 7 • 2 26 + 1, time in seconds.

	+ 1	2 8	2 9	2 10	2 11	2 12	2 13	2 14	2 15	2 16	2 17	2 18
	NTL	0.0065 0.015 0.037 0.090 0.20	0.48	1.1	2.5	5.4	12	26
	Flint	0.0093 0.025 0.065 0.17	0.45	1.2	3.0	8.0	22	63	250
	Alg. 1, Cantor-Zassenhaus, χ = 2	0.010	0.020 0.043 0.092 0.20	0.42	0.90	1.9	4.0	8.6	18
	Alg. 1, Cantor-Zassenhaus, modified 0.009	0.018 0.036 0.091 0.20	0.41	0.96	2.1	4.4	10	22
	Alg. 2, Mignotte-Schnorr	0.063	0.12	0.25	0.51	1.0	2.0	3.9	7.8	15	31	62
	Alg. 4, Moenck	0.025	0.051 0.10	0.21	0.44	0.91	1.9	3.9	8.2	17	35
	Alg. 8, Graeffe, randomized	0.003	0.006 0.021 0.074 0.25	0.70	1.7	3.6	7.5	15	29
	Alg. 11, Graeffe, heuristic	0.001	0.002 0.003 0.007 0.015 0.031 0.065 0.13	0.24	0.59	1.4
	d + 1	2 8	2 9	2 10	2 11	2 12	2 13	2 14	2 15	2 16	2 17	2 18
	Flint	0.032 0.084 0.22	0.58	1.5	3.8	9.5	23	58	150	470
	Alg. 1, Cantor-Zassenhaus, χ = 2	0.035 0.075 0.16	0.37	0.78	1.7	3.8	8.6	19	41	90
	Alg. 1, Cantor-Zassenhaus, modified 0.025 0.055 0.11	0.28	0.63	1.3	3.1	7.0	14	33	76
	Alg. 2,	0.22	0.46	0.96	2.0	4.2	8.7	18	39	83	180	370
	Alg. 4, Moenck	0.89	0.18	0.37	0.76	1.6	3.4	7.2	16	33	72	150
	Alg. 8, Graeffe, randomized	0.011 0.038 0.14	0.56	1.1	2.6	6.8	16	40	92	210
	Alg. 11, Graeffe, heuristic	0.005 0.007 0.016 0.032 0.067 0.14	0.29	0.61	1.3	2.7	5.6

Table 2 .

 2 Randomized root finding in degree d over Fp, with p = 5 • 2 55 + 1, time in seconds.

Acknowledgements

Bruno Grenet was partially supported by a LIX-Qualcomm -Carnot postdoctoral fellowship.