domain: math info

text:
Éric Schost

Catherine St-Pierre

p-adic algorithm for bivariate Gröbner bases

We present a p-adic algorithm to recover the lexicographic Gröbner basis G of an ideal in Q[x, y] with a generating set in Z[x, y], with a complexity that is less than cubic in terms of the dimension of Q[x, y]/⟨G⟩ and softly linear in the height of its coefficients. We observe that previous results of Lazard's that use Hermite normal forms to compute Gröbner bases for ideals with two generators can be generalized to a set of t ∈ N + generators. We use this result to obtain a bound on the height of the coefficients of G, and to control the probability of choosing a good prime p to build the p-adic expansion of G.

Introduction

There exists a rich literature dedicated to the solution of polynomial systems in two variables [START_REF] González-Vega | An improved upper complexity bound for the topology computation of a real algebraic plane curve[END_REF][START_REF] Emiris | Real solving of bivariate polynomial systems[END_REF][START_REF] Diochnos | On the asymptotic and practical complexity of solving bivariate systems over the reals[END_REF][START_REF] Alberti | Topology and arrangement computation of semi-algebraic planar curves[END_REF][START_REF] Rouillier | On solving systems of bivariate polynomials[END_REF][START_REF] Berberich | An elimination method for solving bivariate polynomial systems: Eliminating the usual drawbacks[END_REF][START_REF] Emeliyanenko | On the complexity of solving a bivariate polynomial system[END_REF][START_REF] Bouzidi | Rational univariate representations of bivariate systems and applications[END_REF][START_REF] Lebreton | On the complexity of solving bivariate systems: the case of non-singular solutions[END_REF][START_REF] Bouzidi | Improved algorithm for computing separating linear forms for bivariate systems[END_REF][START_REF] Kobel | Improved complexity bounds for computing with planar algebraic curves[END_REF][START_REF] Mehrabi | A softly optimal Monte Carlo algorithm for solving bivariate polynomial systems over the integers[END_REF][START_REF] Kobel | On the complexity of computing with planar algebraic curves[END_REF][START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF][START_REF] Diatta | Bounds for polynomials on algebraic numbers and application to curve topology[END_REF][START_REF] Dahan | Lexicographic Gröbner bases of bivariate polynomials modulo a univariate one[END_REF], due in part to their numerous applications in real algebraic geometry and computeraided design. Our focus in this paper is on the complexity of computing the lexicographic Gröbner basis of a zero-dimensional ideal in Q[x, y], specifically by means of p-adic techniques based on Newton iteration. An important aspect of this work is to give bit-size bounds for such a Gröbner basis, as well as bounds on the number of primes of bad reduction.

p-adic techniques have been considered in the context of Gröbner basis computations (in an arbitrary number of variables) for decades. In 1983 and 1984, Ebert and Trinks addressed the question of modular algorithms for Gröbner bases [START_REF] Ebert | Some comments on the modular approach to Gröbner bases[END_REF][START_REF] Trinks | On improving approximate results of Buchberger's algorithm by Newton's method[END_REF], specifically for systems without multiple roots; these techniques were also used in geometric resolution algorithms [START_REF] Giusti | When polynomial equation systems can be solved fast?[END_REF][START_REF] Giusti | Straight-line programs in geometric elimination theory[END_REF][START_REF] Giusti | Lower bounds for diophantine approximation[END_REF][START_REF] Giusti | A Gröbner-free alternative for polynomial system solving[END_REF].

1
The absence of multiple roots allows for simple and efficient algorithms; for arbitrary inputs, the question is more involved.

Winkler gave the first p-adic algorithm to construct a Gröbner basis [START_REF] Winkler | A p-adic approach to the computation of Gröbner bases[END_REF] that applies to general inputs; Pauer refined the discussion of good prime numbers [START_REF] Pauer | On lucky ideals for Gröbner basis computations[END_REF], and Arnold revisited, and simplified, these previous constructions in [START_REF] Arnold | Modular algorithms for computing Gröbner bases[END_REF]. No complexity analysis was provided; these p-adic algorithms remain complex (they not only lift the Gröbner basis, but also the transformation matrix that turns the input system into its Gröbner basis), and to our knowledge, achieve linear convergence.

In the specific context of bivariate equations, p-adic techniques have already been put to use in previous work, first in the particular case of non-multiple roots [START_REF] Lebreton | On the complexity of solving bivariate systems: the case of non-singular solutions[END_REF], then to compute a set-theoretic description of all roots, even in the presence of multiplicities [START_REF] Mehrabi | A softly optimal Monte Carlo algorithm for solving bivariate polynomial systems over the integers[END_REF]. However, in this case, the latter algorithm does not reveal the local structure at multiple roots.

In [START_REF] Schost | Newton iteration for lexicographic gröbner bases in two variables[END_REF], we presented a form of Newton iteration specifically tailored to lexicographic Gröbner bases in two variables. It crucially rests on results due to Conca and Valla [START_REF] Conca | Canonical Hilbert-Burch matrices for ideals of k[x, y][END_REF], who gave an explicit parametrization of bivariate ideals with a given initial ideal: our lifting algorithm works specifically with the parameters introduced by Conca and Valla. Our contribution in this paper is to build on [START_REF] Schost | Newton iteration for lexicographic gröbner bases in two variables[END_REF] to give a complete p-adic algorithm: we quantify bad primes, show how to initialize the lifting process, give bounds on the size of the output, and analyze the cost of the whole algorithm.

The following theorem gives a slightly simplified form of our main result, where the probability of success and the number of input polynomials are kept constant (the more precise version is given in the last section). In what follows, the height of a nonzero integer u is log(|u|); if G is a family of polynomials in Q[x, y], we define deg(G) = dim Q Q[x, y]/G and let h(G) be the maximum height of the numerators and denominators of its coefficients.

Theorem 1. Let F = (f 1 , . . . , f t ) be in Z[x, y], with degree at most d, height at most h, and with finitely many common solutions in C 2 . Let G be the lexicographic Gröbner basis of I for the order x ≺ y and write s = |G|, δ = deg(G), b = h(G).

For P > 0, assuming P ∈ O(1) and t ∈ O(1), there is an algorithm that computes G with probability of success at least 1 -1/2 P using a number of bit operations softly linear in

d 2 h + (d ω+1 + δ ω ) log(h) + (d 2 δ + dδ 2 + s 2 δ 2 )(b + log(h)).
With the notation in the theorem, the bitsize of the input is linear in d 2 h, and that of the output is linear in sδb. If all solutions of F have multiplicity one, previous forms of Newton iteration achieve better runtimes, softly linear in the output size [START_REF] Giusti | A Gröbner-free alternative for polynomial system solving[END_REF][START_REF] Schost | Computing parametric geometric resolutions[END_REF][START_REF] Lebreton | On the complexity of solving bivariate systems: the case of non-singular solutions[END_REF][START_REF] Dahan | Lifting techniques for triangular decomposition[END_REF] (but instead of a Gröbner basis, they compute a triangular decomposition of V (F), or change coordinates). Hence, it makes sense to apply our techniques only to multiple solutions.

This is what provides the motivation for our second result, where we compute the Gröbner basis of the ⟨x, y⟩-primary component J of I. As a natural extension, one can consider combining this with [START_REF] Hyun | Change of basis for m-primary ideals in one and two variables[END_REF], which shows how to put an arbitrary primary component of I in correspondence with the ⟨x, y⟩-primary component of a related ideal in K[x, y] for a finite extension

K of Q. Theorem 2. Let F = (f 1 , . . . , f t ) be in Z[x, y]
, with degree at most d, height at most h, with finitely many common solutions in C 2 . Let G 0 be the lexicographic Gröbner basis of the ⟨x, y⟩-primary component of I for the order x ≺ y and write r = |G 0 |, η = deg(G 0 ), c = h(G 0 ). For P > 0, assuming P ∈ O(1) and t ∈ O(1), there is an algorithm that computes G 0 with probability of success at least 1 -1/2 P using a number of bit operations softly linear in

d 2 h + (d ω η + η ω ) log(h) + η 2 c.
Outlook. Inspired by Lazard [START_REF] Lazard | Ideal bases and primary decomposition: case of two variables[END_REF], we prove in Section 2 that the Hermite Normal form of an "extended Sylvester matrix" built from f 1 , . . . , f t gives the coefficients of a detaching basis of the ideal I they generate. We also present a variant of this result, where replacing the Hermite normal form by the Howell normal form yields a Gröbner basis of a localization of I. We use these results in two manners: to compute the initial Gröbner basis modulo p, prior to entering Newton iteration, and to obtain bit-size bounds for the output (over Q) and quantify bad choices of the prime p. The underlying algorithms for the above theorems are in Section 5.

Using matrix normal forms

In this section, we assume I = ⟨f 1 , . . . , f t ⟩ ⊂ K[x, y], for t ≥ 2, and we derive the lexicographic Gröbner basis of I, or its primary component at the origin, from either Hermite or Howell normal forms of matrices over K[x], for an arbitrary field K. These results are direct extensions of previous work of Lazard's [START_REF] Lazard | Ideal bases and primary decomposition: case of two variables[END_REF], who used Hermite forms in the case t = 2.

In what follows, for a subset S ⊂ K[x, y] and n ≥ 0, we let S <(.,n) be the subset of all f in S with deg y (f ) < n; notation such as S ≤(.,n) is defined similarly. In particular, if S is an ideal of K[x, y], S <(.,n) is a free

K[x]-module of rank at most n. For S = K[x, y] itself, K[x, y] <(.,n) is a free K[x]-module of rank n, equal to 0≤i<n K[x]y i .
For such an n, we also let

π n denote the K[x]-module isomorphism K[x, y] <(.,n) → K[x] n , which maps f 0 + • • • + f n-1 y n-1 to the vector [f n-1 • • • f 0 ] ⊤ .

Detaching bases

Let I be an ideal in K[x, y] and let G = (g 0 , . . . , g s ) be its reduced minimal Gröbner basis for the lexicographic order induced by y ≻ x, listed in decreasing order; we write n i = deg y (g i ) for all i (so these exponents are decreasing). We define polynomials A 0 , A 1 , . . . as follows: for 0 ≤ i < n s , A i = 0, and if there exists k in {0, . . . , s} such that n k = i, A i = g k ; otherwise, A i is obtained by starting from yA i-1 , and reducing all its terms of y-degree less than i by G.

For example, if I has a Gröbner basis of the form (y -f (x), g(x)), the polynomials A i are given by A 0 = g, A 1 = y -f and for i ≥ 2, A i = y i -(f i mod g) (see [START_REF] Ayoub | On constructing bases for ideals in polynomial rings over the integers[END_REF] for a previous discussion).

Lemma 1. For i ≥ n s , deg y (A i ) = i.
Proof. This is true for i of the form n k . For i in n k , . . . , n k-1 -1, we proceed by induction, with the remark above establishing the base case (for k = 0, we consider all i ≥ n 0 ). Assume deg y (A i-1 ) = i -1, so that deg y (yA i-1 ) = i. Because we use the lexicographic order x ≺ y, the reduction of the terms of y-degree less than i in yA i-1 does not introduce terms of y-degree i or more.

Lemma 2. For n ≥ n s , the K[x]-module I ≤(.,n) is free of rank n -n s + 1, with basis A ns , . . . , A n .

Proof. The polynomials A ns , . . . , A n are all nonzero, with pairwise distinct y-degrees, so they are K[x]-linearly independent. They all belong to I ≤(.,n) , so it remains to prove that they generate it as a K[x]-module. This is done by induction on n ≥ n s . Take f in I ≤(.,n) , and write it as

f = f n y n + g, with f n in K[x] and g in K[x, y] ≤(.,n-1) . Let h ∈ K[x] be the polynomial coefficient of y n in A n , so that A n = h n y n + B n , with B n in K[x, y] ≤(.,n-1) . Write the Euclidean division f n = qh n + r in K[x], with deg x (r) < deg x (h n ): f = (qh n + r)y n + g = qh n y n + ry n + g = qA n -qB n + ry n + g.
The polynomial -qB n + ry n + g is in I, so its normal form modulo G is zero. The terms -qB n + g have y-degree less than n, so their normal form has y-degree less than n as well; since ry n is already reduced modulo G, it must be zero. It follows that f = qA n + g -qB n , with g -qB n in I ≤(.,n-1) . If n = n s , this latter polynomial must vanish, proving the base case. Else, by induction assumption, it is a K[x]-linear combination of A ns , . . . , A n-1 .

