Jhon J Ipus

Javier Sebastián Blázquez

Sergio Lozano-Perez

Alejandro Conde
email: conde@us.es

Javier S Blázquez

Prof A Conde

Microstructural evolution characterization of Fe-Nb-B ternary systems processed by ball milling

Keywords: amorphous alloys, ball-milling, FIB, nanocrystals, TEM user supplied): amorphous alloys, ball-milling, FIB Ball milling, nanocrystalline alloys, amorphous alloys, FIB, TEM

Introduction

Fe based amorphous and nanocrystalline alloys are very interesting due to their magnetic properties [1], which are strongly linked with their metastable microstructure [2]. Mechanical alloying by ball milling has been shown as a very suitable and versatile technique in the production of both amorphous and nanocrystalline alloys in wide compositional ranges [3]. Generally, mechanical alloying is developed from a starting mixture of pure powders, which progressively become alloyed to form an amorphous and/or a supersaturated solid solution. The study of microstructural evolution is frequently followed by X-ray diffraction (XRD) [3,4,5,6,7,8,9], characterizing the phase evolution, average values of microstructural parameters (grain size and microstrains) and compositional information derived from evolution of the average lattice parameter.

However, transmission electron microscopy (TEM) studies on ball milled powder samples are not so abundant (especially for magnetic materials) [4,8,10,11,12,13,14] due to the difficulties associated with sample preparation. TEM imaging and microanalysis supply local microstructural and compositional information overcoming the limitations of XRD. This is of special interest for the detection of minoritary phases (not only due to small volume fraction but also to low scattering power).

In this study, selected samples of Fe-Nb-B alloys prepared by ball milling were studied by different TEM imaging, diffraction and microanalysis techniques in order to obtain complementary information on the microstructure previously studied by XRD and Mössbauer spectrometry [15]. Scanning electron microscopy (SEM) images obtained for short milling times show an initial mixture of the constituent elements forming heterogeneous powder particles [15].

Experimental

In that previous study, evolution of the powder during milling process indicates that the maximum compositional homogenization was achieved after 50 h under the milling conditions used. However, dispersion on composition persists being constant up to the maximum explored milling time (400 h). The ratio between Fe and Fe+Nb content after 50 h milling is 90 ± 4 % and 95 ± 8 % for Nb10 and Nb5, respectively. Contamination from milling media yields the incorporation of Cr to the powder (1.8 and 0.6 at. % Cr for Nb10 and Nb5, respectively) [15]. Different electron microscopy techniques were used in this study, including transmission (TEM), scanning transmission (STEM) and high resolution electron microscopy (HREM). Conventional bright field (BF) and dark field (DF) were used, as well as high angle annular dark field (HAADF) in STEM mode. In diffraction mode, selected area diffraction (SAD) and convergent beam electron diffraction (CBED) patterns were acquired in order to identify the different phases. Several analytical techniques were used, including energy dispersive X-ray spectroscopy (EDX), using an Oxford Instruments EDS ultra-thin window detector running Inca software, capable of a resolution of 137 eV for the Mn K α peak; electron energy loss spectroscopy (EELS) and energy filtered transmission electron microscopy (EFTEM) using a Gatan Image Filter (GIF) which provided an energy resolution of ~ 1eV. Experiments were performed using three different microscopes: a Philips CM-20 operated at 200 kV, a Philips CM200 operated at 200 kV and a Jeol JEM3000F operated at 297 kV.

Sample preparation

The preparation of TEM specimens was performed using a FIB FEI 200, following the lift out procedure [16,17,18,19], applied to individual powder particles.

The procedure to prepare TEM samples using FIB is schematically shown in figure 1 for a typical powder sample. The main steps followed are:

• A protective Pt layer is deposited on the surface of a selected powder particle to prevent Ga + damage in the region of interest during cutting and thinning processes.

• A slice of approximately 15 µm x 5 µm x 2 µm is cut from the powder particle (figure 1a-b) and welded with the help of an in-situ micromanipulator (figure 1c) to a dedicated Cu grid for TEM examination (figure 1d).

• Finally, the TEM sample is thinned to electron transparency (figure 1e-f).

Results and Discussion

Figure 2 shows BF TEM images for Nb10 alloy after 50 h milling and both alloys after 400 h milling together with a SAD pattern of the former sample. The observed microstructure for the Nb10 alloy after 50 h consists on a majority of nanocrystals with regular shape and average crystal size <D> ~10 nm (bcc-Fe type from SAD patterns in agreement with XRD results [15]). In particular, long shaped nanocrystals (Nb rich as it will be shown later in EDX maps) and larger bright inclusions in BF images (size ~100 nm) are randomly dispersed throughout the sample.

For both alloys after 400 h milling, long shaped nanocrystals do not appear and the α-Fe type nanocrystals present a smaller average crystal size, <D> ~7 nm, than that of Nb10 after 50 h milling (this agrees with XRD results, where the crystal size is just slightly reduced after 50 h milling for both alloys). Moreover, large bright inclusions (smaller than for 50 h milling) are also observed. These inclusions exhibit a very weak contrast with respect to the holes in the samples, looking very bright in BF images (sometimes even saturating the negative when contrast in the nanocrystalline matrix is obtained).

The brightness of these crystals in BF images could be ascribed to a low power scattering of the constituent elements of the crystals, characteristic of light elements as B, which is present in the studied alloys. Diffraction analysis evidences their crystalline nature and, moreover, lattice periodicity can be observed in HREM images. In addition to microstructural characterization, chemical analyses were performed for the studied samples. Figure 4 shows EDX compositional maps for the three studied samples. For all the samples, a HAADF image is shown for identification of the different regions and individual maps of Fe, Nb and Cr (due to contamination from milling media [15]). In the case of Nb5 after 400 h milling image, the Pt layer deposited during FIB sample preparation is observed and EDX map of this element (not shown) evidences the sharp edge between this layer and the region of interest of the sample studied. The elemental maps have been equalized in order to enhance individual elemental contributions. For Nb10 alloy after 50 h milling, Fe and Nb are not completely mixed. The matrix formed by regular shaped nanocrystals is Fe rich and the long shaped crystals are identified as Nb rich particles. Moreover, Fe and Nb (and even Cr) free regions are observed in the EDX maps, corresponding to the bright inclusions in BF images and confirming its mainly pure B character (see corresponding HAADF image). In fact, X-rays emitted by B atoms have very low energy and cannot be efficiently detected by EDX. Besides Fe and Nb and, as expected from previous results [15], EDX maps show the presence of Cr in the alloy due to contamination from grinding media as Cr rich clusters. It is worth mentioning the Cr-contrast between this previous EDX experiments performed in SEM (below 1 at.% after 50 h milling) [15].

Compositional maps for both alloys after 400 h milling do not show independent

