Alin Bostan
email: alin.bostan@inria.fr

Xavier Caruso
email: xavier.caruso@normalesup.org

AND Gilles Christol
email: christol.gilles@gmail.com

Philippe Dumas
email: philippe.dumas@inria.fr

FAST COEFFICIENT COMPUTATION FOR ALGEBRAIC POWER SERIES IN POSITIVE CHARACTERISTIC

We revisit Christol's theorem on algebraic power series in positive characteristic and propose yet another proof for it. This new proof combines several ingredients and advantages of existing proofs, which make it very wellsuited for algorithmic purposes. We apply the construction used in the new proof to the design of a new efficient algorithm for computing the N th coefficient of a given algebraic power series over a perfect field of characteristic p. It has several nice features: it is more general, more natural and more efficient than previous algorithms. Not only the arithmetic complexity of the new algorithm is linear in log N and quasi-linear in p, but its dependency with respect to the degree of the input is much smaller than in the previously best algorithm. Moreover, when the ground field is finite, the new approach yields an even faster algorithm, whose bit complexity is linear in log N and quasi-linear in √ p.

Introduction

Given a perfect field k of characteristic p > 0, we address the following question: how quickly can one compute the N th coefficient f N of an algebraic power series

f (t) = n≥0 f n t n ∈ k[[t]],
where N is assumed to be a large positive integer? This question was recognized as a very important one in complexity theory, as well as in various applications to algorithmic number theory: Atkin-Swinnerton-Dyer congruences, integer factorization, discrete logarithm and point-counting [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF][START_REF] Bostan | Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator[END_REF].

As such, the question is rather vague; both the data structure and the computation model have to be stated more precisely. The algebraic series f will be specified in k [[t]] as some root of a polynomial E(t, y) in k[t, y], of degree d = deg y E ≥ 1 and of height h = deg t E. To do this specification unequivocally, we will make several assumptions. First, we assume that E is separable, that is E and its derivative E y = ∂E/∂y are coprime in k(t) [y]. Second, we assume that E is irreducible 1 in k(t) [y]. Note that both assumptions are satisfied if E is assumed to be the minimal polynomial of f and that irreducibility implies separability as soon as we know that E has at least one root in k [[t]]. The polynomial E might have several roots in k [[t]]. In order to specify uniquely its root f , we further assume that we are given 1 The first assumption is not always implied by the second one, as exemplified by E = y pt ∈ Fp[t, y], and in general by any irreducible polynomial E in k[t, y p ]. 1 a nonnegative integer ρ together with f 0 , . . . , f 2ρ in k such that

E(t, f 0 + f 1 t + • • • + f 2ρ t 2ρ ) ≡ 0 (mod t 2ρ+1 ), E y (t, f 0 + f 1 t + • • • + f ρ t ρ ) ≡ 0 (mod t ρ+1 ).
In other words, the data structure used to represent f is the polynomial E together with the initial coefficients f 0 , . . . , f 2ρ . (Actually ρ+1 coefficients are enough to ensure the uniqueness of f . However 2ρ+1 coefficients are needed to ensure its existence; for this reason, we will always assume the coefficients of f are given up to index 2ρ.) We observe that it is always possible to choose ρ less than or equal to the t-adic valuation of the y-resultant of E and E y , hence a fortiori ρ ≤ (2d-1)h.

Under these assumptions, the classical Newton iteration [START_REF] Kung | All algebraic functions can be computed fast[END_REF] allows the computation of the first N coefficients of f in quasi-linear complexity Õ(N ). Here, and in the whole article (with the notable exception of Section 4), the algorithmic cost is measured by counting the number of basic arithmetic operations (+, -, ×, ÷) and applications of the Frobenius map (x → x p ) and of its inverse (x → x 1/p ) in the ground field k. The soft-O notation Õ(•) indicates that polylogarithmic factors in the argument are omitted. Newton's iteration thus provides a quasi-optimal algorithm to compute f 0 , . . . , f N . A natural and important question is whether faster alternatives exist for computing the coefficient f N alone.

With the exception of the rational case (d = 1), where the N th coefficient can be computed in complexity O(log N ) by binary powering [START_REF] Fiduccia | An efficient formula for linear recurrences[END_REF], the most efficient algorithm currently known to compute f N in characteristic 0 has complexity Õ( √ N ) [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF]. It relies on baby step / giant step techniques, combined with fast multipoint evaluation.

Surprisingly, in positive characteristic p, a radically different approach leads to a spectacular complexity drop to O(log N ). However, the big-O term hides a (potentially exponential) dependency in p. The good behavior of this estimate with respect to the index N results from two facts. First, if the index N is written in radix p as (N ℓ-1 . . . N 1 N 0 ) p , then the coefficient f N is given by the simple formula (1)

f N = [(S N ℓ-1 • • • S N1 S N0 f )(0)] p ℓ ,
where the S r (0 ≤ r < p) are the section operators defined by

(2) S r n≥0

g n t n = n≥0 g 1/p pn+r t n .
Note that for the finite field F p the exponents p ℓ in (1) and 1/p in (2) are useless, since the Frobenius map x → x p is the identity map in this case. Second, by Christol's theorem [START_REF] Christol | Ensembles presque periodiques k-reconnaissables[END_REF][START_REF] Christol | Suites algébriques, automates et substitutions[END_REF][START_REF] Harase | Algebraic elements in formal power series rings[END_REF], the coefficient sequence of an algebraic power series f over a perfect field k of characteristic p > 0 is p-automatic: this means that f generates a finite-dimensional k-vector space under the action of the section operators. Consequently, with respect to a fixed k-basis of this vector space, one can express f as a column vector C, the section operators S r as square matrices A r (0 ≤ r < p), and the evaluation at 0 as a row vector R. Formula (1) then becomes

(3) f N = [RA N ℓ-1 • • • A N1 A N0 C] p ℓ .
Since ℓ is about log N , and since the size of the matrices A r does not depend on N , formula (3) yields an algorithm of complexity O(log N ). This observation (for any p-automatic sequence) is due to Allouche and Shallit [START_REF] Allouche | The ring of k-regular sequences[END_REF]Cor. 4.5]. However, this last assertion hides the need to first find the linear representation (R, (A r ) 0≤r<p , C).

As shown in [START_REF] Bostan | Fast computation of the N th term of an algebraic series over a finite prime field[END_REF]Ex. 5], already in the case of a finite prime field, translating the p-automaticity in terms of linear algebra yields matrices A r whose size can be about d 2 hp 2d . Thus, their precomputation has a huge impact on the cost with respect to the prime number p.

In the particular case of a prime field k = F p , and under the assumption E y (0, f 0 ) = 0, this was improved in [START_REF] Bostan | Fast computation of the N th term of an algebraic series over a finite prime field[END_REF] by building on an idea originally introduced by Christol in [START_REF] Christol | Ensembles presque periodiques k-reconnaissables[END_REF]: one can compute f N in complexity Õ((h + d) 5 hp) + O((h + d) 2 h 2 log N ). So far, this was the best complexity result for this task.