For n ≥ n 0 , the detaching basis of I in degree n is the sequence (A ns , . . . , A n ). Because we take n ≥ n 0 , this is (in general) a non-minimal Gröbner basis of I, and we can recover G from it by discarding redundant entries (that is, all polynomials whose leading term is a multiple of another leading term).

Using Hermite normal forms

Given F = (f 1 , . . . , f t ) in K[x, y], we prove that the Hermite normal form of a certain Sylvester-like matrix associated to them gives a lexicographic detaching basis of the ideal I they generate. In [START_REF] Lazard | Ideal bases and primary decomposition: case of two variables[END_REF], Lazard covered the case t = 2, under an assumption on the leading coefficients (in y) of the f i 's.

We extend his work (in a direct manner) to situations where such assumptions do not hold. First, to polynomials F = (f 1 , . . . , f t ) in K[x, y], we associate an integer ∆(F), defined as follows. Definition 1. Let F = (f 1 , . . . , f t ) ∈ K[x, y] t , with (A ns , . . . , A n 0 ) as detaching basis in degree n 0 , with n 0 and n s the maximal, resp. minimal y-degree of the polynomials in the lexicographic Gröbner basis of ⟨f 1 , . . . , f t ⟩, for the order x ≺ y.

We let ∆(F) be the minimal integer ∆ such that for i = n s , . . . , n 0 , there exist w i,1 , . . . , w i,t in K[x, y] t , all of y-degree less than ∆, and such that

A i = w i,1 f 1 + • • • + w i,t f t .
The following proposition gives the basic result using this definition, allowing us to extract a detaching basis from a Hermite form computation. We use column operations, with Hermite normal forms being lower triangular. The first nonzero entry in a nonzero column is called its pivot, its index being called the pivot index. Pivots are monic in x. Proposition 1. Let F = (f 1 , . . . , f t ) be in K[x, y], for t ≥ 2, of y-degree at most d y , and assume that they generate an ideal I = ⟨f 1 , . . . , f t ⟩ of dimension zero. For i = 1, . . . , t, write

f i = f i,0 + • • • + f i,dy y dy , with all f i,j in K[x].
For D ≥ ∆(F), let c 1 , . . . , c K be the nonzero columns of the Hermite normal form

H of S = [S 1 • • • S t ] ∈ K[x] (dy+D)×tD
, where

S i =         f i,dy . . . . . . f i,0 f i,dy . . . . . . f i,0         ∈ K[x] (dy+D)×D .
Then, there exists K ′ ≤ K such that π -1 dy+D (c K ′ ) is monic in y; with K ′ the largest such integer, π -1 dy+D (c K ), . . . , π -1 dy+D (c K ′ ) is a detaching basis of I.

Thus, while we do not know the y-degrees n i of the elements in the Gröbner basis of I, as long as D ≥ ∆(F), it is enough to consider the last nonzero columns of H, stopping when we find (through π -1 dy+D ) a polynomial that is monic in y.

Proof. Let D ≥ ∆(F) be as in the proposition. Let us index the columns of each block S i by y D-1 , . . . , y, 1, and its rows by y dy+D-1 , . . . , y, 1. Then, S i is the matrix of the map K[x, y] <(.,D) → K[x, y] <(.,dy+D) given by w i → w i f i . The matrix S itself maps a vector (w 1 , . . . , w t ), with all entries of y-degree less than D, to t i=1 w i f i ∈ I <(.,dy+D) . Let G = (g 0 , . . . , g s ) be the lexicographic Gröbner basis of I = ⟨f 1 , . . . , f t ⟩, f i ≻ f i+1 , with deg y (g i ) = n i for all i. Since we assume that I has dimension zero, we have n s = 0, and g 0 is monic in y.

Let A 0 , . . . , A n 0 be the detaching basis of I in degree n 0 . We denote by c 1 , . . . , c K the nonzero columns of the Hermite form H of S, and we let H i = π -1 dy+D (c i ), for i = 0, . . . , n 0 . We will prove that A i = H K-i for i = 0, . . . , n 0 . Since g 0 is the only polynomial in A 0 , . . . , A n 0 which is monic in y, this will establish the proposition, with K ′ = K -n 0 .

Since both A i and H K-i are in I, to prove that they are equal, it is enough to prove that for all i, A i -H K-i is reduced with respect to the Gröbner basis G of I. Because D ≥ ∆(F), we deduce that A 0 , . . . , A n 0 are in the column span of S. Since they have respective y-degrees 0, . . . , n 0 , we see that deg y (H K-i ) = deg y (A i ) = i for all i = 0, . . . , n 0 . In addition, for all such i, we can write A i = i j=0 a i,j H K-j , for some a i,j in K[x]. However, Lemma 2 shows that for the same index i, we can write

H K-i = i j=0 b i,j A j , for some b i,j in K[x].
Because both A i and H K-i have leading y-coefficients that are monic in x, it follows that b i,i = a i,i = 1 for all i. This proves that A i and H K-i have the same coefficient of y-degree i (call it

M i ∈ K[x]
), and thus that A i -H K-i has y-degree less than i.

By the definition of a detaching basis, all terms of y-degree less than i in A i are reduced with respect to G. On the other hand, by the property of Hermite forms, for j < i, the coefficient of y-degree j in H K-i is reduced with respect to M j . Since we saw that M j is also the coefficient of y j in A j , this proves that all terms of y-degree less than i in H K-i are reduced with respect to A 0 , . . . , A i-1 , and thus with respect to G. Altogether, A i -H K-i itself is reduced with respect to G, which is what we set out to prove.

We call HermiteGroebnerBasis(F, D) a procedure that takes as input F = (f 1 , . . . , f t ) and D, and returns the lexicographic Gröbner basis of I = ⟨f 1 , . . . , f t ⟩ obtained by computing the Hermite normal form of S as above, extracting the Gröbner basis of I from its detaching basis. Here, we take for d y the maximum degree of the f i 's, and we assume that we have D ≥ ∆(F) and D ≥ d y .

The assumption that the ideal I has dimension zero implies that it contains a non-zero polynomial in K[x]; as a result, its detaching basis has entries of y-degrees 0, 1, . . . , so that the Hermite form of S is lower triangular with d y + D non-zero diagonal entries. In other words, S has rank d y + D (seen as a matrix over K(x)).

If t = 2 and D = d y , this matrix is square, but in general, it may have more columns than rows (recall that we assume D ≥ d y ). Using the algorithm of [START_REF] Labahn | Rank-sensitive computation of the rank profile of a polynomial matrix[END_REF], we can permute the columns of S to find a (d y + D) × tD matrix S ′ whose leading (d y + D) × (d y + D) minor is nonzero; this takes O˜(tD ω d) operations in K, with d the maximum degree of the f i 's. Let us define the tD × tD square matrix

S sq = S ′ 0 (t-1)D-dy,dy+D I (t-1)D-dy,(t-1)D-dy (1)
together with its Hermite form H sq ; the first d y + D rows of it give us the Hermite form H of S. The Hermite form of S sq is computed in O˜(t ω D ω d) operations in K [START_REF] Labahn | Fast, deterministic computation of the Hermite normal form and determinant of a polynomial matrix[END_REF]. This gives the overall cost of computing the lexicographic Gröbner basis of I, assuming an upper bound on ∆(F) is known. To our knowledge, not much exists in the literature on complete cost analysis for bivariate ideals, apart from Buchberger's algorithm, with cost 3 2 (t + 2(d + 2) 2 ) 4 [START_REF] Buchberger | A note on the complexity of constructing gröbner-bases[END_REF].

The following proposition gives various bounds on ∆(F), whose strength depends on the assumptions we make on F. The first one is a direct extension of Lazard's [32, Lemma 7], and is linear in the y-degree of the input. The others are based on results from [START_REF] Kollar | Sharp effective nullstellensatz[END_REF][START_REF] Dickenstein | The membership problem for unmixed polynomial ideals is solvable in single exponential time[END_REF], which involve total degree considerations.

Proposition 2. Let F = (f 1 , . . . , f t ) be in K[x, y] of degree at most d ≥ 1,
and y-degree at most d y , and let

I = ⟨f 1 , . . . , f t ⟩ ⊂ K[x, y]. Set d ′ = max(d, 3
). Then:

• if there exists i in {1, . . . , t} such that the coefficient of y d in f i is a nonzero constant, ∆(F) ≤ ∆ 1 (d y ) = d y • if t = 2 and I has finitely many zeros over K, ∆(F) ≤ ∆ 2 (d) = 2d ′ 2 + d ′ ∈ O(d 2 )
• if I has finitely many zeros over

K, ∆(F) ≤ ∆ 3 (d) = 16d ′ 4 + 2d ′ 2 + 2d ′ ∈ O(d 4 ) First item.
In what follows, without loss of generality, we assume that the coefficient of y dy in f t is 1. We prove a slightly more general claim: any polynomial f in I <(.,2dy) can be written as ,dy) . This is enough to conclude since all entries A ns , . . . , A n 0 in the detaching basis of I in degree n 0 have y-degree at most d y ≤ 2d y -1 (because we use a lexicographic order with x ≺ y).

f = w 1 f 1 + • • • + w t f t , with all w i in K[x, y] <(.
Let thus f be given in I <(.,2dy) . There is a family w = (w 1 , . . . , w t ) in

K[x, y] such that f = t i=1 w i f i , since f is in I. For such a family w, we define S w = {i | deg y (w i ) ≥ d y }.
For any w such that S w is not empty, we further set ν w = min(S w ) ∈ {1, . . . , t}, and we let ν be the maximal value of these ν w 's. It is well-defined, since there is a vector w for which S w is not empty (we can replace (w t-1 , w t ) by (w t-1 + gf t , w t -gf t-1 ) for any g in K[x, y]).

Let w be such that ν = ν w . We claim that S w ̸ = {t}: otherwise we would have deg y (w t f t ) ≥ 2d y , while deg y (w i f i ) < 2d y for all other i's; this would contradict the assumption deg y (f ) < 2d y . This shows that ν < t. Let us further refine our choice of w, by taking it such that, among all those vectors for which S w is not empty and ν w = ν, the y-degree of w ν is minimal. Let us then write e = deg y (w ν ) (so that e ≥ d y ) and let c ∈ K[x] be the coefficient of y e in w ν . We can use it to rewrite f as

f = t i=1 w i f i + cy e-dy f ν f t -cy e-dy f t f ν .

If we set

w ′ i =     
w ν -cy e-dy f t when i = ν; w t + cy e-dy f ν when i = t; w i otherwise,

we still have f = t i=1 w ′ i f i . By construction, deg y (w ′ i ) = deg y (w i ) < d y for all i < ν, so none of 1, . . . , ν -1 is in S w ′ . If ν is in S w ′ , then the inequality deg y (w ′ ν ) < deg y (w ν
) contradicts the choice of w, so that ν is not in S w ′ . This shows that S w ′ is empty, since otherwise its minimum element would be greater than ν.

For the second and third items, we use results from [START_REF] Dickenstein | The membership problem for unmixed polynomial ideals is solvable in single exponential time[END_REF], for which we need total degree bounds on the input polynomials F = (f 1 , . . . , f t ) and the elements A 0 , . . . , A n 0 in the detaching basis (here n s = 0 since I having finitely many solutions implies that it contains a nonzero polynomial in K[x]). For the inputs f i , we have the degree bound deg(f i ) ≤ d ≤ d ′ . For the A i 's, we have the bounds deg x (A i ) ≤ d 2 (by Bézout's theorem) and deg y (A i ) ≤ d for i ≤ n 0 , so their total degree is at most

D = d ′2 + d ′ .
Second item. When t = 2 and I has dimension zero (that is, has a finite, nonzero number of solutions in K), f 1 , f 2 are in complete intersection, so that we have

A i = w i,1 f 1 + w i,2 f 2 , with deg y (w i,j ) ≤ D + d ′ 2
for all i, j, by Theorem 5.1 in [START_REF] Dickenstein | The membership problem for unmixed polynomial ideals is solvable in single exponential time[END_REF]. Overall, the resulting degree bound is 2d ′ 2 + d ′ . If we assume that I = K[x, y], we know that there are [START_REF] Kollar | Sharp effective nullstellensatz[END_REF]. Multiplying this by A j , for j ≤ n 0 , we obtain the expression (g

g 1 , g 2 in K[x, y] such that g 1 f 1 + g 2 f 2 = 1, with deg(g i ) ≤ d ′ 2
1 A j )f 1 +(g 2 A j )f 2 = A j , with deg y (g i A j ) ≤ d ′ 2 +d in this case.

Third item. [15, Corollary 3.4] gives equalities

A i = w i,1 f 1 + • • • + w i,t f t , with deg y (w i,j ) ≤ D + 16d ′ 4 + d ′ 2 + d ′ for all i, j.

Using the Howell form

We now investigate how using another matrix normal form, the Howell form [START_REF] Howell | Spans in the module Z s m . Linear and Multilinear Algebra[END_REF], yields information about certain primary components of an ideal I as above. Howell forms are defined for matrices with entries in a principal ideal ring A; below, we will take A = K[x]/x k , for an integer k. Again, we consider column operations then an n × m matrix H over A = K[x]/x k is in Howell normal form if the following (taken from [START_REF] Storjohann | Algorithms for matrix canonical forms[END_REF]Chapter 4]) hold:

1. let r ≤ m be the number of nonzero columns in H; then these nonzero columns have indices 1, . . . , r

2. H is in lower echelon form: for i = 1, . . . , r, let j i ∈ {1, . . . , n} be the index of the first nonzero entry in the ith column; then,

j 1 < • • • < j r 3.
all pivots H j i ,i , for i = 1, . . . , r, are of the form x c i 4. for i = 1, . . . , r and k = 1, . . . , i -1, H j i ,k is reduced modulo H j i ,i

5. for i = 0, . . . , r, any column in the column span of H with at least j i leading zeros is an A-linear combination of columns of indices i + 1, . . . , r (here, we set j 0 = 0)

For any M in A n×m , there is a unique H in Howell normal form in A n×m , and a not necessarily unique invertible matrix U in A m×m such that H = M U . The matrix H is the Howell normal form of M .

Given f 1 , . . . , f t as before, we are interested here in computing the lexicographic Gröbner basis of J = ⟨f 1 , . . . , f t , x k ⟩, for a given integer k. In particular, if (0, 0) is in V (f 1 , . . . , f t ), and no other point (0, β) is, for β ̸ = 0, J is the ⟨x, y⟩-primary component of I = ⟨f 1 , . . . , f t ⟩, if k is large enough.

The following proposition shows how to reduce this computation to a Howell normal form calculation. In what follows, the canonical lift of an element in A = K[x]/x k to K[x] is its unique preimage of degree less than k; this carries over to vectors and matrices (and in particular to the output of the Howell form computation). Contrary to what happens for Hermite forms, there is no guarantee that the polynomials extracted from the Howell form are a detaching basis, as we may be missing the first polynomial (that belongs to K[x]) and its multiples. The proposition below restores this by considering a few extra columns, if needed. Proposition 3. Let f 1 , . . . , f t be in K[x, y], for t ≥ 2, of y-degree at most d y , and assume that they generate an ideal of dimension zero. Let k be a positive integer and

A = K[x]/x k .
For D ≥ ∆(f 1 , . . . , f t , x k ), let B ∈ A (dy+D)×tD be the Howell normal form of S = S mod x k , with S as in Proposition 1, and let B lift be its canonical lift to K[x] (dy+D)×tD .

Let h 1 , . . . , h L be the nonzero columns of B lift , and let r ∈ {1, . . . , d y +D} be the pivot index of h L . Set L ′′ = L+d y +D -r and, for i = L+1, . . . , L ′′ let

h i = [0 • • • 0 x k 0 • • • 0] ⊤ , with x k at index r + i -L ∈ {r + 1, . . . , d y + D}.
Then, there exists L ′ ≤ L such that π -1 dy+D (h L ′ ) is monic in y; with L ′ be the largest such integer, π -1 dy+D (h L ′′ ), . . . , π -1 dy+D (h L ′ ) is a detaching basis of ⟨f 1 , . . . , f t , x k ⟩.

Proof. Let Γ = (Γ 0 , . . . , Γ σ ) be the lexicographic Gröbner basis of J = ⟨f 1 , . . . , f t , x k ⟩, listed in decreasing order, with Γ i of y-degree ν i for all i; since x k is in J, ν σ = 0. Then, let C 0 , . . . , C ν 0 be the detaching basis of J in degree ν 0 , with deg y (C i ) = i for all i.

We know that the first polynomials in the detaching basis are of the form

C 0 = x ℓ , C 1 = yx ℓ , . . . , C ν σ-1 -1 = y ν σ-1 -1 x ℓ , for some ℓ ≤ k. If ℓ = k, they all vanish modulo x k , but the next polynomial C ν σ-1 does not. If ℓ < k,
none of them vanishes modulo x k . Thus, we define ρ = ν σ-1 in the former case and ρ = 0 in the latter.

Let further D ≥ ∆(f 1 , . . . , f t , x k ) be as in the proposition. If we consider the extended Sylvester matrix T ∈ K[x] (dy+D)×(t+1)D built from f 1 , . . . , f t , x k , the assumption on D shows that each π dy+D (C i ) is in the column span of T . For i = 0, . . . , ν 0 , we let v i be the column vector π dy+D (C i ) mod x k ∈ A dy+D ; the previous paragraph shows that the nonzero vectors v i are precisely v ρ , . . . , v ν 0 . By reduction modulo x k of the membership relations above, we see that v ρ , . . . , v ν 0 are in the A-span of the columns of S.

Lazard's structure theorem [START_REF] Lazard | Ideal bases and primary decomposition: case of two variables[END_REF]Theorem 1] shows that every polynomial Γ j in the reduced Gröbner basis of J is of the form Γ j = x m j γ j , with γ j monic in y and m j ≤ ℓ (the inequality is strict, except for j = 0). It follows that for i = ρ, . . . , ν 0 , the pivot in v i is also a power of x, at index d y + D -i (precisely, it is x m j , for j the largest integer such that ν j ≤ i).

Let η 1 , . . . , η L be the nonzero columns in the Howell form B of S. By definition of the Howell form, the former observation implies that for i = ρ, . . . , ν 0 , v i is in the A-span of those η j 's starting with at least d y + D -i -1 zeros. For such an i, since the entry at index d y + D -i in v i is nonzero, there exists (exactly) one η j with pivot index d y + D -i. We now prove that the pivot in η L is at index d y + D -ρ. Recall that we write h 1 , . . . , h L for the canonical lifts of η 1 , . . . , η L to K[x] dy+D ; in particular, the pivot index r of h L , as defined in the proposition, is also the pivot index of η L , so our claim is r = d y + D -ρ. Suppose that the pivot in η L is at an index different from d y + D -ρ. By the previous discussion, it can only lie at a larger index, say m > d y + D -ρ; this may happen only if ρ > 0, in which case we saw that ρ = ν σ-1 = deg y (Γ σ-1 ) and Γ σ = x k .

Let H 1 , . . . , H L be obtained by applying π -1 dy+D to h 1 , . . . , h L . It follows that H L has y-degree d y + D -m < ρ = deg y (Γ σ-1 ), and x-degree less than k = deg x (Γ σ ). Thus, H L is reduced with respect to the Gröbner basis Γ of J. On the other hand, because η L is in the column span of S, its canonical lift h L is in the column space of S, up to the addition of a vector with entries in x k K[x]. In other words, H L is in J, so that H L must be zero, a contradiction.

Thus, the pivot index of η L is exactly d y + D -ρ, that is, the same as that of v ρ . Our previous discussion on the pivots in the vectors η i then implies that for i = ρ, . . . , ν 0 , the pivot index of η L+ρ-i is d y + D -i, that is, the same as that of v i . This implies that

v i = i j=ρ α i,j η L+ρ-j , (2)
for some coefficients α i,j in A = K[x]/x k . On the other hand, all polynomials H L , . . . , H L+ρ-ν 0 are in J (by the argument we used for H L ). By Lemma 2, we deduce that for i = ρ, . . . , ν 0 , H L+ρ-i can be written as H L+ρ-i = i j=ρ β i,j C j , for some coefficients β i,j in K[x]. Applying π dy+D and reducing modulo x k , this gives

η L+ρ-i = i j=ρ βi,j v j , (3)
with βi,j = β i,j mod x k for all i, j. We know that the pivots of both v i and η L+ρ-i are powers of x (the latter, by the properties of the Howell form), so Eq. ( 2) and Eq. ( 3) show that the pivots in v i and η L+ρ-i are the same, for i = ρ, . . . , ν 0 . Back in K[x, y], we deduce that C i and H L+ρ-i have the same coefficient in y i , for i = ρ, . . . , ν 0 . As in the proof of Proposition 1, we deduce that we actually have C i = H L+ρ-i for i = ρ, . . . , ν 0 : we observe that their terms of y-degree less than i are reduced with respect to Γ; it follows that C i -H L+ρ-i is both in J and reduced with respect to its lexicographic Gröbner basis, so it vanishes.

Taking i = ν 0 , we deduce in particular that H L+ρ-ν 0 is monic in y (and no H i of larger index has this property), so the index L ′ defined in the proposition is L ′ = L + ρ -ν 0 ; the corresponding polynomials are C ν 0 , . . . , C ρ . Since we saw that r = d y + D -ρ, the integer L ′′ in the proposition is L ′′ = L + ρ, and through π -1 dy+D , the columns h L+1 , . . . , h L+ρ become y ρ-1 x k , . . . , x k (there is no such column if ρ = 0). These are precisely the polynomials C ρ-1 , . . . , C 0 that were missing if ρ > 0.

We call HowellGroebnerBasis(F, k, D) a procedure that takes as input F = (f 1 , . . . , f t ), k and D, and returns the lexicographic Gröbner basis of ⟨f 1 , . . . , f t , x k ⟩ obtained from the Howell form of S, taking for d y the maximum of the degrees of f 1 , . . . , f t , and choosing for D the integer prescribed by Proposition 2. In this case, there is no need to make S square: the algorithm of [START_REF] Storjohann | Algorithms for matrix canonical forms[END_REF]Chapter 4] computes its Howell form using O˜(tD ω k) operations in K.

The main application we will make of Howell form computation is to obtain the Gröbner basis of the ⟨x, y⟩-primary component of an ideal such as I = ⟨f 1 , . . . , f t ⟩. In order to do so, we will assume that we are in "nice" coordinates, in the sense that the projection on the first factor V (F) → K is one-to-one. Lemma 3. Let F = (f 1 , . . . , f t ) be in K[x, y], and suppose that the projection on the first factor V (F) → K is one-to-one. Let further J be the ⟨x, y⟩primary component of I = ⟨f 1 , . . . , f t ⟩, with m the smallest integer such that x m is in J. Then: the smallest power of x in the ideal H = ⟨f 1 , . . . , f t , x k ⟩ is x min(m,k) , and for k ≥ m, H = J.

Proof. First, we establish that J = ⟨f 1 , . . . , f t , x m ⟩. For one direction, all f i 's, as well as x m , are in J by definition. Conversely, the assumption on V (F) implies that we can write ⟨f 1 , . . . , f t ⟩ = JJ ′ , with J ′ having no solution above x = 0; in particular, there exist polynomials u, v with ux m + v = 1 and v in J ′ . From this, we get J = (ux m + v)J, and every element in ux m J is a multiple of x m , while every element in vJ is in ⟨f 1 , . . . , f t ⟩.

Suppose k ≥ m. We now have polynomials u ′ , v ′ with u ′ x k-m +v ′ = 1 and v ′ in J ′ . Multiplying by x m shows that x m is in the ideal H = ⟨f 1 , . . . , f t , x k ⟩, so that H = J (this proves the last claim in the lemma). In this case, the smallest power of x in H is thus x m .

Suppose k ≤ m. In this case, we prove that the minimal power of x in H = ⟨f 1 , . . . , f t , x k ⟩ is x k . First, note that in this case, H = ⟨f 1 , . . . , f t , x m , x k ⟩ = J + ⟨x k ⟩, and let x e be the minimum power of x in H; suppose e < k, so that e < m. It follows that x e is the normal form of a polynomial of the form f x k , modulo the Gröbner basis G of J. However, Lazard's structure theorem [32, Theorem 1] implies that through reduction modulo such a Gröbner basis, no term of x-degree less than k can appear; a contradiction. This allows us to design an algorithm GroebnerBasisAtZero that computes the Gröbner basis of J (under the position assumption in the lemma), even though we do not know m in advance: we call HowellGroebnerBasis with inputs the polynomials (f 1 , . . . , f t , x k ), for k = 2 i , with i = 0, 1, . . . , until the output does not contain x k . Indeed, the lemma shows that if x k is in the Gröbner basis of H = ⟨f 1 , . . . , f t , x k ⟩, we have k ≤ m, while if it is not, we have reached k > m, and the output is the Gröbner basis of J.

Altogether, we do O(log(m)) calls to HowellGroebnerBasis, with k ≤ 2m. With d the maximum degree of f 1 , . . . , f t , the runtime is O˜(tD ω m) operations in K, with D in {∆ 1 (d y ), ∆ 2 (d), ∆ 3 (d)}, depending on our assumptions on f 1 , . . . , f t (recall that d y and d are the maximum y-degree, resp. degree, of the input).

Coefficient size and bad reductions

Our goal now is to give height bounds on the elements in the lexicographic Gröbner basis of some polynomials F = (f 1 , . . . , f t ), working specifically over K = Q. In this section, we assume that the input polynomials have integer coefficients.

The height u ∈ Z -{0} is simply log(|u|). The key quantity H(F), together with a nonzero integer β F ∈ Z, are defined as follows.

Definition 2. Consider polynomials F = (f 1 , . . . , f t ) in Z[x, y], let I be the ideal they generate in Q[x, y], with lexicographic Gröbner basis G = (g 0 , . . . , g s ). We define H(F) as the smallest integer such that there exists β F nonzero in Z for which we have:

• the polynomials β F g 0 , . . . , β F g s are in Z[x, y]
• all coefficients of β F g 0 , . . . , β F g s (which include in particular β F itself ) have height at most H(F)

• for any prime p in Z, if p does not divide β F , G mod p is the lexicographic Gröbner basis of ⟨f 1 mod p, . . . , f t mod p⟩ in F p [x, y].

In order to give upper bounds on H(F), we introduce two functions B(n, d, h) andC(t, d, D, h). The first one is defined by whereas C(t,d,D,h) is defined by

B(n, d, h) = (N + 1)h + N log(N ) + log(n(d + 1)), with N = n 2 d -nd + n,
C(t, d, D, h) = B(tD, d, h) + h + log(2). In particular, B(n, d, h) is in O˜(n 2 dh) and C(t, d, D, h) is in O˜(t 2 D 2 dh). Proposition 4. Let F = (f 1 , . . . , f t ) be in Z[x, y], for t ≥ 2, such that the ideal I = ⟨f 1 , . . . , f t ⟩ ⊂ Q[x,
y] has dimension zero. Suppose that all f i 's have y-degree at most d y , degree at most d, and coefficients of height at most h.

(i) if there exists i in {1, . . . , t} such that the coefficient of

y dy in f i is a nonzero constant, H(F) ≤ C(t, d, ∆ 1 (d y ), h) ∈ O˜(t 2 d 3 h) (ii) if t = 2, H(F) ≤ C(2, d, ∆ 2 (d), h) ∈ O˜(d 5 h) (iii) in general, H(F) ≤ C(t, d, ∆ 3 (d), h) ∈ O˜(t 2 d 9 h).
The proposition will follow from height bounds for Hermite forms of matrices due to Storjohann, with the results in the previous section. We do not have a direct equivalent for Howell forms (we are not aware of previous work about height bounds or primes of bad reduction in this context): if we are interested in the primary component of I at the origin, we may apply the results above to the polynomials f 1 , . . . , f t , x k , for a large enough k.

To our knowledge, no comparable bounds were given in this setting. Several previous results discussed the case of radical ideal with finitely many solutions. If their Gröbner basis G is a triangular set, the results in [START_REF] Dahan | Sharp estimates for triangular sets[END_REF] show that the polynomials in G have coefficients with numerator and denominator of height O˜(d 3 h + d 4 ). Our result does not feature the term d 4 , but this might be due to the proof techniques of [START_REF] Dahan | Sharp estimates for triangular sets[END_REF], which are not limited to systems in two variables. If we keep the radicality assumption, but allow arbitrary leading terms, the best previous bound we are aware of is O˜(d 7 h + d 8 ), from [START_REF] Dahan | Size of coefficients of lexicographical Groöbner bases: The zero-dimensional, radical and bivariate case[END_REF].

We start the proof with a result due to Storjohann.

Proposition 5 ([41, Section 6.2]). Let A be in Z[x] n×n , with nonzero determinant and entries of degree at most d > 0 and height at most h. Let further H be the Hermite normal form of A. Then, there exists α nonzero in Z such that all entries of αH are in Z[x], α and the coefficients of all entries of αH have height at most B(n, d, h), and for any prime p, if p does not divide α, then H mod p is the Hermite normal form of A mod p in

F p [x] n×n .
Let then f 1 , . . . , f t be as in Proposition 4. First, we define integers γ and D through the following case discussion. If we are in case (i), we know that at least one of the f i 's has a coefficient of y-degree d y in Z -{0}; let γ be such a coefficient. We let D = ∆ 1 (d y ) from Proposition 2. In case (ii) or (iii), we let γ = 1, and we take respectively D = ∆ 2 (d) or D = ∆ 3 (d), with notation as above. In any case, we know that ∆(F) ≤ D, so we can apply Proposition 1; it shows that we can recover the (minimal, reduced) lexicographic Gröbner basis of I = ⟨f 1 , . . . , f t ⟩ from the columns of the Hermite form of the Sylvester-like matrix S defined in that proposition.

As in the previous section, there is a (d y + D) × tD matrix S ′ obtained by permuting the columns of S whose leading (d y + D) × (d y + D) minor is nonzero. Consider again the tD × tD square matrix S sq of Eq. ( 1) and its Hermite form H sq ; the first d y + D rows of H sq are the Hermite form H of S.

Since S sq has nonzero determinant, we let α be the nonzero integer associated to it by Proposition 5, and set β = αγ. Then, all entries of βH sq , and thus of βH, are in Z[x], the latter having coefficients of height at most C(t, d, D, h). By Proposition 1, these bounds apply in particular to the Gröbner basis (g 0 , . . . , g s ) of I.

Suppose that p is a prime that does not divide β. Because p does not divide α, Proposition 5 shows that Hsq = H sq mod p is the Hermite normal form of Ssq = S sq mod p. Considering only the first tD rows, we see that H = H mod p is the Hermite normal form of S = S mod p. Now, let us prove that we still have ∆( F) ≤ D.

• If we are in case (i), since p does not divide γ, at least one of the polynomials fi = f i mod p has its coefficient of y-degree d y a nonzero constant in F p .

Since all fi 's have y-degree at most d y , we deduce ∆( F) = d y in this case (first item of Proposition 2)

• If we are in case (ii) or (iii), the discussion above shows that ḡ0 and ḡs are in the ideal ⟨ f1 , . . . , ft ⟩, so that this ideal admits finitely many solutions in an algebraic closure of F p . Using the second and third items of Proposition 2 gives our claim.

We can then apply Proposition 1 to F = ( f1 , . . . , ft ), and deduce that the columns of the Hermite form of S give a detaching basis, and in particular the lexicographic Gröbner basis of ⟨ f1 , . . . , ft ⟩. This proves the proposition.

Applying changes of coordinates

Now, we quantify changes of coordinates that ensure desirable properties. We write γ for a 2 × 2 matrix γ = [γ i,j ] 1≤i,j≤2 with entries in Q, and we identity

M 2 (Q) with Q 4 through γ → [γ 1,1 , γ 1,2 , γ 2,1 , γ 2,2 ]. For γ in GL 2 (Q) as above and f in Q[x, y], we write f γ = f (γ 1,1 x + γ 2,1 y, γ 1,2 x + γ 2,2 y).
For F = (f 1 , . . . , f t ) as in the previous sections, the best degree and height bounds ∆(F) and H(F) apply when the input equations have a particular property: at least one f i has a term of maximal degree that involves y only. Geometrically, this means that the curve V (f i ) ⊂ Q 2 has no vertical asymptote; we also say that it is in Noether position. The following lemma is straightforward.

Lemma 4. Take f in Q[x, y] of degree d. Then there exists a hypersurface Y 1 ⊂ Q 4 of degree at most d such that if γ is in Q 4 -Y 1 , the coefficient of y d in f γ is nonzero.
Another favorable situation, illustrated when we dealt with Howell forms, occurs when the projection V (F) → Q given by (α, β) → α is one-to-one. Again, the proof is standard (see e.g. [START_REF] Rouillier | Solving zero-dimensional systems through the Rational Univariate Representation[END_REF]), once we see that V (F) has cardinality at most d 2 . Lemma 5. Let F = (f 1 , . . . , f t ) be in Q[x, y] of degrees at most d, and suppose that V (F) is finite. Then there exists a hypersurface Y 2 ⊂ Q 4 of degree at most d 4 such that if γ is invertible and in Q 4 -Y 2 , the projection on the first factor V (f γ ) → Q is one-to-one. and we compute either the Gröbner basis G = (g 0 , . . . , g s ) of I = ⟨f 1 , . . . , f t ⟩, or the Gröbner basis G 0 = (g 0 0 , . . . , g 0 r ) of the ⟨x, y⟩-primary component of I using p-adic approximation, for a prime p. In what follows, we give details for the computation of G; we will mention what modifications are needed if we want to compute G 0 .

The algorithm is randomized; it takes a parameter P ≥ 1, our goal being to obtain the correct output with probability at least 1 -1/2 P . Throughout, we assume that f 1 has maximum degree among the f i 's (we write d = deg(f 1 )) and that I has dimension zero. Let δ = deg(I) = dim Q Q[x, y]/I, δ ≤ d 2 and let b be the maximum height of the numerators and denominators of the coefficients in G. Each polynomial in G has at most δ + 1 coefficients, so the total bit-size of the output is O(sδb).

Overview

We start by presenting the main steps of the algorithm, leaving out some details of the analysis for the next subsection. Runtimes are given in terms of bit operations; here, we use the fact that operations (+, ×) modulo a positive integer M take O˜(log(M )) bit operations, as does inversion modulo M if M is prime [START_REF] Zur Gathen | Modern Computer Algebra[END_REF].

Introducing a change of coordinates. We first choose a change of variables γ with coefficients in Z. Applying it to the input equations F gives polynomials H = (h 1 , . . . , h t ), which we do not need to compute explicitly (as they may have large height). We let B = (B 0 , . . . , B σ ) be the lexicographic Gröbner basis of these polynomials in Q[x, y] (as with H, we do not compute it explicitly).

We assume that γ satisfies the assumptions of Lemmas 4 and 5, so that their conclusions hold.

Computing Gröbner bases modulo p. Next, we choose two primes p, p ′ , and compute the Gröbner bases B p of (H mod p), and B p ′ of (H mod p ′ ). We assume that neither p nor p ′ divides the integers β F and β H from Definition 2 applied to F and H, respectively. In particular, all denominators in B are invertible modulo p and p ′ , and B p = B mod p and B p ′ = H mod p ′ .

To compute B p and B p ′ , the algorithm reduces the O(td 2 ) coefficients of F modulo p and p ′ . Then, we apply γ to the results, to obtain H mod p and H mod p ′ . Due to Lemma 4, the coefficient of y d in h 1 is a nonzero constant; if this is still the case modulo p and p ′ , we use HermiteGroebnerBasis with D = d to get B p and B p ′ ; otherwise, we raise an error.

Changing coordinates in B p and B p ′ . Using the Gröbner bases B p and B p ′ of (H mod p) and (H mod p ′ ), we compute the Gröbner bases of (F mod p) and (F mod p ′ ). This is done using the algorithm of [START_REF] Neiger | Computing syzygies in finite dimension using fast linear algebra[END_REF]. Since pp ′ does not divide β F , we deduce that we obtain G 1 = G mod p and G ′ 1 = G mod p ′ . Computing G k . At each step of the main loop, we start from G k/2 = G mod p k/2 , and we compute G k = G mod p k . For this, we first need F mod p k ; then, we use procedure LiftOneStepGroebner from [START_REF] Schost | Newton iteration for lexicographic gröbner bases in two variables[END_REF]Remark 7.3] to obtain G k . Rational reconstruction. We next attempt to recover all rational coefficients of G, given those of G k = G mod p k . For each coefficient α of G k , we attempt to recover a pair (η, θ) in Z × N, with |η| < p k/2 /2 and θ ≤ p k/2 , θ invertible modulo p and α = η/θ mod p k .

By assumption, all nonzero coefficients of G have numerators and denominators of height at most b, it follows that if p k/2 > 2e b , we will succeed and correctly recover the corresponding coefficient in G [START_REF] Zur Gathen | Modern Computer Algebra[END_REF]Theorem 5.26]. For smaller values of k, rational reconstruction may find no solution (in which case we reenter the lifting loop at precision 2k), or may already terminate; in this case, its output G rec may be different from G. Testing for correctness. The final step in the loop is a randomized test, using G ′ 1 = G mod p ′ as a witness to detect those cases where rational reconstruction returned an incorrect result. We attempt to reduce G rec modulo our second prime p ′ (G red ); if this fails (because p ′ divides one of the denominators in it), we reenter the lifting loop at precision 2k. We simply compare G red and G ′ 1 = G mod p ′ . If they coincide, we return G rec , otherwise, we reenter the lifting loop.

Analysis

We assume that choosing a random integer in a set {0, . . . , A} (uniform distribution) uses O˜(log(A)) bit operations. We assume that we have an oracle O, which takes as input an integer C, and returns a prime number in I = [C + 1, . . . , 2C], uniformly distributed within the set of primes in I, using O˜(log(C)) bit operations. Parameters choice. The change of variables γ needs to avoid a hypersurface Y ⊂ Q 4 of degree at most A 1 = d 4 + d. We choose its entries uniformly at random in {0, . . . , 2 P +2 A 1 }; the cost of getting γ is negligible. Then, by the De Millo-Lipton-Schwartz-Zippel lemma, the probability that γ lies on Y is at most 1/2 P +2 . In what follows, we assume that this is the case, so all polynomials H = F γ have coefficients of height at most h ′ = h + d(P + 5 + log(A 1 )) ∈ O˜(h + dP ).

Let β F and β H be the nonzero integers from Definition 2 applied to respectively F and H, and define

Algorithm 1 GroebnerBasis(F) Require: F = (f 1 , . . . , f t ) in Z[x, y], d = max{deg f i }, Ensure: the lexicographic Gröbner basis of F in Q[x, y] 1: choose γ in M 2 (Z)
B i ← HermiteGroebnerBasis(H i , d) 7: G 1 ← ChangeCoordinatesGroebner(B p , γ -1 mod p) 8: G ′ 1 ← ChangeCoordinatesGroebner(B p ′ , γ -1 mod p ′ ) 9: k ← 1 10: repeat 11: k ← 2k 12: G k ← LiftOneStepGroebner(F mod p k , G k/2 ) 13: error, G rec ← RationalReconstruction(G k ) 14:
if not error then error, G red ← G rec mod p ′ 15: until not error and

G red = G ′ 1 16: return G rec C F = C(t, d, ∆ 3 (d), h) ∈ O˜(t 2 d 9 h) C H = C(t, d, ∆ 1 (d), h ′ ) ∈ O˜(t 2 d 4 hP ).
Proposition 4 proves height(β F ) ≤ C F and height(β H ) ≤ C H . In particular, the height bound b on the coefficients of G satisfies b ≤ C F , so b is in O˜(t 2 d 9 h). Let µ 1 be the coefficient of y d in h 1 , which has height at most h ′ . Our first requirement on p and p ′ is that neither of them divides µ = β F β H µ 1 . This is a nonzero integer, with height(µ) ≤ A 2 , where we set

A 2 = C F + C H + h ′ ∈ O˜(t 2 d 9 hP ).
Finally, we want to ensure that in the verification step, if G rec and G differ, their reductions modulo p ′ , called G red and G ′ 1 , differ as well. Below, we let k 0 be the first k which is a power of two and such that, at step k, rational reconstruction correctly computes G rec = G. For this, it suffices that p k/2 > 2e b , and one verifies this implies that k 0 ≤ 8b ∈ O˜(t 2 d 9 h). Since all indices k we go through are powers of two, there are at most log 2 (8b) incorrect indices k.

Suppose then that at step k < k 0 , we have found G rec with rational coefficients; they all have numerators and denominators at most p k/2 ≤ 2e b ; on the other hand, the coefficients of G have numerators and denominators at most e b . If G rec and G differ, there exists a monomial whose coefficients in G rec and G are different; it suffices that p ′ does not divide the numerator of their difference. This number has an absolute value of at most 4e 2b . Taking all k < k 0 into account, our last requirement is that p ′ also not divide a certain nonzero integer µ ′ (that depends on p). This integer µ ′ has height at most log 2 (8b)(2b + log(4)), so that height(µ ′ ) ≤ A 3 , with

A 3 = log 2 (8C F )(2C F + log(4)) ∈ O˜(t 2 d 9 h).
To summarize, once γ avoids Y , it suffices that p does not divide µ and p ′ does not divide µµ ′ to ensure success. We can then finally make our procedure for choosing p and p ′ explicit: bit operations per iteration. Here, n 0 = deg y (g 0 ) and m s = deg x (g s ). Rational reconstruction takes O˜(k log(p)) bit operations per coefficient, for a total of O˜(sδk log(p)). For the test: reduction modulo p ′ takes O˜(b+log(p ′ )) bit operations per coefficient, for a total of O˜(sδ(b + log(p ′ )).

• Let B = 2 P +3 ⌈A 2 ⌉.
Furthermore, both log(p) and log(p ′ ) are in O˜(P +log(tdh)). Besides, the definition of k 0 implies that at all lifting steps, k log(p) is in O˜(b + log(p)), that is O˜(b + P + log(tdh)). After some straightforward simplifications, the runtime becomes softly linear in td 2 h, (t ω d ω+1 + δ ω )(P + log(tdh)) and (s 2 n 0 m s + tδ(d 2 + dm s + sδ))(b + P + log(tdh)).

In order to get a better grasp on this runtime, let us assume that P and the number of equations t are fixed constants, and use the upper bounds n 0 , m s ≤ δ. This gives a total bound softly linear in The first term is the input size, the second describes computations done modulo small primes, and the last one computations are done modulo higher powers of p. The output size O(sδb) bits.

Computing the ⟨x, y⟩-primary component

Finally, we describe how to modify the algorithm if we are only interested in the Gröbner basis G 0 = (G 0 0 , . . . , G 0 r ) of the ⟨x, y⟩-primary component J of I. In what follows, we let η be the degree J, and c be the maximum height of the numerators and denominators of the coefficients of G 0 : the output total size is O(rηc) bits.

As before, we use a change of coordinates γ, and we call B 0 the Gröbner basis of J γ . Then, we use GroebnerBasisAtZero instead of Hermite-GroebnerBasis, modulo p and p ′ . Since we are in generic coordinates, we can use degree D = d, so the runtime is O˜(td ω m log(pp ′ )) bit operations, where m is the maximal x-degree of the polynomials in B 0 . We will use the bound m ≤ η.

Then LiftOneStepPunctualGroebnerBasis from [START_REF] Schost | Newton iteration for lexicographic gröbner bases in two variables[END_REF]Rk 7.3] can be used with an initial cost (bit operations) of O˜(tη ω log(p)) and O˜(tη 2 k log(p)) at the k th iteration. The rest of the algorithm is unchanged except for a slight difference in conditions of success. Now, γ has to avoid a hypersurface Y ′ of degree at most d 4 + d, in order to guarantee that J γ satisfies Lemmas 4 and 5. The primes p and p ′ should divide the denominator of no coefficient in G 0 and B 0 ; besides, these polynomials reduced modulo p (resp. p ′ ) should still define the ⟨x, y⟩primary components of f 1 mod p, . . . , f t mod p and f γ 1 mod p, . . . , f γ t mod p (resp. modulo p ′ ).

The ⟨x, y⟩-primary component of ⟨f 1 , . . . , f t ⟩ is the ideal generated by F ′ = (f 1 , . . . , f t , x d 2 , y d 2 ); similarly for H = (f γ 1 , . . . , f γ t ), giving us polynomials H ′ . It is then sufficient that neither p nor p ′ divides the integers β F ′ β H ′ from Definition 2. Their heights are in O˜(t 2 d 6 h) and O˜(t 2 d 6 h ′ ), where h ′ is the height bound on H.

The rest of the analysis is conducted as before. Given a fixed integer P , we deduce that we can compute the Gröbner basis G 0 , with a probability of success of at least 1 -1/2 P , using O˜(td 2 h + (td ω η + η ω )(P + log(tdh)) + tη 2 (c + log(tdh))) bit operations. Assuming t and P are fixed, this is softly linear in d 2 h + (d ω η + η ω ) log(h) + η 2 c. To wit, the input size is linear in dh and that the output size is in O(rηc) ⊂ O(η 2 c), with r the number of polynomials in G 0 η its degree and c the bit-size of its coefficients.

d 2 h

 2 + (d ω+1 + δ ω ) log(h) + (d 2 δ + dδ 2 + s 2 δ 2 )(b + log(h)).

  2: choose two different primes p, p ′ ; do steps 3-6 for i ∈ {p, p ′ } 3: if γ mod i is not invertible then raise an error 4: H i ← ChangeCoordinates(F mod i, γ mod i) 5: if the coefficient of y d in H i (1) is zero then raise an error 6:

  We use the oracle O to obtain a uniformly sampled prime number in [B + 1, . . . , 2B]. There are at least B/(2 log(B)) primes in this interval, and at most log(µ)/ log(B) of them can divide µ, so the probability that p does is at most 2 log(µ)/B, which is at most 1/2 P +2 .• Let B ′ = 2 P +3 ⌈A 2 + A 3 ⌉. We use the oracle O to pick p ′ in the interval [B ′ + 1, . . . , 2B ′ ],and as a result, the probability that p ′ divides µµ ′ is at most 1/2 P +2 .Altogether, the probability that γ avoids H, p does not divide µ and p ′ does not divide µµ ′ (and thus that the algorithm succeeds) is thus at least 1 -3/2 P +2 ≥ 1 -1/2 P . Complexity. To find B p and B p ′ : reducing the coefficients, changing coordinates and HermiteGroebnerBasis uses O˜(td 2 (log(pp ′ )), O˜(td 2 (h + log(pp ′ ))) [44, Corollary 9.16] and O˜(t ω d ω+1 (log(pp ′ ))) bit operations, respectively. Inverting the γ on B p and B p ′ takes O˜(δ 3 ) operations in F p ′ , which is O˜(δ 3 log(p ′ )) bit operations. To compute G k : coefficients reduction takes O˜(td 2 (h + k log(p))) bit operations. Algorithm LiftOneStepGroebner takes a one-time cost of tδ ω log(p) bit operations, plus for O˜(s 2 n 0 m s + tδ(d 2 + dm s + sδ))k log(p))

Main algorithmsWe can finally present our main algorithms, where we use Newton iteration to compute lexicographic Gröbner bases: we are given F = (f 1 , . . . , f t ) in Z[x, y],

Acknowledgments. We thank Arne Storjohann and Vincent Neiger for answering our questions on Hermite normal form computations. Schost is supported by an NSERC Discovery Grant. St-Pierre thanks NSERC, Alexander Graham Bell Canada Graduate Scholarship, FQRNT and the European Research Council (ERC) under the European Union's Horizon Europe research and innovation programme, grant agreement 101040794 (10000 DIGITS) for their generous support.


Abdennour Abbas

Thomas Dargent

Dominique Croix

Michel Salzet

Bertrand Bocquet
email: bertrand.bocquet@iemn.univ-lille1.fr

Ex-vivo detection of neural events using THz BioMEMS

Keywords: 87, 85, Ox Biomedical instrumentation and transducers, including micro-electro-mechanical systems (MEMS) -41, 20, -q Applied classical electromagnetism

In this paper, we present a Millimeter-wave BioMEMS (Biological MicroElectro-Mechanical System) dedicated to the ex vivo detection of nitric oxide synthase (NOS) activity. The latter is involved in neurodegenerative phenomena. The BioMEMS is fabricated by using polydimethylsiloxane (PDMS) sealed on a glass substrate supporting gold coplanar waveguides (CPWs). The NOS activity detection is performed by Millimeter-wave (MMW) transmission signals through CPWs placed under the lesion site of an immobilised leech nerve cord. Tests are carried out in the frequency range 140-220 GHz, and the results obtained show that MMW transmission spectroscopy combined with microfluidic technology can be used for monitoring biochemical events in aqueous environment and consequently in biological models.

Introduction

Today, the emergence of the microtechnology, combined with the microelectronic process, allows the creation of very sophisticated miniaturized objects for biological analysis. In these integrated circuits called BioMEMS, we can 2001). The microwave and MMW region up to several hundred GigaHertz, could present a very interesting alternative in terms of free-label investigation [START_REF] Facer G R, Notterman | Dielectric spectroscopy for bioanalysis: From 40 Hz to 26.5 GHz in a microfabricated wave guide[END_REF][START_REF] Hefti | Sensitive detection method of dielectric dispersions in aqueous-based, surfacebound macromolecular structures using microwave spectroscopy[END_REF], and more selective detection of radicals and biomolecules [START_REF] Smye | The interaction between TeraHertz radiation and biological tissue[END_REF][START_REF] Siegel | Terahertz technology in biology and medicine[END_REF]. This spectrum could also provide additional information in molecular biology, as for example, the conformational states of proteins [START_REF] Globus | THz characterization of lysozyme at different conformations[END_REF][START_REF] Markelz | THz time domain spectroscopy of biomolecular conformational modes[END_REF]. The results are expressed in terms of dielectric spectroscopy. Alternatively, an approach in far field is less attractive due to the poor spatial resolution of the wavelength. The use of planar waveguides inside integrated circuits makes possible these measurements. Nevertheless, increasing the frequency range up to the TeraHertz spectrum potentially has several interesting advantages.

Here, we want to detect the nitric oxide synthase (NOS) activity in the leech nerve cord after injury. This study is of great interest due to the critical role of the nitric oxide (NO), a gaseous molecule produced by NOS activity, in several biochemical networks, especially in the enhancement of nerve cord repair [START_REF] Chen | Nitric Oxide Influences Injury-Induced Microglial Migration and Accumulation in the Leech CNS[END_REF]. Today, the detection of NOS activity is performed by histochemical staining that requires fluorescent tags, whereas the monitoring of NO generation is generally realized by amperometric measurements [START_REF] Amatore | Monitoring in Real Time with a Microelectrode the Release of Reactive Oxygen and NitrogenSpecies by a Single Macrophage Stimulated by its Membrane Mechanical Depolarization[END_REF]. But these techniques remain difficult and depend strongly on the measurement conditions, such as the distance between the probe and the nerve cord, the molecular diffusion phenomena and the short half-life of NO in aqueous media (~ 6 s).

The idea here is to immobilize a nerve cord inside a MMW microfluidic system (Fig. 1(a)). We have validated this approach by using a mixed technology such as polymer on silicon for the analysis of protein solutions [START_REF] Mille | New technology for high throughput THz BioMEMS[END_REF] or living cells [START_REF] Treizebré | Planar excitation of Goubau transmission lines for THz[END_REF][START_REF] Treizebré | Nanometric metal wire as a guide for THz investigation of living cells[END_REF]. The present work is dedicated to neural tissues investigations.

Experimental

Animal Model

Our biological model is the leech Hirudo medicinalis, which is a very interesting and largely studied invertebrate. A lesion or a cut of its nerve cord, that represents the central nervous system (CNS), causes biochemical events leading to a complete neuron regeneration and restoration of the biological function in approximatively four weeks after damage [START_REF] Modney | Regeneration of a Central Synapse Restores Nonassociative Learning[END_REF]. Note that in mammals, the synapse regeneration is successful in the peripheral nervous system (PNS), but fail very quickly in the CNS [START_REF] Fawcett | Astrocytic and neuronal factors affecting axon regeneration in the damaged central nervous system[END_REF]. After nerve cord injury, the fast production of NO is one of the early events occurring at the lesion [START_REF] Kumar | Nerve Injury Induces a Rapid Efflux of Nitric Oxide (NO) Detected with a Novel NO Microsensor[END_REF]. This production is catalysed by an isoenzymes family called: NO synthases (NOS, EC. 1.14.13.39), according to the reaction (1):

L-arginine + 2 NADPH + 2 O 2 + 2 H + → NO + L-citrulline + 2 NADP + + 2 H 2 O (1)
(NADPH: reduced form of Nicotinamide Adenine Dinucleotide Phosphate)

Moreover, a treatment with the NOS inhibitor L-N G -Nitroarginine methyl ester (L-NAME) blocks the repair mechanism, showing the important role of the NOS to modulate the axon growth [START_REF] Chen | Nitric Oxide Influences Injury-Induced Microglial Migration and Accumulation in the Leech CNS[END_REF][START_REF] Shafer | C L Injury-induced expression of endothelial nitric oxide synthase by glial and microglial cells in the leech central nervous system within minutes after injury[END_REF]. However, the application of exogenous NO with the NO donor Spermine NONOate (SPNO) increases the concentration and leads to a total blocking of regeneration. In mammals, under pathological circumstances, the important increase in NOS activity induces toxicity and/or neurons apoptosis [START_REF] Bonfoco | Apoptosis and necrosis: Two distinct events induced, respectively, by mild and intense insults with N-methyl-D-aspartate or nitric oxide/superoxide in cortical cell cultures[END_REF]. The produced NO is involved in neurodegenerative phenomena, particularly Alzheimer's diseases and Parkinsonism [START_REF] Zhang | Role of nitric oxide in Parkinson's disease Pharmacology & Therapeutics[END_REF]. An important question is: how the precise regulation of the NOS activity in the leech induces neuron regeneration and not cell death? The real time monitoring of NOS activity is a main step towards understanding the role of NO in nerve repair and how man loses his capacity to regenerate his CNS during the neurodegenerative disorders.

BioMEMS Design

Before designing our microsystem, the average dimensions of the connectives and ganglions of the leech Hirudo medicinalis were determined. The microchannel (width: 450 µm, depth: 1 mm, length: 4 cm, volume less than 15 µl) is matched to the nerve cord dimensions for improving its immobilization and limiting the medium evaporation. The MMW BioMEMS design uses a classical technology based on PDMS due to the large dimensions of the microchannel [START_REF] Mcdonald | Fabrication of microfluidic systems in poly(dimethylsiloxane)[END_REF]. The PDMS is prepared by mixing the precursor Sylgard 184 (Dow Corning)

with the crosslinking agent at 10:1 weight ratio. We use a coplanar waveguide (CPW) which is well suited to this measurement [START_REF] Gupta | Microstrip Lines and Slotlines Second édition[END_REF]. It is placed under the probes of the vectorial network analyzer (VNA) (Fig. 1(b)). . We obtain a good correlation between the computed and the measured values on water (Fig. 3).

The microsystem is realized on a 2-inch-square-glass substrate. A thin binding of 200 Å Chromium (Cr) layer was sputtered on the substrate before the deposition of a 0.5 µm Gold (Au) layer. Then, both Gold and Chromium are etched in order to create the CPW lines. Two-step photolithography has been adopted for the substrate processing.

The PDMS channel was obtained by molding in a mechanically etched piece of Teflon ® . To obtain a permanent bond of the device, both glass substrate and PDMS channel were exposed separately to plasma oxidation before joining them together.

Spectral Measurements

The nervous chains were extracted from adult leeches (Hirudo medicinalis) weighing from 2 to 3 mg, and maintained in Ringer's solution (pH 7.4) for their survival.

Results and Discussion

In order to study the transmission spectral characteristics of the CPW, we initially carried out tests separately on well-known deionized water and on the Ringer's solution (115 mM NaCl, 4 mM KCl 2 , and 1.8 mM CaCl 2 ), buffered at pH 7.4 with 10 mM Tris maleate (Fig. 3(a)). As it can be seen with the computed values obtained on pure water, the transmission coefficient decreases with increasing frequency, but remains valuable for a good measurement. The very small volume matches our measurements well despite the high water absorption. Note that we can observe some small rebound phenomena due to the standing waves minimized by the impedance optimisation. The same figure

shows a very small difference between the spectral responses of water and Ringer's solution, which can be explained by the fact that moderate ionic concentrations have negligible effects on transmission spectra [START_REF] Xu | 0.15-3.72 THz absorption of aqueous salts and saline solutions[END_REF]. The measured values were stables and reproducibles and the measurement errors are estimated at ±0.1 dB. This result shows that the Ringer's solution can be a useful survey medium for in vivo/ex vivo studies by MMW spectroscopy.

Therefore, when the nerve cord is immobilized inside the microchannel containing Ringer's solution, we observe an improvement of the transmission coefficient (Fig. 3(a) and (b)). This is mainly due to the displacement of water by nerve tissues.

The second step was the calibration with L-NAME in the millimolar concentration range. Figure 4 shows that the Beer's law is validated for the millimolar concentration, in the frequency range 145-200 GHz. We can see clearly a linear change in MMW transmission with the concentration of L-NAME. The increase in transmission in this case could be explained by the increasing number of bound water molecules (hydration shell) that present low-mobility and consequently lower dielectric constant than free water [START_REF] Mickana | Terahertz spectroscopy of bound water in nano suspensions[END_REF]. In these experiments, the detection limit of the bioMEMS reaches down to 0.01 g/ml. This limit, obtained by varying L-NAME concentrations, concerns only hydrated biomolecules that increase transmission and not the NOS activity products (probably NO) that increase absorption as we can see below. The last step of these experiments is a preliminary comparative study of the leech nerve cord before and after injury. The aim is to block the NOS activity after injury, and observe the change in MMW transmission features. First, nerve cords were injured directly inside the microchannel filled with Ringer's solution without any treatment (Fig. 5). As it can be seen, the lesion causes a decreased transmission coefficient from approximately 0.8±0.1 dB comparatively with the intact cord. Curves (b) and (c) show respectively that after injury, NOS activity reaches immediately a peak of intensity, and decrease over time never returning to preinjury level.

Secondly, nerve cords were exposed to 1 mM L-NAME during 40 min before injury. L-NAME is the antagonist of L-arginine amino-acid. Thus, it inhibits the NO synthase enzyme and then block the NO and L-citruline production.

Here, the L-NAME was used as a reference sample, and any statistically significant differences in the transmission of the intact nerve cord and the nerve cord injured was observed after L-NAME treatment. Similar results obtained by the amperometric method are described in the literature [START_REF] Kumar | Nerve Injury Induces a Rapid Efflux of Nitric Oxide (NO) Detected with a Novel NO Microsensor[END_REF]. Note that the L-NAME concentration at 1 mM have no effect on transmission features according to the limited sensitivity of the present bioMEMS. The transmission change described above is likely due to the highest absorbance of millimeter-waves by the released products.

Conclusions