Fe or Nb rich regions but both elements are homogeneously dispersed throughout the samples except for those regions where Fe and Nb are not detected, ascribed to B crystals. Hence, the milling process used is unable to completely incorporate the B into the Fe(Nb) rich matrix. B crystals were present since the starting mixture, although (in a different crystalline structure, and remain for the length of the milling time explored in this study, although its size decreases from 50 to 400 h. With the milling time, homogeneous Cr content of the matrix increases in agreement with SEM results [15] and Cr cluster are not so well defined as for the sample milled 50 h.

In order to confirm the B composition of the bright inclusions in BF images, EELS and EFTEM experiments were performed. Figure 5 (up) shows the B content corresponding to the line profile drawn on the HAADF image obtained for Nb10 after 400 h milling. As the image provides Z-contrast, the B rich inclusion appears dark and the Fe(Nb) rich matrix appears brighter. In the inclusions, B signal is high and outside the inclusions is negligible. Figure 5 (down) shows EFTEM images for Nb5 after 400 h milling filtered for B and Fe energies, respectively, along with the corresponding BF image for comparison. The B rich inclusion (left hand) and the Fe rich nanocrystalline matrix (right hand) can be clearly distinguished.

The observed B crystalline phase in this study by TEM was not detected using XRD in previous studies on the same samples [15] because of the low scattering factor of B atoms. In fact, several authors studying Fe-X-B (X=Nb, Zr,…) systems achieve a final microstructure by ball milling consisting on amorphous and/or supersaturated solid solution nanocrystals without mentioning any B rich phase [6,7,9]. On the other hand, The evolution of Nb during ball milling is typical for ductile materials [3,[START_REF] Davis | [END_REF], developing long crystals with a large surface, which facilitates the incorporation of this element into the Fe nanocrystalline matrix forming a supersaturated α-Fe(Nb) solid solution. The present results do not disagree with the preferential partitioning of Nb to intercrystalline region derived from Mössbauer results [15]. In fact, the resolution of the microscopic technique used and the sample thickness prevent any conclusion at scales below ~1 nm (intercrystalline region).

The evolution of B during ball milling is typical for brittle materials [3], forming inclusions which, in this case, can not be incorporated into the Fe(Nb) supersaturated solid solution. It is worth mentioning that B compounds as BN and B 4 C are among the harder materials known.

It is worth mentioning that using conventional XRD techniques, stationary state of the milled powder would be achieved after 50 h milling [15]. At such a time, no further evolution was detected for Nb5 alloy and a supersaturated solid solution is formed as the stationary phase. In the case of Nb10 alloy, a progressive amorphization is detected as milling progresses from 50 to 400 h, but without change in the global composition distribution of the powder particles (except for the contamination effects).

Extrapolation of crystalline volume fraction dependence on milling time to very long times leads to a saturation value for crystalline fraction of ~10 %. However, very long milling could imply recrystallization processes [22]. Results obtained in the present work show that even far beyond the time considered for stationary state using global techniques B inclusions are detected.

Conclusions

The combination of different diffraction, imaging and analytical TEM techniques has allowed the characterization of the microstructural evolution of ball milled Fe-Nb-B alloys. While Nb and Fe followed the expected evolution during the milling process, B was found to remain undissolved for long milling times. This behavior, undetected by conventional X-ray diffraction analysis, is crucial to understand the evolution of the microstructure, the level of amorphization of the Fe matrix and some physical properties.

  x-y Nb x B y (x=5, y=10 and x=10, y=15) compositions were prepared from elemental powders by ball milling in a planetary mill Fritsch Pulverisette 4 Vario. from elemental powders (>99% purity). The initial powder mass was 30 g and the ball to powder ratio was 10:1. Hardened steel balls were used with a composition 86 at. % Fe and 14 at. % Cr. The rotational speed of the disk which supports the vials was 150 rpm and that of the vials was 300 rpm in the opposite direction. The opening and closing of the vials were done under argon atmosphere in a Saffron Omega glove box to avoid oxygen and humidity contamination. For simplicity, the studied alloys will be named in the following by their Nb content: Nb10 for Fe 75 Nb 10 B 15 and Nb5 for Fe 85 Nb 5 B 10 .

Figure 3 (

 3 Figure 3 (up right) shows a HREM image of one of these crystalline inclusions

  cluster and the matrix, which indicates a very low Cr content in agreement with

  [11] detected B 2 O 3 phase in Fe(Co,Ni)ZrB system submitted to ball milling. In the present study, EELS experiments do not show any significant oxygen in the B inclusions, discarding this oxide in the systems studied in this paper.

Figure 1 .

 1 Figure captions

Figure 2 .

 2 Figure 2. From left to right SAD image of Nb10 after 50 h milling and BF TEM images

Figure 3 .

 3 Figure 3. (up left): BF image of Nb5 after 400h showing a B inclusion and the Pt layer,

Figure 4 .

 4 Figure 4. HAADF image and EDX maps of Nb10 alloy after 50h milling (left column),

Figure 5 .

 5 Figure 5. (up) HAADF image of Nb10 alloy after 400h and B profile corresponding to

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgements

This work was supported by MEC of the Spanish Government and EU FEDER (Project MAT2007-65227) and the PAI of the regional Government of Andalucía (Project P06-FQM-01823). The TEM investigations were supported by the IP3 project of the 6 th Framework Programme of the European Commission: ESTEEM Contract number 026019. J.J.I. acknowledges a fellowship from the Spanish Ministry of Education and Science. J.S.B. acknowledges a research contract from this Regional Government.


Pierrick Gaudry
email: gaudry@lix.polytechnique.fr

Éric Schost
email: eric.schost@polytechnique.fr

A low-memory parallel version of Matsuo, Chao and Tsujii's algorithm

We present an algorithm based on the birthday paradox, which is a low-memory parallel counterpart to the algorithm of Matsuo, Chao and Tsujii. This algorithm computes the group order of the Jacobian of a genus 2 curve over a finite field for which the characteristic polynomial of the Frobenius endomorphism is known modulo some integer. The main tool is a 2-dimensional pseudo-random walk that allows to heuristically choose random elements in a 2-dimensional space. We analyze the expected running time based on heuristics that we validate by computer experiments. Compared with the original algorithm by Matsuo, Chao and Tsujii, we lose a factor of about 3 in running time, but the memory requirement drops from several GB to almost nothing. Our method is general and can be applied in other contexts to transform a baby-step giant-step approach into a low memory algorithm.

Introduction

Jacobians of small genus curves have now become an important tool for publickey cryptography; computing the Zeta function of such curves remains one of the central problems in this area.

For elliptic curves, the question found a first answer with Schoof's algorithm and the subsequent improvements (see [START_REF] Blake | Elliptic curves in cryptography[END_REF] and the references therein). Furthermore, if the characteristic of the definition field is small, p-adic methods initiated with Satoh's algorithm [START_REF] Satoh | The canonical lift of an ordinary elliptic curve over a finite field and its point counting[END_REF] give a tremendous speed-up.

For higher genus curves, the p-adic methods (based on either Mestre's [START_REF] Mestre | Utilisation de l'AGM pour le calcul de E(F2n )[END_REF] or Kedlaya's [START_REF] Kedlaya | Counting points on hyperelliptic curves using Monsky-Washnitzer cohomology[END_REF] algorithms) also yield very satisfactory solutions in small characteristic, both in theory and in practice. However, in large characteristic, the question remains delicate. From the theoretical point of view, Pila's algorithm [START_REF] Pila | Frobenius maps of abelian varieties and finding roots of unity in finite fields[END_REF] and subsequent improvements [START_REF] Adleman | Counting points on curves and abelian varieties over finite fields[END_REF][START_REF] Huang | Counting points on curves over finite fields[END_REF] give polynomial time solutions, following Schoof's strategy. However, these have been turned into practical algorithms for genus 2 curves only [START_REF] Gaudry | Counting points on hyperelliptic curves over finite fields[END_REF], and it was only very recently that a Jacobian of cryptographic size was counted that way [START_REF] Gaudry | Construction of secure random curves of genus 2 over prime fields[END_REF].

In this paper, we concentrate on the last part of a Schoof-like genus 2 algorithm. We first recall the basics of such algorithms.

Let us denote by χ ∈ Z[T ] the characteristic polynomial of the Frobenius endomorphism, so that χ(1) equals the Jacobian order. A basic task is to compute the reduction of χ modulo ℓ, for sufficiently many small primes or prime powers ℓ. Such information is obtained through the study of the torsion in the Jacobian; if the characteristic is medium, the Cartier-Manin operator [START_REF] Gaudry | Counting points on hyperelliptic curves over finite fields[END_REF][START_REF] Bostan | Linear recurrences with polynomial coefficients and computation of the Cartier-Manin operator on hyperelliptic curves[END_REF] can also be used.

Thanks to Weil's bounds, one can then collect modular information until χ is known. However, it is far better to switch to a baby-step giant-step (BSGS) algorithm before the end, since collecting modular information is costly. For instance, in [START_REF] Gaudry | Construction of secure random curves of genus 2 over prime fields[END_REF], the highest prime ℓ was 19, whereas without the BSGS phase, it would have been 59.

If χ( 1) is known modulo m, then using standard BSGS techniques, the time necessary to compute the Jacobian order varies like 1/ √ m. In 2002, Matsuo, Chao and Tsujii [START_REF] Matsuo | An improved baby step giant step algorithm for point counting of hyperelliptic curves over finite fields[END_REF] showed that if we use not only the value χ(1) modulo m, but rather all coefficients of χ modulo m, it is possible to speed-up the BSGS computation. The runtime of their method (called MCT in the following) varies like 1/m, which is an important improvement. The main drawback is the space complexity: the largest example shown in [START_REF] Matsuo | An improved baby step giant step algorithm for point counting of hyperelliptic curves over finite fields[END_REF] used 12 GB of central memory, whereas the runtime (5 days on a single processor) was quite reasonable.

Standard BSGS techniques have low-memory, parallelizable, probabilistic counterparts, based on the rho or lambda (kangaroo) methods of Pollard's [START_REF] Pollard | Monte Carlo methods for index computation mod p[END_REF]: such techniques were presented in [START_REF] Van Oorschot | Parallel collision search with cryptanalytic applications[END_REF][START_REF] Teske | Speeding up Pollard's rho method for computing discrete logarithms[END_REF][START_REF] Teske | Computing discrete logarithms with the parallelized kangaroo method[END_REF][START_REF] Gaudry | Counting points on hyperelliptic curves over finite fields[END_REF][START_REF] Stein | The parallelized Pollard kangaroo method in real quadratic function fields[END_REF]. In this paper, we apply the ideas of Matsuo et al. to improve the variant of [START_REF] Gaudry | Counting points on hyperelliptic curves over finite fields[END_REF]; we obtain a probabilistic algorithm, with a heuristic complexity analysis, but which requires almost no memory and is immediately parallelizable. The expected running time is about 3 times the running time of the MCT algorithm.

In order to simplify the exposition, we concentrate on genus 2 curves. However, our idea works for more general curves as soon as the characteristic polynomial χ is known modulo some integer m. Another application of our method is the algorithm of [START_REF] Bauer | Point counting on Picard curves in large characteristic[END_REF] for Picard curves, which is a BSGS type algorithm.

The paper is organized as follows. Section 2 defines the necessary notation and recalls the BSGS algorithm of [START_REF] Matsuo | An improved baby step giant step algorithm for point counting of hyperelliptic curves over finite fields[END_REF]. Then, to introduce methods based on the birthday paradox, we start in Section 3 by a special case, where the modular information is complete enough, so we can use the method of [START_REF] Gaudry | Counting points on hyperelliptic curves over finite fields[END_REF]. In Section 4, we deal with the general case, and introduce bidimensional analogues of these techniques. Section 5 finally presents our experimental results.

MCT algorithm for genus hyperelliptic curves

Characteristic polynomial of Frobenius endomorphism

Let C be a genus 2 curve, defined over the finite field F q with q elements. The Jacobian group of C is denoted by J(C) and the characteristic polynomial of the Frobenius endomorphism is denoted by χ(T ). This polynomial has the form

χ(T ) = T 4 -s 1 T 3 + s 2 T 2 -qs 1 T + q 2 ,
where s 1 and s 2 are integers. The group order of the Jacobian is then given by #J(C) = χ(1) = q 2 + 1 -s 1 (q + 1) + s 2 . By Weil's theorem, the roots of χ have absolute value √ q and bounds on s 1 and s 2 follow directly. Better bounds can be found using the fact that the roots of χ come in pairs of complex conjugates [START_REF] Lenstra | A hyperelliptic smoothness test[END_REF]Proposition 7.1]:

|s 1 | ≤ 4 √ q, 2|s 1 | √ q -2q ≤ s 2 ≤ s 2 1 4 + 2q. ( 1
)
The values of (s 1 , s 2 ) satisfying these bounds form the hatched zone in Figure 1.

s2 s1 6q 2q -2q -6q √ q - √ q 4 √ q -4 √ q
Fig. 1. Bounds on s1 and s2

Review of the MCT algorithm

In large characteristic, point counting algorithms work by collecting modular information; this information is then recombined using the Chinese Remainder Theorem and the work is finished using some BSGS strategy. We now describe the MCT algorithm for this last step.

From now on, we assume that the characteristic polynomial χ is known modulo some positive integer m, i.e. s 1 and s 2 are known modulo m. Therefore we introduce new variables for the known and unknown parts of s 1 and s 2 :

s 1 = s 1 + m s 1 and s 2 = s 2 + m s 2 ,
so that our goal is now to find s 1 and s 2 . To this effect, a random divisor D is first picked in J(C): the strategy is to compute the order of D, hoping that it is large enough to be able to conclude (the case when J(C) is highly non-cyclic is rare in practice and easily tackled).

The order of D divides the group order χ(1), therefore we have the equality

(q 2 + 1 -s 1 (q + 1) + s 2 ) • D + (-s 1 (q + 1) + s 2 ) • m • D = 0.
To obtain a BSGS algorithm, it is usual to separate the two unknowns, one on each side of the equation. Here the unknowns s 1 , s 2 lie in intervals that are of different sizes, therefore it is necessary to split again s 2 into two pieces. Let n be a parameter to be fixed later, and let us write s 2 = t 2 + nu 2 with 0 ≤ t 2 < n.

Then χ(1) • D = 0 rewrites

q 2 + 1 -s 1 (q + 1) + s 2 + m(-s 1 (q + 1) + nu 2 ) • D = -t 2 m • D.
The algorithm proceeds as follows: first, all possible values for the right-hand side are computed and stored in a data-structure in which searching is fast. Then the left-hand side is computed for all possible values of s 1 and u 2 , until a match is found with a value of the right-hand side. For each value of s 1 , the bounds of Equation ( 1) are used to find the range of the possible values for u 2 ; a precise study of the area of the space of search leads to an optimal value n ≈ q 3/4 /m, yielding a running time of about 4q 3/4 /m operations in J(C).

3 Preliminaries: the case m ≥ 8 √ q

The main drawback of the previous method is the storage requirement. We now introduce low-memory variants and start by describing a special case, for which most features of the general treatment are already present. Suppose that χ is known modulo m, with m ≥ 8 √ q. Because |s 1 | ≤ 4 √ q ≤ m/2, the coefficient s 1 is known exactly. Corresponding to this value of s 1 , we have bounds for s 2 which yield bounds on s 2 , which is the only remaining unknown. We are thus in the setting of the search of a group order in a bounded arithmetic progression. Among the many variants inspired by Pollard's lambda (or kangaroos) method that yield a low-memory solution to this question, we chose the one based on the birthday paradox described in [START_REF] Gaudry | Counting points on hyperelliptic curves over finite fields[END_REF]. We recall it here briefly; detailed descriptions and analyses of several other variants can be found in [START_REF] Pollard | Monte Carlo methods for index computation mod p[END_REF][START_REF] Van Oorschot | Parallel collision search with cryptanalytic applications[END_REF][START_REF] Teske | Speeding up Pollard's rho method for computing discrete logarithms[END_REF][START_REF] Teske | Computing discrete logarithms with the parallelized kangaroo method[END_REF][START_REF] Stein | The parallelized Pollard kangaroo method in real quadratic function fields[END_REF].

Random search in two intersecting intervals

Let K = q 2 + 1 -s 1 (q + 1) + s 2 . Then the group order can be written #J(C) = K + m s 2 , where K is known and s 2 is unknown. Performing a suitable shift, it is possible to adjust K by some multiple of m so that the bounds on s 2 become symmetric, i.e. such that | s 2 | ≤ B for some integer B.

As before, we pick at random an element D in J(C), of presumably large order. We then define two sets of divisors that we denote W and T :

W = {(K + mσ 2 ) • D ; σ 2 ∈ [-B, B]} , T = {mσ 2 • D ; σ 2 ∈ [-B, B]} .
By definition W and T intersect, since the zero divisor is in both of them, by taking σ 2 = s 2 for W and σ 2 = 0 for T . More precisely, the size of the intersection is #(W ∩ T ) ∈ [B + 1, 2B + 1], depending on how far s 2 is from 0. The algorithm then proceeds as follows: we pick random elements uniformly alternatively in W and T . If the same divisor is obtained as an element of both W and T then (a multiple of) the order of D can be deduced.

Assuming that picking a random element in W or T has unit cost, that all elements are stored in table, and that a search in the table has unit cost, then by the birthday paradox, a collision can be obtained in expected time and space O( √ B). The constant in the O( ) can actually be made more explicit. First, for fixed s 2 , the expected running time grows like (2B + 1)π/γ, where γ = #(W ∩ T )/(2B + 1) = (2B + 1 -| s 2 |)/(2B + 1). Thus, assuming that s 2 is uniformly distributed in [-B, B], the expected number of operations is asymptotic to

1 2B + 1 B -B (2B + 1)π γ( s 2 ) d s 2 = B -B π 2B + 1 -| s 2 | d s 2 , which itself is asymptotic to 2(2 - √ 2) √ 2πB ≈ 2.07 √ 2B.
Though, one should note that even if C is randomly chosen uniformly among all curves, the number of points of the Jacobian is not uniformly distributed, as the concentration is slightly higher in the middle of the Hasse-Weil interval; therefore s 2 is not uniformly distributed.

Pseudo-random walk and distinguished points

We now address the main issues raised in the above algorithm: the generation of random elements in W and T and their storage. The key to the answer is to replace randomness by a deterministic pseudo-random walk.

Let r ≥ 0 and ℓ be parameters to be fixed later. The pseudo-random walk is initialized as follows: for all k in [1, r],

an offset O k ∈ J(C) is precomputed as O k = α k m • D, where α k is a random integer in [0, 2ℓ].
We also need a hash function H that maps elements of J(C) to [1, r]; this hash function should have good statistical properties, but no cryptographically strong property, like onewayness, is required. Typically, H is obtained by taking a few bits in the internal representation of the elements and the integer obtained this way is taken modulo r. Then, starting with an element P in W (resp. in T ) for which we know the corresponding σ 2 , we define another point

Q by Q = P + O H(P ) .
Assuming that ℓ is not too large compared to B, with high probability the point Q is still an element of W (resp. of T ). Furthermore, the value of σ 2 corresponding to Q is obtained by adding α H(P ) to the value of σ 2 for P .

Iterating this process yields a chain of pseudo-randomly chosen elements in W (resp. in T ). However, the chain should not be too long, to keep the probability to go out of the domain moderate. Thus the average length ℓ of an offset must be adapted according to the average number of steps we expect to do in one chain (see below). With this device it is then possible to produce each new pseudo-random element in W or T for one operation in J(C).

We now deal with storage requirements. To this effect, we introduce the concept of distinguished points, originally appeared in [START_REF] Delescaille | How easy is collision search? Application to DES[END_REF]. We say that an element of J(C) is a distinguished point if its image by a second hash function is 0. Again, we do not ask much of this hash function, except that its behavior looks independent of any arithmetic property of the element, and that the probability of being distinguished can be effectively estimated and tuned to a prescribed value. As usual, looking at some bits in the internal representation of the elements is a good way of doing. We denote by p D the probability for an element of being distinguished.

The algorithm then proceeds as follows: starting with a random point alternatively in W and T , we produce pseudo-random elements using the pseudorandom walk, until a distinguished point is hit. Then this point is stored and another chain is started. The length of the chains is 1/p D on average, and the parameter ℓ should be tuned accordingly. If all the parameters are well chosen, then after (say) 1000 chains on average, we have produced enough points to expect a collision. If this occurs at a point that is not distinguished, then the two chains continue on the same track because the pseudo-random walk is deterministic; the two chains will end at the same distinguished point, thus allowing the detection and the solution of the problem.

Many practical experiments were made to test this strategy against the idealized algorithm described above. They are very satisfactory and in [START_REF] Teske | Speeding up Pollard's rho method for computing discrete logarithms[END_REF] it was suggested that taking r = 20 is enough to simulate a random walk in this context.

As a conclusion, in the case when m ≥ 8 √ q, the number of points can be computed using a parallel low-memory algorithm that requires on average O( q/m) operations in J(C). This is worse than the O(q 3/4 /m) complexity announced in [START_REF] Matsuo | An improved baby step giant step algorithm for point counting of hyperelliptic curves over finite fields[END_REF], but their complexity analysis implicitly excludes that case.

The general case

In the case when m < 8 √ q, there are several choices for s 1 and still many more for s 2 . We could loop over all the possible values of s 1 and look for a corresponding s 2 using the algorithm of the previous section, but this approach has worse complexity than the MCT algorithm. The workaround is to take into account the bidimensional nature of the problem.

Recall that we have s 1 = s 1 + m s 1 and s 2 = s 2 + m s 2 , where s 1 and s 2 are known integers in [0, m -1], the goal being to find s 1 and s 2 . From the bounds (1) on s 1 and s 2 , we deduce similar bounds for s 1 and s 2 . In order to make the description more generic, we write these bounds in the following form:

B 1,min ≤ s 1 ≤ B 1,max , B 2,min ≤ s 2 ≤ B 2,max .

Random search in two intersecting rectangles

Let D be a random element of J(C); as previously, we assume that the order of D is large enough compared to the group order. Since χ(1) • D = 0, we have

q 2 + 1 -s 1 (q + 1) + s 2 • D + (-s 1 (q + 1) + s 2 ) • m • D = 0.
Let K be q 2 + 1 -s 1 (q + 1) + s 2 , which is a known integer; then we have to find s 1 and s 2 such that K • D + (-s 1 (q + 1) + s 2 ) m • D = 0. Let R be a rectangle containing the possible values for the pair ( s 1 , s 2 ):

R = {(σ 1 , σ 2 ); σ 1 ∈ [B 1,min , B 1,max ], σ 2 ∈ [B 2,min , B 2,max ]} .
Since B 1,min (resp. B 2,min ) is not necessarily the opposite of B 1,max (resp. B 2,max ), we have to normalize the situation. We define a value K ′ to be used in place of K, so as to center our search:

K ′ = K + m - B 1,min + B 1,max 2 (q + 1) + B 2,min + B 2,max 2
.
We can now define two sets of points:

W = {K ′ • D + (-σ 1 (q + 1) + σ 2 ) • m • D; (σ 1 , σ 2 ) ∈ R} , T = {(-σ 1 (q + 1) + σ 2 ) • m • D; (σ 1 , σ 2 ) ∈ R} .
We assume that these two sets have cardinality exactly #R. This may not hold in general, but is true with the further assumptions that D is of large order and that m is larger than 8. By construction, the sets W and T have a non-trivial intersection. Let D W be in W and D T in T . We write (σ 1W , σ 2W ) the values corresponding to D W and (σ 1T , σ 2T ) the values corresponding to D T . Then, assuming again that the order of D is large enough, D W = D T if and only if

σ 1W -σ 1T = s 1 -⌊(B 1,min + B 1,max )/2⌋ , σ 2W -σ 2T = s 2 -⌊(B 2,min + B 2,max )/2⌋ . ( 2
)
Hence it is easily checked that

N = #(W ∩ T ) ∈ #R 4 , #R .
When B 1,min = -B 1,max and B 2,min = -B 2,max , we get the picture of Figure 2.

A first version of our algorithm now proceeds as follows: random elements of W and T are constructed by picking random elements in R. These elements are stored in a data structure in which it is possible to detect quickly collisions between an element of W and an element of T . Together with these elements, we also store the corresponding pair (σ 1 , σ 2 ). On average, due to the birthday paradox, a collision occurs after having constructed O( √ N ) elements of W and T ; taking the difference of the pairs σ T and σ W then gives the result by Equations [START_REF] Adleman | Counting points on curves and abelian varieties over finite fields[END_REF].

Since the bounds on s 1 and s

2 yield | s 1 | = O( √ q/m) and | s 2 | = O(q/m), we have N = O(q 3/2 /m 2 )
, and therefore the expected number of points to construct is in O(q 3/4 /m). Just as in the unidimensional case, we now give an estimate of the constant hidden in the O( ), using simplifying assumptions. For such estimates, we can assume that the problem is centered, and therefore we put B 1 = B 1,max = -B 1,min and B 2 = B 2,max = -B 2,min ; hence #R = (2B 1 + 1)(2B 2 + 1). We denote by γ ∈ [ 1 4 , 1] the ratio #(W ∩ T )/#R: this parameter is easily computed as:

γ( s 1 , s 2 ) = (2B 1 + 1 -| s 1 |)(2B 2 + 1 -| s 2 |) #R .
From the birthday paradox, we see that the expected number of elements that have to be created before a collision between an element of W and an element of T occurs is asymptotic to π #R/γ. We now assume that s 1 and s 2 are uniformly distributed. Then the average number of points to construct grows like

1 #R B1 e s1=-B1 B2 e s2=-B2 π #R γ( s 1 , s 2 ) d s 1 d s 2 .
This integral is easily computed and shows that the expected number of points to construct is asymptotic to 8

√ π(3 -2 √ 2) √ #R ≈ 2.43 √ #R. Now, bounds (1)
on s 1 and s 2 yield approximate bounds for s 1 and s 2 :

B 1,max -B 1,min = 8 √ q/m, B 2,max -B 2,min = 8q/m, hence #R = 64 q 3/2 /m 2 .
The approximate value of 2.43 then yields a running time of about 19.5 q 3/4 /m operations in J(C). Hence, we see that a constant factor of about 5 is lost compared to the original (memory consuming) MCT algorithm. This difference is partly due to the fact that in the MCT algorithm, the BSGS approach allows to search only in the area described in Figure 1, whereas ours does not take this specificity into account.

Note that this analysis is idealized in several places: first, the pseudo-random walk that will be used below is not purely random and can cause some discrepancies. Next, the assumption that s 1 and s 2 are uniformly distributed in a rectangle is actually wrong. However, we expect that this should still give a good estimate.

A bidimensional pseudo-random walk

The main questions are now the generation of random elements in W and T , and how to avoid storing them. These are the same issues as in Section 3 and similar techniques will be applied to solve them. The definition of distinguished points still makes sense in this new context, and will be used with no modification. However, the previous pseudo-random walk needs to be converted into a bidimensional one.

Let r, ℓ 1 and ℓ 2 be parameters to be fixed later: r controls the number of offsets that we are going to precompute, and ℓ 1 and ℓ 2 are the average lengths of the horizontal and vertical offsets. For each k and k ′ in [1, r], we select a random non-negative integer α k,k ′ uniformly in [0, 2ℓ 1 ] and a random non-negative integer β k,k ′ uniformly in [0, 2ℓ 2 ]. Then for each k and k ′ in [1, r] and b in {0, 1}, we compute and store the offsets

O k,k ′ ,b = (-1) b α k,k ′ (q + 1)mD + β k,k ′ mD ∈ J(C),
where D is the base point whose order is to be computed.

Starting with a point P in W (resp. in T ) for which we know the corresponding pair (σ 1 , σ 2 ), we define another point as follows. We compute k, k ′ , and b as pseudo-random deterministic functions of P , by using some hash functions. Then we define Q = P + O k,k ′ ,b . If ℓ 1 and ℓ 2 are small enough, with high probability Q is still in W (resp. in T ) and the corresponding pair is given by

(σ 1 -(-1) b α k,k ′ , σ 2 + β k,k ′ ).
Iterating this process allows us to produce chains of pseudo-randomly chosen elements in W (resp. in T ). As before, the chains can not be too long, otherwise they go out of W (resp. of T ). The cost of producing one element is one group operation.

Note that we have only used positive values for the offset in the second direction: this is intended to reduce the chance of creating cycles. However, experiments with alternative strategies turned out to yield similar results.

Setting the parameters

Let λ be such that the expected number of points to construct is λ √ #R and let C be the number of chains we expect to construct. Note that C is fixed by the user; it should be large enough so that averaging considerations make sense (say C ≥ 1000), and small enough so that the cost of initializing a chain is negligible compared to the cost of the steps that are done in that chain. Also the number of chains is essentially the number of distinguished points that have to be stored and therefore should be small enough. Knowing an approximation of λ and having fixed C, the probability of being distinguished follows from

p D = C λ √ #R .
Now we fix ℓ 1 and ℓ 2 ; to control the probability of going out of W (or T ), we first evaluate the average length of a chain. There are about 1/p D steps and each step goes on average ℓ 1 in the first direction and ℓ 2 in the second direction. In the second direction, all the offsets are positive and therefore the length of a chain is ℓ 2 /p D on average. We want this to be small enough compared to the size of R in that direction, say one tenth of B 2,max -B 2,min :

ℓ 2 = (B 2,max -B 2,min ) p D 10 .
In the first direction, the situation is different since the offsets can be positive or negative, but still with an average absolute value of ℓ 1 . The central limit theorem applied to the 1-dimensional random walks gives that the average distance to the origin after 1/p D steps is about 2 2/3π ℓ 1 / √ p D : the factor 2/ √ 3 corresponds to the standard deviation of the lengths of the offsets. For convenience, we approximate 2 2/3π by 9/10. Again, imposing that this value is one tenth of the size of R in that direction yields:

ℓ 1 = (B 1,max -B 1,min ) √ p D 9 .
It may happen that ℓ 1 or ℓ 2 is very small, and even smaller than 1. This is especially the case when the bounds on s 1 and s 2 have a different order of magnitude. Unfortunately this is the case in the genus 2 point counting case, where s 2 is on average about √ q times larger than s 1 . A solution would be to enlarge C, but this is not satisfactory since it implies more storage. A better choice is to modify the random walk as follows. Assume that ℓ 1 is small (ℓ 1 and ℓ 2 can not be simultaneously small). Then with probability p, we add either the same offset O k,k ′ ,b as before or a modified offset O k,k ′ = β k,k ′ mD which does not include a progression in the first direction. The probability p is fixed so that the apparent mean value of ℓ 1 is the one we wanted. Note that the decision of adding O k,k ′ ,b or O k,k ′ is a deterministic choice that depends on D.

Reducing the search space

In order to have a better understanding of the distribution of (s 1 , s 2 ), we ran some statistics. Note that similar statistics, supported by theoretical and heuristic considerations were done in [START_REF] Stein | Explicit bounds and heuristics on class numbers in hyperelliptic function fields[END_REF], but with the purpose of finding the mean value of the class number.

Here, for p = 10 6 + 3, we randomly selected 10, 000 monic squarefree polynomials of degree 5 over F p and computed the (s 1 , s 2 ) values for the corresponding curves. As expected, the pairs (s 1 , s 2 ) tend to be not too close to the borders of the domain. In Figure 3, we represented the domain where (s 1 , s 2 ) are allowed to stay according to bounds (1), and inside the domain, the darkness of a point means that the density of pairs (s 1 , s 2 ) is high.

On the picture, we see that there are very few pairs for which s 2 is large, because the "wings" are very thin. In fact, these points correspond to curves which are close to maximal curves, and it is no surprise that they are rare. More precisely, in our tests, the proportion of curves for which s 2 is larger than 3q is about 2.2% and the proportion of curves for which s 2 is larger than 4q is about 0.23% (remember that in theory, s 2 can be as large as 6q). In view of these results, in order to improve our constant 19.5, and remarking that there is no point in spending too much time for large values of s 2 , we will restrict the rectangle R to the following bounds:

R = {(σ 1 , σ 2 ); σ 1 ∈ [-2.5 √ q/m, 2.5 √ q/m], σ 2 ∈ [-2q/m, 3q/m]} .
In our statistics, for more than 97% of the curves, s 2 < 3q/m, and therefore the overlapping factor of W and T is at least 1/4; and for about 99.7% of the curves, s 2 < 4q/m and therefore the overlapping factor is at least 3 25 = 0.12. The overlapping decreases as the point (s 1 , s 2 ) gets closer to the end of the "wings" of the arrow on the picture and if s 2 > 5.5 q, the sets W and T do not overlap.

With this strategy, the area of R is reduced to 25 q 3/2 /m 2 so that the expected runtime is about 12 q 3/4 /m group operations. Therefore we lose "only" a factor of 3 compared to the original MCT algorithm; this strategy is used in the experiments presented below.

We expect that this is very unlikely to get curves with s 1 , s 2 outside of the above rectangle R by random constructions. In case the algorithm does not find an answer after say 10 times the above expected time it could pay to start a classical MCT algorithm to search deterministically in the ends of the wings if we have enough memory for this much smaller subproblem. Otherwise, another method is to start a chain corresponding to these wings in the area outside the rectangle with a small probability, so that we do not perturb the average runtime but we can guarantee that the program finishes.

Practical experiments

We did two kinds of experiments: first with a high level implementation using the Magma computer algebra system [START_REF] Bosma | Handbook of Magma functions[END_REF], we ran some simulations to check the validity of our approach in various situations. Then, we wrote a parallel C++ implementation of our algorithm using the NTL [START_REF] Shoup | NTL: A library for doing number theory[END_REF] and the MPICH [START_REF]MPICH: A portable implementation of MPI[END_REF] libraries, in order to run tests with real sized curves. In our experiments, we used the reduced search space given in Section 4.4, with #R ≈ 25 q 3/2 /m 2 .

Simulations

In order to test the validity of our heuristics on sufficiently many examples of reasonable size, we used the following simulated algorithm. Let p be an integer. We pick two integers s 1 and s 2 at random uniformly in the area of Figure 1. Then we form the integer N = p 2 + 1 -s 1 (p + 1) + s 2 and work in the group Z/N Z instead of J(C). An integer m is chosen, and using only p and s 1 , s 2 modulo m, our goal is to recover N . Thus, we can construct appropriate examples with different values for the parameters at almost no cost.

We chose several pairs (p, m) yielding always #R ≈ 5 × 10 10 , so that p D can be taken to be 2 -8 , and with ℓ 1 varying from 0.01 to 21. For each such pair (p, m), 100 random pairs (s 1 , s 2 ) where tested and the average number of steps is measured and compared to p Our conclusion is that even if the rectangle R is very thin (i.e. ℓ 1 is tiny), then the measured running time is quite close to the predictions. In the extreme, the last case ℓ = 0.01 corresponds to a case where m ≈ 4 √ q, therefore B 1,max -B 1,min is just 2. In that case, the heuristics we used hardly make sense, but this is a good surprise that the average running time is still within a factor 2 of the heuristic analysis.

Several runs on the same curve

We ran our software implementation many times on a given curve, in order to check that the average measured running time is close to the heuristic estimate, for a medium sized problem.

Let p = 5 × 10 24 + 8503491 and f a random monic squarefree polynomial of degree 5 over F p . Using the Schoof-like method described in [START_REF] Gaudry | Construction of secure random curves of genus 2 over prime fields[END_REF] we have deduced the values of s 1 and s 2 modulo m = 44696171520. The curve is such that s 1 / √ p ≈ -0.84 and s 2 /p ≈ 0.38, therefore it is not close to any border of the bounds.

With the MCT algorithm, the computation is feasible and requires about 1.5 GB of memory. We ran our software 100 times on that input, and for each run we measured the number of operations in the Jacobian. The average value is 11.27 p 3/4 /m, which is in accordance with our estimates. The minimal value is 2.69 p 3/4 /m and the maximal value is 31.55 p 3/4 /m. In Figure 4 we give the histogram for the number of runs whose running time is in a given range.

A larger example

In order to test the scalability of our method, we ran a larger example. Let p = 5 × 10 24 + 8503491 as above, and f a random monic squarefree polynomial of degree 5 over F p . We now suppose that s 1 and s 2 are known modulo only m = 1655413760 (note that the modulus 44696171520 mentioned above equals 27 × 1655413760).

With this more restricted information, we cannot conclude using the original MCT algorithm: the computation would require to store about 5 × 10 9 points together with their indices. Even using hash tables, it seems difficult to use less than 8 bytes for each entry, which means at least 40 GB of memory.

We ran our algorithm on that input. The probability of being distinguished was set to p D = 2 -24 . The program was run in parallel on a cluster of 24 Pentium IV at 3 GHz. After 4 hours of computation, 544 distinguished points were computed and a useful collision occurred. About 9 × 10 9 steps were performed: this is a "lucky" run, since this is about 4.3 p 3/4 /m. The memory requirement was about 1.5 MB on each node, mostly for the executable code, not the data. The amount of communication is reduced to a few KB between the nodes and a "master" node.

For that curve, s 1 / √ p ≈ -1.16 and s 2 /p ≈ 1.25, so the curve was not exceptional, in the sense that it was not too close to the border of the area of Figure 1.

Curves of higher genus

The MCT algorithm can be extended to higher genus curves for which the characteristic polynomial of the Frobenius endomorphism is known modulo some integer m. In [START_REF] Izadi | Counting points on an abelian variety over a finite field[END_REF], different variants of this extension are studied and compared. Our extension of the original genus 2 MCT algorithm also applies to all these variants with a few modifications. Indeed, the number k of unknown coefficients in the characteristic polynomial can be larger than 2 and then a k-dimensional random walk should be designed in order to search for a collision in the intersection of two k-dimensional boxes. It is a tedious task to fill in the details, in particular the expected constants hidden in the O( ), but for fixed genus, we certainly lose only a constant factor compared to the memory-costly algorithms.

Another range of application of our method is the BSGS algorithm developed in [START_REF] Bauer | Point counting on Picard curves in large characteristic[END_REF] for counting points of the Jacobians of Picard curves. Without giving many details, let us mention that their algorithm ends by the search for a collision in two arithmetic progressions; therefore our algorithm applies almost directly to that case and should dramatically reduce the space complexity, which is the main drawback of their method.

Conclusion

We have presented a low-memory, parallelizable analogue of the algorithm by Matsuo, Chao and Tsujii for genus 2 hyperelliptic curves; the method works in other cases, such as the BSGS algorithm of [START_REF] Bauer | Point counting on Picard curves in large characteristic[END_REF] for Picard curves.

The main tool we used is a bidimensional pseudo-random walk. As usual with this kind of algorithms, it is impossible to make a rigorous analysis and heuristics and computer experiments are necessary to validate the approach. Our numerical data are positive in that sense, therefore our algorithm can be used when the memory constraint becomes problematic.

Fig. 2 .

 2 Fig. 2. Intersection of W and T

Fig. 3 .

 3 Fig. 3. Statistics on (s1, s2)

Fig. 4 .

 4 Fig. 4. Histogram showing the number of runs whose running time is in the given range (divided by p 3/4 /m).

Table 1 .

 1 3/4 /m. The results are reported in Table1. Simulations with cyclic groups Z/N Z. Each line corresponds to 100 runs.

	p	m	ℓ1	ℓ2	Avg ratio nb jumps / (p 3/4 /m)
	14 118 949 953674 1.2 × 10 17 144675925 5.7 × 10 7 9.2 × 10 8 1.5 × 10 10 1.5 × 10 14 1.9 × 10 18 1157407407 0.05 3, 242, 542 21 8, 040 10 15, 250 5 30, 350 0.5 304, 000 0.1 1, 621, 271 2.5 × 10 21 250000000000 0.01 19, 455, 252	12.7 12.7 12.5 12.6 14.3 14.5 22.2

Acknowledgements

We thank Gérard Guillerm for giving us access to his PC cluster, Gilles Schaeffer for sharing his knowledge of the central limit theorem, and Edlyn Teske for many helpful comments on a first version of this paper. Part of the computations were done at LIX, on a cluster paid by the ACI Cryptologie funding.


Alin Bostan
email: alin.bostan@stix.polytechnique.fr

Pierrick Gaudry
email: gaudry@lix.polytechnique.fr

Éric Schost
email: eric.schost@polytechnique.fr

Linear recurrences with polynomial coefficients and computation of the Cartier-Manin operator on hyperelliptic curves

We improve an algorithm originally due to Chudnovsky and Chudnovsky for computing one selected term in a linear recurrent sequence with polynomial coefficients. Using baby-steps / giant-steps techniques, the nth term in such a sequence can be computed in time proportional to √ n, instead of n for a naive approach. As an intermediate result, we give a fast algorithm for computing the values taken by an univariate polynomial P on an arithmetic progression, taking as input the values of P on a translate on this progression. We apply these results to the computation of the Cartier-Manin operator of a hyperelliptic curve. If the base field has characteristic p, this enables us to reduce the complexity of this computation by a factor of order √ p. We treat a practical example, where the base field is an extension of degree 3 of the prime field with p = 2 32 -5 elements.

Introduction

In this paper, we investigate some complexity questions related to linear recurrent sequences. Specifically, we concentrate on recurrences with polynomial coefficients; our main focus is on the complexity of computing one selected term in such a recurrence.

A well-known particular case is that of recurrences with constant coefficients, where the nth term can be computed with a complexity that is logarithmic in n, using binary powering techniques.

In the general case, there is a significant gap, as for the time being no algorithm with complexity polynomial in log(n) is known. Yet, in [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF], Chudnovsky and Chudnovsky proposed an algorithm that allows to compute one selected term in such a sequence without computing all intermediate ones. This algorithm appears as a generalization of those of Pollard [START_REF] Pollard | Theorems on factorization and primality testing[END_REF] and Strassen [START_REF] Strassen | Einige Resultate über Berechnungskomplexität[END_REF] for integer factorization; using baby-steps / giant-steps techniques, it requires a number of operations which is roughly linear in √ n to compute the nth term in the sequence.

Our main contribution is an improvement of the algorithm of [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF]; for simplicity, we only give the details in the case when all coefficients are polynomials of degree 1, as the study in the general case would follow in the same manner. The complexity of our algorithm is still (roughly) linear in √ n; Chudnovsky and Chudnovsky actually suggested that this bound might be essentially optimal. We improve the time and space complexities by factors that are logarithmic in n; in practice, this is far from negligible, since in the application detailed below, n has order 2 32 . A precise comparison with Chudnovsky and Chudnovsky's algorithm is made in Section 3.

Along the way, we also consider a question of basic polynomial arithmetic: given the values taken by a univariate polynomial P on a set of points, how fast can we compute the values taken by P on a translate of this set of points? An obvious solution is to make use of fast interpolation and evaluation techniques, but we show that one can do better when the evaluation points form an arithmetic sequence.

Computing the Cartier-Manin operator. Our initial motivation is an application to point-counting procedures in hyperelliptic curve cryptography, related to the computation of the Cartier-Manin operator of curves over finite fields. We now present these matters in more detail.

The Cartier-Manin operator of a curve defined over a finite field, together with the Hasse-Witt matrix, are useful tools to study the arithmetic properties of the Jacobian of that curve. Indeed, the supersingularity, and more generally the p-rank, can be read from the invariants of the Hasse-Witt matrix. In the case of hyperelliptic curves, this matrix was used in [START_REF] Gaudry | Counting points on hyperelliptic curves over finite fields[END_REF][START_REF] Matsuo | An improved baby step giant step algorithm for point counting of hyperelliptic curves over finite fields[END_REF] as part of a point-counting procedure for cryptographic-oriented applications.

Indeed, thanks to a result of Manin, computing the Cartier-Manin operator gives the coefficients of the Zeta function modulo p; this partial information can then be completed by some other algorithms. However, in [START_REF] Gaudry | Counting points on hyperelliptic curves over finite fields[END_REF] and [START_REF] Matsuo | An improved baby step giant step algorithm for point counting of hyperelliptic curves over finite fields[END_REF], the method used to compute the Hasse-Witt matrix has a complexity which is essentially linear in p.

It turns out that one can do better. The entries of the Hasse-Witt matrix of a hyperelliptic curve y 2 = f (x) defined over a finite field of characteristic p are coefficients of the polynomial h = f (p-1)/2 , so they satisfy a linear recurrence with rational function coefficients. Using our results on linear recurrences, this remark yields an algorithm to compute the Hasse-Witt matrix whose complexity now grows like √ p, up to logarithmic factors, instead of p.

We demonstrate the interest of these techniques by a point-counting example, for a curve of genus 2 defined over a finite field whose characteristic just fits in one 32-bit machine word; this kind of fields have an interest for efficiency reasons [START_REF] Bailey | Optimal extension fields for fast arithmetic in public-key algorithms[END_REF].

Note finally that other point-counting algorithms, such as the p-adic methods used in Kedlaya's algorithm [START_REF] Kedlaya | Countimg points on hyperelliptic curves using Monsky-Washnitzer[END_REF], also provide efficient point-counting procedures in small characteristic, but their complexity remains at least linear in p [START_REF] Gaudry | Counting points in medium characteristic using Kedlaya's algorithm[END_REF]. On the other hand, Kedlaya's algorithm outputs the whole Zeta function and should be preferred if available. Therefore, the range of application of our algorithm is when the characteristic is too large for Kedlaya's algorithm to be run.

Organization of the paper. We start in Section 2 with our algorithm for shifting a polynomial given by its values on some evaluation points. This building block is used in Section 3 to describe our improvement on Chudnovsky and Chudnovsky's algorithm. In Section 4 we apply these results to the computation of the Cartier-Manin operator of a hyperelliptic curve. We conclude in Section 5 with a numerical example.

Notation.

In what follows, we give complexity estimates in terms of number of base ring operations (additions, subtractions, multiplications and inversions of unit elements) and of storage requirements; this last quantity is measured in terms of number of elements in the ring. We pay particular attention to polynomial and matrix multiplications and use the following notation.

-Let R be a commutative ring; we suppose that R is unitary, its unit element being denoted by 1 R , or simply 1. Let ϕ be the map We also assume that the product of two degree d polynomials can be computed in space O(d); this is the case for all classical algorithms, such as naive, Karatsuba and Schönhage-Strassen multiplications.

N → R sending n to n • 1 R = 1 R + • • • + 1 R (n
-We let ω be a real number such that for every commutative ring R, all n × n matrices over R can be multiplied within O(n ω ) operations in R. The classical multiplication algorithm gives ω = 3. Using Strassen's algorithm [START_REF] Strassen | Gaussian elimination is not optimal[END_REF], we can take ω = log 2 (7) 2.81. We assume that the product of two n × n matrices can be computed in space O(n 2 ), which is the case for classical as well as Strassen's multiplications.

In the sequel, we need the following classical result on polynomial arithmetic over R. The earliest references we are aware of are [START_REF] Moenck | Fast modular transforms via division[END_REF][START_REF] Borodin | Fast modular transforms[END_REF], see [START_REF] Zur Gathen | Modern computer algebra[END_REF] for a detailed account. We also refer to [START_REF] Bostan | Tellegen's principle into practice[END_REF] for a solution that is in the same complexity class, but where the constant hidden in the O( ) notation is actually smaller than that in [START_REF] Zur Gathen | Modern computer algebra[END_REF].

Multipoint evaluation. If P is a polynomial of degree d in R[X] and r 0 , . . . , r d are points in R, then the values P (r 0 ), . . . , P (r d ) can be computed using O(M(d) log(d)) operations in R and O(d log(d)) space.

Shifting evaluation values

In this section, we address a particular case of the question of shifting evaluation values of a polynomial. The question reads as follows: Let P be a polynomial of degree

d in R[X],
where R is a commutative unitary ring. Let a and r 0 , . . . , r d be in R. Given P (r 0 ), . . . , P (r d ), how fast can we compute P (r 0 + a), . . . , P (r d + a)?

A reasonable condition for this question to make sense is that all differences r i -r j , i = j, are units in R; otherwise, uniqueness of the answer might be lost. Under this assumption, using fast interpolation and fast multipoint evaluation, the problem can be answered within O(M(d) log(d)) operations in R. We now show that the cost reduces to M(2d)+O(d) operations in R, in the particular case when r 0 , . . . , r d are in arithmetic progression, so we gain a logarithmic factor.

Our solution reduces to the multiplication of two suitable polynomials of degree at most 2d; O(d) additional operations come from additional pre-and post-processing operations. As mentioned in Section 1, all operations made below on integer values actually take place in R.

The algorithm underlying Proposition 1 is given in Figure 1; we use the notation coeff(Q, k) to denote the coefficient of degree k of a polynomial Q. We stress the fact that the polynomial P is not part of the input of our algorithm.

Input P (0), . . . , P (d) and a in R Output P (a), . . . , P (a + d)

-Compute δ(0, d) = d j=1 (-j), δ(i, d) = i i -d -1 δ(i -1, d) i = 1, . . . , d ∆(a, 0, d) = d j=0 (a -j), ∆(a, k, d) = a + k a + k -d -1 ∆(a, k -1, d) k = 1, . . . , d -Let P = d i=0 P (i) δ(i, d) X i , S = 2d i=0 1 a + i -d X i , Q = P S.
-Return the sequence ∆(a, 0, Proof. Our assumption on R enables to write the Lagrange interpolation formula:

d)•coeff(Q, d), . . . , ∆(a, d, d)•coeff(Q, 2d).
P = d i=0 P (i) d j=0,j =i (X -j) d j=0,j =i (i -j)
.

From now on, we denote by δ(i, d) the denominator d j=0,j =i (i -j) and by P i the ratio P (i)/δ(i, d).

First note that all δ(i, d), i = 0, . . . , d, can be computed in O(d) operations in R. Indeed, computing the first value δ(0, d) = d j=1 (-j) takes d multiplications. Then for i = 1, . . . ,d,δ(i,d) can be deduced from δ(i -1, d) for two ring operations using the formula

δ(i, d) = i i -d -1 δ(i -1, d),
so their inductive computation requires O(d) multiplications as well. Thus the sequence P i , i = 0, . . . , d, can be computed in admissible time and space O(d) from the input sequence P (i). Accordingly, we rewrite the above formula as

P = d i=0 P i d j=0,j =i (X -j).
For k in 0, . . . , d, let us evaluate P at a + k:

P (a + k) = d i=0 P i d j=0,j =i (a + k -j).
Using our assumption on a, we can complete each product by the missing factor a + k -i:

P (a + k) = d i=0 P i d j=0 (a + k -j) a + k -i =   d j=0 (a + k -j)   • d i=0 P i 1 a + k -i . (1)
Just as we introduced the sequence δ(i, d) above, we now introduce the sequence

∆(a, k, d) defined by ∆(a, k, d) = d j=0 (a + k -j).
In a parallel manner, we deduce that all ∆(a, k, d), k = 0, . . . , d can be computed in time and space O(d), using the formulas:

∆(a, 0, d) = d j=0 (a -j), ∆(a, k, d) = a + k a + k -d -1 ∆(a, k -1, d).
Let us denote Q k = P (a + k)/∆(a, k, d). We now show that knowing P i , i = 0, . . . , d, we can compute Q k , k = 0, . . . , d in M(2d) base ring operations and space O(d); this is enough to conclude, by the above reasoning.

Using the coefficients ∆(a, k, d), Equation ( 1) reads

Q k = d i=0 P i 1 a + k -i . (2)
Let P and S be the polynomials:

P = d i=0 P i X i , S = 2d i=0 1 a + i -d X i ;
then by Equation ( 2), for k = 0, . . . , d, Q k is the coefficient of degree k + d in the product P S. This concludes the proof.

We will conclude this section by an immediate corollary of this proposition; we first give a few comments.

-An alternative O(M(d)) algorithm which does not require any inversibility hypotheses can be designed in the special case when a = d + 1. The key fact is that for any degree d polynomial P , the sequence P (0), P (1), . . . is linearly recurrent, of characteristic polynomial Q(X) = (1 -X) d+1 . Thus, if the first terms P (0), . . . , P (d) are known, the next d+1 terms P (d+1), . . . , P (2d+1) can be recovered in O(M(d)) using the algorithm in [27, Theorem 3.1]. -The general case when the evaluation points form an arbitrary arithmetic progression reduces to the case treated in the above proposition. Indeed, suppose that r 0 , . . . , r d form an arithmetic progression of difference δ, that P (r 0 ), . . . , P (r d ) are known and that we want to compute the values P (r 0 + a), . . . , P (r d + a), where a ∈ R is divisible by δ. Introducing the polynomial Q(X) = P (δX + r 0 ), we are under the hypotheses of the above proposition, and it suffices to determine the shifted evaluation values of Q by a/δ. -The reader may note the similarity of our problem with the question of computing the Taylor expansion of a given polynomial P at a given point in R. The algorithm of [START_REF] Aho | Evaluating polynomials at fixed sets of points[END_REF] solves this question with a complexity of M(d) + O(d) operations in R and space O(d). The complexity results are thus quite similar; it turns out that analogous generating series techniques are used in that algorithm.

-In [START_REF] Hanrot | The middle product algorithm, I. Speeding up the division and square root of power series[END_REF], an operation called middle product is defined: Given a ring R, and A, B in R[X] of respective degrees d and 2d, write Proof. For any s ∈ N, let us denote by P s (X) the polynomial P (2 s X). We prove by induction that all values P s (0), . . . , P s (d) can be computed in time s M(2d) + O(sd) and space O(d), which is enough to conclude. The case s = 0 is obvious, as there is nothing to compute. Suppose then that P s (0), . . . , P s (d) can be computed in time s M(2d)+O(sd) and using O(d) temporary space allocation. Under our assumption on R, Proposition 1 shows that the values P s (d + 1), . . . , P s (2d + 1) can be computed in time M(2d) + O(d), using again O(d) temporary space allocation. The values P s (0), P s (2), . . . , P s (2d) coincide with P s+1 (0), P s+1 (1), . . . , P s+1 (d), so the corollary is proved.

AB = C 0 + C 1 X d+1 + C 2 X 2d+2 ,

Computing one selected term of a linear sequence

In this section, we recall and improve the complexity of an algorithm due to Chudnovsky and Chudnovsky [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF] for computing selected terms of linear recurrent sequences with polynomial coefficients. The results of the previous section are used as a basic subroutine for these questions.

As in the previous section, R is a commutative ring with unity. Let A be a n × n matrix of polynomials in R[X]. For simplicity, in what follows, we only treat the case of degree 1 polynomials, since this is what is needed in the sequel. Nevertheless, all results extend mutatis mutandis to arbitrary degree.

For r in R, we denote by A(r) the matrix over R obtained by specializing all coefficients of A at r. In particular, for k in N, A(k • 1 R ) is simply denoted by A(k), following the convention used up to now. Given a vector of initial conditions U 0 = [u 1 , . . . , u n ] t ∈ R n and given k in N, we consider the question of computing the kth term of the linear sequence defined by the relation

U i = A(i)U i-1 for i > 0, that is, the product U k = A(k)A(k -1) • • • A(1)U 0 .
For simplicity, we write

U k = k i=1 A(i) U 0 ,
performing all successive matrix products, i = 1, . . . , k, on the left side. We use this convention hereafter.

In the particular case when A is a matrix of constant polynomials, and taking only the dependence on k into account, the binary powering method gives a time complexity of order O(log(k)) base ring operations.

In the general case, the naive solution consists in evaluating all matrices A(i) and performing all products. With respect to k only, the complexity of this approach is of order O(k) base ring operations. In [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF], Chudnovsky and Chudnovsky propose an algorithm that reduces this cost to essentially O √ k . We first recall the main lines of this algorithm; we then present some improvements in both time and space complexities.

The algorithm of Chudnovsky and Chudnovsky. The original algorithm uses baby-step / giant-step techniques, so for simplicity we assume that k is a square in N. Let C be the n × n matrix over R[X] defined by

C = √ k i=1 A(X + i),
where A(X + i) denotes the matrix A with all polynomials evaluated at X + i.

By assumption on A, the entries of C have degree at most √ k. For r in R, we denote by C(r) the matrix C with all entries evaluated at r. Then the requested output U k can be obtained by the equation

U k =   √ k-1 j=0 C(j √ k)   U 0 . (3)
Here are the main steps of the algorithm underlying Equation (3), originally due to [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF].

Baby steps. The "baby steps" part of the algorithm consists in computing the polynomial matrix C. In [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF], this is done within Summing all the above costs gives an overall complexity bound of

O(n ω M( √ k
O n ω M √ k + n 2 M √ k log(k)
base ring operations for computing a selected term of a linear sequence. Due to the use of fast evaluation algorithms in degree √ k, the space complexity is

O(n 2 √ k + √ k log(k)
). In the particular case when A is the 1 × 1 matrix [X], the question reduces to the computation of k j=1 j in the ring R. For this specific problem, note that the ideas presented above were already used in [START_REF] Pollard | Theorems on factorization and primality testing[END_REF][START_REF] Strassen | Einige Resultate über Berechnungskomplexität[END_REF], for the purpose of factoring integers.

Avoiding multiplications of polynomial matrices.

In what follows, we show how to avoid the multiplication of polynomial matrices, and reduce the cost of the above algorithm to

O(n ω √ k + n 2 M( √ k) log(k)) base ring operations, storing only O(n 2 √
k) elements of R. Our improvements are obtained through a modification of the baby steps phase; the underlying idea is to work with the values taken by the polynomial matrices instead of their representation on the monomial basis. This idea is encapsulated in the following proposition.

Proposition 2 Let A be a n×n matrix with entries in R[X], of degree at most 1.

Let N ≥ 1 be an integer and let C be the n × n matrix over R[X] defined by

C = N i=1 A(X + i).
Then one can compute all scalar matrices C(0), C(1), . . . , C(N ) within O(n ω N ) operations in R and with a memory requirement of O(n 2 N ) elements in R.

Proof. We first compute the scalar matrices [A(1), A(2), . . . , A(2N )] . Since all entries of A are linear in X, the complexity of this preliminary step is O(n 2 N ), both in time and space.

Then, we construct the matrices (C j ) 0≤j≤N and (C j ) 0≤j≤N , which are defined as follows: we let C 0 and C 0 equal the identity matrix I n and we recursively define

C j = A(N + j)C j-1 for 1 ≤ j ≤ N, C j = C j-1 A(N -j + 1) for 1 ≤ j ≤ N.
Explicitly, for 0 ≤ j ≤ N , we have

C j = A(N + j) • • • A(N + 1)
and

C j = A(N ) • • • A(N -j + 1), thus C N -j = A(N ) • • • A(j + 1).
Computing all the scalar matrices (C j ) and (C j ) requires 2N matrix multiplications with entries in R; their cost is bounded by O(n ω N ) in time and by O(n 2 N ) in space. Lastly, the formula From this proposition, we deduce the following corollary, which shows how to compute the scalar matrices used in the giant steps.

C(j) = A(N + j) • • • A(N + 1)A(N ) • • • A(j + 1) = C j C N -j , 0 ≤ j ≤ N enables to recover C(0), C (1)
Corollary 2 Let A and C be polynomial matrices as in Proposition 2. If the elements 1, . . . , 2N + 1 are units in R, then for any integer s ≥ 1, the sequence

C(0), C(2 s ), . . . , C(2 s (N -1)) can be computed using O(n ω N +n 2 s M(N )) operations in R and O(n 2 N ) memory space.
Proof. This is an immediate consequence of Proposition 2 and Corollary 1.

The above corollary enables us to perform the "giant steps" phase of Chudnovsky and Chudnovsky's algorithm in the special case when N = 2 s ; this yields the 4 s th term in the recurrent sequence. Using this intermediate result, the following theorem shows how to compute the kth term, for arbitrary k, using the 4-adic expansion of k.

Theorem 1 Let A be a n × n matrix with linear entries in R[X] and let U 0 be in R n . Suppose that (U i ) is the sequence of elements in R n defined by the linear recurrence U i+1 = A(i + 1)U i , for all i ≥ 0.

Let k > 0 be an integer and suppose that 1, . . . , 2 √ k + 1 are units in R. Then the vector U k can be computed within

O n ω √ k + n 2 M √ k log(k) operations in R and using memory space O(n 2 √ k).
The proof of Theorem 1 is divided in two steps. We begin by proving the proposition in the particular case when k is a power of 4, then we treat the general case.

The case k is a power of 4. Let us suppose that N = 2 s and k = N 2 , so that k = 4 s . With this choice of k, Corollary 2 shows that the values C(0), C(N ), . . . , C((N -1)N ) can be computed within the required time and space complexities. Then we go on to the giant step phase described at the beginning of the section, and summarized in Equation (3). It consists in performing √ k successive matrix-vector products, which has a cost in both time and space of O(n 2 √ k).

The general case. We now consider the general case. Let k = s i=0 k i 4 i be the 4-adic expansion of k, with k i ∈ {0, 1, 2, 3} for all i. Given any t, we will denote by k t the integer t-1 i=0 4 i k i . Using this notation, we define a sequence (V t ) 0≤t≤s as follows: we let V 0 = U 0 and, for 0 ≤ t ≤ s we set

V t+1 = A( k t + 4 t k t ) • • • A( k t + 1)V t . (4)
It is easy to verify that V s+1 = U k . Therefore, it suffices to compute the sequence (V t ) within the desired complexities. Supposing that the term V t has been determined, we estimate the cost of computing the next term V t+1 . If k t is zero, we have nothing to do. Otherwise, we let V (0) t+1 = V t , and, for 1 ≤ j ≤ k t , we let A (j) (X) = A (X + k t + 4 t (j -1)). Then we define V (j)

t+1 by V (j) t+1 = A (j) (4 t ) • • • A (j) (1)V (j-1) t+1 , j = 1, . . . , k t .
By Equation ( 4), we have V kt t+1 = V t+1 . Thus, passing from V t to V t+1 amounts to computing k t selected terms of a linear recurrence of the special form treated in the previous paragraph. Using the complexity result therein and the fact that all k t are bounded by 3, the total cost of the general case is thus

O s t=0 n ω 2 t + n 2 t M(2 t ) = O n ω 2 s + n 2 s s t=0 M(2 t ) .
Using the fact that 2 s ≤ √ k ≤ 2 s+1 and the assumptions on the function M, we easily deduce that the whole complexity fits into the bound O n ω √ k+ n 2 M √ k log(k) , as claimed. Similar considerations also yield the bound concerning the memory requirements. This concludes the proof of Theorem 1.

Comments. The question of a lower time bound for computing U k is still open. The simpler question of reducing the cost to O n ω √ k + n 2 M √ k base ring operations, that is gaining a logarithmic factor, already raises challenging problems.

As the above paragraphs reveal, this improvement could be obtained by answering the following question: Let P be a polynomial of degree d in R[X]. Given r in R, how fast can we compute P (0), P (r), . . . , P (rd) from the data of P (0), P (1), . . . , P (d)? A complexity of order O(M(d)) would immediately give the improved bound mentioned above. We leave it as an open question.

The Cartier-Manin operator on hyperelliptic curves

We finally show how to apply the above results to the computation of the Cartier-Manin operator, and start by reviewing some known facts on this operator.

Let C be a hyperelliptic curve of genus g defined over the finite field F p d with p d elements, where p is the characteristic of F p d . We suppose that p > 2 and that the equation of C is of the form y 2 = f (x), where f ∈ F p d [X] is a monic squarefree polynomial of degree 2g +1. The generalization to hyperelliptic curves of the Hasse invariant for elliptic curves is the so-called Hasse-Witt matrix, which is defined as follows:

Definition 1 Let h k be the coefficient of degree k in the polynomial f (p-1)/2 . The Hasse-Witt matrix is the g × g matrix with coefficients in F p d given by H = (h ip-j ) 1≤i,j≤g .

This matrix was introduced in [START_REF] Hasse | Zyklische unverzweigte Erweiterungskörper vom primzahlgrade p über einem algebraischen Funktionenkörper der Charakteristik p[END_REF]; in a suitable basis, it represents the operator on differential forms that was introduced by Cartier in [START_REF] Cartier | Une nouvelle opération sur les formes différentielles[END_REF]. Manin then showed in [START_REF] Manin | The Hasse-Witt matrix of an algebraic curve[END_REF] that this matrix is strongly related to the action of the Frobenius endomorphism on the p-torsion part of the Jacobian of C. The article [START_REF] Yui | On the Jacobian varietes of hyperelliptic curves over fields of characteristic p > 2[END_REF] provides a complete survey about those facts; they can be summarized by the following theorem:

Theorem 2 (Manin) Let C be a hyperelliptic curve of genus g defined over F p d . Let H be the Hasse-Witt matrix of C and let

H π = HH (p) • • • H (p d-1 )
, where the notation H (q) means element-wise raising to the power q. Let κ(t) be the characteristic polynomial of the matrix H π and let χ(t) be the characteristic polynomial of the Frobenius endomorphism of the Jacobian of C. Then χ(t) ≡ (-1) g t g κ(t) mod p.

This result provides a quick method to compute the characteristic polynomial of the Frobenius endomorphism and hence the group order of the Jacobian of C modulo p, when p is not too large. Combined with a Schoof-like algorithm and / or a baby-step / giant-step algorithm, it can lead to a full point-counting algorithm, in particular for genus 2 curves, as was demonstrated in [START_REF] Gaudry | Counting points on hyperelliptic curves over finite fields[END_REF][START_REF] Matsuo | An improved baby step giant step algorithm for point counting of hyperelliptic curves over finite fields[END_REF].

The obvious solution consists in expanding the product f (p-1)/2 . Using balanced multiplications, and taking all products modulo X gp this can be done in O(M(gp)) base field operations, whence a time complexity within O(M(p)), if g is kept constant. In what follows, regarding the dependence in p only, we show how to obtain a complexity of O(M( √ p) log(p)) base field operations, using the results of the previous sections. We will make the assumption that the constant term of f is not zero. Note that if it is zero, the problem is actually simpler:

writing f = Xf 1 , the coefficient of degree ip-j in f (p-1)/2 is the coefficient of degree ip-j -(p-1)/2 in f (p-1)/2 1 .
Hence we can work with a polynomial of degree 2g instead of 2g + 1 and the required degrees are slightly less.

Furthermore, for technical reasons, we assume that g < p. This is not a true restriction since for g ≥ p, all the coefficients of f (p-1)/2 up to degree g(p -1) are needed to fill in the matrix H.

Introduction of a linear recurrent sequence. In [START_REF] Flajolet | The SIGSAM challenges: Symbolic asymptotics in practice[END_REF], Flajolet and Salvy already treat the question of computing a selected coefficient in a high power of some given polynomial, as an answer to a SIGSAM challenge. The key point of their approach is that h = f (p-1)/2 satisfies the following first-order linear differential equation

f h - p -1 2 f h = 0.
From this, we deduce that the coefficients of h satisfy a linear recurrence of order 2g + 1, with coefficients that are rational functions of degree 1.

Explicitly, let us denote by h k the coefficient of degree k of the polynomial h, and for convenience, set h k = 0 for k < 0. Similarly, the coefficient of degree k of f is denoted by f k . From the above differential equation, for all k in Z, we deduce that

2g+1 i=0 k + 1 - (p + 1)i 2 f i h k+1-i = 0. We set U k = [h k-2g , h k-2g+1 , . . . , h k ] t
, and let A(k) be the (2g + 1) × (2g + 1) companion matrix:

A(k) =         0 1 0 • • • 0 0 0 1 • • • 0 . . . . . . . . . . . . . . . 0 0 0 . . . 1 f2g+1((2g+1)(p-1)/2-(k-2g-1)) f0k • • • • • • • • • f1((p-1)/2-k+1) f0k        
.

The initial vector U 0 = [0, . . . , 0, f

(p-1)/2 0
] t can be computed using binary powering techniques in O(log(p)) base field operations; then for k ≥ 0, we have U k+1 = A(k + 1)U k . Thus, to answer our specific question, it suffices to note that the vector U ip-j gives the coefficients h ip-j for j = 1, . . . , g that form the ith row of the Hasse-Witt matrix of C.

Yet, Theorem 1 cannot be directly applied to this sequence, because A(k) has entries that are rational functions, not polynomials. Though the algorithm could be adapted to handle the case of rational functions, we rather use the very specific form of the matrix A(k), so only a small modification is necessary. Let us define a new sequence V k by the relation

V k = f k 0 k!U k .
Then, this sequence is linearly generated and we have V k+1 = B(k + 1)V k , where

B(k) = f 0 kA(k).
Therefore, the entries of the matrix B(k) are polynomials of degree at most 1. Note also that the denominators f k 0 k! satisfy the recurrence relation

f k+1 0 (k + 1)! = (f 0 (k + 1)) • (f k 0 k!).
Thus, we will compute separately, first V p-1 , V 2p-1 , . . . V gp-1 and then the denominators f p-1 0 (p -1)!, . . . , f gp-1 0 (gp -1)!. To this effect, we proceed iteratively. Let us for instance detail the computation of the sequence V p-1 , V 2p-1 , . . . V gp-1 . Knowing V 0 , we compute V p-1 using Theorem 1. Then we shift all entries of B by p, so another application of Theorem 1 yields V 2p-1 . Iterating g times, we obtain V p-1 , V 2p-1 , . . . V gp-1 as requested; the same techniques are used to compute f p-1 0 (p-1)!, . . . , f gp-1 0 (gp-1)!. Then the vectors U p-1 , U 2p-1 , . . . U gp-1 are deduced from

U k = 1 f k 0 k! V k .
Lifting to characteristic zero. A difficulty arises from the fact that the characteristic is too small compared to the degrees we are aiming to, so p! is zero in F p d . The workaround is to do computations in the unramified extension K of Q p of degree d, whose residue class field is F p d . The ring of integers of K will be denoted by O K ; any element of O K can be reduced modulo p to give an element of F p d . On the other hand, K has characteristic 0, so p is invertible in K.

We consider an arbitrary lift of f to O K [X]. The reformulation in terms of linear recurrent sequence made in the above paragraph can be performed over K; the coefficients of f (p-1)/2 are computed as elements of K and then projected back onto F p d . This is possible, as they all belong to O K .

Using the iteration described above, we separately compute the values in K of the vectors V ip-1 and the denominators f ip-1 0 (ip-1)!, for i = 1, . . . , g. To this effect, we apply g times the result given in Theorem 1; this requires to perform

O g ω+1 √ p + g 3 M( √ p) log(p) ,
operations in K and to store O(g 2 √ p) elements of K.

Computing at fixed precision. Of course, we do not want to compute in the field K at arbitrary precision: for our purposes, it suffices to truncate all computations modulo a suitable power of p. To evaluate the required precision of the computation, we need to check when the algorithm operates a division by p.

To compute the vectors V ip-1 and the denominators f ip-1 0 (ip -1)!, for i = 1, . . . , g, we use Theorem 1. This requires that all integers up to 2 √ p + 1 are invertible, which holds as soon as p ≥ 11.

Then, for all i = 1, . . . , g, to deduce U ip-1 from V ip-1 , we need to divide by f ip-1 0 (ip -1)!. The element f 0 is a unit in O K , so the only problem comes from the factorial term. With our assumption that g < p, we have i < p and then the p-adic valuation of (ip -1)! is exactly i -1. Therefore the worst case is i = g, for which we have to divide by p g-1 . Hence computing the vectors V ip-1 modulo p g is enough to know the vectors U ip-1 modulo p, and then to deduce the Hasse-Witt matrix. The matrix H by itself gives some information on the curve C, for instance H is invertible if and only if the Jacobian of C is ordinary [START_REF] Yui | On the Jacobian varietes of hyperelliptic curves over fields of characteristic p > 2[END_REF]Corollary 2.3]. However, as stated in Theorem 2, the matrix H π and in particular its characteristic polynomial χ(t) tell much more and are required if the final goal is point-counting. Thus, we finally concentrate on the cost of computing the characteristic polynomial of H π .

The matrix H π is the "norm" of H and as such can be computed with a binary powering algorithm. For simplicity, we assume that d is a power of 2, then denoting

H π,i = HH (p) • • • H p 2 i -1 . we have H π,i+1 = H π,i • (H π,i ) p 2 i
.

Hence the computation of H π,i+1 from H π,i costs one matrix multiplication and 2 i matrix conjugations. A matrix conjugation consists in raising all the entries to the power p, therefore it costs O(g The cost of the characteristic polynomial computation is bounded by the cost of a matrix multiplication [START_REF] Keller-Gehrig | Fast algorithms for the characteristic polynomial[END_REF] and is therefore negligible compared to the other costs.

If we are interested only in the complexity in p and d, i.e. if we assume that the genus is fixed, we get a time complexity for computing χ(t) mod p in

O ((M( √ p) + d) M(d log(p)) log(p)) .
Case of large genus. In case of large genus, the algorithm of Theorem 1 is asymptotically not the fastest. In this paragraph, we assume that the function M is essentially linear and we do not take into account the logarithmic factors; adding appropriate epsilons in the exponents would yield a rigorous analysis.

The cost in bit-operations of Theorem 3 is at least g 4 √ pd whereas the cost of the naive algorithm is linear in gpd. If g > p 1/6 , then g 4 √ p > gp, and therefore the naive algorithm is faster.

Point-counting numerical example

We have implemented our algorithm using Shoup's NTL C++ library [START_REF] Shoup | NTL: A library for doing number theory[END_REF]. NTL does not provide any arithmetic of local fields or rings, but allows to work in finite extensions of rings of the form Z/p g Z, as long as no division by p occur; the divisions by p are well isolated in the algorithm, so we could handle them separately. Furthermore, NTL multiplies polynomials defined over this kind of structure using an asymptotically fast FFT-based algorithm.

To illustrate that our method can be used as a tool in point-counting algorithms, we have computed the Zeta function of a (randomly chosen) genus 2 curve defined over F p 3 , with p = 2 32 -5. Such a Jacobian has therefore about 2 192 elements and should be suitable for cryptographic use if the group order has a large prime factor. Note that previous computations were limited to p of order 2 23 [START_REF] Matsuo | An improved baby step giant step algorithm for point counting of hyperelliptic curves over finite fields[END_REF].

The characteristic polynomial χ of the Frobenius endomorphism was computed modulo p in 3 hours and 41 minutes, using 1 GB of memory, on an AMD Athlon MP 2200+. Then we used the Schoof-like algorithms of [START_REF] Gaudry | Counting points on hyperelliptic curves over finite fields[END_REF] and [START_REF] Gaudry | Cardinality of a genus 2 hyperelliptic curve over GF(5 • 10 24 + 41). e-mail to the NMBRTHRY mailing list[END_REF] to compute χ modulo 128 × 9 × 5 × 7, and finally we used the modified baby-step / giant-step algorithm of [START_REF] Matsuo | An improved baby step giant step algorithm for point counting of hyperelliptic curves over finite fields[END_REF] to finish the computation. These other parts were implemented in Magma [START_REF] Bosma | The Magma algebra system. I. The user language[END_REF] and were performed in about 15 days of computation on an Alpha EV67 at 667 MHz. We stress that this computation was meant as an illustration of the possible use of our method, so little time was spent optimizing our code. In particular, the Schoof-like part and the final baby-step / giant-step computations are done using a generic code that is not optimized for extension fields.

Numerical data. The irreducible polynomial P (t) that was used to define F p 3 as F p [t]/(P (t)) is

t 3 + 1346614179t 2 + 3515519304t + 3426487663.
The curve C has equation y 2 = f (x) where f is given by f (x) = x 5 + (2697017539t 2 + 1482222818t + 3214703725)x 3 + (676673546t 2 + 3607548185t + 1833957986)x 2 + (1596634951t 2 + 3203023469t + 2440208439)x+ 2994361233t 2 + 3327339023t + 862341251.

Then the polynomial characteristic χ(T ) of the Frobenius endomorphism is given by T 4 -s 1 T 3 + s 2 T 2 -p 3 s 1 T + p 6 , where

s 1 = 332906835893875, s 2 = 142011235215638946167187570235.
The group order of the Jacobian is then This number has a large prime factor of size 2 158 , therefore that curve is cryptographically secure.

Measure of the complexity in p. To check the practical asymptotic behaviour of our algorithm, we ran our implementation on a genus 2 curve defined over F p 3 with p = 2 34 -41. We performed only the Cartier-Manin step, and not the full point-counting algorithm. As the characteristic is about 4 times larger than in the previous example, a complexity linear in √ p means a runtime multiplied by about 2. On the same computer, the runtime is 8 hours and 48 minutes.

Hence the ratio of the runtimes is about 2.39. The defect of linearity can be explained by taking into account the logarithmic factors. Assuming that M(n) is O(n log(n) log(log(n))), and neglecting the multi-logarithmic factors, the complexity announced in Theorem 3 is in O( √ p(log(p)) 3 ). With this estimate, the expected ratio between the runtimes becomes about 2.40, that is very close to the measure. This validates our analysis.

Conclusion

In this paper, we have presented an improvement of an algorithm by Chudnovsky and Chudnovsky to compute selected terms in a linear sequence with polynomial coefficients. This algorithm is then applied to the computation of the Cartier-Manin operator of hyperelliptic curves, thus leading to improvements in the point-counting problems that occur in cryptography. This strategy extends readily to curves of the form y r = f (x) with r > 2, for which the Hasse-Witt matrix has a similar form. For more general curves, Mike Zieve pointed to us the work of Stöhr and Voloch [START_REF] Stöhr | A formula for the Cartier operator on plane algebraic curves[END_REF] that gives formulas that still fit in our context in some cases.

Finally, Mike Zieve pointed out to us the work of Wan [START_REF] Wan | Computing zeta functions over finite fields[END_REF] that relates Niederreiter's polynomial factorization algorithm to the computation of the Cartier-Manin operator of some variety. The link with our work is not immediate, as that variety has dimension zero. Nevertheless, this remains intriguing, especially if we think of Pollard-Strassen's integer factoring algorithm as a particular case of Chudnovsky and Chudnovsky's algorithm.

  times); the map ϕ extends to a map Z → R. When the context is clear, we simply denote the ring element ϕ(n) by n. -We denote by M : N → N a function that represents the complexity of univariate polynomial multiplication, i.e. such that over any ring R, the product of two degree d polynomials can be computed within M(d) base ring operations. Using the algorithms of [25, 24, 8], M(d) can be taken in O(d log(d) log(log(d))). We suppose that the function M verifies the inequality M(d 1 ) + M(d 2 ) ≤ M(d 1 + d 2 ) for all positive integers d 1 and d 2 ; in particular, the inequality M(d) ≤ 1 2 M(2d) holds for all d ≥ 1. On the other hand, we make the (natural) hypothesis that M(cd) ∈ O(M(d)) for all c ≥ 1.

Fig. 1 .

 1 Fig. 1. Shifting evaluation values

  )) base ring operations, as products of polynomial matrices with entries of degree O( √ k) are required. Giant steps. In the second part the matrix C is evaluated on the arithmetic progression 0, √ k, 2 √ k, . . . , ( √ k -1) √ k and the value of U k is obtained using Equation (3). Using fast evaluation techniques, all evaluations are done within O(n 2 M( √ k) log(k)) base ring operations, while performing the √ k successive matrix-vector products in Equation (3) adds a negligible cost of O(n 2 √ k) operations in R.

  , . . . , C(N ) in time O(n ω N ) and space O(n 2 N ).

Theorem 3

 3 Overall complexity. Storing an element of O K /p g O K requires O(dg log(p)) bits, and multiplying two such elements can be done with O(M(dg log(p))) bitoperations. From the results of Section 3, we then deduce the following theorem on the complexity of computing the Hasse-Witt matrix. Let p a prime, d ≥ 1 and C a hyperelliptic curve defined over F p d by the equation y 2 = f (x), with f of degree 2g + 1. Then, assuming g < p, one can compute the Hasse-Witt matrix of C with a complexity of O g ω+1 √ p + g 3 M( √ p) log(p) M(dg log(p)) bit-operations and O dg 3 √ p log(p) storage.

6277101691541605395917785080771825883860189465813625993977 = 3 3 ×

 3 13 × 67 × 639679 × 417268068727536370810010172344236025455933953139.

  with all C i of degree at most d; then the middle product of A and B is the polynomial C 1 . This is precisely what is needed in the above algorithm. Up to considering the reciprocal polynomial of A, the middle product by A can be seen as the transpose of the map of multiplication by A. General program transformation techniques[START_REF] Bürgisser | Algebraic complexity theory[END_REF][START_REF] Hanrot | The middle product algorithm, I. Speeding up the division and square root of power series[END_REF] then show that it can be computed in time M(d) + O(d), but with a possible loss in space complexity. In[START_REF] Bostan | Tellegen's principle into practice[END_REF], it is shown how to keep the same space complexity, at the cost of a constant increase in time complexity. Managing both requirements remains an open question, already stated in[START_REF] Kaltofen | Challenges of symbolic computation: my favorite open problems[END_REF] Problem 6].

Corollary 1 Let R be a commutative ring with unity, and d ∈ N such that 1, . . . , 2d + 1 are units in R. Let P be a degree d polynomial in R[X] such that the sequence P (0), . . . , P (d) is known. For any s in N, the sequence P (0), P (2 s ), . . . , P (2 s d) can be computed in time s M(2d) + O(sd) ∈ O(s M(d)) and space O(d).

  p d . The general case where d is not a power of 2 is handled by adjusting the recursive step according to the binary expansion of d and yields the same complexity up to a constant factor.

2 log(p)) operations in F p d . The matrix we need to compute is H π = H π,log 2 (d) . Hence the cost of computing H π is O dg 2 log(p) + g ω log(d) operations in F

Acknowledgements

Part of our numerical experiments were done with the computer of the Medicis center [1]. We thank Takakazu Satoh and Mike Zieve for their interest in our work and their insightful comments and also Bruno Salvy for his numerous helpful suggestions.