Contributions. We further improve the complexity result from [START_REF] Bostan | Fast computation of the N th term of an algebraic series over a finite prime field[END_REF]

down to Õ(d 2 hp + d ω h) + O(d 2 h 2 log N ) (Theorem 3.4, Section 3.2).
Here ω is the exponent of matrix multiplication. In the case where k is a finite field, we propose an even faster algorithm, with bit complexity linear in log N and quasi-linear in √ p (Theorem 4.1, Section 4). It is obtained by blending the approach in Section 3.2 with ideas and techniques imported from the characteristic zero case [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF]. All these successive algorithmic improvements are consequences of our main theoretical result (Theorem 2.2, Section 2.2), which can be thought of as an effective version of Christol's Theorem (and in particular reproves it).

Effective version of Christol's theorem

We keep the notation of the introduction. Christol's theorem is stated as follows.

Theorem 2.1 (Christol).

Let f (t) in k[[t]
] be a formal power series that is algebraic over k(t), where k is a perfect field with positive characteristic. Then there exists a finite-dimensional k-vector space containing f (t) and stable by the section operators.

The aim of this section is to state and to prove an effective version of Theorem 2.1, on which our forthcoming algorithms will be built. Our approach follows the initial treatment by Christol [START_REF] Christol | Ensembles presque periodiques k-reconnaissables[END_REF], which is based on Furstenberg's theorem [START_REF] Furstenberg | Algebraic functions over finite fields[END_REF]Thm. 2]. For the application we have in mind, it turns out that the initial version of Furstenberg's theorem will be inefficient; hence we will first need to strengthen it, considering residues around the moving point f (t) instead of residues at 0. Another new input we shall use is a globalization argument allowing us to compare section operators at 0 and at f (t). This argument is formalized throught Frobenius operators and is closely related to the Cartier operator used in a beautiful geometric proof of Christol's theorem due to Deligne [START_REF] Deligne | Intégration sur un cycle évanescent[END_REF] and Speyer [START_REF] Speyer | Christol's theorem and the Cartier operator[END_REF], and further studied by Bridy [START_REF] Bridy | Automatic sequences and curves over finite fields[END_REF].

2.1. Frobenius and sections. Recall that the ground field k is assumed to be a perfect field of prime characteristic p, for example a finite field F q , where q = p s . Let K = k(t) be the field of rational functions over k and let

L = K[y]/(E).
Since k is a perfect field, the Frobenius endomorphism F : k → k defined by x → x p is an automorphism of k. It extends to a ring homomorphism, still denoted by F, from L[t 1/p ] to L which raises an element of L[t 1/p ] = L 1/p to the power p. This homomorphism is an isomorphism and its inverse writes (4)

F -1 = p-1 r=0 t r/p S r ,
where each S r , with 0 ≤ r < p, maps L onto itself.

The use in (4) of the same notation as in Formula ( 2) is not a mere coincidence. The algebraic series f provides an embedding of L into the field of Laurent series k((t)), which is the evaluation of an element P (y) of L at the point y = f (t). We will call eval f : L → k((t)) the corresponding map, which sends P (y) to P (f (t)). The Frobenius operator extends from L to k((t)), and the same holds for the sections S r (0 ≤ r < p). These extensions are exactly those of Eq. ( 2). The S r 's in Eq. ( 4) then appear as global variants of the S r 's in Eq. [START_REF] Bostan | Fast computation of the N th term of an algebraic series over a finite prime field[END_REF]. Moreover, global and local operators are compatible, in the sense that they satisfy (5)

F • eval f = eval f • F, S r • eval f = eval f • S r .
As for rational functions, the Frobenius operator and the section operators induce, respectively, a ring isomorphism F from K[t 1/p ] onto K and maps σ r (0 ≤ r < p) from K onto K such that F -1 = p-1 r=0 t r/p σ r . The operators F and S r (0 ≤ r < p) are not K-linear but only k-linear. More precisely, for any λ in K[t 1/p ], µ in K, and z in L,

(6) F(λz) = F(λ) F(z) and S r (µz) = p-1 s=0 t ⌊ r+s p ⌋ σ s (µ) S r-s (z).
In other words both F and F -1 are actually semi-linear.

≡ Q E y (mod E).
The rest of this subsection is devoted to the proof of Theorem 2.2. Although mainly algebraic, the proof is based on the rather analytic remark that any algebraic function in k(t)[f ] can be obtained as the residue at T = f of some rational function in k(t, T ) (see Lemma 2.3). This idea was already used in Furstenberg [START_REF] Furstenberg | Algebraic functions over finite fields[END_REF], whose work has been inspiring for us. The main new insight of our proof is the following: we replace several small branches around zero by a single branch around a moving point. In order to make the argument work, we shall need further to relate the behavior of the section operators around 0 and around the aforementioned moving point. This is where the reinterpretation of the S r 's in terms of Frobenius operators will be useful.

We consider the ring H = k((t))[[T ]] of power series over k((t)). Its fraction field is the field K = k((t))((T )) of Laurent series over k((t)). There is an embedding k((t))[y] → H taking a polynomial in y to its Taylor expansion around f . Formally, it is simply obtained by mapping the variable y to f +T . It extends to a field extension k((t))(y) → K. We will often write P (t, f +T ) for the image of P (t, y) ∈ k((t))(y) in K. The field K is moreover endowed with a residue map res : K → k((t)), defined by res

∞ i=v a i T i = a -1 (by convention, a -1 = 0 if v > -1). It is clearly k((t))-linear.
Lemma 2.3. For any polynomial P ∈ k((t))[y], the following equality holds:

res P (t, f +T ) E(t, f +T ) = P (t, f ) E y (t, f ) .
Proof. Since f is a simple root of E, the series E(t, f +T ) has a simple zero at T = 0. This means that it can be written E(t, f +T ) = T • q(T ) with q ∈ H, q(0) = 0. Taking the logarithmic derivative with respect to T gives

E y (t, f +T ) E(t, f +T ) = 1 T + q ′ (T ) q(T ) ,
akin to [START_REF] Furstenberg | Algebraic functions over finite fields[END_REF]Formula (15), p. 276], from which we derive

P (t, f +T ) E(t, f +T ) = g(T ) T + g(T ) q ′ (T ) q(T ) ,
where g(T ) = P (t, f +T )/E y (t, f +T ). Since E y (t, f +T ) does not vanish at T = 0, the series g(T ) has no pole at 0. Therefore, the residue of g(T )/T is nothing but g(0). Besides the residue of the second summand g(T ) q ′ (T )/q(T ) vanishes. All in all, the residue of

P (t, f +T )/E(t, f +T ) is g(0) + 0 = P (t, f )/E y (t, f ).
We now introduce analogues of section operators over K. For this, we first observe that the Frobenius operator x → x p defines an isomorphism

F : K[t 1/p , T 1/p ] → K. Moreover K[t 1/p , T 1/p ] is a field extension of K of degree p 2 . A basis of K[t 1/p , T 1/p ] over K is of course (t r/p T s/p ) 0≤r,s<p
, but it will be more convenient for our purposes to use a different one. It is given by the next lemma.