Studies of the interaction between millimeter-waves and reactive molecules in liquid phase would lead to a better understanding of the biological activity. This work demonstrates that the combination between millimeter-waves and microfluidic circuits inside a BioMEMS allows us to perform ex vivo measurements. The main interest is the integration on a micrometric scale to obtain fine spatial resolution and selective detection. The second advantage is the very small volume which lead to a huge increase in the local concentration of the reaction products. We have specifically designed a microsystem for ex vivo analysis of leech nerve cords and the preliminary results show two main considerations: Firstly, our measurements show clearly that it is possible to use millimeter wave spectroscopy for the detection of biochemical reactions in aqueous solution. Secondly, we have characterized the NOS activity around a leech nerve cord. Future experimental tests will be dedicated to the investigation of analytical parameters and the real time monitoring of NOS activity. We will obtain quantitative measurements of the reaction products for studying the effects of several drugs on NO production, and then on nerve repair mechanisms.

  ird-00436065, version 1 -25 Nov 2009 Author manuscript, published in "Medical Science Monitor 15, 9 (2009) MT121-125" Manuscript published in the journal Medical Science Monitor Med Sci Monit 2009; 15(9):MT121-125 Manuscript ID: 878163 2 mix electronic and microfluidic functions. The large majority of biosensors use the electrochemical principle or optical detection. But these techniques present some disadvantages. The principal one of them is probably the alteration of the biological activity by chemical reaction or by the fluorescent tags bonded on molecules (Facer

  Manuscript published in the journal Medical Science Monitor Med Sci Monit 2009; 15(9):MT121-125 Manuscript ID: 878163 This planar waveguide is constituted by three conductors and designed on a glass substrate (Fig. 2(a)). The latter has low losses in the MMW propagation and its transparency allows further use of contrast phase microscopy for the observation of the nerve cord under measurement. We have optimized the propagation condition on the real microsystem composed by a tank filled with liquid water. Many geometric configurations in terms of slot width to central conductor width ratio are possible in order to obtain the wanted impedance of 50 Ω, which is recommended for the connection at the VNA. Simulations are performed with a full 3D software named Microwave Studio ® from Computer Simulation Technology (CST) and based on a finite element method (FEM) model. The true sizes take into account the propagation of the fundamental non-dispersive mode and a small excursion of the electric field for reducing the radiated losses. An example of a simulated structure is shown on the Figure 2(b)

  They are manually placed and immobilized inside the open microchannel after PDMS-glass bonding. Dissection and lesions performed during the experiment are carried out in aseptic conditions using Patscheff scissors. L-NAME hydrochloride was purchased from Cayman Chemical. Solutions with different concentrations were prepared by dissolving L-NAME powder in the Ringer's solution, and injections were performed using a Hamilton Syringe, gauge 33 (internal diameter: 0.21 mm) obtained at NHamilton Bio. The transmission measurements were carried out at room temperature, in aerobic conditions. The VNA is an Anritsu 37147C associated with mixers of reference V05VNA2-T/R from OML (Oleson Microwave Laboratories) working in a bandwidth of 140-220 GHz. We use Line-Reflect-Match (LRM) calibration with a calibration kit Manuscript published in the journal Medical Science Monitor Med Sci Monit 2009; 15(9):MT121-125 Manuscript ID: 878163 reference 101-190B of Cascade Microtech. The analyzer is gauged in the plane of the measurement probes. Here, we have exploited only the transmission modulus relative to the wave absorption.

Figure 1 Figure 1 .Figure 2 Figure 2 .

 1122 Figure 1

Figure 3 Figure 3 .Figure 4 Figure 4 .

 3344 Figure3

Figure 5 Figure 5 .

 55 Figure 5


Abdou Kâ Diongue
email: abdou.diongue@gmail.com

Dominique Guegan
email: dguegan@univ-paris1.fr

Bertrand Vignal
email: bertrand.vignal@edf.fr

Abdou Kâ Diongue

Dominique Guégan

Forecasting electricity spot market prices with a k-factor GIGARCH process

Keywords: Conditional mean, Conditional variance, Electricity prices, Forecast, GIGARCH process

In this article, we investigate conditional mean and conditional variance forecasts using a dynamic model following a k-factor GIGARCH process. Particularly, we provide the analytical expression of the conditional variance of the prediction error. We apply this method to the German electricity price market for the period August 15, 2000 -December 31, 2002 and we test spot prices forecasts until one month ahead forecast. The forecasting performance of the model is compared with a SARIMA-GARCH benchmark model using the year 2003 as the outof-sample. The proposed model outperforms clearly the benchmark model. We conclude that the k-factor GIGARCH process is a suitable tool to forecast spot prices, using the classical RMSE criteria.

Introduction

In finance, commodity price forecasting is a crucial issue, in particular for electricity. Indeed, price forecasts can help evaluating bilateral contracts.

For such a commodity, price forecasts can be based on a spot price analysis.

Since electricity spot prices have a behavior which presents specific features, the price forecasts problem can be complex. In another hand, electricity price volatility has a marked variability in time. We observe both high and low periods of price reaction. Recently spot price volatility has been studied and modelled using a volatility model depending on time. Using such a model, Benini, Maracini, Pelacchi and Venturini (2002) investigate several markets. In addition, electricity spot prices exhibit long memory behavior combined with a periodic behavior. Recent works have taken into account these last features using related ARFIMA models, see [START_REF] Koopman | Periodic seasonal reg-ARFIMA-GARCH models for daily electricity spot prices[END_REF] or Diongue and Guégan (2008).

After modelling these electricity prices, the forecasting problem arises. In the literature, two approaches have been considered: parametric models using AR, ARX, AR-GARCH, ARX-GARCH and Regime switching models, Misiorek, Truek and Weron (2006); and non parametric methods like the neuronal nets for instance, Ramsay and Wang (1997) and Conejo, Contreras, Espinola and Plazas (2005).

In this paper, in order to provide robust forecasts for spot electricity prices, we propose a new approach based on the k-factor GIGARCH process, (Guégan, 2000), which allows taking into account a lot of stylized facts observed on the electricity spot prices, in particular stochastic volatility, long memory and periodic behaviors. The investigation of this model is done in [START_REF] Guégan | A Prospective Study of the k-factor Gegenbauer Process with Heteroscedastic Errors and an Application to Inflation Rates[END_REF]. Diongue and Guégan (2004) introduced the parameter estimation of the k-factor GIGARCH process. Here, we provide new developments which concern the expression of the forecasts using the k-factor GIGARCH process and we give their properties. We apply these results on the German electricity prices market providing forecasting prices up until a one month-ahead.

These results are totally new in the sense that, in most published papers, the previsions concern mostly the one day-ahead horizon and here we are interested by long term prediction. For comparison purpose, we compare our approach with a benchmark model in terms of forecasting, using the RMSE criteria.

This paper is organized as follows. The next Section presents the data and contains the main empirical findings. In Section 3, we specify some notations, define the k-factor GIGARCH process along with new theoretical results on forecasting in mean and in variance. In the fourth Section, we provide forecasts for the German spot prices data set. The last section is dedicated to conclusions.

The data set

We consider the hourly series of the EEX spot prices, denoted (S t ) t=1,••• ,N , from August 15, 2000 to December 31, 2003, which yields N = 29616 hourly observations. For this electricity spot market, the prices are fixed each day, for each of the 24 hours of the next day. This spot market has two particular interests. On the one hand, it allows a buyer (seller) to supply their bid (load) on an hourly basis (physical market) and on the other hand, the spot market could correspond to the price reference for the contracts. In order to make the series (S t ) t=1,••• ,N stationary, we take its logarithmic form and define (X t ) t such that, ∀t = 1, • • • , N , X t = log (S t ). We build two subsets from the original series. The first one, from August 15, 2000 to December 31, 2002 (figure 1) and containing T = 20856 observations, is used for the parameter estimation procedure and the second one, from January 1, 2003 to December 31, 2003, is used for validation. In figure 1, we observe that there are periods of low volatility (Figure 1, Regions 1) following by periods of high volatility (Figure 1, Regions 2). Figure 2 provides a zoom on the evolution of the spot prices, in its logarithmic scale, from October 14, 2000 to October 27, 2000 (two weeks). We observe that the series exhibit two kinds of seasonality: daily and weekly seasonalities.

Insert Figures 1 and2 Table 1 gives the summary statistics of the log electricity spot prices. The mean and the standard deviation are quite small, while the estimated measure of skewness, -2.6, is significantly negative, indicating that the log electricity spot prices has non-symmetric distribution. Furthermore, the large value of the kurtosis statistic, 29.67, suggests that the underlying data are leptokurtic, or fat-tailed and sharply peaked around the mean when compared with the standard Normal distribution (Figure 3). The results of the non-Normality test agree with previous literature concerning this kind of behavior for these log electricity spot prices data, Misiorek, Trueck and Weron (2003). The Box-Pierce Q-tests of up to twenty-fourth order serial correlation for the levels and squares of the mean-corrected log electricity spot prices (Q (24) and Q 2 (24)) are both significant. In summary, the diagnostics suggest that a GARCH-class model would be appropriate to explain the evolution of the log electricity spot prices, along with an error distribution that allows for greater kurtosis than the Gaussian distribution.

Insert Table 1 and Figure 3 The autocorrelation function of the EEX spot prices (Figure 4) exhibits a slow decay pattern at two seasonality lags corresponding to daily and weekly seasonalities. The autocorrelation function does not converge exponentially to zero in covariance sense, [START_REF] Guégan | How can we define the concept of long memory? An econometric survey[END_REF], for discussions on this concept), thus we assume that the series has a long range dependence behavior. Moreover, the spectral density represented by the periodogram (Figure 5), is unbounded at three frequencies corresponding respectively, to the daily, weekly and half daily seasonalities.

Insert Figures 4 and 5

Thus, this empirical and statistical analysis highly justifies the use of the k-factor Gegenbauer process with time varying conditional variance errors to model the data set (X t ) t .

3 The k-factor GIGARCH process

Notations and model

Assume that (ξ t ) t∈Z is a white noise process with unit variance and mean zero. Let (X t ) t∈Z be a k-factor GIGARCH process defined by the following equations:

Φ (B) k i=1 I -2ν i B + B 2 d i (X t -µ) = Θ (B) ε t , (3.1)
where

ε t = h t ξ t with h t = a 0 + r i=1 a i ε 2 t-i + s j=1 b j h t-j . (3.2)
In the relationship (3.1)-(3.2), µ represents the mean of the process

(X t ) t∈Z , d = (d 1 , • • • , d k ), the memory parameters and ν = (ν 1 , • • • , ν k ) the frequency location parameters with, 0 < d i < 1 2 if |ν i | < 1, and 0 < d i < 1 4 if |ν i | = 1 for i = 1, • • • , k, with k a nonnegative integer. Φ (B)
and Θ (B) are the autoregressive and the moving average polynomials in

B (B n X t = X t-n , n ≥ 0), defined by Φ (B) = I -p i=1 φ i B i and Θ (B) = I -q j=1 θ j B j
, where p and q are nonnegative integers. We suppose that all the roots of the polynomials Φ (B) and Θ (B) are outside the unit circle. We assume that

a 0 > 0, a 1 , . . . , a r , b 1 , . . . , b s ≥ 0, and r i=1 a i + s j=1 b j < 1. The k-factor GIGARCH process (3.1)-(3.
2) under the previous assumptions is well defined. We can use it to forecast electricity spot prices. The interest of such a heteroskedastic model is that it allows emphasizing the optimal predictor for conditional variance. A basic property for the k-factor GARMA process with time varying innovations following a GARCH(r, s) model is that the long memory part provides a way to model the conditional mean whereas the heteroscedastic part provides a way to model the conditional variance. From a forecasting point of view, the k-factor GIGARCH model [START_REF] Guégan | A New Model: The k-factor GIGARCH Process[END_REF] and the k-factor GARMA model (Giraitis and leipus, 1995) will give the same optimal predictor for the conditional mean and the difference between their respective forecasts appears mainly in the building of the confidence interval.

Forecasting equation for the k-factor GIGARCH process

In this section, we provide new results concerning the predictions based on a strong white noise (ξ t ) t∈Z .

In a first step, we provide the expression of the forecasts in mean. Observing 2) have been estimated using Whittle method, see Diongue and Guégan (2004). Under the regular conditions specified in section 3.1, this model can be rewritten in its infinite moving average form:

X 1 , • • • , X T ,
X t = ∞ j=0 βj d, ν, γ ε t-j , (3.3)
where the weights βj d, ν, γ j∈Z follow the expression:

βj d, ν, γ = ψj d, ν + min(j,p) i=1 φi βj-i d, ν - min(j,q) i=1 θi ψj-i d, ν, γ , (3.4)
and

ψj d, ν = 0≤l 1 ,••• ,l k ≤j, l 1 +•••+l k =j C l 1 d1 , ν1 • • • C l k dk , νk . (3.5)
The weights C l i d, ν

l i ∈Z
are the Gegenbauer polynomials. Denoting the

vectorial parameters d = d1 , • • • , dk , γ = φ1 , • • • , φp , θ1 , • • • , θq , and ω = â0 , â1 , • • • , âr , b1 , • • •
, bs , we can derive the h-time-step-ahead least square predictor, Xt (h), given by the following expression:

Xt (h) = ∞ j≥h βj d, ν, γ ε t+h-j , (3.6)
which can be rewritten as

Xt (h) = ∞ j≥0 βj+h d, ν, γ ε t-j . (3.7)
Unfortunately, in practice, the process (ε t ) t∈Z cannot be observed, thus this latter expression cannot be used to obtain the forecasts in mean from the process. Therefore, we need to investigate an autoregressive form of the kfactor GIGARCH process. Following Proposition 4 in Ferrara and Guégan (2001), it is straightforward to show that a useful and easy formula to update, for computing the least-squares predictor Xt (h), is given by

Xt (h) = - h-1 j=1 αj d, ν, γ Xt (h -j) - j≥0 αj+h d, ν, γ X t-j , (3.8
êt,h = X t+h -Xt (h) = h-1 j=0 βj d, ν, γ ε t+h-j . (3.10)
The conditional mean squared error for this predictor has the following expression:

E ê2 t,h | F t = h-1 j=0 β2 j d, ν, γ E ε 2 t+h-j | F t , (3.11)
where F t is the σ-algebra generated by ε 0 , • • • , ε t , (the noise which appears in

(3.1) -(3.2)
). We remark that the conditional mean squared error depends on the future conditional variance of the innovation process ε 2 t t∈Z . Its expression, denoted var t (ê t,h ), is given by

var t (ê t,h ) = h i=1 β2 h-i d, ν, γ ωi + h k=1 β2 h-k d, ν, γ s-1 i=0 δi,k h t-i + m-1 i=0 ρi,k ε 2 t-i , (3.12)
with m = max (r, s), ωh = e

′ 1 I + Γ + • • • + Γ h-1 e 1 â0 , δj,h = -e ′ 1 Γ h e m+j+1 if j = 0, • • • , s -1, ρj,h = e ′ 1 Γ h (e j+1 + e m+j+1 ) if j = 0, • • • , s -1, ρj,h = e ′ 1 Γ h e j+1 if j = s, • • • , m -1.
Here {e j } m+s j=1 is the standard base of R m+s and the matrix Γ is defined such that its first row is equal to â1 + b1 , • • • , âm + bm , -b1 , • • • , -bs , and the other elements are all equal to zero, except the first one under diagonal elements which are all equal to one. We use the expression (3.12) to build the confidence interval around Xt (h).

A real case study

In this section, we investigate the German Energy EXchange market (EEX) corresponding to the Wholesale spot electricity market and the associated spot prices. The German spot prices have several key features such as: long memory persistence with seasonality and non-constant volatility. These stylized facts strongly suggest to use a GIGARCH process in order to make forecast. Thus, we adjust several k-factor GIGARCH processes on this data set and compare their forecasting performances with a benchmark model, a SARIMA-GARCH model.

Estimated models

In this subsection, we compare competitively three different modelling approaches in order to forecast this hourly series for different forecast horizons, starting from the first of January 2003. The first model is the SARIMA model of [START_REF] Box | Time series: Theory and Methods[END_REF] with conditional heteroskedastic noise. We denote it M 1 . The second model is a 3-factor GIGARCH model denoted M 2 . Lastly, the third one is a 1-factor GIGARCH process estimated after removing the weekly seasonality by applying the filter 1 -B 168 to the data set. We denote it M 3 1 . Below, the values given in brackets correspond to the standard deviation of the estimated parameters.

The model M 1

The autocorrelation function of the EEX spot prices (Figure 4) exhibits an exponential rate at the daily lag. Then, using the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC), we investigate, in that sense, the best SARIMA(p, d, q) × (P, D, Q) 24 -GARCH(1, 1) model. 1 The results of this model are available from authors upon request.

The selected model M 1 is defined by:

(I -B) I -B 24 φ (1) (B) Φ (1) B 24 (X t -2.97) = θ (1) (B) Θ (1) B 24 ε t , (4.1)
where,

φ (1) (B) = I - φ(1) 1 B,
with φ(1) 1 = 0.648(7.425 × 10 -4 ), and

Φ (1) B 24 = I - φ(1) 24 B 24 - φ(1) 48 B 48 , φ (1
) 24 = 0.1221(5.756 × 10 -5 ), φ(1) 48 = -0.0579(5.5707 × 10 -5 ). We also get

θ (1) (B) I - θ(1) 1 B, with θ (1
) 1 = 0.7479(5.6393 × 10 -4 ), and

Θ (1) B 24 = I - θ (
t = â(1) 0 + â(1) 1 ε 2 t-1 + b(1) 1 h (1) t-1 , (4.2)
with â(1) 0 = 0.00872(6.6668 × 10 -9 ), â(1) 1 = 0.3799(1.8223 × 10 -5 ), and b(1) 1 = 0.4867(1.9028 × 10 -5 ).

Model M 2

The periodogram of the EEX spot prices is unbounded at three frequencies, thus we consider that the series (X t ) t∈Z can be modelled using a 3-factor GIGARCH process. The observed prices (X t ) t , in their logarithmic scale, have the following representation:

Φ (2) (B) 3 i=1 I -2ν (2) i B + B 2 d(2) i (X t -2.97) = Θ (2) (B) ζ t , (4.3)
For each model, forecasts are computed for h = 1, • • • , 720 hours (from one hour to one month). As we use the prices in a logarithmic scale, we need to compute the exponentials of these log-forecastings. In order to get these forecasts, we assume the Gaussianity of the white noise (ξ t ) t∈Z in equation (3.2). Then, the forecasts for the spot prices at time T can be computed using the following expression

ŜT (h) = exp XT (h) + 1 2 var [ê T,h ] . (4.6)
Using two particular weeks, on Figures 67 2.

Insert Figures 6 and 7

In Table 2, we provide errors for the thirty days ahead, for models (M 1 ), (M 2 ) and (M 3 ), for some specific hours. We notice that the model (M 2 ) outperforms the benchmark model (M 1 ), even for short-term prediction. In addition, we empirically prove that, with a well specified nonlinear longmemory model, the predictions are more accurate than with a misspecified linear long-memory model. This kind of misspecification and its disastrous consequences for long range forecasting have already been pointed out by [START_REF] Ray | Modelling long-memory processes for optimal long-range prediction[END_REF] through a simulation study. We also perform forecast-accuracy comparisons in order to evaluate if the accuracy of the SARIMA-GARCH and 3-factor GIGARCH models are equal. Visual inspection of the sample autocorrelation function of the loss differential series in Figure 8 is suggestive of long-range dependence behaviour. Moreover, the Box-Pierce Q-test of up to twenty five order, Q (25) = 5.6558 × 10 3 , is highly significant relative to its asymptotic null distribution of χ2 25 . According to this latter result, the simple F test, Morgan-Granger-Newbold test and Meese-Rogoof test are not applicable.Therefore, we perform the comparing predictive accuracy test proposed [START_REF] Diebold | Comparing predictive accuracy[END_REF]. For the SARIMA-GARCH and the 3-factor GIGARCH models' forecasts, the estimated S 1 statistic 2 , S 1 = -2.3546 implying a p value of 0.0185, is well above the critical value at the 5% level of significance, thereby rejecting the null hypothesis of equal forecast roots expected square error. However, when considering the 3-factor GIGARCH and 1-factor GIGARCH models, we obtain S 1 = -1.6216, implying a p value of 0.1049. In this latter case, we do not reject at 5% levels the null hypothesis of equal roots expected square error. Thus, these comparative predictive accuracy tests proposed by [START_REF] Diebold | Comparing predictive accuracy[END_REF] confirm the supremacy of model M 3 for forecasting this electricity data set.

Insert Table 2 5 Conclusion

In this paper we have investigated forecasting's method using a stochastic model such as a k-factor GIGARCH process. We derive a least square predictor and its properties and we characterize the conditional variance error of this predictor. The results provided in section 3.2 are new and permit to obtain close form expressions for predictions using a k-factor GIGARCH process.

This forecasting method is applied to the German hourly electricity spot market prices. We adjust three different models on this data set: a SARIMA-GARCH model as a benchmark, a 1-factor GIGARCH model and a 3-factor GIGARCH model. The forecasting results for the year 2003 with the estimated models are highly convening in the sense of RMSE criteria. The model M 2 (3-factor GIGARCH process) provides better forecasts in the sense of RMSE criteria than model M 1 (SARIMA-GARCH process) and model M 3

(1-factor GIGARCH process), when modelling EEX prices on the period under study. In addition, the comparing predictive accuracy test from [START_REF] Diebold | Comparing predictive accuracy[END_REF] suggest that the expected square error for the models forecasts are not equal, confirming the fact that this new modelling improves the forecasts for this kind of data set.

Thanks to this study, we have detected presence of seasonal long memory and heteroskedascity in electricity spot prices and provide interesting modelling to take them into account. Note that these features appear frequently inside the European electricity market.

  using the stationary and invertible k-factor GIGARCH process defined by equations (3.1)-(3.2), we compute the h-time-step-ahead least square predictor, denoted Xt (h) , h > 0. The predictor is obtained by minimizing the mean squared error between forecasts and observations and our approach follows the work of[START_REF] Ferrara | Forecasting with k-factor Gegenbauer processes: Theory and Applications[END_REF]. We assume that the parameters of the model (3.1)-(3.

  ) where the weights αj d, ν, γ j∈Z are defined by αj d, ν, γ = πj d, ν -min(j,p) i=1 φi πj-i d, ν + min(j,q) i=1 θi αj-i d, ν, γ , (3.9) and πj d, ν = ψjd, ν . In practice, we use the expression (3.8) to get the forecasts. Now, we provide the expression for the variance's errors for the forecasts. The h-time-step-ahead prediction error for the predictor Xt (h) given in (3.8) is denoted e t,h for h > 0, and is equal to:

  , we illustrate the seven step ahead forecasting performance of the M 2 model (3-factor GIGARCH model) and the benchmark M 1 model (SARIMA-GARCH). The first week goes from 16, 2003 to March 22, 2003 while the second goes from September 8, 2003 to September 14, 2003. We notice that the forecasts from the 3-factor GIGARCH model (model M 2 ) are closer to the observed prices values than those obtained from the benchmark model (model M 1 ), with respect to the RMSE, see Table

Figure 1 :

 1 Figure 1: Evolution of the EEX spot market prices in its logarithmic scale from 15/08/2000 to 31/12/2002.

Figure 2 :

 2 Figure 2: Two weeks zoom concerning the evolution of the EEX spot market prices in its logarithmic scale from october 14, 2000 to october 27, 2000.

Figure 3 :Figure 4 :

 34 Figure 3: Non-parametric density of the log electricity spot prices (solid line) and probability density function of the Normal distribution (dotted line)

Figure 5 :

 5 Figure 5: Periodogram of the EEX spot market prices in its logarithmic scale from 15/08/2000 to 31/12/2002.

Figure 6 :

 6 Figure 6: Seven days ahead spot prices forecast for EEX spot market, for both the model M1 and M2, from March 16, 2003 to March 22, 2003.

Figure 7 :

 7 Figure 7: Seven days ahead spot prices forecast for EEX spot market, for both the model M1 and M2, from September 8, 2003 to September 14, 2003.

Figure 8 :

 8 Figure 8: Autocorrelation function for the loss-differential series.

Table 1 :

 1 Statistics of daily log returns of the S&P 500 stock market index.The Jarque-Bera test critical value at significance level of 5% is 5.85423.Q and Q 2 are the Box Pierce statistics for the levels and square of the log electricity spot prices, respectively, using 24 lags. The critical value at level 5% is 36.4150.

	Number of observations 20856	Skewness	-2.6
	Mean	2.97	Kurtosis	29.67
	Standard deviation	0.6241	Jarque-Bera test 7.3098×10 5
	Minimum	-4.6052	Q (24)	8.4597×10 4
	Maximum	6.5387	Q 2 (24)	3.7156×10 4

Table 2 :

 2 RMSE of the two models associated to the spot prices at some hours each day from one day-ahead to one month-ahead spot prices forecasts (the bold numbers are RMSE provided by model M2 while the italic are from model M3).

The estimated statistic S1 for testing the null hypothesis of equal forecast accuracy is d/ 2π fd (0) /T , where d is the sample mean loss-differential and fd (0) is a constant estimate of f d (0) given by 2π fd (0) = T -1 τ =-(T -1) 1 (τ /S (d)) γd (τ ), with γd is the sample autocovariance of d, 1 (τ /S (d)) is the uniform lag window, S (T ) is the lag truncation and T the sample size of the loss differential series.

Acknowledgments: We wisk to thank the two referees as well as the editor for helpful comments and suggestions, which have permittedto improve the presentation of this paper.

where d(2) 1 = 0.135(0.009), d(2) 2 = 0.214(0.006) and d(2) 3 = 0.141(0.015). The polynomials Φ (2) (B) and Θ (2) (B) are respectively equal to:

where φ24 = 0.863(0.179), φ168 = 0.951(0.221) and

where â0

(2) = 0.0141(3.39.10 -4 ) and â1

(2) = 0.756(0.365). Here, we have retained three different G-frequencies for which the spectral density function is unbounded:

where M is the size of the out-of-sample forecast period and Xt (i) is the predicted value of X t+i . We use this criteria because we are interested in the measure of the variation between the observed data and the forecast data.