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

On the complexity of the LickteigRoy subresultant algorithm

In their 1996 article, Lickteig and Roy introduced a fast divide and conquer variant of the subresultant algorithm which avoids coecient growth in defective cases. The present article concerns the complexity analysis of their algorithm over eective rings endowed with the partially dened division routine. This leads to new convenient complexity bounds for gcds, especially when coecients are in abstract polynomial rings where evaluation/interpolation schemes are not supposed to be available.

Introduction

Euclidean polynomial remainder sequences are a cornerstone of computer algebra for gcd, lcm, modular inversion, BerlekampMassey algorithm, Padé approximant, etc. This article brings new complexity results for computing subresultant polynomials over commutative rings endowed with the partially dened division routine. We design a divide and conquer algorithm with a bounded coecient growth even in defective cases. In particular, this leads to new deterministic complexity bounds for bivariate gcds.

Notations and denitions

Until the end of the article, A is a commutative ring with unity, and Frac(A) is its total ring of fractions (also called the total quotient ring), namely S ¡1 A where S is the set of the nonzero divisors in A. In other words, elements of Frac(A) are of the form a/s, with a 2 A and s 2 S. When A is integral, that is an integral domain, then Frac(A) is its eld of fractions.

Divisions. For algorithmic purposes, A is assumed to be endowed with the partially dened division routine: precisely, if a and b are two elements of A such that b divides a, then this routine returns a/b. Operations in A at our disposal are: addition, subtraction, multiplication, and this partially dened division. Let A and B be two polynomials in A [x]. If B = / 0, we say that the division of A by B is well dened when there exist Q and R in A[x] such that A = Q B + R and deg R < deg G. These polynomials Q and R are respectively written quo(A; B) and rem(A; B).

The leading coecient of a polynomial A is written lc(A). The pseudo-division of A by B = / 0 in A[x] is the division of lc(B) deg A¡deg B+1 A by B: the remainder (resp. quotient), written prem(A; B) (resp. pquo(A; B)), is called the pseudo-remainder (resp. pseudoquotient). Pseudo-divisions have the advantage to be well dened and easily computable without divisions in A.

Subresultant and Euclidean sequences. Our two input polynomials are written F = P i=0 n 0 f i x i and G = P i=0 n 1 g i x i , and are of respective degrees n 0 and n 1 . Throughout the article, we assume that n 0 > n 1 . For 0 6 k < n 1 , the kth subresultant coecient of F and G is written s k 2 A, and the associated kth subresultant polynomial S k 2 A[x] (usual denitions are recalled in Section 2.1). The subresultant polynomial S k is said to be defective when its degree is strictly less than k. Of course, when n 0 < n 1 , without loss of generality, we may swap the two input polynomials since the subresultant sequences coincide up to signs.

The extended Euclidean algorithm consists in computing the remainder sequence, recursively dened by R 0 = F , C 0 = 1, D 0 = 0, R 1 = G, C 1 = 0, D 1 = 1, and

R i+1 = R i¡1 ¡ E i R i , C i+1 = C i¡1 ¡ E i C i , D i+1 = D i¡1 ¡ E i D i
, where E i = quo(R i¡1 ; R i ). Consequently we have R i+1 = rem(R i¡1 ; R i ), and R i = C i F + D i G for all i > 0. The extended Euclidean sequence (R i ; C i ; D i ) i is well dened over Frac(A) whenever the leading coecients of the nonzero R i are invertible in Frac(A). In this case, the sequence ends after w ¡ 1 division steps with R i = / 0 for all 0 6 i 6 w, and R w+1 = 0 the last nonzero polynomial R w is gcd(F ; G) whenever A is an integral domain.

For all 0 6 i 6 w + 1, we let n i = deg R i , and call (n i ) i the Euclidean degree sequence. This sequence is said to be normal when n i+1 = n i ¡ 1 for all 1 6 i < w. When n i+1 < n i ¡ 1 for some i we say that a degree gap occurs at step i. In addition, it is easy to verify that deg C i < n 1 ¡ n i and deg D i < n 0 ¡ n i by induction on 1 < i 6 w + 1.

Complexity model. For complexities, we shall use computation trees over A with the total complexity point of view. This means that complexity estimates charge a constant cost for each arithmetic operation in A (addition, subtraction, multiplication, and division in our framework) and the equality test. All constants in A are though to be freely at our disposal. See denitions in [START_REF] Bürgisser | Algebraic complexity theory[END_REF]Chapter 4].

A univariate polynomial of degree n is represented by the vector of its n + 1 coecients. We write M: N ! Z for a function that bounds the cost of a polynomial product algorithm in terms of the number of ring operations performed independently of the coecient ring, assuming a unity is available. In other words, two polynomials of degrees at most n over such a ring A may be multiplied with M(n) arithmetic operations in A. The fastest known algorithm, due to Cantor and Kaltofen [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF], provides us with M(n) = O(n log n log log n) = O ~(n). Here, the soft-Oh notation f (n) 2 O ~(g(n)) means that f (n) = g(n) log O (1) g(n) (we refer the reader to [START_REF] Zur Gathen | Modern computer algebra[END_REF]Chapter 25,Section 7] for technical details). In order to simplify cost analyses, we make the customary assumption that n 7 ! M(n) / n is non decreasing, which implies the super-additivity of M, namely M(n 1 ) + M(n 2 ) 6 M(n 1 + n 2 ) for all n 1 > 0 and n 2 > 0.

Over concrete rings or elds, explicitly presented over a nite eld, or Z, we shall use Turing machines endowed with suciently many tapes, in order to benet from standard algorithms. Integers are though to be represented by their binary expansion. Additions take linear time, and we write I(n) for a function that bounds the cost of an integer product in size n, with the same customary assumptions as for M.

Related work

The Euclidean algorithm has a long history in computational mathematics, which goes back to Euclid. Nowadays it is widely used in computer algebra systems for gcds of integers and polynomials, with softly linear time in most common situations. The naive algorithm has quadratic cost in the input size (in the Turing model for integers, and in the computation tree model for univariate polynomials over abstract elds). The key ideas of the fast divide and conquer algorithm are due to Lehmer, for integers [START_REF] Lehmer | Euclid's algorithm for large numbers[END_REF]: Euclidean quotients only depend on higher bits, and their total bit size does not exceed the input size. The rst softly linear cost has been achieved by Knuth [START_REF] Knuth | The analysis of algorithms[END_REF], namely O(I(n) log 4 n) for input size n.

At present time, the best known asymptotic complexity bound is owned by Schönhage [START_REF] Schönhage | Schnelle Berechnung von Kettenbruchentwicklungen[END_REF], namely O(I(n) log n). For polynomials over elds, these algorithms have been adapted rst to normal sequences by Moenck [START_REF] Moenck | Fast computation of GCDs[END_REF], and Aho, Hopcroft, Ullman [START_REF] Aho | The Design and Analysis of Computer Algorithms[END_REF], and then completed by Brent, Gustavson, and Yun [START_REF] Brent | Fast solution of Toeplitz systems of equations and computation of Padé approximants[END_REF], who reached the analogue cost O(M(n) log n) for general sequences and input degree n. We shall refer to these fast variants of the Euclidean algorithms as the half-gcd algorithm. For modern pleasant presentations of this algorithm we refer the reader to classical books [START_REF] Bürgisser | Algebraic complexity theory[END_REF][START_REF] Cohen | A Course in Computational Algebraic Number Theory[END_REF][START_REF] Zur Gathen | Modern computer algebra[END_REF][START_REF] Yap | Fundamental Problems in Algorithmic Algebra[END_REF].

In the polynomial case, the half-gcd algorithm works ne when coecient sizes do not grow during the computation (typically over nite elds). Over integers or polynomial rings, it is well known that intermediate sizes grow very quickly, and straightforward implementations are not practical even in medium sizes. The subresultant theory provides a nice solution to this issue. First, it gives a simple condition to decide if a given degree occurs in the Euclidean degree sequence. Second, it oers polynomial expressions for the coecients of these remainders in terms of the coecients of the input polynomials (up to a suitable renormalisation). These polynomial expressions even turn out to be irreducible over Z [START_REF] Kahoui | An elementary approach to subresultants theory[END_REF]Lemma 4.2]. Informally speaking, this means that subresultant polynomials have generically optimal sizes, and are thus convenient representative of Euclidean remainders.

Subresultant coefficients were introduced by Sylvester in 1840 as determinants of matrices nowadays called Sylvester matrices. The terminology subresultant was coined latter by Bôcher in 1907 [START_REF] Bôcher | Introduction to higher algebra[END_REF] for the subresultant coecients. Straightforward naive computations of the subresultant polynomials by means of their dening determinants lead essentially to a cubic cost, which is far from ecient in general. The rst algorithm with quadratic cost goes back to Habicht [START_REF] Habicht | Eine Verallgemeinerung des Sturmschen Wurzelzählverfahrens[END_REF]: he showed that only the subresultant polynomials S n i and S n i ¡1 are nonzero, and that they may be computed recursively as follows, by means of pseudo-divisions:

S n i+1 = (s n i ¡1;n i n i ¡n i+1 ¡1 S n i ¡1 )/s n i n i ¡n i+1 ¡1
;

(1)

S n i+1 ¡1 = prem(S n i ; S n i ¡1 )/s n i n i ¡n i+1 ; ( 2
)
where s n i ¡1;n i represents the coecient of degree n i in S n i ¡1 . In particular, if all the nonzero subresultant coecients s k are invertible in Frac(A), then each nonzero subresultant polynomial is proportional over Frac(A) to the Euclidean remainder of the same degree (in particular the Euclidean sequence is well dened). The rst use of subresultants in computer algebra is due to Collins [START_REF] Collins | Subresultants and reduced polynomial remainder sequences[END_REF], who coined the subresultant polynomial terminology, and reported on the practical impact. He also showed lower bounds for the expression swell in the Euclidean algorithm. Then, Brown and Traub [START_REF] Brown | On Euclid's algorithm and the theory of subresultants[END_REF], followed by Loos [START_REF] Loos | Generalized polynomial remainder sequences[END_REF], extended Collins' algorithm over unique factorization domains, with easier proofs for the proportionality to the Euclidean remainder sequence. Overall, this led to major practical algorithms. First, over a eld, one may compute whatever polynomial remainder sequence (including using the half-gcd algorithm) in order to deduce subresultants by suitable a posteriori renormalisations. Second, over rings for which the coecient size is an issue, it is often possible to use the specialization property of the subresultant in order to reduce the bulk of computations over nite elds, thanks to the multi-modular approach which was initiated by Collins [START_REF] Collins | The calculation of multivariate polynomial resultants[END_REF] and Brown [START_REF] Brown | On Euclid's algorithm and the computation of polynomial greatest common divisors[END_REF].

For non normal sequences, and when multi-modular techniques are not available, the coecient growth becomes an issue for large gaps in the degree sequence (see our Example 18). In order to make the discussion precise on coecient growth, we consider that A is a polynomial ring K[t], where K is an eective eld, and that the degrees in t of F and G are 6d. Then the degrees in t of the subresultant polynomials are 6 (n 0 + n 1 ) d. The rst part of the solution for formula (1) is due to Lazard: in a manu-script remained unpublished, he showed that s n i ¡1;n i e / s n i e¡1 is in A for 1 6 e 6 n i ¡ n i+1 , and that its degree in t remains 6(n 0 + n 1 ) d (see our Lemma 7). Ducos [START_REF] Ducos | Algorithme de Bareiss, algorithme des sous-résultants[END_REF] found the second part of the solution, for the coecient growth involved in (2): he designed the rst algorithm with a quadratic number of operations in the ground ring, and with a degree in t not exceeding twice the degree bound in t of the subresultants, namely 2 (n 0 + n 1 ) d. At present time, it is not known how to apply the divide and conquer paradigm to his algorithm.

In their 1996 article [START_REF] Th | Cauchy index computation[END_REF], Lickteig and Roy designed an alternative solution for the coefcient growth involved by degree gaps, by limiting this growth to a factor of 3. Their main improvement is an exact polynomial division scheme in A[x] in replacement of pseudodivisions. They appealed to Lazard's optimization, and replaced formula (2) by

(¡1) n i ¡n i+1 ¡1 s n i ¡1;n i+1 s n i+1 S n i = Q n i+1 S n i ¡1 + s n i 2 S n i+1 ¡1 ; (3) with Q n i+1 2 A[x]
. This means that the division of

s n i ¡1;n i+1 s n i+1 S n i by S n i ¡1 is well dened in A[x]
, and this allows to obtain S n i+1 ¡1 as rem((¡1)

n i ¡n i+1 ¡1 s n i ¡1;n i+1 s n i+1 S n i ; S n i ¡1 )/s n i 2 .
We shall refer to formula (3) as the LickteigRoy division. Roughly speaking, the coecients may grow by a factor of 3, but Lickteig and Roy showed how this division process may benet from the half-gcd strategy. However they did not explain how to perform this division eciently over an abstract ring A. They contented themselves to rely on multi-modular techniques for the concrete ring A = Z (see [37, p. 335]).

In his 2001 article [START_REF] Ducos | Optimizations of the subresultant algorithm[END_REF], Ducos proposed an algorithm for performing the LickteigRoy division with a coecient growth of only 2, but with a total quadratic cost. Lombardi, Roy, and Safey El Din [START_REF] Lombardi | New structure theorem for subresultants[END_REF], achieved the same growth and quadratic cost with an alternative algorithm.

When A is a multivariate polynomial ring D[t 1 ; :::; t r ], the most ecient techniques for subresultants rely on fast multi-point evaluation and interpolation algorithms (see [START_REF] Arnold | Faster sparse multivariate polynomial interpolation of straight-line programs[END_REF][START_REF] Van Der Hoeven | Sparse polynomial interpolation in practice[END_REF] for instance for recent advances). These algorithms require D to have suciently many elements, which is not very restrictive in practice. In fact, if necessary, we may often perform computations over an algebraic extension of D. Nevertheless for an abstract domain D, the complexity of subresultant computations is of theoretical interest. If d now represents a bound on the partial degrees in t 1 ; :::; t r of F and G, then the coecient size of the k-th subresultant polynomial is O((n 0 + n 1 ¡ 2 k) r d r ). A growth of a factor of in the partial degrees implies a growth of a factor of r in the coecient sizes. Consequently the algorithms presented so far behave well only when n 0 is suciently larger than d r . Otherwise one may appeal to suitable linear algebra techniques as those designed in [START_REF] Abdeljaoued | Minors of Bezout matrices, subresultants and the parameterization of the degree of the polynomial greatest common divisor[END_REF][START_REF] Abdeljaoued | Bezout matrices, subresultant polynomials and parameters[END_REF][START_REF] Akritas | A new method for computing polynomial greatest common divisors and polynomial remainder sequences[END_REF][START_REF] Akritas | Sylvester's forgotten form of the resultant[END_REF][START_REF] Akritas | Matrix computation of subresultant polynomial remainder sequences in integral domains[END_REF][START_REF] Akritas | Fast matrix computation of subresultant polynomial remainder sequences[END_REF][START_REF] Bini | Fast fraction-free triangularization of Bezoutians with applications to sub-resultant chain computation[END_REF][START_REF] Kerber | Division-free computation of subresultants using Bezout matrices[END_REF]. We shall not investigate this situation in the present article.

For a modern use of Habicht's original ideas, the reader might consult Reischert's article [START_REF] Reischert | Asymptotically fast computation of subresultants[END_REF], which also contains a divide and conquer variant of the subresultant algorithm based on formulas (1) and (2) (thus discarding the coecient growth in defective cases). In an other more recent article in this vein, El Kahoui [START_REF] Kahoui | An elementary approach to subresultants theory[END_REF] proved the LickteigRoy division formula over any commutative ring. Finally, let us mention the 2003 article by von zur Gathen and Lücking [START_REF] Gathen | Subresultants revisited[END_REF], which contains a comprehensive history of the Euclidean algorithm, lower bounds for the coecient swell, and also comparisons between performances of usual implementations with quadratic costs.

Our contributions

Our rst contribution concerns the cost analysis of the LickteigRoy algorithm. In fact we propose to perform the LickteigRoy divisions by a divide and conquer algorithm over A. We show that the quotient Q n i+1 in formula (3) may be obtained with O(M(n i ¡ n i+1 ) log(n i ¡ n i+1 )) operations in A. This is a logarithmic factor higher than by using Newton's iteration when A is a eld. At rst sight, one might thus fear a total cost of O(M(n 0 ) log 2 n 0 ) for the complete subresultant algorithm, but it is not so, because the total contribution of polynomial divisions is O(

P i M(n i ¡ n i+1 ) log (n i ¡ n i+1 )) = O(M(
n 0 ) log n 0 ). Our Section 3.1 concerns fast divisions in A[x]: the presented algorithm is certainly elementary, but we need to recall it for properly analyzing coecient growths.

In Section 2 we prove the LickteigRoy formula from scratch. Our reasons for repeating this proof are twice. First, in their article [START_REF] Th | SylvesterHabicht sequences and fast Cauchy index computation[END_REF], Lickteig and Roy consider the Habicht remainder sequence over an integral domain (this sequence coincides to the subresultant sequence up to signs). Second, they assume n 1 = n 0 ¡ 1, and therefore we have to detail the two rst division steps in the general case. We could not rely on [START_REF] Kahoui | An elementary approach to subresultants theory[END_REF] neither, mostly because of the coecient growth in the second division step (see part 5 of Theorem 9). Algorithm 4 and Theorem 15 constitute our main results for any A, whence our second contribution.

Our third contribution concerns rened complexity results when A is a univariate polynomial ring of the form B[t], where B is a commutative ring with unity endowed with its partially dened division routine. Complexity results are stated in Section 4.4. Then Section 4.5 contains corollaries for the deterministic complexity of bivariate gcds. Our new complexity bounds improve on previously known ones from the asymptotic point of view. Unfortunately, these new bounds do not turn out to be relevant to practice: for computations that last several minutes, Ducos' algorithm or evaluation/interpolation strategies are faster. Nevertheless, for testing purposes, we included an open source implementation of our main algorithm in Mathemagix [START_REF] Van Der Hoeven | Mathemagix User Guide[END_REF][START_REF] Van Der Hoeven | [END_REF].

Subresultants

This section recalls formulas needed for the fast subresultant algorithm presented in the next section. Proofs are established from scratch, on the top of basic linear algebra statements. Recall that our input polynomials are F = P i=0 n 0 f i x i , with degree n 0 , and G = P i=0 n 1 g i x i , with degree n 1 6 n 0 . The A-module of polynomials of degrees <n is written

A[x] <n .

Denitions and main properties

For all 0 6 k < n 1 , the k-th Sylvester map of F and G is dened as:

A[x] <n 1 ¡k A[x] <n 0 ¡k ! A[x] <n 0 +n 1 ¡k (U ; V ) 7 ! U F + V G:
In the canonical monomial basis 1; x; x 2 ; :::, its matrix is

Syl k = 0 B B B B B B B B B B B B @ f 0 g 0 f 0 g n 1 g 0 f n 0 f n 0 g n 1 1 C C C C C C C C C C C C A : (4)
The coecients of F occupy the n 1 ¡ k rst columns, and those of G the n 0 ¡ k last ones. The determinant of the submatrix of Syl k obtained by discarding the k-th rst rows, is called the k-th subresultant coecient of F and G, and is written s k . In other words, s k is the determinant of the map

A[x] <n 1 ¡k A[x] <n 0 ¡k ! A[x] <n 0 +n 1 ¡2k (5) (U ; V ) 7 ! (U F + V G) quo x k :
The determinant s 0 is the classical resultant of F and G, written Res(F ; G). We dene the k-th subresultant polynomial, written S k , by the following determinant over A[x]:

S k = F x n 1 ¡k¡1 F G x n 0 ¡k¡1 G f k+1 f 2k¡n 1 +2 g k+1 g 2k¡n 0 +2 g n 1 f n 0 f n 0 g n 1 ; ( 6
)
where the coecients of F occupy the n 1 ¡ k rst columns, and those of G the n 0 ¡ k last ones. The coecient of degree l in S k , written s k;l , may be obtained as the following determinant over A:

s k;l = f l f l+k¡n 1 +1 g l g l+k¡n 0 +1 f k+1 f 2k¡n 1 +2 g k+1 g 2k¡n 0 +2 g n 1 f n 0 f n 0 g n 1 :
In this way, we see that S k has degree 6k, and that s k is the coecient s k;k of degree k in S k . When s k is zero, the subresultant polynomial S k is said to be defective. Then we introduce the two polynomials U k and V k , called the cofactors of S k :

U k = 1 x n 1 ¡k¡1 0 0 f k+1 f 2k¡n 1 +2 g k+1 g 2k¡n 0 +2 g n 1 f n 0 f n 0 g n 1 ; (7) V k = 0 0 1 x n 0 ¡k¡1 f k+1 f 2k¡n 1 +2 g k+1 g 2k¡n 0 +2 g n 1 f n 0 f n 0 g n 1 : (8)
The coecient of degree l in U k (resp. in V k ) is written u k;l (resp. v k;l ). From these denitions we straightforwardly obtain the Bézout relation

S k = U k F + V k G;
which turns out to be unique whenever s k is invertible in Frac(A), as detailed in the following proposition.

Proposition 1. If s k is invertible in Frac(A), then S k , U k and V k are the unique polyno- mials S ; U ; V in Frac(A)[x] satisfying the following conditions: S = U F + V G, deg S = k, deg U < n 1 ¡ k, deg V < n 0 ¡ k, and s k = lc(S).
Proof. The uniqueness follows from the injectivity of the map dened in ( 5), for which U k ; V k is the unique preimage of (s k ; 0; :::; 0).

We nally introduce the transition matrix

T k = U k V k U k ¡1 V k ¡1 , for k < n 1 , so we have S k S k ¡1 = T k F G .

Highest subresultants and conventions

When n 0 > n 1 , it makes sense to consider the case k = n 1 , for which the Sylvester map is

A[x] <0 A[x] <n 0 ¡n 1 ! A[x] <n 0 (U ; V ) 7 ! U F + V G;
and to use the following conventions:

s n 1 = g n 1 n 0 ¡n 1 , S n 1 = g n 1 n 0 ¡n 1 ¡1 G, U n 1 = 0 and V n 1 = g n 1 n 0 ¡n 1 ¡1 . When n 0 = n 1 , it is sometimes convenient to use the same values for s n 1 ,
S n 1 , U n 1 and V n 1 notwithstanding that S n 1 and V n 1 are not dened over A but over Frac(A), as long as g n 1 is not a zero divisor in A.

When n 1 > 1 and k = n 1 ¡ 1, the Sylvester map,

A[x] <1 A[x] <n 0 ¡n 1 +1 ! A[x] <n 0 +1 (U ; V ) 7 ! U F + V G; has for matrix 0 B B B B B B B B @ f 0 g 0 g 0 g n 1 f n 0 g n 1 1 C C C C C C C C A :
Elementary calculations then yield

S n 1 ¡1 = prem(F ; G): (9)
We introduce

Q n 1 = pquo(F ; G), so we have n 2 = deg S n 1 ¡1 < n 1 .
We shall dene the other quotients Q n i latter on. Notice that we may also consider the rst transition matrix T n 1

T n 1 = 0 s n 1 / g n 1 g n 1 n 0 ¡n 1 +1 ¡Q n 1 ! ;
but it is not dened over A in general, although det T n 1 = ¡s n 1 2 actually belongs to A.

Proportionality to the Euclidean sequence

If the Euclidean sequence is well dened over Frac(A), if R i has degree n i > 0 for some i > 2, and if s n i is invertible in Frac(A), then R i is proportional to S n i over Frac(A), by Proposition 1. The following lemma is the key step to establish the converse: every subresultant polynomial is proportional to a remainder in the Euclidean sequence.

Lemma 2. If n 1 > 0 and if g n 1 is invertible in Frac(A)
, then R 2 is well dened over Frac(A), and the following properties hold:

1. When n 1 > 1, we have S n 1 ¡1 = prem(F ; G) = g n 1 n 0 ¡n 1 +1 R 2 ; 2. When n 2 < k < n 1 ¡ 1, we have S k = 0; 3. When n 2 > 0, we have S n 2 = lc(R 2 ) n 1 ¡n 2 ¡1 g n 1 n 0 ¡n 2 R 2 ; 4. When 0 6 k < n 2 , we have S k = (¡1) (n 1 ¡k)(n 2 ¡k) g n 1 n 0 ¡n 2 S ^k, where S ^k is the k-th polynomial subresultant of G and R 2 , over Frac(A).
Proof. Let us write r 2;j for the coecient of degree j in R 2 . When n 2 6 k < n 1 , interpreting the division of F by G in terms of column operations over Frac(A), we observe that

S k = R 2 x R 2 x k ¡n2 R 2 x k ¡n2+1 R 2 x n1¡k ¡1 R 2 G x n0¡k ¡1 G 0 0 r 2;n2 r 2;2k ¡n1+2 g k+1 g 2k ¡n0+2 0 0 r 2;n2 0 0 0 g n1 g n1
;

For k = n 1 ¡ 1, we recover the formula S n 1 ¡1 = g n 1 n 0 ¡n 1 +1 R 2 , already seen in equation ( 9).

Then we obtain S k = 0 when n 2 < k < n 1 ¡ 1, and

S n 2 = r 2;n 2 n 1 ¡n 2 ¡1 g n 1 n 0 ¡n 2 R 2 . For k < n 2 , the matrix Syl k is column equivalent to 0 B B B B B B B B B B B B B B B B B B B @ r 2;0 g 0 n 0 ¡ n 2 ! g 0 g 0 r 2;n2 r 2;0 g 0 g n1 r 2;n2 g n1 g n1 ¡ n 1 ¡ k ¡! ¡ n 2 ¡ k ¡! g n1 1 C C C C C C C C C C C C C C C C C C C A ; which implies that S k = (¡1) (n 1 ¡k)(n 2 ¡k) g n 1 n 0 ¡n 2 S ^k.
Corollary 3. If g n 1 and all the nonzero s k are invertible in Frac(A), then the Euclidean remainder sequence is well dened over Frac(A), and the following properties hold: 1. S n i ¡1 and S n i+1 are proportional to R i+1 over Frac(A), when 1 6 i < w;

2. S k = 0, when n i+1 < k < n i ¡ 1 and 1 6 i < w;

3. S k = 0, when k < n w .
Proof. The proof is done by induction on w. The case w = 1 means that R 2 = 0, and the previous lemma implies that all the subresultant polynomials are zero. Now assume w > 2.

From the previous lemma, we know that: R 2 is well dened over Frac(A), R 2 is proportional to S n 1 ¡1 and S n 2 over Frac(A), the leading coecient of R 2 is invertible in Frac(A), S k = 0 for all n 2 < k < n 1 ¡ 1, and the other S k are proportional to the subresultant polynomials of G and R 2 . Since the length of the Euclidean sequence of G and R 2 is w ¡ 1, we may use the corollary by induction.

Corollary 4. If g n 1 and all the nonzero s k are invertible in Frac(A), then, for 2 6 i < w the transition matrix

T n i is the unique matrix T = T 1;1 T 1;2 T 2;1 T 2;2 ! over Frac(A)[x] satisfying
the following properties:

S n i S n i ¡1 ! = T F G , deg T 1;1 < n 1 ¡ n i , deg T 1;2 < n 0 ¡ n i , deg T 2;1 < n 1 ¡ n i+1 , deg T 1;2 < n 0 ¡ n i+1 .
Proof. Proposition 1 already implies T 1;1 = U n i and T 1;2 = V n i . By the previous corollary, S n i ¡1 is proportional to S n i+1 , and the conclusion follows again from Proposition 1 since deg

U n i ¡1 < n 1 ¡ (n i ¡ 1) 6 n 1 ¡ n i+1 and deg V n i ¡1 < n 0 ¡ (n i ¡ 1) 6 n 0 ¡ n i+1 .
About the uniqueness of the cofactors U n w ¡1 and V n w ¡1 of S n w ¡1 = 0, we begin with the following simple lemma: Lemma 5. For all 1 6 k 6 n 1 , the following identities hold:

u k¡1;n 1 ¡k = (¡1) n 1 ¡k g n 1 s k ; v k¡1;n 0 ¡k = (¡1) n 1 ¡k¡1 f n 0 s k :
Proof. These formulas are straightforward, from expanding determinants of equations [START_REF] Akritas | Fast matrix computation of subresultant polynomial remainder sequences[END_REF] and [START_REF] Arnold | Faster sparse multivariate polynomial interpolation of straight-line programs[END_REF] along their rst rows. Corollary 6. Assume w > 2, and that f n 0 , g n 1 and all the nonzero s k are invertible in Frac(A). Then, the cofactors U n w ¡1 and V n w ¡1 are the unique polynomials A and B of Frac(A)[x] satisfying the following properties:

A F + B G = 0, deg A = n 1 ¡ n k , deg B = n 0 ¡ n k , and lc(A) = (¡1) n 1 ¡n w g n 1 s n w .
Proof. Let a j (resp. b j ) represent the coecient of degree j in A (resp. in B). We consider the relation

(A ¡ U n w ¡1 ) F + (B ¡ V n w ¡1 ) G = 0: By Lemma 5 we have deg(A ¡ U n w ¡1 ) < n 1 ¡ n k . On the other hand, using a n 1 ¡n k f n 0 + b n 0 ¡n k g n 1 = 0, we observe that (¡1) n 1 ¡n w ¡1 b n 0 ¡n w V n w ¡1 /(f n 0 s n w ) = (¡1) n 1 ¡n w a n 1 ¡n w V n w ¡1 /(g n 1 s n w ) = V n w ¡1 : Using Lemma 5 again, we obtain that deg(B ¡ V n w ¡1 ) < n 0 ¡ n k . Since the map A[x] <n 1 ¡n w A[x] <n 0 ¡n w ! A[x] <n 0 +n 1 ¡2n w (U ; V ) 7 ! (U F + V G) quo x n w
is injective, we deduce that A = U n w ¡1 and B = V n w ¡1 .

Structure theorem

We are now ready to prove Lazard's lemma and the LickteigRoy division (equation ( 3)).

Lemma 7. Assume that f n 0 , g n 1 , and all the nonzero s k are invertible in Frac(A), and let 1 6 i < w. Then, the following properties hold:

1. When n i+1 6 k < n i , the ratio s n i ¡1;k n i ¡k /s n i n i ¡k¡1 equals a minor of size n 0 + n 1 ¡ 2 k of Syl k , hence it belongs to A (Lazard's lemma); 2. S n i+1 = (s n i ¡1;n i+1 n i ¡n i+1 ¡1 S n i ¡1 )/s n i n i ¡n i+1 ¡1
.

Proof. From Corollary 3, we know that S n i ¡1 has degree n i+1 . Thanks to Lemma 5, over Frac(A) we may write

S ~ni ¡1 := (¡1) n 1 ¡n i ¡1 S n i ¡1 f n 0 s n i = U ~ni ¡1 F + V ~ni ¡1 G; ( 10
)
where

U ~ni ¡1 = (¡1) n 1 ¡n i ¡1 U n i ¡1 / (f n 0 s n i ), and V ~ni ¡1 = (¡1) n 1 ¡n i ¡1 V n i ¡1 / (f n 0 s n i ), so V ~ni ¡1 is monic of degree n 0 ¡ n i . Assume n i+1 6 k 6 n i .
The Bézout relation [START_REF] Berthomieu | Reduction of bivariate polynomials from convex-dense to dense, with application to factorizations[END_REF] rephrases in terms of the following column operations on the matrix (4) dening Syl k : for l from n i ¡ k ¡ 1 down to 0, the column x n 0 ¡n i +l G(x) may be replaced by x l S ~ni ¡1 (x) by adding to it

X j =0 n 1 ¡n i u ~ni ¡1;j x j +l F (x) + X j =0 n 0 ¡n i ¡1 v ~ni ¡1;j x j+l G(x);
where s ~ni ¡1;j , u ~ni ¡1;j , v ~ni ¡1;j are the coecients of degree j in S ~ni ¡1 , U ~ni ¡1 , V ~ni ¡1 respectively.

The Sylvester matrix Syl k is thus column equivalent over Frac(A) to

f 0 n i ¡ k ! g 0 s ~ni¡1;0 s ~ni¡1;0 f ni+1¡1 s ~ni¡1;ni+1¡1 f ni+1 g ni+1 s ~ni¡1;ni+1 g 0 f ni+1+ni¡k ¡1 g ni+1+ni¡k ¡1 s ~ni¡1;ni+1 f ni+1+ni¡k g ni+1+ni¡k f 0 f ni¡1 f 2ni¡n1 g ni¡1 g 2ni¡n0 f ni f 2ni¡n1+1 g ni g 2ni¡n0+1 g n1 f n0 f n0 g n1 f n0 n 1 ¡ n i ! f n0 n 0 ¡ n i ! n i ¡ k ! 1 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C A :
When n i+1 6 k < n i , we discard the n i+1 rst grayed rows and the other k ¡ n i+1 grayed rows corresponding to the monomials x n i+1 +n i ¡k ; :::; x n i ¡1 , and obtain the determinant

f n 0 n i ¡k s n i s ~ni ¡1;n i+1 n i ¡k = f n 0 n i ¡k s n i s n i ¡1;n i+1 f n 0 s n i n i ¡k = s n i ¡1;n i+1 n i ¡k s n i n i ¡k¡1 ;
which equals the minor of Syl k obtained by deleting the same rows. This concludes part 1.

When k = n i+1 this minor coincides to s n i+1 = s n i ¡1;n i+1 n i ¡n i+1 /s n i n i ¡n i+1 ¡1 , whence part 2, since S n i+1 is proportional to S n i ¡1 by Corollary 3.

Proposition 8. For all 2 6 i 6 w, we have

det T n i = U n i V n i U n i ¡1 V n i ¡1 = (¡1) n 1 ¡n i ¡1 s n i 2 ; (¡1) n 1 ¡n i ¡1 s n i 2 T n i ¡1 = V n i ¡1 ¡V n i ¡U n i ¡1 U n i :
Proof. By left multiplying booth sides of

S n i S n i ¡1 = U n i V n i U n i ¡1 V n i ¡1 F G by ¡ V n i ¡1 ¡V n i , we obtain V n i ¡1 S n i ¡ V n i S n i ¡1 = F det T n i . By Proposition 1, we have deg(V n i S n i ¡1 ) 6 n 0 ¡ 2, and deg(V n i ¡1 S n i ) 6 n 0 . Consequently det T n i is in A, and we have v n i ¡1;n 0 ¡n i s n i = f n 0 det T n i .
The conclusion follows from Lemma 5. Theorem 9. Let A be a commutative ring, and let F, G be polynomials of respective degrees n 0 > n 1 > 0. We assume that f n 0 , g n 1 , and all the nonzero subresultant coecients are invertible in Frac(A). Then we have:

1.

S n 1 ¡1 = prem(F ; G), Q n 1 = pquo(F ; G), when n 1 > 1;
2. S k = 0 when n i+1 < k < n i ¡ 1 and 1 6 i 6 w; ¡2 2 A, and when n i+1 < n i ¡ 1;

3. s n i ¡1;n i+1 n i ¡k¡1 /s n i n i ¡k¡2 equals a minor of size n 0 + n 1 ¡ 2 k of Syl k , when n i+1 6 k < n i ¡ 1; 4. s n i S n i+1 = n i S n i ¡1 , where n i = s n i ¡1;n i+1 n i ¡n i+1 ¡1 /s n i n i ¡n i+1

The division of s

n 1 ¡1;n 2 s n 2 G by S n 1 ¡1
is well dened over A, and we have

(¡1) n 1 ¡n 2 ¡1 s n 1 ¡1;n 2 s n 2 G = Q n 2 S n 1 ¡1 + g n 1 s n 1 S n 2 ¡1 ; with Q n 2 2 A[x];
6. When i > 2, the division of s n i ¡1;n i+1 s n i+1 S n i by S n i ¡1 is well dened over A, and we have

(¡1) n i ¡n i+1 ¡1 s n i ¡1;n i+1 s n i+1 S n i = Q n i+1 S n i ¡1 + s n i 2 S n i+1 ¡1 ; with Q n i+1 2 A[x]:
Proof. Parts 1 to 4 have already been proved. As for part 5, from Corollary 3 and part 2 of Lemma 7, there exist ' invertible in Frac(A), and

2 Frac(A)[x] of degree n 1 ¡ n 2 , such that S n 2 S n 2 ¡1 = 0 n 1 /s n 1 ' (¡1) n 1 ¡n 2 G S n 1 ¡1
;

which implies S n 2 S n 2 ¡1 = 0 n 1 /s n 1 ' (¡1) n 1 ¡n 2 0 1 g n 1 s n 1 ¡Q n 1 F G = g n 1 n 1 ¡ n 1 Q n 1 /s n 1 (¡1) n 1 ¡n 2 g n 1 s n 1 ' ¡ (¡1) n 1 ¡n 2 Q n 1 F G by using S n 1 ¡1 = prem(F ; G) = g n 1 s n 1 F ¡ Q n 1 G. Assume w > 2. Since deg(g n 1 n 1 ) = 0 < n 1 ¡ n 2 , deg Q n 1 = n 0 ¡ n 1 < n 0 ¡ n 2 , deg < n 1 ¡ (n 2 ¡ 1), deg (Q n 1 ) = n 0 ¡ n 2 < n 0 ¡ (n 2 ¡ 1)
, we may apply Corollary 4 to obtain

T n 2 = g n 1 n 1 ¡ n 1 Q n 1 /s n 1 (¡1) n 1 ¡n 2 g n 1 s n 1 ' ¡ (¡1) n 1 ¡n 2 Q n 1 : (11)
In particular, Proposition 8 gives det

T n 2 = (¡1) n 1 ¡n 2 ¡1 s n 2 2 = ' g n 1 n 1 , which leads to ' = (¡1) n 1 ¡n 2 ¡1 s n 2 2 / (g n 1 n 1 ) = (¡1) n 1 ¡n 2 ¡1 s n 1 ¡1;n 2 s n 2 / (g n 1 s n 1 )
. On the other hand g n 1 s n 1 belongs to A[x], which gives part 5, with

Q n 2 = (¡1) n 1 ¡n 2 +1 g n 1 s n 1 .
If w = 2, then S n 2 ¡1 = 0, and we may divide ' and by lc() without loss of generality. Corollary 6 implies (¡1)

n 1 ¡n 2 g n 1 s n 1 = U n 2 ¡1 and ' ¡ (¡1) n 1 ¡n 2 Q n 1 = V n 2 ¡1
, which again yields [START_REF] Bini | Fast fraction-free triangularization of Bezoutians with applications to sub-resultant chain computation[END_REF], and the conclusion follows as in the case w > 2.

Now for part 6, from Corollary 3 and part 2 of Lemma 7, again, there exist invertible in Frac(A), and 2 Frac(A)[x] of degree n i ¡ n i+1 , such that

S n i+1 S n i+1 ¡1 ! = 0 n i /s n i (¡1) n 1 ¡n i+1 ! S n i S n i ¡1
; which implies

S n i+1 S n i+1 ¡1 ! = 0 n i /s n i (¡1) n 1 ¡n i+1 ! T n i F G = n i U n i ¡1 /s n i n i V n i ¡1 /s n i U n i + (¡1) n 1 ¡n i+1 U n i ¡1 V n i + (¡1) n 1 ¡n i+1 V n i ¡1 ! F G : Assume w > i + 2. Since deg U n i ¡1 < n 1 ¡ (n i ¡ 1) 6 n 1 ¡ n i+1 , deg V n i ¡1 < n 0 ¡ (n i ¡ 1) 6 n 0 ¡ n i+1 , deg( U n i + U n i ¡1 ) < max (n 1 ¡ n i ; n 1 ¡ (n i+1 ¡ 1)) 6 n 1 ¡ (n i+1 ¡ 1); and deg( V n i + V n i ¡1 ) < max (n 0 ¡ n i ; n 0 ¡ (n i+1 ¡ 1)) 6 n 0 ¡ (n i+1 ¡ 1);
we may apply Corollary 4 again to obtain

T n i+1 = n i U n i ¡1 /s n i n i V n i ¡1 /s n i U n i + (¡1) n 1 ¡n i+1 U n i ¡1 V n i + (¡1) n 1 ¡n i+1 V n i ¡1 ! : (12)
In particular Proposition 8 gives us

det T n i+1 = (¡1) n 1 ¡n i+1 ¡1 s n i+1 2 = ¡( n i /s n i ) det T n i = (¡1) n 1 ¡n i s n i 2 s n i+1 /s n i ¡1;n i+1 ; which leads to = (¡1) n i ¡n i+1 ¡1 s n i ¡1;n i+1 s n i+1 /s n i 2 . From ( ) T n i = ¡ U n i+1 ¡1 V n i+1 ¡1 ; ( 13
U n i + (¡1) n 1 ¡n i+1 U n i ¡1 = c U n i+1 ¡1 and V n i + (¡1) n 1 ¡n i+1 V n i ¡1 = c V n i+1 ¡1
for some element c invertible in Frac(A). Without loss of generality we may divide and by c so that (12) still holds. The conclusion then follows as in the case w > i + 2.

Atomic transition matrices

We introduce the following atomic transition matrices M n i over Frac(A)[x]:

M n 1 = 0 s n 1 / g n 1 g n 1 s n 1 ¡Q n 1 ! ; M n i+1 = 0 n i s n i (¡1) n i ¡n i+1 ¡1 s n i ¡1;n i+1 s n i+1 ¡Q n i+1 ! /s n i 2 ; when 1 6 i < w; so S n 1 S n 1 ¡1 ! = M n 1 F G
, and by Theorem 9 we have

S n i+1 S n i+1 ¡1 ! = M n i+1 S n i S n i ¡1
; when 1 6 i < w:

Notice that s n i 2 M n i+1 has all its entries in A[x], while the entry s n 1 / g n 1 in M n 1 does not necessarily belong to A. In addition, we have the following useful formulas:

det M n 1 = ¡s n 1 2 ; det M n i+1 = (¡1) n i ¡n i+1 s n i+1 2 /s n i 2 ; T n i = M n i M n 1 ; T n 2 = g n 1 n 1 ¡( n 1 Q n 1 )/s n 1 ¡Q n 2 ((¡1) n 1 ¡n 2 +1 s n 1 ¡1;n 2 s n 2 + Q n 1 Q n 2 )/(g n 1 s n 1 ) ! : (14)
By construction, T n 2 has all its entries in A[x].

Lemma 10. When 1 6 i < j 6 w, the product

s n i 2 M n j M n i+1 = T n j (s n i 2 T n i ¡1 ) = (¡1) n 1 ¡n i ¡1 U n j V n j U n j ¡1 V n j ¡1 ! V n i ¡1 ¡V n i ¡U n i ¡1 U n i
has all its entries in A[x], the matrix of its degrees is

n i+1 ¡ n j ¡1 n i ¡ n j ¡1 n i+1 ¡ n j n i ¡ n j !
, with the convention that a negative degree (in the top left entry) means the zero polynomial.

Proof. Since T n j = M n j M n 1 = M n j M n i+1 T n i has all its entries in A[x], so has s n i 2 M n j M n i+1 , by Proposition 8. When j = i + 1 the degree matrix of M n i+1 is <0 0 0 n i ¡ n i+1 . Then the degree matrix of M n i+2 M n i+1 is 0 n i ¡ n i+1 n i+1 ¡ n i+2 n i ¡ n i+2
, and the one of

M n i+3 M n i+2 M n i+1 is n i+1 ¡ n i+2 n i ¡ n i+2 n i+1 ¡ n i+3 n i ¡ n i+3
. The conclusion follows easily by induction.

Half subresultant algorithm

The key idea of the half-gcd algorithm is the computation of the atomic transition matrices by a divide and conquer approach. In our framework over a ring endowed with its partial division routine, we rst need to describe how well dened divisions in A[x] may be performed fast. Let A = P i>0 a i x i be a polynomial in A[x], and let n > 0, we shall use the upper and lower truncations written dAe n := P i=0 n¡1 a i x i and bAc n := P i>0 a i+n x i .

Fast polynomial division

Let A and B be two polynomials in A[x] such that the division of A by B is well dened.

In the same vein as the Barrett or SievekingKung division algorithms, we rst compute quotients and then remainders. Quotients are obtained from jet computations. Precisely, if and ¡ are two jets in AJxK/(x n ), the division /¡ is said to be well dened whenever ¡ = / 0 and there exists in AJxK/(x n ) such that = ¡ . If a > 0 is a real number, then bac represents the largest integer 6a, and dae is the smallest integer >a.

Algorithm 1

Input. = P i=0 n¡1 ' i x i and ¡ = P

i=0 n¡1 i x i in AJxK/(x n ).
Output. /¡. Assumption: the division of by ¡ is well-dened, and 0 = / 0.

1. If n = 1 then return ' 0 / 0 . 2. Let h = dn/2e. 3.
Call recursively the algorithm with de h and d¡e h to obtain l = de h / d¡e h in AJxK/(x h ).

Call recursively the algorithm with bc h ¡ b l d¡e h c h ¡ l b¡c h and d¡e h in

AJxK/(x n¡h ), and write h the jet in return.

Return

l + h x h , seen in AJxK/(x n ). Proposition 11. Algorithm 2 is correct and takes O(M(n) log n) operations in A.
Proof. Let represent the quotient / ¡ in AJxK/ (x n ). First, it is clear that the division in step 3 is well-dened and that l = d e h . Then, using 2 h > n, we obtain

= d e h d¡e h + (d e h b¡c h + d¡e h b c h ) x h + O(x n )
which yields the correctness of step 4:

h = bc h ¡ bd e h d¡e h c h ¡ d e h b¡c h d¡e h + O(x n¡h ) = b c h + O(x n¡h ):
The correctness of the algorithm follows by strong induction on n. Its cost function

C(n) satises C(n) = C(h) + C(n ¡ h) + O(M(n)), which classically leads to C(n) = O(M(n) log n).
For a polynomial A 2 A[x] of degree 6n, we write rev(A; n) for the reverse polynomial x n A(1/x). The latter algorithm classically allows to compute polynomial quotients.

Proof. Let R = rem(A; B), Q = quo(A; B), Q ~= rev(Q; l), and R ~= rev(R; n). From A = Q B + R we deduce A ~= Q ~B ~+ R ~, and then Q ~= A ~/ B ~+ O(x l+1
). This proves the correctness. The complexity simply follows from the latter proposition.

Once the quotient Q of A by B is computed, then R = A ¡ Q B may be obtained with M(n) operations in A.

Main divide and conquer routine

Theorem 9 gives formulas to compute

S n 1 ¡1 , Q n 1 , S n 2 , Q n 2 ,
Q n 3 +deg Q n 4 ) and bS n 2 ¡1 c n 4 ¡(deg Q n 3 +deg Q n 4 )
. By induction, we thus prove that Q n 3 ; :::; Q n i may be obtained from bS

n 2 c n i ¡(deg Q n 3 ++deg Q n i ) and bS n 2 ¡1 c n i ¡(deg Q n 3 ++deg Q n i ) . By using n i = n 2 ¡ (n 2 ¡ n 3 ) ¡ ¡ (n i¡1 ¡ n i ) = n 2 ¡ (deg Q n 3 + + deg Q n i ); the low truncation order n i ¡ (deg Q n 3 + + deg Q n i ) rewrites into n 2 ¡ 2 (deg Q n 3 + + deg Q n i ).
If we x an integer l > 0, in order to compute Q n 3 ; :::; Q n j with j maximal such that deg Q n 3 + + deg Q n j 6 l, it suces to low truncate S n 2 and S n 2 ¡1 at order n 2 ¡ 2 l. In order to handle negative truncation order, it could be convenient to use Laurent polynomials from a theoretical point of view, but for practice we prefer to keep computations in A[x]. In fact the construction we have just sketched simply works ne by setting bAc n = A for all n < 0, thanks to the following lemma: Proof. If l < 0 then bB c l = B and the lemma is correct. Now assume l > 0. We need to prove that the coecient P p+q=i;q>0 a p b q+l of degree i in A bB c l equals the coecient of degree i + l in A B for all i > deg A. This assertion is correct because p 6 deg A implies q = i ¡ p > 0.

Let i > 2 and 0 6 l 6 n i . If we are given bS n i c and bS n i ¡1 c , where = n i ¡ 2 l, then the same reasoning shows that we may compute Q n i+1 ; :::; Q n j with j maximal such that n i ¡ n j 6 l. In fact, we shall better compute the subresultants s n i ; :::; s n j ¡1 and the numerators N n i+1 ; :::; N n j of the atomic transition matrices made from Q n i+1 ; :::; Q n j respectively, namely

N n i+1 = s n i 2 M n i+1 .
At the same time we shall return the composite transition matrix M n i+1 ;n j = M n j M n i+1 , or more precisely its numerator N n i+1 ;n j = s n i 2 M n i+1 ;n j , which has coecients in A[x] according to Lemma 10. We may then recover S n j and S n j ¡1 as

S n j S n j ¡1 ! = N n i+1 ;n j S n i S n i ¡1 /s n i 2 :
Upon this strategy we are now able to present the adaptation of the half-gcd algorithm to subresultants. But before, notice that given l and bAc l one may deduce deg A whenever bAc l = / 0: this is simply deg bAc l if l < 0, or l + deg bAc l otherwise.

Algorithm 3

Input. n 0 ; :::; n i , an integer l 2 f0; :::; n i g, and bS n i c and bS n i ¡1 c , where = n i ¡ 2 l.

Output. n i+1 ; :::; n j , s n i+1 ; :::; s n j , N n i+1 ; :::; N n j , and N n i+1 ;n j with j > i maximal such that n i ¡ n j 6 l.

Assumptions: i > 2, f n 0 , g n 1 ,
k ¡1 c = / 0. 5. If n k+1 < n k ¡ 1 then compute n k = s n k ¡1;n k+1 n k ¡n k+1 ¡1 / s n k n k ¡n k+1 ¡2
, and bS

n k+1 c = ( n k bS n k ¡1 c )/s n k , otherwise set n k = s k . 6. Compute Q n k+1 from (¡1) n k ¡n k+1 ¡1 s n k ¡1;n k+1 s n k+1 bS n k c and bS n k ¡1 c , by means of Algorithm 2. 7. Build N n k +1 = 0 @ 0 n k s n k (¡1) n k ¡n k+1 ¡1 s n k ¡1;n k+1 s n k+1 ¡Q n k+1 1 A . 8. Let = + n k ¡ n k+1 , and deduce bS n k+1 ¡1 c = j (¡ 1) n k ¡n k+1 ¡1 s n k ¡1;n k+1 s n k+1 S n k ¡ Q n k+1 S n k ¡1 k /s n k
2 by means of Lemma 13. 9. Call recursively the algorithm with n 0 ; :::; n k+1 , l ¡ (n i ¡ n k+1 ), bS n k+1 c and bS n k+1 ¡1 c . Let n k+1 ; :::; n j , s n k+2 ; :::; s n j , N n k+2 ; :::; N n j , and N n k+2 ;n j be the data obtained in return.

10. Return n i+1 ; :::; n j , s n i+1 ; :::; s n j , N n i+1 ; :::; N n j and (N n k+2 ;n j ((

N n k +1 N n i+1 ;n k ) / s n k 2 ))/s n k+1 2 .
Proposition 14. Algorithm 3 is correct and takes O(M(l) log l) operations in A.

Proof. The proof is done by strong induction on l. If l = 0 then = n i , bS n i ¡1 c = 0, and the output is correct. Now assume l > 1 and that the algorithm is correct up to l ¡ 1.

In step 1 the case bS n i ¡1 c = 0 means that > n i+1 , whence n i ¡ n i+1 > 2 l > l, and the output is correct. Otherwise we obtain n i+1 , and the output is again correct whenever l < n i ¡ n i+1 .

In step 2 we necessarily have 0 6 h < l, which implies that the recursive call is valid and returns a correct result by induction. By Lemma 10 the degrees of the entries of N n i+1 ;n k are 6n i ¡ n k . Lemma 13 ensures that we may safely obtain bS n k c and bS n k ¡1 c in step 3.

If bS n k ¡1 c = 0 in step 4, then this means that > n k+1 , whence

n i ¡ n k+1 > n i ¡ = n i ¡ ¡ (n i ¡ n k ) = 2 l ¡ (n i ¡ n k ) > l
and the output is correct. If bS n k ¡1 c = / 0 then we may determine n k+1 so the output is also correct when l < n i ¡ n k+1 .

The computation in step 5 follows from part 4 of Theorem 9. In step 6, we have n i ¡ n k+1 6 l, and therefore, from

= n i ¡ 2 l, we obtain = + n i ¡ n k = 2 n i ¡ n k ¡ 2 l 6 n k+1 ¡ (n k ¡ n k+1
) so we may safely obtain Q n k+1 by part 6 of Theorem 9.

In step 9, we verify the requested conditions to the recursive call: = n k+1 ¡ 2 (l ¡ (n i ¡ n k+1 )), and l ¡ (n i ¡ n k+1 ) < l.

In step 10,

s n i 2 M n i+1 ;n k +1 = s n i 2 (N n k +1 / s n k 2 ) (N n i+1 ;n k / s n i 2 ) = (N n k +1 N n i+1 ;n k ) / s n k 2
has its entries in A[x] by Lemma 10. Repeating this argument we obtain

s n i 2 M n i+1 ;n j = (N n k+2 ;n j ((N n k +1 N n i+1 ;n k )/s n k 2 ))/s n k+1 2
. We are done with the correctness. Let C(l) represent the cost of the algorithm called with l. Steps 1, 4, 5 and 7 involve O(l) operations in A.

Step 2 costs C(h). Steps 3, 8 and 10 perform O(M(l)) operations in A. In step 6, we appeal to Proposition 12 to get a cost

O(M(n k ¡ n k+1 ) log(n k ¡ n k+1 )). Step 9 amounts to C(l ¡ (n i ¡ n k+1 )): notice that n i ¡ n k+1 > h + 1 implies l ¡ (n i ¡ n k+1 ) 6 l ¡ (h + 1) < l /2. Overall there exists a constant c > 0 such that C(l) 6 C(h) + C(l ¡ (n i ¡ n k+1 )) + c (M(l) + M(n k ¡ n k+1 ) log(n k ¡ n k+1 )): Since h and l ¡ (n i ¡ n k+1 ) are 6l /2, it is classical to deduce C(l) = O M(l) log l + X k=i j ¡1 M(n k ¡ n k+1 ) log(n k ¡ n k+1 ) ! = O(M(l) log l):

Top level algorithm

We are now ready to present the main algorithm for computing all the atomic transition matrices.

Algorithm 4

Input. F and G in A[x] of respective degrees n 0 > n 1 .

Output. n 0 ; :::; n w , s n 1 ; :::; s n w , N n 3 ; :::

; N n w , S n 1 ¡1 , Q n 1 , n 1 , Q n 2 .
Assumption: f n 0 , g n 1 , s n 1 ; :::; s n w are nonzero divisors in A.

1. If F = 0 then w = ¡1 and return nothing none of the quantities are dened.

2. If G = 0 then w = 0 and return n 0 the other quantities are not dened.

Compute S

n 1 ¡1 = prem(F ; Q), Q n 1 = pquo(F ; G), and s n 1 = g n 1 n 0 ¡n 1 . 4. If S n 1 ¡1 = 0 then w = 1 and return n 0 ; n 1 , s n 1 , S n 1 ¡1 , Q n 1 the other quantities are not dened. 5. Let n 2 = deg S n 1 ¡1 . If n 2 < n 1 ¡ 1 then compute S n 2 = ( n 1 S n i ¡1 ) / s n 1 , where n 1 = s n 1 ¡1;n 2 n 1 ¡n 2 ¡1 /s n 1 n 1 ¡n 2 ¡2
, otherwise let n 1 = s n 1 .

Perform the division

(¡1) n 1 ¡n 2 ¡1 s n 1 ¡1;n 2 s n 2 G = Q n 2 S n 1 ¡1 + g n 1 s n 1 S n 2 ¡1 , in order to obtain Q n 2 and S n 2 ¡1 . 7. If S n 2 ¡1 = 0 then w = 2 and return n 0 ; n 1 ; n 2 , s n 1 ; s n 2 = lc(S n 2 ), S n 1 ¡1 , Q n 1 , n 1 ,
Q n 2 the other quantities are not dened. 8. Call Algorithm 3 with n 0 ; n 1 ; n 2 , l = n 2 , bS n 2 c and bS n 2 ¡1 c , where = ¡n 2 . Let n 3 ; :::; n j , s n 3 ; :::; s n j , N n 3 ; :::; N n j , and N n 3 ;n j represent the data obtained in return.

9. Return n 0 ; :::; n j , s n 1 ; :::; s n j , N n 3 ; :::

; N n j , S n 1 ¡1 , Q n 1 , n 1 , Q n 2 .
Theorem 15. Algorithm 4 is correct and takes O(M(n 0 ) log n 0 ) operations in A.

Proof. The correctness follows from Theorem 9 and Proposition 14 after noticing that j necessarily coincides to w in step 8. The divisions in steps 2 and 5 take O(M(n 0 ) log n 0 ) by Proposition 12. The cost of step 8 is given in Proposition 14.

Example 16. Let us briey illustrate the algorithm with

A = Z, n 0 = 5, n 1 = 4, F = ¡3 x 5 ¡5 x 4 +3 x 3 ¡2 x 2 +4 x+2; G = ¡5 x 4 + 4 x 3 ¡2 x 2 +3 x ¡1: We obtain S n 1 ¡1 = ¡43 x 3 ¡ 21 x 2 + 4 x + 87, n 2 = 3,

Bivariate case

In this section we study the complexity of Algorithm 4 when A is a polynomial ring

B[t],
where B is a commutative ring with unity and endowed with its partially dened division routine. First, we illustrate the coecient growth issue.

Coecient growth

The sizes of the coecients of the subresultant polynomials may be bounded from their dening determinant. Proof. This follows from the denitions in equations ( 6), [START_REF] Akritas | Fast matrix computation of subresultant polynomial remainder sequences[END_REF], and (8).

Example 18.

Let A = Z[t]. We consider the following family of polynomials parametrized by the integer m > 2:

F = t x 3m ¡ t x 2m + x m ¡ 1; G = t x 3m + t x 2m + x m :
We have n 0 = n 1 = 3 m. The convention for S n 1 and formula (9) lead to:

S n 1 = t ¡1 G = x 3m + x 2m + t ¡1 x m 2 t ¡1 Z[t]; S n 1 ¡1 = prem(F ; G) = t F ¡ t G = ¡ 2 t 2 x 2m ¡ t:
We obtain n 2 = 2 m, and formulas (1) and (2) give:

S n 2 = ( n 1 S n 1 ¡1 )/s n 1 = (¡2 t 2 ) m¡1 S n 1 ¡1 = (¡2 t 2 ) m x 2m ¡ (¡2) m¡1 t 2m¡1 ; S n 2 ¡1 = prem(S n 1 ; ¡S n 1 ¡1 )/s n 1 m+1 = (¡2) m t 2m+1 x m ¡ t (¡2 t 2 ) m :
Then we get n 3 = m, and

S n 3 = (s n 2 ¡1;n 3 m¡1 S n 2 ¡1 ) / s n 2 m¡1 = ((¡2 t 2 ) m t / (¡2 t 2 ) m ) m¡1 S n 2 = (¡2) m t 3m (x m ¡ 1); S n 3 ¡1 = prem(S n 2 ; ¡S n 2 ¡1 )/(¡2 t 2 ) m(m+1) = (¡2) m t 3m+1 ¡ (¡2) m¡1 t 3m :
Finally n 4 = 0, and

S n 4 = ¡(1 + 2 t) m¡1 (2 t + 1) t 3m :
This example shows that the intermediate sizes grow linearly with the gap sizes, when computing subresultants by means of formulas ( 1) and [START_REF] Abdeljaoued | Bezout matrices, subresultant polynomials and parameters[END_REF]. Of course, for normal sequences of subresultants, these formulas simplify to S n i+1 = prem(S n i ; S n i ¡1 )/s n i 2 = rem(s n i¡1 2 S n i ; S n i ¡1 ), and intermediate coecient sizes just increase by a factor of three.

In the rest of this section, we show that the combination of Lazard's lemma and of the LickteigRoy division allows Algorithm 4 to preserve a coecient growth bounded by a constant factor only.

Computation of the n i

Here we analyze the computation of S n i+1 from S n i ¡1 when n i+1 < n i ¡

Polynomial divisions

We need now to revisit the costs of the division algorithms of Section 3. Proof. The quotient Q may be obtained with the claimed complexity thanks to the latter proposition by using Algorithm 2. Then the remainder is obtained as

A ¡ Q B with O(M(d n)) additional operations in B. Lemma 22. Assume A = B[t]. Let A; B 2 A[x] be of respective degrees n > m > 0. Then the degree in t of pquo(A; B) is 6deg t A + (n ¡ m) deg t B.
Proof. We follow the naive pseudo-division algorithm of A by B:

1. Set R = A and Q = 0; 2. For l from n down to m replace Q by b m Q + r l x l¡m and R by b m R ¡ r l x l¡m B.
At the end, Q and R are respectively the pseudo-quotient and pseudo-remainder. At each step in the loop the degrees in t of R and Q increase by at most deg t B.

Main complexity bounds

In order to apply Proposition 21, we need to bound the degrees in t of all the quotients Q n i arising in the subresultant algorithm. Some care is necessary since A is not assumed to be integral. Our sole assumption is that f n 0 , g n

O(M((d 0 + d 1 n 0 ) n 0 ) log n 0 + n 0 M(d 0 + d 1 n 0 ) log(d 0 + d 1 n 0 )) + O(M(D n 1 ) log n 1 + n 1 M(D) log D log n 1 ):
operations in B. In addition, the underlined expressions may be discarded when B is a eld.

Proof. We rst analyse the cost C(l) of Algorithm 3 in terms of operations in B:

C(l) 6 C(h) + C(l ¡ (n i ¡ n k+1 )) + c (M(D l) + l M(D) log D +M(D (n k ¡ n k+1 )) log(n k ¡ n k+1 ) + (n k ¡ n k+1 ) M(D) log D)): Since h and l ¡ (n i ¡ n k+1 ) are <l /2, it is classical to deduce C(l) = O(M(D l) log l + l M(D) log D log l):
Then we analyze the cost of Algorithm 4 in terms of operations in B:

steps 1, 2 take O(d 0 n 0 + d 1 n 1 ), step 3 takes O(M((d 0 + d 1 n 0 ) n 0 ) log n 0 + n 0 M(d 0 + d 1 n 0 ) log(d 0 + d 1 n 0 ))
O(M((d 0 + d 1 n 0 ) n 0 ) log n 0 + n 0 M(d 0 + d 1 n 0 ) log(d 0 + d 1 n 0 )) + O(M(D n 1 ) log n 1 + n 1 M(D) log D log n 1 ): operations in B. Corollary 26. Assume A = B[t]. Let F ; G 2 A[x] be of respective degrees n 0 > n 1 , let d 0 = deg t F, d 1 = deg t G, and D = d 1 n 0 + d 0 n 1 . If f n 0 , g n 1 ,
O(M(D n 0 ) + M((d 0 + d 1 n 0 ) n 0 ) log n 0 + n 0 M(D) log(D)) + O(M(D n 1 ) log n 1 + n 1 M(D) log D log n 1 )
operations in B. The sole computation of one subresultant without its associated cofactors takes

O(M((d 0 + d 1 n 0 ) n 0 ) log n 0 + n 0 M(d 0 + d 1 n 0 ) log(d 0 + d 1 n 0 )) + O(M(D n 1 ) log n 1 + n 1 M(D) log D log n 1 ):
In addition, the underlined expressions may be discarded when B is a eld. Let us mention that complexities in terms of convex hulls of the supports of A and B may further be derived thanks to the algorithm in [START_REF] Berthomieu | Reduction of bivariate polynomials from convex-dense to dense, with application to factorizations[END_REF]. Previously such deterministic costs for bivariate gcds were involving hypotheses on the cardinality of K in order to use fast evaluation/interpolation strategies, as in [START_REF] Zur Gathen | Modern computer algebra[END_REF]Chapter 10].

The next corollary concerns the cost of the bivariate multi-gcd problem, which is for instance useful for computing separable decompositions [START_REF] Lecerf | Fast separable factorization and applications[END_REF]. If A is a eld, written K, which contains suciently many elements, then may be interpolated from n + 1 values of t. For each value a of t the specialized resultant Res x (A(x); a ¡ B(x)) takes O(M(n) log n) operations in K. This leads to a total cost O(n M(n) log n), which is smaller than the general bound of the latter corollary. However, if K = F 2 , we need to perform this evaluation/interpolation procedure over L = F 2 with being the rst integer such that 2 > n + 1: the resultant may be thus obtained with O(n M(n) M(log n) log n) operations in F 2 . With M(n) = n log n log log n, this cost rewrites O(n 2 log 3 n (log log n) 2 log log log n); which is higher than the one of the latter corollary. Notice that there also exist algorithms with subquadratic costs based on power projections and the NewtonGirard formula [26, Section 2], but they also require hypotheses on K (such as characteristic zero, or to be a nite eld, etc).

Implementation.

The algorithms presented in this article have been implemented in the C++ library called Algebramix of Mathemagix [START_REF] Van Der Hoeven | [END_REF] from 2009. The source code is available from our SVN server http://gforge.inria.fr/projects/mmx/. For the design of the C++ libraries see [START_REF] Van Der Hoeven | Modular SIMD arithmetic in Mathemagix[END_REF]Section 4].

When A is a eld, the simplest algorithm with quadratic cost is adapted from the naive version of the Euclidean algorithm: it is available from polynomial_naive.hpp. The halfgcd is implemented in polynomial_dicho.hpp.

Over rings, the le polynomial_ring_naive.hpp contains the following implementations of subresultants: over any commutative ring we appeal to Berkowitz' algorithm [START_REF] Berkowitz | On computing the determinant in small parallel time using a small number of processors[END_REF] to compute the dening determinants without division; when a partial division routine is available in the ground ring, we use Ducos' algorithm [START_REF] Ducos | Optimizations of the subresultant algorithm[END_REF].

Algorithm 4 can be found in polynomial_ring_dicho.hpp. The latter le includes the fast polynomial division routines. These C++ implementations are rather intricate because of several optimizations for various contexts, taking into account for instance which subresultants or cofactors are actually requested. For the sake of convenience, in revision 10523, we added a simple implementation of Algorithm 4 in the Mathemagix language [START_REF] Van Der Hoeven | Mathemagix User Guide[END_REF]: see gregorix/mmx/subresultant_lickteig_roy.mmx, and gregorix/mmx/ subresultant_test.mmx for the tests.

Concerning eciencies, unfortunately our implementations of Algorithm 4 did not reveal to be competitive to Ducos' algorithm even in large sizes for bivariate polynomials. In the special case of Corollary 29, the evaluation/interpolation strategy is in general observed to be faster than the direct use of Algorithm 4.

  x i and B = P i=0 m b i x i in A[x] of respective degrees n and m > 0. Output. quo(A; B). Assumption: the division of A by B is well dened. 1. Let l = n ¡ m, A ~= rev(A; n) and B ~= rev(B; m).

2 .

 2 Compute C ~= A ~/B ~+ O(x l+1 ) with Algorithm 1. 3. Return rev(C ~; l), where C ~is seen as a polynomial of degree l. Proposition 12. Algorithm 2 is correct and takes O(M(n ¡ m) log(n ¡ m)) operations in A.

Lemma 13 .

 13 Let A = / 0; B 2 A[x], and let l 2 Z. If l < 0 then we have bA B c l+deg A = bA bB c l c l+deg A , otherwise we have bA B c l+deg A = bA bB c l c deg A .

and S n 2 =

 2 S n 1 ¡1 . Then we have S n 2 ¡1 = ¡415 x 2 ¡482 x+890, and we enter Algorithm 3 with S n 2 , S n 2 ¡1 and l = 3. In a recursive call with bS n 2 c 1 , bS n 2 ¡1 c 1 and l = 1 we obtain Q n 3 , which has degree 1. Then we deduce S n 3 ¡1 = ¡11348 x + 13885, and then Q n 3 = 4709420 x+11232011. Finally we have n 4 = 1, Q n 4 = 2724234924 x ¡3333274755, and S n 4 ¡1 = ¡240063.

Lemma 17 .

 17 Assume A = B[t]. If deg t F 6 d 0 and deg t G 6 d 1 then we have deg t s n 1 6 (n 0 ¡ n 1 ) d 1 , and, for 0 6 k < n 1 , the partial degrees of S k , U k and V k are 6(n 1 ¡ k) d 0 + (n 0 ¡ k) d 1 .

Proof.

  The correctness is a consequence of part 3 of Theorem 9. When A = B[t], we use the degree bound provided by Lemma 17, and part 3 of Theorem 9 also gives that deg t a = O(D) holds during the execution. Then the cost analysis is rather standard: each product takes O(M(D)), and each division costs O(M(D) log D) by Proposition 12 (the factor log D may be discarded when B is a eld).

steps 1 and 4

 4 involve O(D l) operations in B, step 2 takes C(h), step 3 takes O(M(D l) + l M(D) log D), step 5 takes O(M(D) log D log l + l M(D) log D) by using Lemma 19, step 6 takes O(M(D (n k ¡ n k+1 )) log(n k ¡ n k+1 ) + (n k ¡ n k+1 ) M(D) log D) by combining Lemma 23 to Proposition 21, step 7, 8 and 10 take O(M(D l) + l M(D) log D) by using Lemma 24, step 9 takes C(l ¡ (n i ¡ n k+1 )). Consequently, there exists a constant c > 0 such that

  by combining Lemma 22 to Proposition 21, step 4 just performs O((d 0 + d 1 n 0 ) n 0 ) operations in B, step 5 takes O(M(D) log D log n 1 + n 1 M(D) log D by using Lemma 19, step 6 takes O(M(D n 1 ) log n 1 + n 1 M(D) log D) by combining Lemma 23 to Proposition 21, steps 7 takes O(D n 1 ), step 8 takes C(n 2 ). The total cost of Algorithm 4 is therefore

Corollary 28 .

 28 Let K be a eld, and let A and B 1 ; :::; B s be polynomials inK[t; x]. Let n = deg x A, d = deg t A, n i = deg x B i , d i = deg t B i , and assume n 1 + + n s = O(n), s = O(n), and d 1 + + d s = O(d).Then, the gcds of A and B 1 , A and B 2 ,:::, A and B s may be computed withO(M(d n 2 ) log n + n M(d n) log d) operations in K.Proof. First we compute the content a and the primitive partA ^of A seen in K[t][x] with O(n M(d) log d) operations in K. Now let 1 6 i 6 s. We compute the content b i and the primitive part B i of B i with O(n i M(d i ) log d i ) operations in K.Then we deduce the gcd of a and b i with O(M(d) log d) operations in K. By Corollary 26, we may compute the last nonzero subresultant of A ^and B i with O(M((d+ d i n) n) log n + M(D i n i ) log n i ) operations in K, where D i = d i n + d n i = O(d n).It has degree 6n i in x and O(D i ) in t, hence its primitive part may be obtained withO(n i M(D i ) log(D i )) = O(n i M(d n) log(d n)).The total cost for computing all the gcds is O(nM(d) log d + P i=1 s [n i M(d i ) log d i + M((d + d i n) n) log n + M(d n n i ) log n i + n i M(d n) log(d n)]) = O(M(d n 2 ) log n + n M(d n) log d).An other application of Corollary 26 is the computation of characteristic polynomials in a A-algebra of the form A[x]/(A(x)).Corollary 29. Let A be an integral domain endowed with its partially dened division, and let A and B be two polynomials in A[x] of degrees 6n. Then, the resultant (t) = Res x (A(x); t ¡ B(x)) may be computed with O(M(n 2 ) log n + n M(n) log 2 n) operations in A. If A is a eld then the cost simplies to O(M(n 2 ) log n).

Proof.

  It suces to use Corollary 26 with d 0 = 0 and d 1 = 1.

  and S n 2 ¡1 by means of the above division algorithm. Let us assume these polynomials already computed, and let us examine how to obtain the other quotients Q n i+1 for i > 2. The key observation is that Q n 3 only depends on the n 2 ¡ n 3 + 1 highest coecients of S n 2 and S n 2 ¡1 , so it can be obtained from bS n 2 c n 3 ¡deg Q n 3 and bS n 2 ¡1 c n 3 ¡deg Q n 3 . Then Q n 4 only depends on the n 3 ¡ n 4 + 1 highest coecients of S n 3 and S n 3 ¡1 . But the coecient of degree j in S n 3 ¡1 only depends on the coecients of S n 2 and S n 2 ¡1 of degrees >j ¡ deg Q n 3 . Consequently Q n 3 and Q n 4 may be obtained from bS n 2 c n 4 ¡(deg

  s n 1 ; :::; s n w are nonzero divisors in A.1. If bS ni ¡1 c = 0 or l < n i ¡ n i+1 thenreturn nothing. Notice that n i+1 is determined from bS n i ¡1 c and when bS n i ¡1 c = / 0. 2. Let h = bl /2c and call the algorithm recursively with n 0 ; :::; n i , h, and bS n i c n i ¡2h , bS n i ¡1 c n i ¡2h . Let n i+1 ; :::; n k , s n i+1 ; :::; s n k , N n i+1 ; :::; N n k , and N n i+1 ;n k be the data obtained in return. 3. Let = + n i ¡ n k , and compute If bS n k ¡1 c = 0 or l < n i ¡ n k+1 then return n i+1 ; :::; n k , s n i+1 ; :::; s n k , N n i+1 ; :::; N n k , and N n i+1 ;n k . Notice that n k+1 is determined from bS n k ¡1 c and when bS n

	!			!
	bS n k c bS n k ¡1 c	from N n i+1 ;n k	bS n i c bS n i ¡1 c	/s n i 2 .
	4.

  Then we detail the computations of the n i and analyze the cost of the divisions in A[x]. We are interested in deterministic algorithms that do not rely on fast evaluation/interpolation schemes, which would require specic assumptions on B.We recall that the product of two polynomials in B[t][x] of partial degrees 6d in t and 6n in x may be achieved with O(M(d n)) operations in B by means of the classical Kronecker substitution[START_REF] Zur Gathen | Modern computer algebra[END_REF] Section 8.4].

  1. Example 18 shows that n i should not be computed as the division of s n i ¡1;n i+1 n i ¡n i+1 ¡1 by s n i n i ¡n i+1 ¡2 . Instead, we appeal to the divide and conquer algorithm based on Lazard's lemma (part 3 of Theorem 9), as described by Ducos in [21, p. 338], and which is reminiscent of the classical binary powering algorithm. Input. s n i and s n i ¡1;n i+1 in A, for 1 6 i 6 w, and an integer 1 6 l 6 n i ¡ n i+1 ¡ 1. If l is even then return b else return (s n i ¡1;n i+1 b)/s n i . Lemma 19. Algorithm 5 is correct and takes O(log (n i ¡ n i+1 )) operations in A. If A = B[t], deg t F 6 d 0 , and deg t G 6 d 1 , then the algorithm costs O(M(D) log(D) log(n i ¡ n i+1 )) operations in B, where D = n 1 d 0 + n 0 d 1 . If B is a eld then the latter cost simplies to O(M(D) log(n i ¡ n i+1 )).

	Algorithm 5
	Output. s n i ¡1;n i+1 l Assumption: n i+1 < n i ¡ 1. /s n i l¡1 2 A. 1. If l = 1 then return s n i ¡1;n i+1 . 2. Let h = bl /2c. 3. Recursively compute a = s n i ¡1;n i+1 h	/s n i h¡1 , and then b = a 2 /s n i .
	4.

  1 when A = B[t]. We use the notation of Algorithm 1, and write C(n) for its cost. A straightforward induction shows that the degrees in t of the values of ¡ and of the output are always at most d throughout the recursive calls. This implies that the degrees in t of the values of always remain bounded by 2 d throughout the recursive calls. Consequently, there exists a constant c > 0 such that Proposition 21. Assume A = B[t]. Let A; B 2 A[x] be of respective degrees n > m > 0, and such that the division of A by B is well dened, and let d bound the degrees in t of A, B, and of their quotient Q. Then, Q may be obtained with O(M(d (n ¡ m)) log(n ¡ m) + (n ¡ m) M(d) log d) operations in B. If B is a eld then the latter cost simplies to O(M(d (n ¡ m)) log(n ¡ m)). Then the remainder may be deduced from Q with O(M(d n)) additional operations in B.

	C(n) 6 C(h) + C(n ¡ h) + c M(d n):
	In general we have C(1) = O(M(d) log d) by Proposition 12, which simplies to C(1) =
	O(M(d)) when B is a field. The conclusion classically follows by unrolling the latter
	inequality.

Proposition 20. Assume A = B[t]. Let ; ¡ 2 AJxK/(x n ) be such that the division of by ¡ is well dened, and let d bound the degrees in t of , ¡, and . Then, Algorithm 1 takes O(M(d n) log n + n M(d) log d) operations in B. If B is a eld then the latter cost simplies to O(M(d n) log n).

Proof.

  The bound for Q n 1 follows from Lemma 22. By equation (14) the degree bound for Q n 2 follows from the one for U n 2 ¡1 given in Lemma 17. For i > 3, we combine equation (13), Proposition 8, and Lemma 17. Lemma 24. Assume A = B[t], deg t F 6 d 0 , deg t G 6 d 1 , and let D = n 1 d 0 + n 0 d 1 . Then, for all 1 6 i < j 6 w the degrees in t of the entries of s n i 2 M n j M n i+1 are O(D). Theorem 25. Assume A = B[t], deg t F 6 d 0 , deg t G 6 d 1 , and let D = n 1 d 0 + n 0 d 1 . Then, Algorithm 4 takes

	Proof. This is a consequence of Lemmas 10 and 17.

1 , s n 1 ; :::; s n w are nonzero divisors in A.

Lemma 23. Assume A = B[t], deg t F 6 d 0 , deg t G 6 d 1 , and let D = n 1 d 0 + n 0 d 1 . Then, deg Q n 1 6 d 0 + n 0 d 1 ,

and for all 2 6 i 6 w the degree in t of Q n i is O(D).

Proof.

  Corollary 27. Let K be a eld, and let A and B be two polynomials in K[t; x] of degrees 6d in t and 6n in x. Then, the gcd of A and B may be computed withO(M(d n 2 ) log n + n M(d n) log d) operations in K, which simplies to O(M((d n) 3/2 ) log d) whenever n 6 d.Proof. First we compute the contents and primitive parts ofA and B seen in K[t][x] with O(n M(d) log d) operations in K. Then we deduce the gcd of the contents with O(M(d) log d) operations in K. By Corollary 26, we may compute the last nonzero subresultant of the primitive parts with O(M(d n 2 ) log n) operations in K. It has degree 6n in x and O(d n) in t, hence its primitive part may be obtained with O(n M(d n) log(d n)) more operations.

On the complexity of the LickteigRoy subresultant algorithm

On the complexity of the LickteigRoy subresultant algorithm

Proof. After calling Algorithm 4, we have the atomic transition matrices at our disposal, and we may easily build T n 2 via formula [START_REF] Brown | On Euclid's algorithm and the computation of polynomial greatest common divisors[END_REF]. Now assume we want to obtain S n i , S n i ¡1 and T n i for 2 6 i 6 w. For this purpose we compute T n i rst, and then use

M n i M n 3 by the following divide and conquer algorithm that computes s n i 2 M n j M n i+1 for i < j as follows:

1. If j = i + 1 then return s n i 2 M n i+1 (available as N n i+1 from the output of Algorithm 4);

2. Let h = b(i + j) / 2c, and recursively compute

The correctness is ensured by Lemma 10, which also provides us with the degree bound

In addition the degree in t of the latter matrix is O(D) by Lemma 24. Therefore the computation of

If the cofactors are not needed then we may truncate

If B is a eld that contains suciently many elements (for instance when the characteristic is zero), then we may use fast evaluation and interpolation algorithms [START_REF] Zur Gathen | Modern computer algebra[END_REF]Chapter 10].

To simplify the situation we reformulate the costs in terms of d = max (d 0 ; d 1 ). In this way, a subresultant of degree k may be interpolated from D + 1 values of t with O(k M(d n 0 ) log(d n 0 )). All the specializations of F and G may be obtained with

and each specialized subresultant requires O(M(n 0 ) log n 0 ). Thus the total cost is

If k ' n 0 , then the dominant term is k M(d n 0 ) log(d n 0 ), which is asymptotically higher than the cost O(M(d n 0

2

For any k, the ratio of the cost of Corollary 26 over the one of (15) when M(n) = n log n log log n becomes bounded by

Consequently, in all case the cost of Corollary 26 is never asymptotically higher (up to a constant factor) than with the evaluation/interpolation strategy.

Applications

Gcd is a very classical application of subresultants. Our Corollary 26 leads to the following deterministic result.


A Plantec
email: plantec@univ-brest.fr

V Ribaud
email: ribaud@univ-brest.fr

Syseca & Libr

Using and Re-using Application Generators

Application generators translate specifications into products (programs, documentation). They parse specification statements into data structures (called dictionaries or meta-models), from which desired products can be derived. An application generator builder offers a way to define specification languages and associated parsers, to describe and traverse the meta-models structure and to specify the derivation on this structure.

STEP is an ISO 10303 standard developed to facilitate product information sharing by specifying sufficient semantic content for data and their usage. STEP offers very useful technology, such as an object oriented modelling language EXPRESS and a data access interface, and can be applied for the design and the implementation of application generators.

EUGENE is a STEP-based framework, that is intended for the building of application generators. After a concise presentation of EUGENE, this paper presents different application generators built with this framework, mainly data management tools, a browsing tool and conversion tools. Then, we discuss the main benefits induced from our approach.

Introduction

Application generators translate source specifications into target products as programs or documentation. They use data structures called dictionaries or meta-models to store the source specifications and intermediate representations. As in compilers, there are two main functions in a generator: parsing of the source specifications and the target code generation.

Typically, an application generator builder offers a means to define specification languages and associated parsers, to describe and traverse the meta-model structure and to specify the derivation on this structure. In most builders, derivation is specified with templates (or skeletons) of code. Templates contain a mixture of commands operating on the meta-models and "real" code directly inserted into the products. With classical tools and methods, specifying and implementing a generator is very much like specifying and implementing a language.

In this paper, we describe a way of building application generators directly from the object-oriented specification of meta-models. We used this possibility in a tool, called EUGENE [START_REF] Plantec | [END_REF]. Meta-models and code templates are programmed with the object-oriented data modelling langage, EXPRESS.

We argue that building generators from the specification of meta-models can simplify their design and implementation and makes their use and re-use in industrial projects easier. As a consequence, generators can be implemented for specific problems and generated realizations can be well integrated.

EUGENE is briefly presented in chapter 2. In chapter 3, we show different generators built with EUGENE. We discuss the main benefits of EUGENE in chapter 4.

Building application generators with EUGENE

An application generator built with EUGENE uses meta-data in order to generate a target textual representation. Meta-data come from a source specification analysis and are stored in different kind of meta-models.

Derivations are programmed with imperative functions, called translation functions. They are made of fixed parts that are mainly either meta-model traversal routines or string constants directly put into the target and made of variable parts that are values fetched from the meta-models.

Figure 1 shows that an application generator built with EUGENE is only a process that consumes meta-data and produces some realization. The meta-data are themselves produced by other processes or tools. Metadata consist of a representation of the source specification called intermediate representation. With EUGENE, any meta-data producer can be used. Most application generators builders produce the generators while exploiting the definition of the syntax and the semantic of the source specification language [2,3,4]. EUGENE differs from this way, because the generator is not produced from a grammar-like specification, but from the structure and the organization of the meta-data (globally referred to as the meta-models).

STEP description and implementation methods

The EXPRESS language [1] is an object-oriented modelling language. The application data are described in schemata. A schema has the type definitions and the object descriptions of the application called Entities. An entity is made up of attributes and constraint descriptions.

The constraints expressed in an entity definition can be of four kinds: (1) the unique constraint allows entity attributes to be constrained to be unique either solely or jointly, (2) the derive clause is used to represent computed attributes, (3) the where clause of an entity constraints each instance of an entity individually and (4) the inverse clause is used to specify the inverse cardinality constraints. Entities may inherit attributes and constraints from their supertypes.

The STEP physical file format defines an exchange structure using a clear text encoding of product data for which a conceptual model is specified in the EXPRESS language. The mapping from the EXPRESS language to the syntax of the exchange structure is specified in [2].

The Standard Data Access Interface (SDAI) [3] defines an access protocol for EXPRESS-modelled databases and is defined independently from any particular system and language. The representation of this functional interface in a particular programming language is referred to as a language binding in the standard. As an example, ISO 10303-23 is the STEP part describing the C++ SDAI binding [4].

The five main goals of the SDAI are: (1) to access and manipulate data which are described using the EXPRESS language, (2) to allow access to multiple data repositories by a single application at the same time, (3) to allow commit and rollback on a set of SDAI operations, (4) to allow access to the EXPRESS definition of all data elements that can be manipulated by an application process, and (5) to allow the validation of the constraints defined in EXPRESS.

An SDAI can be implemented as an interpretor of EXPRESS schemata or as a specialized data interface. The interpretor implementation is referred to in the standard [3] as the SDAI late binding. An SDAI late binding is generic in nature. The specialized implementation is referred to in the standard as the SDAI early binding.

Figure 1 The building and the using of an application generator

Thus the tool which feeds the meta-models with meta-data can be a classical parser (generally automatically

META-LEVEL

Meta-data (intermediate representation)

Application generator

Eugene

CASE tool Parser

Meta-models & Translation functions

Source specification

Realization

STEP toolkit

STEP file

built in other environments than EUGENE) but another tool not relying on a precise and formal concrete syntax. Let us take the example of the classical work done by a screen generator of a Relational Data Base Management System (RDBMS). The generator uses the SQL table definitions (usually stored in a dictionary) and produces a graphical form (a window with fields, buttons,etc.) for a given user interface toolkit (Windows, X, etc.).

The meta-models

The specification of the generator being built consists in several meta-models. A precise explanation of the nature of these meta-models can be found in [START_REF] Plantec | [END_REF]. To simplify the example, we will use one source metamodel.

The source language meta-model consists in a set of EXPRESS schemata that describes the source language data constructs. The main components of a source language meta-model are types and entities, describing concepts that can be used with the source language. Entities provide buckets to store meta-data while global and local EXPRESS constraints are used to ensure metadata soundness.

Considering the classical example of building form screens from SQL table definitions, the source language is SQL and figure 2

The translation functions

The translation functions are also written in EXPRESS and are specified in the translation schema. The specification of translation functions is a programming activity in which EXPRESS is used as an imperative language. A typical translation function returns a string and is parameterized with types that are entities defined in the source language meta-model. The resulting string represents part of the target textual representation. Because of the nature of parameter types, this activity is often called meta-programming [1,[START_REF] Ladd | A*: A Language for Implementing Language Processors[END_REF].

Figure 3 shows two skeletons of translation functions related to the paper example. column_to_field and table_to_form can be used in order to produce the definitions of widgets for, respectively, a column and a

Figure 3 An example of translation functions

Applications

As most application generator builders, EUGENE can be used to produce software components from source specifications of a high abstraction level.

Two conditions are a prerequisite, namely :  the intermediate representation of the source specification can be described with an EXPRESS meta-schema;  the realization can be produced with imperative processing applied to meta-data.

Data management tools

EUGENE was primarily intended to build data management components generators. The semantic of the domain is well-known (roughly query, insert, update and delete data). Data types are the variable part, they parametrize the realization. The way the implementation is carried out depends on data types and on the storage system (file, database, …) called a repository. Hence, generators for the data management domain are widely used.

Specifications for such generators are essentially data structures as records, classes, schemas. Realizations, built from these specifications, manage repository in a specific language such as C, SQL, C++. The SDAI is an abstraction of the data management domain. From the conceptual description (in EXPRESS) of the domain entities, the SDAI implementation is produced (the functionalities of the SDAI are shown in the box STEP given above). Another kind of generator, called IO, provides the basic access interface to SQL tables through a host language; embedded SQL for query and update is generated from a data structure description.

As depicted in figure 1, a complete generator is made up of a parser and an application generator. The table below presents the realization time (in man-month) of several parsers (we re-used empty yacc-parsers):

Browsing tool

EUGENE was used to build an instance browser generator (see figure 4). Given a data schema as specification source, a data management application in Java is generated. The generated application enables the user to browse the instances of a repository, to read and write the values associated with the instances and to navigate along their relationships. This instance browser generator produces two components : a graphical user interface and a SDAI to manage the instances. The second component is produced by the Java SDAI generator described in a previous paragraph. A threelayered generator produces the graphical user interface, each layer consuming meta-data produced by the former layer.

Figure 4 Java data browser generation

Conversion tools

A conversion tool transforms a source system Meta-data representation into a target system representation. Unlike former tools, a conversion tool has no leverage effect, i.e. the complexity of generated products is similar to the source representation. Documentation tools are kinds of conversion tools. These tools use a meta-model which describes information extracted from source documents. This metamodel is used to reproduce the document according to the target representation.

STEP file

This architecture was used in two tools: ADOC, a program documentation generator and an indexation and navigation tool for blind people.

ADOC produces Latex program documentation from the analysis of structured comments contained in programs sources. The comments should be structured according to a BNF grammar. ADOC uses a Lex/Yacc parser and a generator built with EUGENE. The metaschema, derived from the BNF grammar, is populated by the parser. It is then used by the generator to produce Latex documents.

The second tool aims to make it easier for blind people to navigate within Web documents. Such documents are originally written using a word processor (e.g. Latex) and then generated in HTML using a converter (e.g. Latex2HTML). The tool consists in a HTML scanner and a HTML generator built with EUGENE. The tree of Web pages is scanned to construct a word index and the hierarchical organization in chapter, section, and so on. The scan result, in the form of a STEP file exchange, is consumed by the generator to produce a table of contents, a plan and a word index, all hyperlinked to the original tree. Rather than exploring in depth the Web tree, blind people consult these meta-structures prior to navigating in the Web pages.

The joint use of both tools provides meta-structures to navigate in program documentation. As depicted in figure 5, inter-operability between these two tools can be achieved in two ways : either by the Latex2HTML converter, or by a gateway built with EUGENE which populates the meta-schema used by the HTML generator from the meta-data resulting from the comments analysis.

Figure 5 Conversion tools

Discussion

Usability

EUGENE uses an unique language, EXPRESS, for the meta-model specification and translation functions. Using EUGENE requires essentially data design and imperative programming. These activities are familiar to software engineers.

Parsers are not built with EUGENE, rather with wellknown compiler-compiler such as Yacc. Meta-model implementation is automatically derived from specification. Parsers that have to be developed can use this implementation in order to create an intermediate representation.

Re-use

As mentioned in §3.1, EUGENE was used to build three SDAI-subset generators, for C++, Java and Smalltalk-80. While the specification source language is the same, the three SDAI-subset generators re-use the same parser and source meta-model. While the target languages (C++, Java, Smalltalk-80) are object-oriented, the translation functions of the C++ generator serve as a good example for the translation functions of the two The HTML scanner of §3.3 was developed in a separate project. It was re-used without changes, adding only a feature intended to flush the scan results in a STEP file, while using the SDAI built for the meta-model of the HTML generator.

A design and validation method of data models is defined in the STEP standard [START_REF] Palmer | Guidelines for the development and approval of STEP application p r o t o c o l s . T e c h n i c a l r e p o r t[END_REF]. EUGENE applies a part of this method in order to develop and re-use the metamodels.

Interoperability

Thanks to the STEP standard, meta-data can be automatically encoded in neutral ASCII files. It enables generators to read meta-data independently from their producers. Then a generator built with EUGENE is easily interoperable.

The instances browser depicted in § 3.2 is made up of two components, sharing the same meta-model (the same specification language), thus making it possible for the two generators to inter-operate. The overall project was implemented by different people over a period of three years, and thanks to the integration method no major problem has arisen.

In § 3.3, figure 5 shows that the hypertext index generator consumes meta-data produced by an HTML parser. In order to build HTML indexes for programs comments, a first solution consists in using ADOC generator and Latex2HTML. In order to re-use the hypertext index generator, a second solution could be to produce meta-data needed by this generator directly from the meta-data produced by the comments parser.

Conclusion

This paper has briefly described the EUGENE environment, a STEP-based application generator builder. A generator is automatically built from the specification of the intermediate representation manipulated by the generator (the meta-schema) and from the specification of translation functions.

Different generators are presented: data management tools, browsing tool, conversion tools.

Using EUGENE offers three main benefits : 1. easy use: half of projects had be done within the context of final-year coursework (postgraduate students) and the other half in a software compagny. All the projects succeeded in producing valid prototypes or operational generators. 2. re-use: several meta-schemata and software components (parsers or converters) were re-used between the different projects and are suitable for future re-use. 3. inter-operability: all the projects demonstrated that software integration through meta-data exchange is a suitable characteristic for both the integration of the generator into its environment and its re-use.

  shows a simplified SQL meta-model. The table definition is related to a list of columns.

	END_ENTITY;
	ENTITY column;
	name : STRING;
	domain : simple_type;
	END_ENTITY; ...
	END_SCHEMA ;
	Figure 2 An example of source meta-model: a simple
	SQL dictionary

SCHEMA sql_dictionary; ENTITY simple_type SUPERTYPE OF (ONEOF(real_type, integer_type, string_type)); END_ENTITY; ... ENTITY table; name : STRING; columns : LIST [1:?] OF column; ...

  table.

	SCHEMA sql_traduction;
	USE FROM sql_dictionary ;
	END_LOCAL ;
	(* build statements for the management of a
	* column in a graphical form*)
	RETURN (result) ;
	END_FUNCTION ;
	FUNCTION table_to_form( t : table ) : STRING ;

FUNCTION column_to_field ( c : column) : STRING ; LOCAL result : STRING := ''; LOCAL result : STRING := ''; END_LOCAL ; (* build statements for the management of a * table in a graphical form *) REPEAT no := LOINDEX(t.columns) TO IINDEX(t.columns) ; …. END_REPEAT ; RETURN (result) ; END_FUNCTION ; END_SCHEMA ;