Lemma 2.4. The family (t r/p (f +T ) s/p ) 0≤r,s<p is a basis of K[t 1/p , T 1/p ] over K.

Proof. For simplicity, we set y = f +T ∈ K. We have:

1 y 1/p • • • y (p-1)/p = 1 T 1/p • • • T (p-1)/p • U
where U is the square matrix whose (i, j) entry (for 0 ≤ i, j < p) is j i f i/p . In particular, U is upper triangular and all its diagonal entries are equal to 1. Thus U is invertible and the conclusion follows.

For r and s in {0, 1, . . . , p-1}, we define the section operators S r,s : K → K by

F -1 = p-1 r=0 p-1 s=0 t r/p (f +T ) s/p S r,s .
(These operations look like those used in [2, §3.2], but they are not exactly the same.) Clearly S r,0 extends the operator S r : k((t)) → k((t)) defined by Eq. ( 2) and S r,s (g p 1 g 2 ) = g 1 S r,s (g 2 ) for all g 1 , g 2 ∈ K. We observe moreover that the S r,s 's stabilize the subrings k((t))[y] and k[t, y], since y corresponds to f +T .

Proposition 2.5. The following commutation relation holds over K:

S r • res = res • S r,p-1 . Proof. Let us write g ∈ K as g = ∞ i=v a i T i with v ∈ Z and a i ∈ k((t)) for all i ≥ v.
Its image under F -1 can be expressed in two different ways as follows:

F -1 (g) = ∞ i=v F -1 (a i ) T i/p = p-1 r=0 p-1 s=0 t r/p (f +T ) s/p S r,s (g).
We identify the coefficient in T -1/p . For doing so, we observe that the terms obtained with s < p -1 do not contribute, while the contribution of the term t r/p (f +T ) (p-1)/p S r,p-1 (g) is the residue of t r/p S r,p-1 (g). We then get

F -1 (a -1 ) = p-1 r=0 res • S r,p-1 (g) • t r/p .
Going back to the definition of S r , we derive S r (a -1 ) = res • S r,p-1 (g), from which the lemma follows.

Proof of Theorem 2.2. Let P ∈ k[t, y] and 0 ≤ r < p. We set Q = S r,p-1 (P E p-1 ) ∈ k[t, y]. Combining Lemma 2.3 and Proposition 2.5, we derive the following equalities:

S r P (t, f ) E y (t, f ) = S r • res P (t, f +T ) E(t, f +T ) = res • S r,p-1 P (t, f +T ) E(t, f +T ) = res Q(t, f +T ) E(t, f +T ) = Q(t, f ) E y (t, f ) (compare
S r P (t, f ) E y (t, f ) = Q(t, f ) E y (t, f ) with deg t Q ≤ h and deg y P < d. Moreover deg t Q < h provided that r > 0.
Another remark in this direction is the following: if P has degree at most d-2, the section S r,p-1 (P E p-1 ) has degree at most d-2 for any r ∈ {0, 1, . . . , p-1}. Indeed, P E p-1 has degree at most pd -2 < p(d-1) + p -1. In other words, the subspace k[t, y] <h,≤d-2 is stable by the section operators S r (0 ≤ r < p).

Application to algorithmics

Theorem 2.2 exhibits an easy-to-handle finite dimensional vector space which is stable under the section operators. In this section, we derive from it two efficient algorithms that compute the N th term of f in linear time in log N . The first is less efficient, but easier to understand; we present it mainly for pedagogical purposes.

3.1. First algorithm: modular computations. The first algorithm we will design follows rather straightforwardly from Theorem 2.2. It consists of the following steps:

(1) we compute the matrix giving the action of the Frobenius F with respect to the "modified monomial basis" B = (y j /E y ) 0≤j≤d-1 ;

(2) we deduce the matrix of F -1 with respect to B;

(3) we extract from it the matrices of the section operators S r ; (4) we compute the N th coefficient of f using Formula (1).

FAST COEFFICIENT COMPUTATION FOR ALGEBRAIC POWER SERIES IN CHAR. > 0 7 M =      
1 + 2t 4 + 4t 5 + 3t 6 + 2t 7 + 2t 8 + 2t 12 + 4t 13 + 3t 14 + t 15 + 2t 16 • • • t + 2t 2 + 3t 3 + 4t 5 + 4t 6 + 3t 7 + 3t 8 + 3t 9 + t 10 + 4t 11 + t 13 + 2t 15 + t 16 • • • 1 + 2t + 3t 2 + 3t 5 + 2t 6 + 2t 7 + t 8 + t 9 + 3t 10 + t 11 + 3t 12 + 3t 14 + 4t 15 + 3t 16 •

• • 0 • • •       (mod t 17 ) M -1 =       1 + 3t
4 + 2t + • • • + 2t 9 4t 12 + 2t + • • • + +2t 9 + 4t 12 1 + 2t + • • • + 3t 13 + t 16 2 + 4t + • • • + 4t 5 + 2t 8 0 0 0 1 + 4t + • • • + 4t 13 + t 16       Figure 1.
Frobenius and its inverse in the "modified monomial basis".

Let us be a bit more precise (though we will not give full details because we will design in §3.2 below an even faster algorithm). Let M be the matrix of F in the basis B; its jth column contains the coordinates of the vector F( y j Ey ) =

y pj E p
y in the basis B, which are also the coordinates of y pj /E p-1 y in the monomial basis (1, y, . . . , y d-1 ). It is easily seen that the matrix of F -1 with respect to B is

F -1 (M -1
), which is, by definition, the matrix obtained by applying F -1 to each entry of M -1 . We now discuss the complexity of the computation of M -1 . Thanks to Theorem 2.2 and Eq. ( 4), we know that its entries are polynomials of degree at most h(p-1). However, this bound is not valid for the entries of M . Indeed, in full generality, the latter are rational fractions whose numerators and denominators have degrees of magnitude dhp. In order to save the extra factor d, we rely on modular techniques: we choose a polynomial B of degree h(p-1) + 1 and perform all computations modulo B. To make the method work, B must be chosen in such a way that both M and M -1 make sense modulo B, i.e. B must be coprime with the denominators of the entries of M . The latter condition discards a small number of choices, so that a random polynomial B will be convenient with high probability.

Using fast polynomial and matrix algorithms, the computation of M modulo B can be achieved within Õ(d 2 hp) operations in k, while the inversion of M modulo B requires Õ(d ω hp) operations in k, where ω ∈ [START_REF] Bostan | Fast computation of the N th term of an algebraic series over a finite prime field[END_REF][START_REF] Bostan | Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator[END_REF] is the matrix multiplication exponent. Since we count an application of F -1 : k → k as a unique operation, the cost of the first two steps is Õ(d ω hp) as well. The third step is free as it only consists in reorganizing coefficients. As for the evaluation of the formula (1), each application of S r has a cost of O(d 2 h 2 ) operations in k. The total complexity of our algorithm is then Õ(d ω hp) + O(d 2 h 2 log N ) operations in k.

Remark 3.1. We do not need actually to apply the Frobenius inverse F -1 : k → k since, at the end of the computation, we raise the last intermediate result at the power p ℓ . The complexity Õ(d ω hp) + O(d 2 h 2 log N ) can then be reached even if we do not count an application of F -1 as a single operation.

A detailed example. Consider k = F 5 and the polynomial

E = (t 4 + t + 1)y 4 + y 2 + y -t 4 ∈ k[t, y]. It admits a unique root f in k[[t]] which is congruent to 0 modulo t.
The matrix M of the Frobenius F with respect to the basis B = ( 1 Ey , y Ey , y 2 Ey , y 3 Ey ) writes D -1 • M , where D and the largest entry of M have degrees 55 and 39, respectively. However, by Theorem 2.2, we know that M -1 has polynomial entries of degree at most 16. Noticing that 0 is not a root of the resultant of E and E y , we can compute M and its inverse modulo B(t) = t 17 . The result of this computation is displayed partly on Figure 1. We observe that the maximal degree of the entries of M -1 is 16 and reaches our bound h(p-1) (which is then tight for this example). We furthermore observe that M is block triangular, as expected after Remark 2.6.

Let us now compute the images of y ∈ L under the section operators. For this, we write y = E -1 y • 4t 4 + 2y + 3y 2 in L. We then have to compute the product M -1 • (4t 4 2 3 0) T . As a result, we obtain: 16 + 4t 17 + 4t 20 t + 3t 4 + 2t 5 + t 8 + 3t 9 + 4t 10 + 3t 12 + 3t 13 + 4t 16 1 + 2t 2 + 3t 3 + 4t 4 + t 5 + t 6 + 4t 7 + 4t 8 + 3t 9 + 2t 10 + 2t 13 + 4t 16 0

     t 4 + 4t 8 + 2t 12 + 4t
     .
Rearranging the terms, we finally find that

S 0 (y) = E -1 y • 4t 4 + (2t + 4t 2 )y + (1 + t + 2t 2 )y 2 S 1 (y) = E -1 y • 4t 3 + (1 + 4t 3 )y + (t + 4t 3 )y 2 S 2 (y) = E -1 y • (2t 2 + 4t 3 ) + 3t 2 y + (2 + 4t)y 2 S 3 (y) = E -1 y • 4t + (t + 3t 2 )y + (3 + 4t + 2t 2 )y 2 S 4 (y) = E -1 y • 1 + (3 + 3t)y + (4 + 3t)y 2 .
To conclude this example, suppose that we want to compute the 70th coefficient of f . Applying Eq. ( 1), we find that it is equal to the constant coefficient of S 2 S 4 S 0 f . Therefore we have to compute S 2 S 4 S 0 y. Repeating twice what we have done before, we end up with S 2 S 4 S 0 y = E -1 y • (2 + t 2 ) + (4 + 3t + 3t 3 )y + (2 + 4t 2 + 2t 3 )y 2 . Plugging y = f in the above equality, we get S 2 S 4 S 0 f = 2 + O(t), from which we conclude that f 70 = 2. Remark 3.2. In the above example, only the constant coefficient of f was needed to carry out the whole computation. This is related to the fact that E y (f (t)) has t-adic valuation 0. More generally if E y (f (t)) has t-adic valuation ρ, we will need the first ρ+1 coefficients of f since the final division by E y will induce a loss of t-adic precision of ρ "digits". This does not change the complexity bound, since ρ ≤ deg t Res y (E, E y ) ∈ O(dh).

3.2.

Second algorithm: Hermite-Padé approximation. For obvious reasons related to the size of the computed objects, we cannot hope to achieve a complexity lower than linear with respect to p using the approach of Section 3.1. However, the exponent on d still can be improved. In order to achieve this, we return to Theorem 2.2. The key idea is to leap efficiently from the polynomial P to the polynomial Q in Formula [START_REF] Christol | Suites algébriques, automates et substitutions[END_REF].

Let P = d-1 i=0 a i (t)y i in k[t, y] <h,<d and 0 ≤ r < p. By Theorem 2.2, there exists

Q = d-1 i=0 b i (t)y i in k[t, y] <h,<d such that S r (P/E y ) ≡ Q/E y (mod E), or, equivalently, (8) S r d-1 i=0 a i (t) f (t) i E y (t, f (t)) = d-1 j=0 b j (t) f (t) j E y (t, f (t))
.

The algorithmic question is to recover efficiently the b i 's starting from the a i 's.

Identifying coefficients in Eq. ( 8) yields a linear system over k in the coefficients of the unknown polynomials b i . This system has hd unknowns and an infinite number of linear equations. The point is that the following truncated version of Eq. ( 8)

(9) S r d-1 i=0 a i (t) f (t) i E y (t, f (t)) ≡ d-1 j=0 b j (t) f (t) j E y (t, f (t)) (mod t 2dh )
is sufficient to uniquely determine Q. This is a direct consequence of the following.

Lemma 3.3. If Q in k[t, y] <h,<d satisfies Q Ey (t, f (t)) ≡ 0 (mod t 2dh ), then Q = 0.
Proof. The resultant r(t) of E(t, y) and Q(t, y) with respect to y is a polynomial of degree at most d(h-1) + h(d-1). On the other hand, we have a Bézout relation

E(t, y) u(t, y) + Q(t, y) v(t, y) = r(t),
where u(t, y) and v(t, y) are bivariate polynomials in k[t, y]. By evaluating the previous equality at y = f (t) it follows that

r(t) ≡ Q(t, f (t)) v(t, f (t)) ≡ 0 (mod t 2dh )
holds in k((t)), and therefore r = 0. Thus E and Q have a non-trivial common factor; since

E is irreducible, it must divide Q. But deg y Q < deg y E, so Q = 0.
Solving Eq. ( 9) amounts to solving a Hermite-Padé approximation problem. In terms of linear algebra, it translates into solving a linear system over k in the coefficients of the unknown polynomials b i . This system has dh unknowns and N = 2dh linear equations. Moreover, it has a very special shape: it has a quasi-Toeplitz structure, with displacement rank ∆ = O(d). Therefore, it can be solved using fast algorithms for structured matrices [START_REF] Pan | Structured matrices and polynomials[END_REF][START_REF] Bostan | Solving structured linear systems with large displacement rank[END_REF] in Õ(∆ ω-1 N ) = Õ(d ω h) operations in k. These algorithms first compute a (quasi)-inverse of the matrix encoding the homogenous part of the system, using a compact data-structure called displacement generators (or, ΣLU representation); then, they apply it to the vector encoding the inhomogeneous part. The first step has complexity Õ(∆ ω-1 N ) = Õ(d ω h), the second step has complexity Õ(∆N ) = Õ(d 2 h).

In our setting, we will need to solve log N systems of this type, each corresponding to the current digit of N in radix p. An important feature is that these systems share the same homogeneous part, which only depends on the coefficients of the power series s j (t) = f j Ey(t,f (t)) occurring on the right-hand side of [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF]. Only the inhomogeneous parts vary: they depend on the linear combination

d-1 i=0 a i (t)s i (t).
Putting these facts together yields Algorithm 1 and the following complexity result. Theorem 3.4. Let k be a perfect field with characteristic p > 0. Let E(t, y) be an irreducible polynomial in k[t, y] of height h and degree d. We assume that we are given a nonnegative integer ρ and a polynomial f (t) such that E(t, f (t)) ≡ 0 (mod t 2ρ+1 ) and E y (t, f (t)) ≡ 0 (mod t ρ+1 ).

There there exists a unique series f (t) congruent to f (t) modulo t ρ+1 for which

E(t, f (t)) = 0. Moreover, Algorithm 1 computes the N th coefficient of f for a cost of Õ(d 2 hp + d ω h) + O(d 2 h 2 log N ) operations in k.
Proof. The first assertion is Hensel's Lemma [START_REF] Eisenbud | Commutative Algebra with a View Toward Algebraic Geometry[END_REF]Th. 7.3].

Algorithm Nth coefficient via Hermite-Padé. Input: A polynomial E(t, y) = e d (t)y d + • • • + e 0 (t) and a truncation

g = f 0 + • • • + O(t ρ+1 ) of a series f such that E(t, g) = O(t ρ+1
). Output: The N th coefficient f N of the series f . 1. Precompute the first 2pdh coefficients of the series expansions s j of f (t) j /E y (t, f ), 0 ≤ j < d. . 5. For i = 0, 1, . . . , ℓ -1, (1) write g = P (t, f )/E y (t, f ) as a linear combination of the s j 's, (2) compute the section S Ni (g) at precision O(t 2dh ), ( 3) recover Q such that S Ni (g) = Q/E y by Hermite-Padé, (4) redefine g as Q/E y 6. Replace y by f (t) in g and call ḡ(t) the obtained result. 7. Expand ḡ(t) at precision O(t). 8. Set ḡ0 to the constant coefficient of ḡ(t). 9. Return ḡp ℓ 0 .

Algorithm 1

The precomputation of s j (t) = f (t) j Ey(t,f (t)) modulo t 2dhp for 0 ≤ j < d can be performed using Newton iteration, for a total cost of Õ(d2 hp) operations in k. As explained above, this is enough to set up the homogeneous part of the quasi-Toeplitz system; its inversion has cost Õ(d ω h).

Let us turn to the main body of the computation, which depends on the index N . For each p-digit r = N i of N , we first construct the inhomogeneous part of the system. For this, we extract the coefficients of t pj+r in d-1 i=0 a i (t)s i (t), for 0 ≤ j < d, for a total cost of O(d 2 h 2 ) operations in k. We then apply the inverse of the system to it, for a cost of O(d 2 h 2 ) (using a naive matrix vector multiplication 2 ). This is done ℓ ≈ log N times. The other steps of the algorithm have negligible cost.

Improving the complexity with respect to p

As shown in Theorem 3.4, Algorithm 1 has a nice complexity with respect to the parameters d, h and log N : it is polynomial with small exponents. However, the complexity with respect to p is not that good, as it is exponential in log p, which is the relevant parameter. Thus, when p is large (say > 10 5 ), Algorithm 1 runs slowly and is no longer usable.

For this reason, it is important to improve the complexity with respect to p. In this section, we introduce some ideas to achieve this. More precisely, our aim is to design an algorithm whose complexity with respect to p and N is Õ( √ p) • log N , and remains polynomial in all other relevant parameters. In the current state of knowledge, it seems difficult to decrease further the exponent on p; indeed, the question addressed in this paper is related to other intensively studied questions (e.g., counting points via p-adic cohomologies) for which the barrier Õ( √ p) has not been overcome yet.

Notations and assumptions. We keep the notation of previous sections. We make one additional hypothesis: the ground field k is a finite field. We assume that k is represented as (Z/pZ)[X]/π(X) where π is an irreducible monic polynomial over Z/pZ of degree s. We choose a monic polynomial π ∈ Z[X] of degree s lifting π.

We set W = Z p [X]/π(X) where Z p is the ring of p-adic integers.

The algorithm we are going to design is not algebraic in the sense that it does not only perform algebraic operations in the ground field k, but will sometimes work over W (or, more exactly, over finite quotients of W ). For this reason, throughout this section, we will use bit complexity instead of algebraic complexity.

We use the notation poly(n) to indicate a quantity whose growth is at most polynomial in n. The precise result we will prove reads as follows. If p is bounded by a (fixed) polynomial in d and h, then Theorem 4.1 has been proved already. In the sequel, we will then always assume that p ≫ d, h.

Overview of the strategy. We reuse the structure of Algorithm 1 but speed up the computation of the S Ni (g)'s. Precisely, in Algorithm 1, the drawback was the computation of the f j Ey(t,f ) 's or, almost equivalently, the computation of g = P (t,f )

Ey(t,f )
at sufficient precision. However, only a few (precisely 2dh) coefficients of g are needed, since we are only interested in one of its sections. A classical method for avoiding this overhead is to find a (small) recurrence on the coefficients on g = ∞ n=0 g i t i of the form: [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF] b r (i)

g i+r + b r-1 (i)g i+r-1 + • • • + b 1 (i)g i+1 + b 0 (i)g i = 0.
We then unroll it using matrix factorials (for which fast algorithms are available in the literature [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF] Unrolling the recurrence is straightforward as soon as the leading coefficient b r (i) does not vanish. On the contrary, when b r (i) = 0, the value of g i+r cannot be deduced from the previous ones. Unfortunately, it turns out that b r (i) does sometimes vanish in our setting. We tackle this issue by lifting everything over W and performing all computations over this ring. Divisions by p then become possible but induce losses of precision. We then need to control the p-adic valuation of the denominators, that are the p-adic valuations of the b r (i)'s. We cannot expect to have a good control on them in full generality; even worse, we can build examples where b r (i) vanishes in W for some i. There exists nevertheless a good situation-the so-called ordinary casewhere we can say a lot on the b r (i)'s. With this extra input, we are able to lead our strategy to its end.

The general case reduces to the ordinary one using a change of origin, i.e. replacing t by u+α for some α ∈ k. This change of origin does not seem to be harmless a priori. Indeed the Taylor expansion of g around α (the one we shall compute) has in general nothing to do with the Taylor expansion of g around 0 (the one we are interested in). The sections are nevertheless closely related (see Proposition 4.3). This "miracle" is quite similar to what we have already observed in Proposition 2.5 and again can be thought of as an avatar of the Cartier operator. 4.1. From algebraic equations to recurrences. We consider a bivariate polynomial P (t, y) ∈ k[t, y] with deg t P < h and deg y P < d. We fix moreover an integer r is the range [0, p-1]. Our aim is to compute S r P (t,f )

Ey(t,f ) at precision O(t 2dh ). Set g = P (t,f ) Ey(t,f ) and write g = ∞ i=0 g i t i . By definition S r (g) = ∞ j=0 g
1/p r+pj t j , so that we have to compute the coefficients g r+pj for j < 2dh.

We let L be the leading coefficient of E(t, y) and R be the resultant of E and E y . To begin with, we make the following assumption (which will be relaxed in §4.3):

(H1): Both L and R have t-adic valuation 0.

As explained above, we now lift the situation over W . We choose a polynomial Ê ∈ W [t, f ] of bidegree (h, d) lifting E. We define Êt = ∂ Ê ∂t , Êy = ∂ Ê ∂y . The assumption (H1) implies that the series f lifts uniquely to a series f ∈ W [[t]] such that Ê(t, f ) = 0. We define L as the leading coefficient of Ê(t, y) and set R = Res( Ê, Êy ). We introduce the ring

W K = W [t, ( L R) -1 ]. By (H1), W K embeds canonically into W [[t]
]. We pick a polynomial P ∈ W [t, y] lifting P such that deg t P < h and deg y P < d. We set ĝ = P (t, f ).

We now compute a linear differential equation satisfied by ĝ. For this, we observe that the derivation ∂ ∂t :

W [[t]] → W [[t]] stabilizes the subring W L = W K [ f ]
. Indeed from the relation Ê(t, f ) = 0, we deduce that ∂ f ∂t = -Êt(t, f ) Êy(t, f ) . Thus ∂ f ∂t ∈ W L because Êy (t, f ) is invertible in W L thanks to (H1). Using additivity and the Leibniz relation, we finally deduce that ∂ ∂t takes W L to itself. In particular, all the successive derivatives of ĝ lie in W L . On the other hand, we notice that W L is free of rank d over W K with basis (1, f , . . . , f d-1 ). Let M be the d × d matrix whose jth column (for 0 ≤ j < d) contains the coordinates of ∂ j ĝ ∂t j with respect to the above basis. Similarly let C be the column vector whose entries are the coordinates of

∂ d ĝ ∂t d . Let ∆ d =
det M . We solve the system M X = C using Cramer's formulae and find this way a linear differential equation of the form:

∆ d ∂ d ĝ ∂t d + ∆ d-1 ∂ d-1 ĝ ∂t d-1 + • • • + ∆ 1 ∂ĝ ∂t + ∆ 0 ĝ = 0,
where the other ∆ i 's are defined as determinants as well. In particular, they all lie in W K . Multiplying by the appropriate power of L R, we end up with a differential equation of the form:

(11) âd ∂ d ĝ ∂t d + âd-1 ∂ d-1 ĝ ∂t d-1 + • • • + â1 ∂ĝ ∂t + â0 ĝ = 0
where the âi 's are polynomials in t. We can even be more precise. Indeed, following the above constructions, we find that all entries of M and C are rational functions whose degrees (of numerators and denominators) stay within poly(dh). We then deduce that the degrees of the ∆i 's and âi 's are in poly(dh) as well. Furthermore, they can be computed for a cost of poly(dh) operations in k, that is poly(dh) Õ(s log p) bit operations (recall that s denotes the degree of k over F p )

We write ĝ = ∞ i=0 gi t i i! . The differential equation [START_REF] Deligne | Intégration sur un cycle évanescent[END_REF] translates to a recurrence relation on the gi 's of the form: [START_REF] Eisenbud | Commutative Algebra with a View Toward Algebraic Geometry[END_REF] ∀n ≥ r, b0 (n

)g n + b1 (n)g n-1 + b2 (n)g n-2 + • • • + br (n)g n-r = 0
where the bi 's are polynomials in n over W whose degrees are in poly(dh). Moreover r is at most d + max i deg âi . In particular, r ∈ poly(dh). Finally it is easy to write down explicitly b0 : it is the constant polynomial with value âd (0).

4.2.

The ordinary case. In order to take advantage of Eq. ( 12), we make the following extra assumption, corresponding to the so-called ordinary case:

(H2): âd (0) does not vanish modulo p.

Under (H2), b0 (n) = âd (0) is invertible in W and there is no obstruction to unrolling the recurrence [START_REF] Eisenbud | Commutative Algebra with a View Toward Algebraic Geometry[END_REF]. Let us be more precise. We recall that we want to compute the values of g r+pj for j up to 2dh. Clearly g n is the reduction modulo p of gn n! . In order to get g r+pj , we need to compute gr+pj modulo p v+1 where v is the p-adic valuation of (r + 2dhp)!. Under our assumption that p is large enough compared to d and h, we get v = 2dh. We will then work over the finite ring W ′ = W/p 2dh+1 W .

We first compute the r first coefficients of f modulo p 2dh+1 by solving the equation Ê(t, f ) = 0 (using a Newton iteration for example). Since r ∈ poly(dh), this computation can be achieved for a cost of poly(dh) operations in W ′ , that is poly(dh) Õ(s log p) bit operations. We then build the companion matrix:

M (n) =       1 . . . 1 -br (n) âd (0) -br-1(n) âd (0) • • • -b1(n) âd (0)       ∈ (W ′ [n]) r×r . Obviously, gn-r+1 gn-r+2 • • • gn T = M (n) • M (n-1) • • • M (r) • g0 g1 • • • gr-1 T ,
and computing gn reduces to evaluating the matrix factorial M

(n)•M (n-1) • • • M (r).
Using [START_REF] Chudnovsky | Approximations and complex multiplication according to Ramanujan[END_REF], the latter can be computed within poly(dh) Õ( √ n) operations in W ′ , that is poly(dh) Õ( √ n • s log p) bit operations. All in all, we find that the g r+pj 's (0 ≤ j < 2dh) can all be computed for a cost of poly(dh) Õ(s √ p) bit operations. Plugging this input in Algorithm 1, we end up with an algorithm of bit complexity poly(dh) Õ(s √ p) log N . Theorem 4.1 is thus proved under the extra assumptions (H1) and (H2).

4.3.

Reduction to the ordinary case. We finally explain how (H1) and (H2) can be relaxed. The rough idea is to translate the origin at some point where these two hypotheses hold simultaneously.

The case of complete vanishing. Before proceeding, we need to deal with the case where the whole polynomial âd vanishes modulo p. This case is actually very special; this is shown by the next lemma, whose proof relies on the fact that for a generic g, the minimal-order (homogeneous) linear differential equation over k(t) satisfied by g has order exactly d [START_REF] Chudnovsky | On expansion of algebraic functions in power and Puiseux series[END_REF]. Lemma 4.2. For a generic g ∈ L = k(t)[y]/E(t, y), the reduction of âd modulo p does not vanish.

We say that an element g ∈ L is good if the corresponding âd does not vanish modulo p. Lemma 4.2 ensures that goodness holds generically. It then holds with high probability since we have assumed that the ground field k has a large cardinality. Consequently, even if we were unlucky and g was not good, we can produce with high probability a decomposition g = g 1 + g 2 where g 1 and g 2 are both good (just by sampling g 1 at random). Since moreover the section S r is additive, we can recover S r (g) as S r (g 1 ) + S r (g 2 ).

For this reason, in what follows we will assume safely that g is good.

Change of origin. Let α ∈ W be such that L(α) ≡ 0 (mod p), R(α) ≡ 0 (mod p), âd (α) ≡ 0 (mod p). Such an element exists (since k is assumed to be large) and can be found for a cost of poly(dh) operations in k (e.g., by enumerating its elements). We denote by α ∈ k the reduction of α modulo p and assume that α = 0 (otherwise, we are in the ordinary case). We perform the change of variable τ α : t → u+α. Note that τ α induces isomorphisms k(t) → k(u) and k(t)[y]/E(t, y) → k(u)[y]/E(u-α, y). Furthermore, the polynomial E(α, y) = 0 has d simple roots in an algebraic closure of k. Let f α,0 be one of them. By construction, f α,0 lies in a finite extension ℓ of k of degree at most d. Moreover, by Hensel's Lemma, f α,0 lifts uniquely to a solution: Here S r and S r,u refer to the section operators defined in the usual way. We observe that they stabilize the subfields k(t)[y] E(t,y) and k(u)[y] E(u+α,y) , respectively, since they can alternatively be defined by the relations: [START_REF] Fiduccia | An efficient formula for linear recurrences[END_REF] over k(t)[y] E(t,y) :

f α = f α,0 + f α,1 u + • • • + f α,i u i + • • • ∈ ℓ[[u]]
F -1 =
p-1 r=0 t r/p S r over k(u) [y] E(u-α,y) :

F -1 =
p-1 r=0 u r/p S r,u where F is the Frobenius map (see also Eq. ( 4)). We emphasize that the other section operators S r,⋆ (with r < p-1) do not commute with τ α : the above phenomenon is specific to the index p-1. However, we can relate S r and S p-1,u as follows.

Corollary 4.4. For all g ∈ k(t) [y] E(t,y) , we have S r (g) = τ -1 α • S p-1,u • τ α (t p-1-r g).

Proof. This follows from Proposition 4.3 and from S r (g) = S p-1 (t p-1-r g).

A modified recurrence. In order to use Corollary 4.4, we need to check that τ α (t p-1-r g) fits the ordinary case. Recall the differential equation satisfied by ĝ, âd

∂ d ĝ ∂t d + âd-1 ∂ d-1 ĝ ∂t d-1 + • • • + â1 ∂ĝ ∂t + â0 ĝ = 0.
We set r ′ = p -1 -r and Ĝ = t r ′ ĝ. Applying Leibniz formula to ĝ = t -r ′ Ĝ, we get:

∂ j ĝ ∂t j = j i=0
(-1) i j i r ′ (r ′ + 1) • • • (r ′ + i -1)t -r ′ -i ∂ j Ĝ ∂t j-i , from which we derive the following differential equation satisfied by Ĝ:

0≤i≤j≤d (-1) i âj j i r ′ (r ′ + 1) • • • (r ′ + i -1)t -r ′ -i ∂ j-i Ĝ ∂t j-i .
Reorganizing the terms and multiplying by t r ′ +d , we end up with: , it now just remains to apply τ -1 α to get S r (g). This last operation can be performed for a cost of poly(dh) operations in k as well. All in all, we are able to compute S r (g) for a total bit complexity of poly(dh) Õ(s √ p). Repeating this process log N times, we obtain the complexity announced in Theorem 4.1.

2 .

 2 Precompute the quasi-inverse of the Toeplitz matrix corresponding to the Hermite-Padé approximation problem. 3. Expand N = (N ℓ-1 . . . N 0 ) p with respect to the radix p. 4. Set g = y ∈ L written as E -1 y • -de 0 -(d-1)e 1 y -• • • -e d-1 y d-1

Theorem 4 . 1 .

 41 Under the assumptions of Theorem 3.4 and the above extra assumptions, there exists an algorithm of bit complexity poly(dh) Õ(s √ p) log N that computes the N th coefficient of f .

  to the equation E(u-α, y) = 0. We emphasize that the morphism k(t)[y]/E(t, y) → k(u)[y]/E(u-α, y) does not extend to a mapping k((t)) → ℓ((u)) sending f to f α . The next diagram summarizes the previous discussion:

Proposition 4 . 3 .

 43 The commutation S p-1,u • τ α = τ α • S p-1 holds over k(t)[y] E(t,y) . Proof. Clearly τ α commutes with the Frobenius because it is a ring homomorphism. From the relations (13), we then derivep-1 r=0 u r/p S r,u •τ α = p-1 r=0 (u+α) r/p τ α •S r .Identifying the coefficients in u p-1 p , we get the announced result.

(- 1 )

 1 i âi+j i + j i r ′ (r ′ + 1) • • • (r ′ + i -1)t d-i ∂ j Ĝ ∂t j .Set W L,u = W [u, y]/ Ê(u+α, y) and define the ring homomorphism τ α : W L → W L,u , t → u+ α, y → y. Clearly τ α lifts τ α . Applying τ α to Eq. (14) and noticing that ∂ ∂t = ∂ ∂u , we obtain:i τα (a i+j ) i + j i r ′ (r ′ + 1) • • • (r ′ + i -1)(u+α) d-i ∂ j τ α( Ĝ) ∂u j .Conclusion.The leading term of the latest differential equation (obtained only with j = d and i = 0) is τ α(â d ) (u+ α) d . Its value at u = 0 is then âd (α) αd , which is not congruent to 0 modulo p by assumption. Moreover the other coefficients are polynomials in u whose degrees stay within poly(dh). Therefore, we can apply the techniques of §4.2 and compute S p-1,u (τ α G) at precision O(u 2dh ) for a cost of poly(dh) Õ(s √ p) bit operations. As explained in §3.2, we can reconstruct S p-1,u (τ α G) as an element of k[u, y]/E(u+α, y) for a cost of poly(dh) operations in k using Hermite-Padé approximations. Thanks to Corollary 4.4

  2.2. The key theorem. Let k[t, y] <h,<d be the set of polynomials P ∈ k[t, y] such that deg t P < h and deg y P < d.

	Theorem 2.2. For P ∈ k[t, y] <h,<d and for 0 ≤ r < p, there exists a (unique) polynomial Q in k[t, y] <h,<d such that
	(7)	S r	P E y

  4 + t 8 + t 12 + t 13 + t16

		t 4 + 2t 8	t 8 + t 12	t 12
	2 + 2t + • • • + t 9 + 3t 12	3 + 2t + • • • + 2t 13 + t 16	3 + 4t + t 5 + t 8	1 + 4t 4 + 3t 5 + 2t 9 + 2t 12

One can actually achieve this step for a cost of Õ(d 2 h) operations in k using the quasi-Toeplitz structure; however this is not that useful since the cost of the previous step was already O(d 2 h 2 ).


Celia Lopez-Ongil

Mario Garcia-Valderas

Marta Portela-Garcia

Luis Entrena

Techniques for Fast Transient Fault Grading Based on Autonomous Emulation

come

I. INTRODUCTION

With the advent of very deep submicron technologies fault tolerance (FT) is becoming a concern for an increasing number of applications, since they can be affected by radiation even in the earth surface. During the hardening process of a circuit FT evaluation is a key factor. The injection of real faults artificially generated in a circuit prototype [START_REF] Karlsson | Application of Three Physical Fault Injection Techniques to the Experimental Assessment of the MARS Architecture[END_REF] is a hardware approach for this task. It produces very realistic results, but requires expensive equipment. Also, identification of weak areas is difficult in the circuit prototype, due to the limited observability at the chip pins. So, re-design cost increases notably.

Software approach for FT evaluation is modelling the behaviour of transient faults. Fault model is injected in a circuit description and validation of the faulty circuit response, by means of simulation, will indicate the robustness of the design. In fault simulation re-design for hardening could be easily done; early location of weak areas will imply a significant reduction in re-design cost and time. However the flexibility achieved in the measuring of the FT implies a great cost in terms of CPU time and resources.

Simulation of large designs is currently being substituted by hardware emulation in platform FPGAs that provides a considerable saving of computational resources. Emulation is controlled by a host computer that interacts with the FPGA hardware. This way, software flexibility is maintained while hardware speed is profited. The application of emulation to FT evaluation has been proved to be a very cost effective solution. It is important to note that in fault emulation only functional response is being tested, not physical implementation response.

In this paper, we proposed a new approach for performing FT evaluation by means of autonomous emulation of the circuit in platform FPGAs. This solution makes profit of the hardware resources, executing most of the task in the FPGA instead of in the host, in order to minimise the bottleneck times in the communication between software and hardware. Three different techniques for autonomous emulation are presented and compared, stating that best technique depends on the characteristics of the circuit. This work was partially supported by research project nb. 07T/0052/2003 2.

II. AUTONOMOUS EMULATION SYSTEM

Hardware emulation in FT evaluation consists in testing the response of a faulty circuit, which has been modified in order to behave according to a fault model. Commonly, bitflip is the fault model adopted for SEU effects. Therefore, only memory elements are affected.

The fastest existing emulation solutions for FT evaluation are based on the modification of the circuit under evaluation adding extra hardware to allow fault injection from a host computer [START_REF] Civera | Exploiting FPGA for accelerating Fault Injection Experiments[END_REF]. The main drawback of this approach is the host-emulator communication bottleneck, because it is required for multiple tasks, like stimuli application, fault injection and output values check.

The proposed autonomous emulation system solves this problem by performing all the required tasks in the FPGA. The system includes a controller to manage all the emulation tasks and a RAM to store inputs, outputs and fault classification. The evaluation process can be performed in a single hardware run, reducing emulator-host communication to the beginning and to the end of the whole process.

Three fault injection techniques have been implemented within the proposed emulation system. The mask-scan technique departs from [START_REF] Civera | Exploiting FPGA for accelerating Fault Injection Experiments[END_REF], adding the required components for the system to become autonomous. The state-scan technique and the time-multiplexed technique have been

Techniques for Fast Transient Fault Grading Based on Autonomous Emulation

Celia López-Ongil, Mario García-Valderas, Marta Portela-García, Luis Entrena-Arrontes Microelectronics Group. Electronic Technology Department University Carlos III of Madrid. Spain celia, mgvalder, mportela, entrena@ing.uc3m.es originally developed in this work.

In the mask-scan technique, an additional flip-flop is added to every circuit flip-flop to be used as a fault mask. Mask flip-flops indicate where a fault is to be injected.

In the state-scan technique, circuit flip-flops form a scan chain that allows the insertion of the whole faulty circuit state. These states are stored in RAM and the emulation controller can retrieve and insert them in the circuit.

Finally, the time-multiplexed technique replaces every circuit flip-flop by the structure shown in Figure 1. In this structure, two flip-flops are used to run the faulty circuit and the golden run alternatively, thus performing time multiplexing. A third flip-flop is used as a fault mask, and a state flip-flop is also used to avoid restarting the emulation from the beginning every time. This technique is quite faster than the previous ones, because it allows detecting fault effects disappearing without executing the whole testbench.

III. EXPERIMENTAL RESULTS AND CONCLUSION

The autonomous emulation system has been implemented in a Celoxica RC1000 board, with a Xilinx Virtex-2000E and 8Mb of onboard RAM.

The b14 circuit from the ITC'99 benchmark suite [START_REF] Corno | RT-Level ITC'99 benchmarks and first ATPG results[END_REF], the Viper processor, has been used. The circuit has 32 inputs, 54 outputs and 215 flip-flops. The experiments have been carried out with a set of 160 stimulus vectors and the complete set of single faults. Synthesis results for the original circuit, the modified versions and the complete emulation circuit are shown in Emulation time in the State-Scan technique is longer due to the need of scanning-in the circuit state for each fault. This is because in this circuit (b14) there are a large number of flip-flops (215) and a short number of test bench cycles (160). This method improves when the number of cycles is higher than the flip-flop number. Time-Multiplexed technique is always the fastest.

Finally as a conclusion, this paper presents a new solution for improving the performance of transient fault emulation in platform FPGAs. An autonomous transient fault emulation system is proposed, which executes most of the tasks involved in the fault injection campaign in the FPGA. Experimental results have shown that the proposed system achieves execution times several orders of magnitude faster than existing solutions, with a reasonable cost in terms of extra hardware thanks to the use of platform FPGAs. A comparison has been made between three techniques that apply autonomous emulation, with respect to emulation time and area overhead.

Figure 1 .

 1 Figure 1. Instrument for the time multiplexed technique

table 1 (

 1 Leonardo Spectrum 2003).

		Board	Modified circuit	Emulator System
		/FPGA	(overhead)	(overhead)
		RAM	LUTs	FFs	LUTs	FFs
	b14 original	-	1,172	215	-	-
	Mask Scan	33 / 13.4 kbits	1,657 (41%)	434 (102%)	2,040 (74%)	670 (211%)
	State Scan	7,289 / 13.4 kbits	1,644 (40%)	433 (101%)	1,728 (47%)	518 (140%)
	Time	67 / 5.3	3,836	859	4,162	1,032
	Multiplex.	kbits	(227%)	(300%)	(255%)	(380%)

Table 1 . Synthesis results for the b14 circuit

 1 The overhead due to circuit modification is proportional to the flip-flop number. State-scan technique also uses a higher amount of memory to store the different faulty states to evaluate. Control block overhead depends on the flip-flop number, test bench cycles and circuit inputs and outputs.The set of 34,400 single faults have been classified into a 49.2% failure, 4.4% latent and 46.4% silent faults.Emulation times are shown in table2. With a clock frequency of 25 MHz the average speed obtained is some orders of magnitude better than fault simulation (1300 us/fault) and emulation in[START_REF] Civera | Exploiting FPGA for accelerating Fault Injection Experiments[END_REF](100 us/fault).

	Autonomous	Emulation	Average speed
	System	time (ms)	(us/fault)
	MaskScan	141.11	4.1
	State Scan	386.40	11.2
	Time Mux.	19.95	0.58

Table 2 . Time results for b14 circuit

 2

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE