Alin Bostan
email: alin.bostan@inria.fr

Muhammad Chowdhury

Joris Van Der Hoeven
email: joris.vanderhoeven@math.u-psud.fr

Éric Schost
email: eschost@uwo.ca

Homotopy techniques for multiplication modulo triangular sets

Keywords: Triangular sets, multiplication, complexity

We study the cost of multiplication modulo triangular families of polynomials. Following previous work by Li, Moreno Maza and Schost, we propose an algorithm that relies on homotopy and fast evaluation-interpolation techniques. We obtain a quasi-linear time complexity for substantial families of examples, for which no such result was known before. Applications are given to notably addition of algebraic numbers in small characteristic.

Introduction

Triangular families of polynomials are a versatile data structure, well adapted to encode geometric problems with some form of symmetry [3,[START_REF] Kogan | Computation of canonical forms for ternary cubics[END_REF][START_REF] Foursov | On computer-assisted classification of coupled integrable equations[END_REF][START_REF] Gaudry | Construction of secure random curves of genus 2 over prime fields[END_REF]. However, in spite of this, many complexity questions are still not answered in a satisfying manner.

A high-level question is to provide sharp estimates on the cost of solving polynomial systems by means of triangular representations. This problem has a geometric nature; it itself relies on several difficult lower-level questions, such as the cost of basic operations with triangular sets. In this paper, we address one such question: the arithmetic cost of multiplication of polynomials modulo a triangular set. This justifiably stands as a central question, since many higher-level routines are built on top of it, such as inversion [START_REF] Langemyr | Algorithms for a multiple algebraic extension[END_REF][START_REF] Van Hoeij | A modular GCD algorithm over number fields presented with multiple extensions[END_REF][START_REF] Dahan | On the complexity of the D5 principle[END_REF][START_REF] Li | Fast arithmetic for triangular sets: from theory to practice[END_REF], lifting techniques [START_REF] Dahan | Lifting techniques for triangular decompositions[END_REF] for modular algorithms, solving systems of equations [START_REF] Li | High-performance symbolic computation in a hybrid compiled-interpreted programming environment[END_REF], etc.

Problem statement, overview of our results. We work in R[X 1 , . . . , X n ], where R is a ring, and we are given a set of relations of the form T T n (X 1 , . . . , X n ) . . .

T 2 (X 1 , X 2 )
T 1 (X 1 ).

The polynomials T form a triangular set: for all i, T i is in R[X 1 , . . . , X i ], is monic in X i and reduced modulo T 1 , . . . , T i-1 , in the sense that deg(T i , X j ) < deg(T j , X j ) for j < i. As an aside, note that the case where T i is not monic but with a leading coefficient invertible modulo T 1 , . . . , T i-1 reduces in principle to the monic case; however, inversion modulo T 1 , . . . , T i-1 remains a difficult question [START_REF] Dahan | On the complexity of the D5 principle[END_REF], of cost higher than that of multiplication.

As input, we consider two polynomials A, B reduced modulo T . The direct approach to multiply them modulo T is to perform a polynomial multiplication, followed by the reduction modulo T , by a generalization of Euclidean division. As far as complexity is concerned, when the number of variables grows, this kind of approach cannot give linear time algorithms. Consider for instance the case where all T i have degree 2 in their main variables X i . Then, A and B both have 2 n monomials, but their product before reduction has 3 n monomials; after reduction, the number of monomials is 2 n again. If we let δ = 2 n be a measure of the input and output size, the cost of such an algorithm is at least 3 n = δ log 2 (3) .

In this paper, we show that a different approach can lead to a quasi-linear time algorithm, in cases where the monomial support of T is sparse, or when the polynomials in T have a low total degree. This will for example be the case for systems of the form

X 2 n -2X n-1 . . . X 2 2 -2X 1 X 2 1 or X 2 n -X n-1 . . . X 2 2 -X 1 X 2 1 , (1) 2
whose applications are described later on. Our result also applies to the following construction: start from F ∈ R[X], say F = X 3 -X 2 + X -3, and define the so-called "Cauchy modules" [START_REF] Rennert | Calcul de résolvantes avec les modules de Cauchy[END_REF] F 1 , F 2 , F 3 , which are used in effective Galois theory [START_REF] Rennert | Calcul de résolvantes avec les modules de Cauchy[END_REF][START_REF] Abdeljaouad | Computation of the decomposition group of a triangular ideal[END_REF][START_REF] Renault | A modular algorithm for computing the splitting field of a polynomial[END_REF]:

F 1 (X 1 ) = F (X 1 ) = X 3 1 -X 2 1 + X 1 -3 F 2 (X 1 , X 2 ) = F1(X1)-F1(X2) X1-X2 = X 2 2 + X 2 X 1 -X 2 + X 2 1 -X 1 + 1 F 3 (X 1 , X 2 , X 3 ) = F2(X1,X2)-F2(X1,X3) X2-X3 = X 3 + X 2 + X 1 -1. (2)
For examples [START_REF] Abdeljaouad | Computation of the decomposition group of a triangular ideal[END_REF] and [START_REF] Aho | Evaluating polynomials at fixed sets of points[END_REF], our algorithms give the following results:

• for T as in [START_REF] Abdeljaouad | Computation of the decomposition group of a triangular ideal[END_REF], multiplication modulo T can be performed in quasi-linear time O˜(δ), where δ = 2 n is the input and output size, and where O˜(δ) stands for O(δ(log(δ)) O (1) ).

• for T as in [START_REF] Aho | Evaluating polynomials at fixed sets of points[END_REF], with n = deg(F ), multiplication modulo T can be performed in quasi-linear time O˜(δ), where δ = n! is the input and output size.

No previous algorithm was known featuring such complexity estimates.

Our approach. To obtain this quasi-linear cost, we have to avoid multiplying A and B as polynomials. Our solution is to use evaluation and interpolation techniques, just as FFT multiplication of univariate polynomials is multiplication modulo X n -1.

Fast evaluation and interpolation may not be possible directly, if T does not have roots in R (as in the previous examples). However, they become possible using deformation techniques: we construct a new triangular set U with all roots in R, and multiply A and B modulo S = ηT+(1-η)U, where η is a new variable. The triangular set S has roots in R [[η]], by Hensel's lemma, so one can use evaluation-interpolation techniques over R [[η]].

This idea was introduced in [START_REF] Li | Fast arithmetic for triangular sets: from theory to practice[END_REF], but was limited to the case where all polynomials in T are univariate: T i was restricted to depend on X i only, so this did not apply to the examples above. Here, we extend this idea to cover such examples; our main technical contribution is a study of precision-related issues involved in the power series computations, and how they relate to the monomial support of T.

Previous work. It is only recently that fast algorithms for triangular representations have been throughly investigated; thus, previous results on efficient multiplication algorithms are scarce. All natural approaches introduce in their cost estimate an overhead of the form k n , for some constant k.

The main challenge (still open) is to get rid of this exponential factor unconditionally: we want algorithms of cost O˜(δ), where δ = d 1 • • • d n is the number of monomials in A, B and their product modulo T . For instance, with T i = X di i , the first complexity result of the form O(δ 1+ε ), for any ε > 0, was in [START_REF] Schost | Multivariate power series multiplication[END_REF].

The previous work [START_REF] Li | Fast arithmetic for triangular sets: from theory to practice[END_REF] gives a general algorithm of cost O˜(4 n δ). That algorithm uses fast Euclidean division; for polynomials of low degree (e.g., for d i = 2), the naive approach for Euclidean division can actually give better results. Previous mentions of such complexity estimates (with a constant higher than 4) are in [START_REF] Langemyr | Algorithms for a multiple algebraic extension[END_REF].

As said above, in [START_REF] Li | Fast arithmetic for triangular sets: from theory to practice[END_REF], one also finds the precursor of the algorithm presented here; the algorithm of [START_REF] Li | Fast arithmetic for triangular sets: from theory to practice[END_REF] applies to families of polynomials having T i ∈ R[X i ], and achieves the cost O(δ 1+ε ) for any ε > 0. In that case, the analysis of the precision in power series computation was immediate. Our main contribution here is to perform this study in the general case, and to show that we can still achieve similar costs for much larger families of examples.

Basic notation. Let d = (d 1 , . . . , d n ) be a vector of positive integers. In what follows, these will represent the main degrees of the polynomials in our triangular sets; without loss of generality, we will thus always suppose d i ≥ 2 for all i.

Recall that R is our base ring and that X 1 , . . . , X n are indeterminates over R. We let M d be the set of monomials

M d = X e1 1 • • • X en n | 0 ≤ e i < d i for all i . We denote by Span(M d ) the free R-submodule generated by M d in R[X 1 , . . . , X n ]: Span(M d ) = A = m∈M d a m m | a m ∈ R for all m ∈ M d .
This is thus the set of polynomials A in R[X 1 , . . . , X n ], such that deg(A i , X i ) < d i holds for all i. Finally, we let δ d be the product

δ d = d 1 • • • d n ; this is the cardinality of M d .
Remark that since all d i are at least 2, we have the bounds

2 n ≤ δ d and i≤n d 1 • • • d i ≤ 2δ d .
The former plainly follows from the inequality 2 ≤ d i ; the latter comes from observing that

d 1 • • • d i 2 n-i ≤ d 1 • • • d n = δ d ; this yields d 1 • • • d i ≤ δ d /2 n-i ,
from which the claim follows by summation.

The multi-degree of a triangular set

T = (T 1 , . . . , T n ) is the n-tuple d = (d 1 , . . . , d n ), with d i = deg(T i , X i ) 1≤i≤n . In this case, R[X 1 , . . . , X n ]/ T is a free R-module isomorphic to Span(M d ). We say that a polynomial A ∈ R[X 1 , . . . , X n ] is reduced with respect to T if deg(A, X i ) < d i holds for all i. For any A ∈ R[X 1 , . . . , X n ], there exists a unique A ∈ R[X 1 , . . . , X n ],
reduced with respect to T, and such that A -A is in the ideal T . We call it the normal form of A and write A = A mod T .

Outlook of the paper. In Section 2, we introduce some basic complexity notation. The next section presents basic evaluation-interpolation algorithms for so-called equiprojectable sets, which are extensions of algorithms known for univariate polynomials. We deduce our multiplication algorithm in Section 4; examples, applications and experimental results are in Sections 5 and 6.

Preliminaries

Big-O notation is delicate to use in our situation, since our estimates may depend on several (possibly an unbounded number of) parameters (typically, the multi-degree of our triangular sets). Hence, whenever we use a big-O inequality such as f ∈ O(g), it is implied that there exists a universal constant λ such that f (v 1 , . . . , v s ) ≤ λg(v 1 , . . . , v s ) holds for all possible values of the arguments. When needed, we use explicit inequalities. Finally, the notation f ∈ O˜(g) means that there exists a constant α such that f ∈ O(g log(g) α ), where the big-O is to be understood as above.

Our complexity estimates count additions, multiplications, and inversions, when they are possible. We denote by M : N → N a function such that over any ring, polynomials of degree less than d can be multiplied in M(d) operations, and which satisfies the superlinearity conditions of [START_REF] Zur Gathen | Modern Computer Algebra[END_REF]Chapter 8]. Using the algorithm of Cantor-Kaltofen [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF], one can take

M(d) ∈ O(d lg(d) lg lg(d)), with lg(d) = log 2 (max(d, 2)).
We next let C 0 : N → N be a function such that d → C 0 (d)/d is non-decreasing, C 0 (d) ≥ d holds for all d and such that we have, over any ring R:

(1) for any x 1 , . . . , x d in R and any polynomial A ∈ R[X] of degree less than d, one can compute all values A(x i ) in C 0 (d) additions and multiplications in R;

(2) for any x 1 , . . . , x d in R, one can compute the coefficients of the polynomial (X -

x 1 ) • • • (X -x d ) in C 0 (d)

additions and multiplications in R;

(3) for any x 1 , . . . , x d in R, with x i -x j a unit for i = j, and any values v 1 , . . . , v d in R, one can compute the unique polynomial A ∈ R[X] of degree less than d such that A(x i ) = v i holds for all i in C 0 (d) operations in R.

By the results of [START_REF] Zur Gathen | Modern Computer Algebra[END_REF]Chapter 10], one can take C 0 (d) ∈ O(M(d) log(d)). We continue with the well-known fact that the function C 0 also enables us to estimate the cost of lifting power series roots of a bivariate polynomial by Newton iteration. In the following lemma, η is a new variable over R.

Lemma 1. For any polynomial T in R[η, X], monic in X and with deg(T, X) = d, if the roots a 1 , . . . , a d of T (0, X) are known and mutually distinct, one can compute the roots of

T (η, X) in R[η]/ η in O(C 0 (d)M( )) operations in R.
Proof. The algorithm consists in lifting all roots of T in parallel using Newton iteration, using fast evaluation to compute the needed values of T and ∂T /∂X; it is given in Figure 1, where we use a subroutine called EvalUnivariate to do the evaluation. To obtain simpler estimates, we let C(d) = ΛC 0 (d), where Λ ≥ 1 is the constant implied in the big-O estimate in the former lemma. Hence, problems (1), ( 2) and (3) above can be dealt with in C(d) operations, and the lifting problem of the previous lemma can be solved in C(d)M( ) operations. Finally, we introduce another short-hand notation: for a multi-degree d = (d 1 , . . . , d n ), we write

a i ← a i -v i /w i mod η 2.4 ← 2 3 return [a i mod X | 1 ≤ i ≤ d]
L(d) = i≤n C(d i ) d i ≤ n C(d) d , (3)

Evaluation and interpolation at equiprojectable sets

In this section, we recall from [3] the definition of equiprojectable sets. We prove that one can perform evaluation and interpolation at, and construct the vanishing ideal of, equiprojectable sets in linear time, up to logarithmic factors. We deduce an algorithm for multiplication modulo the vanishing ideal of such sets with a similar complexity. These results extend those given in [START_REF] Pan | Simple multivariate polynomial multiplication[END_REF][START_REF] Li | Fast arithmetic for triangular sets: from theory to practice[END_REF], which dealt with the case of points on a regular grid. The extension to our more general context is rather straightforward, but to our knowledge, it has not appeared in print before.

In all this section, R is a ring; we study subsets of R n and their successive projections on the subspaces R i , for i ≤ n. For definiteness, we let R 0 be a one-point set. Then, for 1 ≤ j ≤ i ≤ n, we let π i,j be the projection

π i,j : R i → R j (x 1 , . . . , x i ) → (x 1 , . . . , x j );
if j = 0, we adapt this definition by letting π i,0 be the constant map R i → R 0 . Finally, since this is the projection we use most, we simply write π = π n,n-1 for the projection

R n → R n-1 .
If V is a subset of R n , for β in π(V ), we let V β be the fiber V ∩ π -1 (β). Hence, if β has coordinates (β 1 , . . . , β n-1 ), the points in V β have the form (β 1 , . . . , β n-1 , a), for some values a in R. In all that follows, a finite set is by convention non-empty.

Equiprojectable sets

Let V be a finite set in R n . Equiprojectability is a property of V that describes a combinatorial regularity in the successive projections of V . For n = 0, we say that the unique non-empty subset of R 0 is equiprojectable. Then, for n > 0, V ⊂ R n is equiprojectable if the following holds:

• the projection π(V ) is equiprojectable in R n-1 , and

• there exists an integer d n such that for all β in π(V ), the fiber V β has cardinality d n .

The vector d = (d 1 , . . . , d n ) is called the multi-degree of V . Remark that for n = 1, any finite V ⊂ R is equiprojectable. One easily sees that if V is equiprojectable, its cardinality equals

δ d = d 1 • • • d n ; more generally, π n,i (V ) ⊂ R i is equiprojectable of cardinality d 1 • • • d i .
When R is a perfect field, it is proved in [3] that equiprojectable sets are exactly the zero-sets of triangular sets that generate radical ideals in R[X 1 , . . . , X n ]; we will discuss this in more detail in Subsection 3.4.

We give first a slightly more precise notation for the fibers V β : if V is equiprojectable, then for all β = (β 1 , . . . , β n-1 ) ∈ π(V ), there exist exactly d n pairwise distinct values v β = [a β,1 , . . . , a β,dn ] in R such that

V β = [(β 1 , . . . , β n-1 , a β,i ) | a β,i ∈ v β ]
and thus

V = [(β 1 , . . . , β n-1 , a β,i ) | β = (β 1 , . . . , β n-1 ) ∈ π(V ), a β,i ∈ v β , 1 ≤ i ≤ d n ].
For instance, n-dimensional grids are special cases of equiprojectable sets, where v β is independent of β. Remark also that for some special choices of v β , improvements in the algorithms below are possible (e.g., if v β is a translation of a set of points in geometric progression, using the algorithm of [START_REF] Aho | Evaluating polynomials at fixed sets of points[END_REF]).

Evaluation

Let V be an equiprojectable set of multi-degree d = (d 1 , . . . , d n ), and let M d , δ d be as in Section 1. We denote by Eval V the evaluation map

Eval V : Span(M d ) → R δ d F → [F (α) | α ∈ V ].
We let C Eval be a function such that for any V equiprojectable of multi-degree d, the map Eval V can be evaluated in C Eval (d) operations. In one variable, with n = 1 and d = (d 1 ), C Eval (d) simply describes the cost of evaluating a polynomial of degree less than d 1 at d 1 points of R, so we can take C Eval (d) = C(d 1 ). More generally, we have the following quasi-linear time estimate.

Proposition 1. One can take C Eval (d) ≤ δ d L(d).
Proof. We will use a straightforward recursion over the variables X n , . . . , X 1 . Let W = π(V ), let e = (d 1 , . . . , d n-1 ) be the multi-degree of W and let A(X 1 , . . . , X n ) ∈ Span(M d ) be the polynomial to evaluate. We write

A = i<dn A i (X 1 , . . . , X n-1 )X i n ,
with A i in Span(M e ), and, for β in R n-1 , we define

A β = i<dn A i (β)X i n ∈ R[X n ]. Hence, for β = (β 1 , . . . , β n-1 ) in R n-1 and x in R, A(β 1 , . . . , β n-1 , x) = A β (x).
As a consequence, to evaluate A at V , we start by evaluating all A i at all points β ∈ W . This gives all polynomials A β , which we evaluate at the fibers v β . The algorithm is given in Figure 2. From this, we deduce that we can take C Eval satisfying the recurrence This implies

C Eval (d 1 , . . . , d n ) ≤ C Eval (d 1 , . . . , d n-1 ) d n + d 1 • • • d n-1 C(d n ). Eval(A, V ) 1 if n = 0 return [A] 2 W ← π(V ) 3 for i = 0, . . . , d n -1 do 3.1 A i ← coeff(A, X n , i) 3.2 val[i] ← Eval(A i , W ) / * val[i] has the form [A i (β) | β ∈ W ] * / 4 for β in W do 4.1 A β ← i<dn A i (β)X i n 5 return [EvalUnivariate(A β , v β ) | β ∈ W ]
C Eval (d 1 , . . . , d n ) ≤ i≤n δ d C(d i ) d i ,
which proves the proposition. 2

Interpolation

Using the same notation as above, the inverse of the evaluation map is interpolation at V :

Interp V : R δ d → Span(M d ) [F (α) | α ∈ V ] → F.
For this map to be well-defined, we impose a natural condition on the points of

V . Let W = π(V ) ∈ R n-1 . We say that V supports interpolation if • if n > 1, W supports interpolation, and
• for all β in W and all x, x in v β , x -x is a unit;

if the base ring is a field, this condition is vacuous. We will see in the following proposition that if V supports interpolation, then the map Interp V is well-defined. Moreover, we let C Interp be such that, for V equiprojectable of multi-degree d, if V supports interpolation, then the map Interp V can be evaluated in C Interp (d) operations (including inversions).

Proposition 2. If V supports interpolation, the map Interp V is well-defined. Besides, one can take C Interp (d) ≤ δ d L(d).
Proof. If n = 0, we do nothing; otherwise, we let W = π(V ). The set of values to interpolate at V has the shape [

f α | α ∈ V ] ∈ R δ d ; we can thus rewrite it as [f β | β ∈ W ],
where each

f β is in R dn . Since V supports interpolation, for β in W , there exists a unique polynomial A β ∈ R[X n ] of degree less than d n , such that Eval(A β , v β ) = f β .
Applying the algorithm recursively on the coefficients of the polynomials A β , we can find a polynomial A such that A(β, X n ) = A β (X n ) holds for all β ∈ W . Then, the polynomial A satisfies our constraints. This provides a right-inverse, and thus a two-sided inverse for the map Eval. The

Interp(f, V ) 1 if n = 0 return [f ] 2 W ← π(V ) 3 for β in W do 3.1 A β ← InterpUnivariate(f β , v β ) 4 for i = 0, . . . , d n -1 do 4.1 c i ← [coeff(A β , X n , i) | β ∈ W ] 4.2 A i ← Interp(c i , W ) 5 return i<dn A i X i n Fig. 3. Interpolation algorithm.
algorithm is given in Figure 3; we use a subroutine called InterpUnivariate for univariate interpolation. As for evaluation, we deduce that we can take C Interp satisfying

C Interp (d 1 , . . . , d n ) ≤ C Interp (d 1 , . . . , d n-1 ) d n + d 1 • • • d n-1 C(d n ),
which gives our claim, as in the case of evaluation. 2

Associated triangular set

Next, we associate to an equiprojectable set V ⊂ R n of multi-degree d = (d 1 , . . . , d n ) a triangular set T = (T 1 , . . . , T n ) of the same multi-degree, which vanishes on V . As soon as V supports interpolation, the existence of T is guaranteed (and is established in the proof of the next proposition). Uniqueness holds as well: if (T 1 , . . . , T n ) and (T 1 , . . . , T n ) both vanish on V and have multi-degree d, then for all i, T i -T i vanishes at V as well and is in Span(M d ); hence, it is zero. We call T the associated triangular set; if R is a field, T is a lexicographic Gröbner basis of the vanishing ideal of V . Proof. We proceed inductively, and suppose that we already have computed T 1 , . . . , T n-1 as the associated triangular set of W = π(V ). We will write d = (d 1 , . . . , d n ) and e = (d 1 , . . . , d n-1 ).

For β in W , let T β be the polynomial a∈v β (X n -a) ∈ R[X n ]. For j < d n , let further T j,n be the polynomial in Span(M e ) that interpolates the jth coefficient of the polynomials T β at W ; for j = d n , we take T dn,n = 1. We then write T n = j≤dn T j,n X j n : this polynomial is in R[X 1 , . . . , X n ], monic of degree d n in X n , has degree less than d i in X i , for i < n, and vanishes on V . Thus, the polynomials T = (T 1 , . . . , T n ) form the triangular set we are looking for. The algorithm is in Figure 4; we use a function PolyFromRoots to compute the polynomials T β .

For a given β in W , the function

PolyFromRoots computes T β in C(d n ) base ring opera- tions; this implies that given T 1 , . . . , T n-1 , one can construct T n using d 1 • • • d n-1 C(d n ) + C Interp (d 1 , • • • , d n-1 )d n operations. The total cost for constructing all T i is thus at most i≤n d 1 • • • d i-1 C(d i ) + i≤n C Interp (d 1 , • • • , d i-1 )d i . AssociatedTriangularSet(V, n) 1 if n = 0 return [] 2 W ← π(V ) 3 (T 1 , . . . , T n-1 ) ← AssociatedTriangularSet(W, n) 4 for β in W do 4.1 T β ← PolyFromRoots(v β ) 5 for j = 0, . . . , d n -1 do 5.1 T j,n ← Interp([coeff(T β , X n , j) | β ∈ W ], W ) 6 return j<dn T j,n X j n + X dn n
δ d i≤n C(d i ) d i + i≤n d 1 • • • d i j≤i-1 C(d j ) d j ≤ δ d i≤n C(d i ) d i + i≤n d 1 • • • d i j≤n C(d j ) d j .
Using the upper bound

i≤n d 1 • • • d i ≤ 2δ d ,
we finally obtain the estimate 3δ d L(d). 2

Multiplication

Using our evaluation and interpolation algorithms, it becomes immediate to perform multiplication modulo a triangular set T associated to an equiprojectable set. Proof. The algorithm is the same as in [22, Section 2.2], except that we now use the more general evaluation and interpolation algorithms presented here. Let A and B be reduced modulo T , and let C = AB mod T . Then for all α in V , C(α) = A(α)B(α). Since C is reduced modulo T , it suffices to interpolate the values A(α)B(α) to obtain C. The cost is thus that of two evaluations, one interpolation, and of all pairwise products; the bounds of Propositions 1 and 2 conclude the proof. 2

Mul(A, B, V ) 1 Val A ← Eval(A, V ) 2 Val B ← Eval(B, V ) 3 Val C ← [ Val A (α)Val B (α) | α ∈ V ] 4 return Interp(Val C , V )

Homotopy techniques for multiplication

Let T be a triangular set in R[X 1 , . . . , X n ]. We saw in the previous section that if T has all its roots in R, and if V (T) supports interpolation, then multiplication modulo T can be done in quasi-linear time. In this section, we extend this approach to an arbitrary T by setting up an homotopy between T and a new, more convenient, triangular set U. This extends the approach of [22, Section 2.2], which dealt with the case where

T i is in R[X i ] for all i.
Let d be the multi-degree of T and assume that there exists an equiprojectable set V in R n which supports interpolation and has multi-degree d. Let U be the triangular set associated to V and let η be a new variable. We then define the set

S in R[[η]][X 1 , . . . , X n ] by S i = ηT i + (1 -η)U i , 1 ≤ i ≤ n.
Since U and T have the same multi-degree d, this set S is triangular, with multi-degree d.

In Subsection 4.1, we prove that S has all its roots in R [[η]]. Thus, we can use evaluation-interpolation techniques to do multiplication modulo S ; this will in turn be used to perform multiplication modulo T .

The algorithm involves computing with power series; the quantity that will determine the cost of the algorithm will be the required precision in η. For e = (e 1 , . . . , e n ) in N n , we define

H 0 (e 1 , . . . , e n ) = deg(X e1 1 • • • X en n mod S ,

Computing the roots of S

We show here that S has all its roots in R[[η]], by a straightforward application of Hensel's lemma.

First, we need some notation. Given positive integers k, and a subset

A ⊂ R[[η]] k , A mod η denotes the set [a mod η | a ∈ A].
Besides, we usually denote objects over R[ [η]] with a superscript, to distinguish them from their counterparts over R. Finally, we extend the notation π to denote the following projection

π : R[[η]] n → R[[η]] n-1 (α 1 , . . . , α n ) → (α 1 , . . . , α n-1 ).
Recall in what follows that V ⊂ R n is equiprojectable of multi-degree d, that its associated triangular set is U, and that S = ηT + (1 -η)U.

Proposition 5. There exists a unique set

V in R[[η]
] n such that the following holds:

• V = V mod η;
• V is equiprojectable of multi-degree d;

• V supports interpolation;

• S is the triangular set associated to V .

Proof. We first claim that for i ≤ n and α = (α 1 , . . . , α n ) in V , the partial derivative

∂S i ∂X i (η = 0, α 1 , . . . , α n ) = ∂U i ∂X i (α 1 , . . . , α n ) = ∂U i ∂X i (α 1 , . . . , α i ) is non zero. Let indeed α = (α 1 , . . . , α i-1 ) ∈ R i-1
. Then, we have by construction

U i (α 1 , . . . , α i-1 , X i ) = a∈v α (X i -a),
so that the previous partial derivative equals

∂U i ∂X i (α 1 , . . . , α i-1 , α i ) = a∈v α ,a =αi (α i -a).
Since V supports interpolation, this quantity is a product of units, so it is a unit as well, establishing our claim.

Since the system S is triangular, its Jacobian determinant is the product of the partial derivatives ∂S i /∂X i . By the previous remark, all these derivatives are units at α, so the Jacobian itself is a unit at α. As a consequence, by Hensel's lemma, for all α in V , there exists a unique α in R[[η]] n such that α = α mod η and S(α ) = 0. We thus let V ⊂ R n be the set of all such α ; hence V = V mod η and S vanishes at V .

Next, we prove that V is equiprojectable of multi-degree d. By induction, we can assume that we have proved that π(V ) is equiprojectable of multi-degree (d 1 , . . . , d n-1 ); it suffices to prove that for each β in π(V ), the fiber V β has cardinality d n .

Let thus α be in V β . We prove that for all γ in V , π(α ) = π(γ ) if and only if π(α) = π(γ), with α = α mod η and γ = γ mod η. To prove our claim, remark first that if π(α ) = π(γ ) then π(α) = π(γ), by reduction modulo η. Conversely, suppose that π(α) = π(γ). Since the system S is triangular, and since α and γ are obtained by lifting α and γ using this system, we deduce that π(α ) = π(γ ), as requested. Thus, V is equiprojectable of multi-degree d.

Finally, we prove that V supports interpolation. This is again done by induction: assume that the projection π(V ) supports interpolation, let β be in π(V ), and let a and a be in v β . By assumption on V , a -a mod η is a unit in R; thus, by Hensel's lemma, a -a is a unit in R[[η]], as requested.

This proves the existence of V with the requested properties. Uniqueness follows in a straightforward manner from the uniqueness property of Hensel's lemma. 2

We continue with complexity estimates: we prove that the roots of S can be computed in quasi-linear time. Proposition 6. Given T, V and > 0, one can compute V mod η in time O(δ d L(d)M( )).

Proof. As before, we proceed inductively: we suppose that the projection W = π(V ) is known modulo η , and show how to deduce V mod η . To do so, we evaluate all coefficients of S n at all points of W modulo η . Then, for each β in W , it suffices to use Hensel's lemma to lift the roots of S n (β , X n ) at precision . The pseudo-code is in Figure 6; for simplicity, we write there W instead of W mod η .

LiftRootsMultivariate(V, S, )

1 n = |S| 2 if n = 0 return [] 3 W ← LiftRootsMultivariate(π(V ), (S 1 , . . . , S n-1 ), ) 4 for i = 0, . . . , d n -1 do 4.1 val i ← Eval(coeff(S n , X n , i), W )
/ * all computations are done modulo η * / 5 for β in W do 5.1 Lemma 1 shows that we can lift the power series roots of a bivariate polynomial of degree d at precision in time C(d)M( ). As a consequence, the overall cost C Lift of the lifting process satisfies

S β ← i<dn val i,β X i n + X dn n 5.2 v β ← LiftRoots(S β , v β , ) 6 return [v β | β ∈ W ]
C Lift (d 1 , . . . , d n , ) ≤ C Lift (d 1 , . . . , d n-1 , ) + C Eval (d 1 , . . . , d n-1 )d n M( ) +d 1 • • • d n-1 C(d n )M( );
the middle term gives the cost of evaluating the coefficients of S n at W mod η (so we apply our evaluation algorithm with power series coefficients); and the right-hand term gives the cost of lifting the roots of S n . This gives

C Lift (d 1 , . . . , d n , ) ≤ i≤n C Eval (d 1 , . . . , d i-1 )d i M( ) + i≤n d 1 • • • d i-1 C(d i )M( ).
As in the proof of Proposition 3, one deduces that the overall sum is bounded by

3δ d i≤n C(d i ) d i M( ) = 3δ d L(d)M( ),
which concludes the proof. 2

Multiplication by homotopy

We continue with the same notation as before. To multiply two polynomials A, B ∈ Span(M d ) modulo T , we may multiply them modulo S over R[η] and let η = 1 in the result. The results of the multiplication modulo S over R[η] and over R[[η]] are the same; when working over R[[η]], we may use the evaluation-interpolation techniques from Subsection 3.5. Indeed, by Proposition 5, S is associated to a subset

V of R[[η]] n that supports interpolation.
Of course, when multiplying A and B modulo S over R[[η]], we cannot compute with (infinite) power series, but rather with their truncations at a suitable order. On the one hand, this order should be larger than the largest degree of a coefficient of the multiplication of A and B modulo S over R[η]. On the other hand, this order will determine the cost of the multiplication algorithm, so it should be kept to a minimum. For e = (e 1 , . . . , e n ) in N n , recall that we defined

H 0 (e 1 , . . . , e n ) = deg(X e1 1 • • • X en n mod S , η)
and H(e 1 , . . . , e n ) = max e 1 ≤e1, ..., e n ≤en H 0 (e 1 , . . . , e n ).

The following proposition relates the cost of our algorithm to the function H; the behavior of this function is studied in the next subsection. Proof. The algorithm is simple: we compute U and use it to obtain V at a high enough precision. In R[X 1 , . . . , X n ], the product AB satisfies deg(AB, X i ) ≤ 2d i -2 for all i ≤ n; since the multiplication algorithm does not perform any division by η, it suffices to apply it with coefficients in R[η]/ η r+1 , with r = H(2d 1 -2, . . . , 2d n -2). The resulting algorithm is given in Figure 7; as before, we write V for simplicity, whereas we should write V mod η r+1 .

Mul(A, B, T, V ) 0. U ← AssociatedTriangularSet(V, n) 1. S ← ηT + (1 -η)U 2. V ← LiftRootsMultivariate(V, S, r + 1) 3. C η ← Mul(A, B, V ) 4. return C η (1, X 1 , . . . , X n ) / * C η is seen in R[η][X 1 , . . . , X n ] * /

Precision analysis

In this section, we study the functions H 0 and H introduced in the previous section, and show how they are related to the monomial support of the polynomials in the triangular set T. Before giving a general bound for these functions, we discuss an example, to motivate our general result.

A worked example

For this example, we assume that n = 2 and that the triangular set S = (S 1 , S 2 ) is as follows:

S 1 = X 3 1 -ηX 2 1 -ηX 1 -η S 2 = X 3 2 -ηX 2 2 -ηX 2 1 .
Here, the multi-degree of S is d = (3, 3). Our question of estimating the functions H 0 and H then essentially boils down to the following: given a monomial m in X 1 , X 2 , what will be the degree in η of the normal form of m modulo S ? Taking m = X 4 1 X 6 2 , we are going to follow the steps of this reduction to show the growth of the degrees.

Reducing m = X 4 1 X 6 2 with respect to S 2 once, we get ηX

4 1 X 5 2 + ηX 6 1 X 3 2 . If we reduce m 1 = ηX 4 1 X 5 2 with respect to S 2 , we get η 4 X 4 1 X 2 2 + • • • + η 4 X 6 1 ; reducing η 4 X 6
1 with respect to S 1 increases the degree in η by 4, reaching a degree 8 (this is maximum we obtain). If we reduce m 2 = ηX 6 1 X 3 2 with respect to S 2 , we obtain η 2 X 6 1 X 2 2 + η 2 X 8 1 ; the reduction of X 8 1 with respect to S 1 increases the degree in η by 6, and the resulting degree will be 8 as well.

The reductions with respect to S 2 amount to replace X 3 2 by ηX 2 2 + ηX 2 1 . They can be described graphically using a tree: the initial monomial m is the root, and the children of a node correspond to the two monomials that appear after one reduction step. This enables us to dispense with writing the η factors: the degree in η of any monomial is the length of the path from the root to the corresponding node. Figure 8 shows this tree for the initial monomial X 4 1 X 6 2 .

X 4 1 X 6 2 X 4 1 X 5 2 X 4 1 X 4 2 X 4 1 X 3 2 X 4 1 X 2 2 X 6 1 X 6 1 X 2 X 6 1 X 2 2 X 6 1 X 3 2 X 8 1
Fig. 8. Tree structure of the reduction by S2.

In this tree, one step to the left reduces the degree in X 2 by 1; one step to the right reduces the degree in X 2 by 3 and increases the degree in X 1 by 2. We perform reductions as long as the degree in X 2 remains at least 3. When we reach a leaf, with degree r in X 1 , we increase the degree in η by r -2 through the reduction by S 1 .

With these rules, our goal is to find an upper bound on the degree in η in the result. We introduce two variables f 1 and f 2 , that respectively represent the number of steps to the left and to the right. After (f 1 , f 2 ) steps, the degree in X 1 is 4 + 2f 2 and the degree in X 2 is 6 -f 1 -3f 2 , so that the end-points of our walk satisfy 0 ≤ 6 -f 1 -3f 2 ≤ 2. Figure 9 depicts this situation in the plane of coordinates (f 1 , f 2 ); the lower oblique line has equation 6 -f 1 -3f 2 = 2 and the higher one has equation 6 -f 1 -3f 2 = 0. Also, at each point with integer coordinates (f 1 , f 2 ), we indicate the degrees in (X 1 , X 2 ); the origin (f 1 , f 2 ) = (0, 0) is in the lower corner on the left. The points in bold face can be reached by our reduction process. This walk only describes reduction with respect to S 2 ; we still need to reduce with respect to S 1 . This is done by remembering that the reduction of a monomial X r1 1 X r2 2 with respect to S 1 increases the degree in η by r 1 -2. Thus, starting from X 4 1 X 6 2 , the degree in η we obtain after reduction is bounded by the maximal value of the sum (f 1 + f 2 ) + (2 + 2f 2 ) = 2 + f 1 + 3f 2 over all (f 1 , f 2 ) that lie in between the two oblique lines: the term (f 1 + f 2 ) gives the degree increase due to the number of steps (reduction with respect to S 2 ), and the term 2 + 2f 2 = 4 + 2f 2 -2 gives the degree increase due to the reduction with respect to S 1 (remember that the degree in X 1 is 4 + 2f 2 ).

To simplify the maximization, we can forget the lower constraint 6 -f 1 -3f 2 ≤ 2, and maximize in the larger area defined by f 1 ≥ 0, f 2 ≥ 0, 6 -f 1 -3f 2 ≥ 0; we can also remove the condition that f 1 and f 2 are integers. Indeed, both simplifications increase the size of the search space, so any upper-bound obtained after simplifications will be an upper bound for the initial problem as well.

After these simplifications, we are left to maximize a linear function over a simplex; the maximum will thus be obtained at one of the vertices. The vertex (f 1 , f 2 ) = (0, 0) is not a maximum (it is the minimum); the other vertices are (f 1 , f 2 ) = (0, 2) and (f 1 , f 2 ) = (6, 0), and the maximum of 2 + f 1 + 3f 2 is 8, obtained at both vertices (in this case, we obtain the exact maximum for our original problem; in most cases, this won't be guaranteed).

General analysis

We will now formalize the former considerations, showing how the monomial structure of the polynomials in S affects the cost of the algorithm. For i ≤ n and ν = (ν 1 , . . . , ν i ) in N i , we will use the notation X ν i = X ν1 1 • • • X νi i and we write the monomial expansion of S i as

S i = X di i + ν∈Ei s ν X ν i , (4)
where E i is the set of exponents that appear in S i , the exponents ν are in N i , and s ν is linear in η. Let us further introduce the coefficients h i defined by h 0 = 0 and for i ≥ 1,

h i = max ν∈Ei h 1 ν 1 + • • • + h i-1 ν i-1 + 1 d i -ν i . ( 5
)
One easily checks that all h i are positive. The following proposition shows that through the coefficients h i , the support E i determines the cost of our algorithm.

Proposition 8. The inequality

H(e 1 , . . . , e n ) ≤ h 1 e 1 + • • • + h n e n
holds for all (e 1 , . . . , e n ) ∈ N n .

Using Proposition 7, this proposition gives as an easy corollary the following statement, where we take e i = 2d i -2 ≤ 2d i ; we still use the previous notation T and V .

Corollary 1. Given A, B, T and V , one can compute AB mod T in time O(δ d L(d)M(r)) ⊂ O˜(δ d r), with r ≤ 2(h 1 d 1 + • • • + h n d n ).
Hence, the lower the h i the better. However, without putting extra assumptions on the monomial supports E i , Corollary 1 only yields estimates of little interest. Even in sparse cases, it remains difficult to simplify the recurrence giving the coefficients h i . Still, several examples in the next section will show that for some useful families of monomial supports, significantly sharper bounds can be derived. The rest of this section is devoted to prove Proposition 8. In all that follows, the multi-degree d = (d 1 , . . . , d n ) and the supports E i are fixed. We also let E i be the set of modified exponents

E i = {ν -(0, . . . , 0, d i ) ∈ Z i | ν ∈ E i }
, so that for all ν = (ν 1 , . . . , ν i ) in E i , ν j ≥ 0 for j < i and ν i < 0. Hence, Equation (5) takes the (slightly more handy) form

h i = max ν∈E i h 1 ν 1 + • • • + h i-1 ν i-1 + 1 -ν i . ( 6
)
Since the proof is rather technical, we will first give a roadmap of it.

Recall that the function H 0 was defined with domain N n ; in what follows, we also see it as a function over N i , for 1 ≤ i ≤ n, by defining H 0 (e 1 , . . . , e i ) = H 0 (e 1 , . . . , e i , 0, . . . , 0), where the right-hand expression contains n -i zeros; for completeness, we write H 0 () = 0 for i = 0. Then, for e = (e 1 , . . . , e i ) in N i , e = (e 1 , . . . , e i-1 ) in N i-1 and i ≥ 1, Lemma 2 proves that

H 0 (e) = H 0 (e ) if e i < d i , H 0 (e) ≤ 1 + max ν∈E i H 0 (e + ν) otherwise;
this translates the fact that if X e is reduced with respect to S i , we can drop the last variable X i ; otherwise, we do one reduction step by S i . Iterating this process, Lemma 3 proves that

H 0 (e) ≤ max (fν ) ν∈E i non-negative integers such that 0 ≤ ei + ν∈E i fν νi ≤ di -1 H 0 (e + ν∈E i f ν ν) + ν∈E i f ν .
This expression generalizes the example seen before, introducing integer-valued variables that count the reduction steps we do; the linear constraints we obtain are of the same nature as the ones seen in the example. Finally, Lemma 4 will "solve" this optimization problem, by observing (as in the example) that we can find a reasonable bound on the maximum by inspecting the vertices of a simplex.

We will now prove the inequalities showed before. The following recurrence relation enables us to control the growth of H 0 .

Lemma 2. For i ≥ 1, let e = (e 1 , . . . , e i ) be in N i and let e = (e 1 , . . . , e i-1 ) in N i-1 . Then the following (in)equalities hold:

H 0 (e) = H 0 (e ) if e i < d i , H 0 (e) ≤ 1 + max ν∈E i H 0 (e + ν) otherwise.
Proof. Let us first suppose e i < d i ; then,

X e1 1 • • • X ei i mod S = (X e1 1 • • • X ei-1
i-1 mod S )X ei i , since the latter product is reduced modulo S . Both sides have thus the same degree in η, and our first claim follows.

We can now focus on the case e i ≥ d i , for which we write f i = e i -d i , so that f i ≥ 0. From Equation (4), we deduce

X fi i S i = X fi+di i + ν∈Ei s ν X fi i X ν i ,
and thus we get

X fi i S i = X ei i + ν∈E i s ν X ei i X ν i ,
by the definition of E i . In our notation, we have

X e1 1 • • • X ei i = X e i .
Thus, after multiplication by

X e1 1 • • • X ei-1
i-1 and term reorganization, the former equality implies that

X e i -X e1 1 • • • X ei-1 i-1 X fi i S i = - ν∈E i s ν X e+ν i .
As a consequence, we deduce that deg(X e i mod S , η) ≤ max

ν∈E i deg(s ν X e+ν i mod S , η). Since for ν in E i , we have deg(s ν X e+ν i mod S , η) = deg(s ν , η) + deg(X e+ν i mod S , η) = 1 + H 0 (e + ν),
the conclusion follows. 2

Iterating the process of the previous lemma, we obtain the following bound. In the next lemma, (f ν ) ν∈E i is a family of integer valued variables. Lemma 3. Let e = (e 1 , . . . , e i ) be in N i . Then the following inequality holds:

H 0 (e) ≤ max (fν ) ν∈E i non-negative integers such that 0 ≤ ei + ν∈E i fν νi ≤ di -1 H 0 (e + ν∈E i f ν ν) + ν∈E i f ν .
Proof. We prove the claim by induction on e i . For e i ≤ d i -1, the family (f ν = 0) ν∈E i satisfies the constraint 0 ≤ e i + ν∈E i f ν ν i ≤ d i -1; for this choice, the value of the function we maximize is precisely H 0 (e), so our claim holds. Suppose now that e i ≥ d i .

Then, the previous lemma gives

H 0 (e) ≤ 1 + max ν∈E i H 0 (e + ν). (7)
Let us fix ν in E i ; then e + ν has non-negative integer coordinates, and its ith coordinate is less than e i . Thus, we can apply the induction assumption, obtaining

H 0 (e+ν) ≤ max (f ν ) ν ∈E i non-negative integers such that 0 ≤ ei + νi + ν ∈E i f ν ν i ≤ di -1 H 0 (e+ν+ ν ∈E i f ν ν )+ ν ∈E i f ν .
To any set of non-negative integers (f ν )

ν ∈E i with 0 ≤ e i + ν i + ν ∈E i f ν ν i ≤ d i -1
appearing in the previous maximum, we associate the non-negative integers (f ν ) ν ∈E i , with f ν = f ν + 1 and f ν = f ν otherwise. These new integers satisfy

0 ≤ e i + ν ∈E i f ν ν i ≤ d i -1
and

H 0 (e + ν + ν ∈E i f ν ν ) + ν ∈E i f ν = H 0 (e + ν ∈E i f ν ν ) + ν ∈E i f ν -1.
Taking maxima, we deduce from the previous inequality

H 0 (e + ν) ≤ max (f ν ) ν ∈E i non-negative integers such that 0 ≤ ei + ν ∈E i f ν ν i ≤ di -1 H 0 (e + ν ∈E i f ν ν ) + ν ∈E i f ν -1.
Substituting in Equation ( 7) and taking the maximum over ν in E i concludes the proof. 2

For i ≤ n, let L i be the linear form (e 1 , . . . , e i ) → h 1 e 1 + • • • + h i e i ,
where the h i are as in Equation ( 5). The following lemma concludes the proof of Proposition 8; as we did for H 0 , for i ≤ n, we extend H to N i , by writing H(e 1 , . . . , e i ) = H(e 1 , . . . , e i , 0, . . . , 0). Lemma 4. For i ≤ n and e = (e 1 , . . . , e i ) in N i , the inequality H(e) ≤ L i (e) holds.

Proof. It is sufficient to prove that H 0 (e) ≤ L i (e) holds; since all coefficients of L i are non-negative, L i is non-decreasing with respect to all of its variables, which implies the thesis.

We prove our inequalities by induction on i ≥ 0. For i = 0, we have H 0 () = L 0 () = 0; hence, our claim vacuously holds at this index. For i ≥ 1, we now prove that if our inequality holds at index i -1, it will also hold at index i. Lemma 3 shows that for any e ∈ N i , we have the inequality

H 0 (e) ≤ max (fν ) ν∈E i non-negative integers such that 0 ≤ ei + ν∈E i fν νi ≤ di -1 H 0 (e + ν∈E i f ν ν) + ν∈E i f ν .
Let ϕ be the natural projection N i → N i-1 , let (f ν ) ν∈E i be non-negative integers that satisfy the conditions in the previous inequality. Since e + ν∈E i f ν ν has degree in X i less than d i , the first point of Lemma 2 shows that

H 0 (e + ν∈E i f ν ν) = H 0 (ϕ(e + ν∈E i f ν ν));
the induction assumption implies that this quantity is bounded from above by

L i-1 (ϕ(e + ν∈E i f ν ν)).
As a consequence, H 0 (e) admits the upper bound max

(fν ) ν∈E i non-negative integers such that 0 ≤ ei + ν∈E i fν νi ≤ di -1 L i-1 (ϕ(e + ν∈E i f ν ν)) + ν∈E i f ν .
This quantity itself is upper-bounded by a similar expression, where we allow the f ν to be non-negative reals numbers; this gives

H 0 (e) ≤ max (fν ) ν∈E i non-negative real numbers such that 0 ≤ ei + ν∈E i fν νi ≤ di -1 L i-1 (ϕ(e + ν∈E i f ν ν)) + ν∈E i f ν .
Since all h i and all ν 1 , . . . , ν i-1 are non-negative, the function of (f ν ) ν∈E i we want to maximize is affine with non-negative coefficients. The domain where we maximize it is defined by the conditions

f ν ≥ 0 for all ν ∈ E i , 0 ≤ e i + ν∈E i f ν ν i ≤ d i -1,
and it is contained in the domain D defined by the conditions

f ν ≥ 0 for all ν ∈ E i , 0 ≤ e i + ν∈E i f ν ν i .
Since all unknowns f ν are non-negative, while the coefficients ν i are negative, the domain D is convex and bounded. Hence, the maximal value we look for is upper-bounded by the maximal value at the end-vertices of D, distinct from the origin; these vertices are

E ν = {f ν = 0 for ν = ν, f ν = - e i ν i }, for ν ∈ E i .
At the point E ν , the objective function takes the value

L i-1 (ϕ(e - e i ν i ν)) - e i ν i .
By the linearity of L i-1 and ϕ, this can be rewritten as

L i-1 (ϕ(e)) -L i-1 (ϕ( e i ν i ν)) - e i ν i = L i-1 (ϕ(e)) - L i-1 (ϕ(ν)) + 1 ν i e i .
As a consequence, we obtain the upper bound

H 0 (e) ≤ L i-1 (ϕ(e)) + max ν∈E i L i-1 (ϕ(ν)) + 1 -ν i e i .
To simplify this further, note that the term L i-1 (ϕ(e)) rewrites as

h 1 e 1 + • • • + h i-1 e i-1 . Similarly, L i-1 (ϕ(ν)) + 1 equals h 1 ν 1 + • • • + h i-1 ν i-1 .
We deduce the inequality

H 0 (e) ≤ h 1 e 1 + • • • + h i-1 e i-1 + max ν∈E i h 1 ν 1 + • • • + h i-1 ν i-1 + 1 -ν i e i ,
which we can finally rewrite as

H 0 (e) ≤ h 1 e 1 + • • • + h i-1 e i-1 + h i e i ,
as requested. 2

Examples

Main family of examples

We give explicit estimates for the coefficients h i of the previous section on the following family of examples. We consider triangular sets T = (T 1 , . . . , T n ) such that T i has the form

T i = X di i + ν∈Di t ν X ν i , t ν ∈ R, (8)
where all X ν i are monomials in X 1 , . . . , X i of total degree at most λ i , for some λ i ∈ N. We let d = max i≤n d i , and we suppose that R contains at least d pairwise distinct values x 1 , . . . , x d , with x i -x j a unit for i = j.

The following proposition illustrates three different situations. The first two cases display a cost quasi-linear in dδ d , which is satisfying, especially for small d; the last one shows that small changes in the assumptions can induce large overheads. We will see in the next subsection cases where d i is constant equal to d, or

d i = n + 1 -i; in such cases, d is logarithmic in δ d and the cost O˜(dδ d ) is thus O˜(δ d ),
which is what we were aiming at.

Proposition 9. With assumptions as above, multiplication modulo T can be performed with the following complexities:

O n δ d C(d) d M(nd) ⊂ O˜(dδ d ) if λ i = d i -1 for all i, O n δ d C(d) d M(n 2 d) ⊂ O˜(dδ d ) if λ i = d i for all i, O n δ d C(d) d M(2 n d) ⊂ O˜(2 n dδ d ) if λ i = d i + 1 for all i.
Proof. First, we construct V : we simply choose the grid

V = [x 1 , . . . , x d1 ] × • • • × [x 1 , . . . , x dn ]. (9)
Thus, we have

U i = (X i -x 1 ) • • • (X i -x di )
; as before we let S = ηT + (1 -η)U. Thus, the monomial support E i associated with S i is contained in

D i = D i ∪ {(0, . . . , 0, ν i ) | 0 ≤ ν i < d i }.
Since each monomial in D i has an exponent of the form (ν 1 , . . . , ν i ), with ν 1 +• • •+ν i ≤ λ i and ν i < d i , we deduce from Equation ( 5) that

h i ≤ max ν∈D i h 1 ν 1 + • • • + h i-1 ν i-1 + 1 d i -ν i ≤ max max ν∈Di h 1 ν 1 + • • • + h i-1 ν i-1 + 1 d i -ν i , 1 .
Let h i = max(h 1 , . . . , h i ), so that

h i ≤ max max ν∈Di h i-1 (ν 1 + • • • + ν i-1 ) + 1 d i -ν i , 1 ≤ max max ν∈Di h i-1 (λ i -ν i ) + 1 d i -ν i , 1 .
(10) Knowing the distribution of the d i and λ i , the former relation makes it possible to analyze the growth of the coefficients h i , and thus of 2(d

1 h 1 + • • • + d n h n ).
Case 1. Suppose first that λ i = d i -1. Then, the former inequality implies h i ≤ 1 for all i, so that 2(d

1 h 1 + • • • + d n h n ) ≤ 2nd.
Case 2. If λ i = d i , then (10) becomes h i ≤ h i-1 + 1, so that h i ≤ i for all i, and thus 2(d

1 h 1 + • • • + d n h n ) ≤ n(n + 1)d.
Case 3. If finally λ i = d i + 1, then (10) becomes h i ≤ 2h i-1 + 1, so that h i ≤ 2 i -1.

In this case, we get 2(d

1 h 1 + • • • + d n h n ) ≤ 2 n+2 d.
To conclude the proof, we simply plug the previous estimates in the cost estimate

O(δ d L(d)M(r)) of Corollary 1, with r ≤ 2(d 1 h 1 + • • • + d n h n ),

Cauchy modules

Cauchy modules [START_REF] Rennert | Calcul de résolvantes avec les modules de Cauchy[END_REF] are a basic construction in Galois theory and invariant theory [START_REF] Sturmfels | Algorithms in invariant theory[END_REF][START_REF] Rennert | Calcul de résolvantes avec les modules de Cauchy[END_REF][START_REF] Abdeljaouad | Computation of the decomposition group of a triangular ideal[END_REF][START_REF] Renault | A modular algorithm for computing the splitting field of a polynomial[END_REF]. Starting from a monic polynomial F ∈ R[X] of degree d, we define a triangular set F 1 , . . . , F d by letting F 1 (X 1 ) = F (X 1 ) and taking iterated divided differences:

F i+1 (X 1 , . . . , X i+1 ) = F i (X 1 , . . . , X i-1 , X i ) -F i (X 1 , . . . , X i-1 , X i+1 ) X i -X i+1 1 ≤ i < d.
The polynomials F 1 , . . . , F d form a triangular set of multi-degree d = (d, d -1, . . . , 1), so that δ d = d!; their interest stems from the fact that they form a system of generators of the ideal (σ i -(-1) i f d-i ) 1≤i≤d , where σ i is the ith elementary symmetric polynomial in X 1 , . . . , X d and f i is the coefficient of X i in F . One easily checks that F i has total degree at most d + 1 -i. Hence, assuming that 0, . . . , d-1 are units in R, we are under the assumptions of Subsection 5.1, with λ i = d i = d + 1 -i for all i and (x 1 , . . . , x d ) = (0, . . . , d -1). As a consequence, Proposition 9 shows that multiplication modulo F 1 , . . . , F d can be done using O d! C(d) M(d 3 ) operations in R, that is, in quasi-linear time O˜(d!). This improves for instance the results given in [START_REF] Gaudry | Evaluation properties of symmetric polynomials[END_REF] on the evaluation properties of symmetric polynomials.

Polynomial multiplication

We show now how to derive quasi-linear time algorithms for univariate multiplication in R[X] from our previous multivariate construction. Unfortunately, our algorithm does not improve on the complexity of Cantor-Kaltofen's algorithm [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF]; however, we believe it is worth mentioning. Precisely, given n ≥ 1, we give here an algorithm to perform truncated multiplication in R[X]/ X 2 n . We introduce variables X 1 , . . . , X n ; computing in

A = R[X]/ X 2 n is equivalent to computing in B = R[X 1 , . . . , X n ]/ V 1 , . . . , V n , with V = (V 1 , . . . , V n ) given by X 1 -X 2 n-1 n . . . X n-1 -X 2 n X 2 n
n , since the dummy variables X 1 , . . . , X n-1 play no role in this representation. However, changing the order of the variables, we see that the ideal V 1 , . . . , V n is also equal to the ideal T 1 , . . . , T n given by

X 2 n -X n-1 . . . X 2 2 -X 1 X 2 1 .
The R-basis of B corresponding to V is (X i n ) i<2 n ; the basis corresponding to T is M d (notation defined in the introduction), with d = (2, . . . , 2). Besides, the change of basis does not use any arithmetic operation, since it amounts to rewrite the exponents i in base 2, and conversely.

Hence, we can apply our multivariate multiplication algorithm modulo T . Remark that the triangular set T satisfies the assumptions of Subsection 5.1 (for any R), with [START_REF] Abdeljaouad | Computation of the decomposition group of a triangular ideal[END_REF]. By Proposition 9, we deduce that the cost of a multiplication in B, and thus in A, is O(2 n nM(n)). Since one can multiply univariate polynomials of degree 2 n using two multiplications in A, this gives the recurrence

d 1 = • • • = d n = d = 2, δ d = 2 n , λ 1 = • • • = λ n = 1 and (x 1 , x 2 ) = (0,
M(2 n ) ≤ k2 n nM(n) and thus M(d) ≤ k d log(d)M(log(d))
for some constants k, k . Unrolling the recursion 1, 2, . . . , times, and taking M(n) ∈ O(n 2 ) to end the recursion, we obtain quasi-linear estimates of the form

M(d) ∈ O(d log(d) 3 ) or M(d) ∈ O(d log(d) 2 log(log(d)) 3 ), . . .
The main noteworthy feature of this multiplication algorithm is that no root of unity is present, though our multivariate evaluation-interpolation routine is somewhat similar to a multivariate Fourier Transform. In particular, the case when 2 is a zero-divisor in R requires no special treatment, contrary to [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF].

Exponential generating series multiplication

We continue with a question somehow similar to the one in the previous subsection. Given two sequences a 0 , . . . , a d and b 0 , . . . , b d in R, we want to compute the sequence c 0 , . . . , c d such that

c k = i+j=k k i a i b j , (11)
where the binomial coefficients are the coefficients of the expansion of (1 + X)

c i i! X i = i≤d a i i! X i i≤d b i i! X i mod X d+1 , (12)
so we can achieve a cost O(M(d)). Under some much milder assumptions on R, we are going to see how to achieve a similar cost through multivariate computations.

We will suppose that there exists a prime p such that for a ∈ N, if gcd(a, p) = 1, then a is a unit in R (this is the case e.g. for R = Z/p k Z). Let n be such that d + 1 ≤ p n , and introduce the triangular set T = (T 1 , . . . , T n ) defined by

X p n -pX n-1 . . . X p 2 -pX 1 X p 1 .
In what follows, for i ≥ 0, (i 0 , i 1 , . . . ) denotes the sequence of its coefficients in base p; thus, for i ≤ d, only i 0 , . . . , i n-1 can be non-zero. Besides, we let f : N → N be defined by f (i) = i!/p v (i!) , where v(i!) is the p-adic valuation of i!. In particular, f (i) is a unit in R.

Proposition 10. Let A = i≤d a i f (i) X i0 n • • • X in-1 1 , B = i≤d b i f (i) X i0 n • • • X in-1 1 , C = i≤d c i f (i) X i0 n • • • X in-1 1
.

Then C = AB mod T .

Proof. Let i, j ≤ d, with k = i + j ≤ d. We start by the obvious remark that

k i = f (k) f (i)f (j) p v(( k i )) (13)
holds in R. Besides, the normal form of the product ai f

(i) X i0 n • • • X in-1 1 by bj f (j) X j0 n • • • X jn-1 1 modulo T is a i b j f (i)f (j) p c(i,j) X k0 n • • • X kn-1 1
, where c i,j is the number of carries held in the addition of i and j in base p. From [31, Eq. (1.6)], c i,j is exactly the valuation of the binomial coefficient k i . Thus, by [START_REF] Foursov | On computer-assisted classification of coupled integrable equations[END_REF], the former product equals k i

a i b j f (k) X k0 n • • • X kn-1 1
.

Summing over all i, j gives our claim. 2

As in the previous subsection, we can apply our multivariate multiplication algorithm modulo T . Remark that the triangular set T satisfies the assumptions of Subsection 5.1, with

d 1 = • • • = d n = d = p, δ d = p n , λ 1 = • • • = λ n =

Application: computing with algebraic numbers

We finally present an application of the previous constructions to computation with algebraic numbers, and give timings of our implementation.

Presentation of the problem

Let k be a field and let f and g be monic polynomials in k[T ], of degrees m and n respectively. We are interested in computing their composed sum h = f ⊕ g. This is the polynomial of degree d = mn defined by

f ⊕ g = α,β (T -α -β),
the product running over all the roots α of f and β of g, counted with multiplicities, in an algebraic closure k of k.

A natural approach consists in computing h(T ) as the resultant of f (T -U ) and g(U ) in U . However, the fastest algorithm for resultants [START_REF] Reischert | Asymptotically fast computation of subresultants[END_REF] has a complexity of order O˜(d 1.5 ) for m = n. To do better, Dvornicich and Traverso [START_REF] Dvornicich | Newton symmetric functions and the arithmetic of algebraically closed fields[END_REF] suggested to compute the power sums

a i = f (α)=0 α i , b i = g(β)=0 β i
of respectively f and g, and deduce the power sums c i of h, by means of Equation [START_REF] Dahan | On the complexity of the D5 principle[END_REF].

In [5], this approach is showed to take time O(M(d)), over fields of characteristic zero or larger than d. Indeed, computing (a i ) i≤d and (b i ) i≤d can be done in O(M(d)) operations, over any field, using Newton iteration for power series division [START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF]. Then, by our assumption on the characteristic, one can compute (c i ) i≤d in quasi-linear time using Equation ( 12), for another M(d) + O(d) operations. Finally, knowing (c i ) i≤d , one can then recover h in time O(M(d)) as well, using fast exponential computation [START_REF] Brent | Multiple-precision zero-finding methods and the complexity of elementary function evaluation[END_REF][START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF][START_REF] Van Der Hoeven | Newton's method and FFT trading[END_REF]7]; this step relies as well on the assumption on the characteristic. If k has positive characteristic less than d, two issues arise: Equation ( 12) makes no sense anymore and (c i ) i≤d are actually not enough to recover h. To our knowledge, no general solution better than the resultant method was known up to now (partial answers are in [5,[START_REF] Schost | Multivariate power series multiplication[END_REF] under restrictive conditions). We propose here a solution that works over finite fields, following an idea introduced in [START_REF] González-Vega | Computing with Newton sums in small characteristic[END_REF].

For simplicity, we consider only k = F p . Since our algorithm actually does computations over rings of the form Z/p α Z, measuring its complexity in F p -operations as we did up to now is not appropriate. Instead, we count bit operations. Thus, we let M Z be such that integers of bit-length can be multiplied using M Z ( ) bit operations; quasi-linear estimates are known as well for M Z , the best to date being Fürer's log( )2 O(log * ( )) [START_REF] Fürer | Faster integer multiplication[END_REF]. Proof. Let Z p be the ring of p-adic integers and let F and G be monic lifts of f and g in Z p [T ], of degrees m and n. Defining H = F ⊕ G ∈ Z p [T ], we have that h = H mod p. Let further (A i ) i≥0 , (B i ) i≥0 and (C i ) i≥0 be the power sums of respectively F , G and H. For any α ≥ 0, the reductions A i mod p α , B i mod p α , and C i mod p α satisfy Equation ( 11), so we can apply the results of Subsection 5.4 to deduce (C i mod p α ) i≤d from (A i mod p α ) i≤d and (B i mod p α ) i≤d .

Besides, taking α = log p (d) + 1, it is proved in [6] that given (C i mod p α ) i≤d , one can compute h in quasi-linear time O(M(d)M Z (log p (d))) bit operations. Remark that this step is non trivial: recovering a polynomial of degree d from its Newton sums requires divisions by 1, . . . , d, and not all these numbers are units in small characteristic.

In the algorithm, the function Lift simply lifts its argument from F p [T ] = Z/pZ[T ] to Z/p α Z[T ]; the function PowerSums computes the first d power sums of its arguments by the algorithm of [START_REF] Schönhage | The fundamental theorem of algebra in terms of computational complexity[END_REF]. Step 7 applies the algorithm of Subsection 5.4, and the last step uses the algorithm presented in [6] to recover h.

Our choice of α implies that log(p α ) = O(log(d)). Thus, operations (+, ×) modulo p α take O(M Z (log(d))) bit operations [START_REF] Zur Gathen | Modern Computer Algebra[END_REF]Chapter 9]. Using Newton iteration, inversions modulo p α take N(p, d) = O(M Z (log(p)) log(log(p)) + M Z (log(d))) bit operations, where the first term stands for the cost computing the inverse modulo p, and the second one for lifting it modulo p α .

The cost of computing (A i ) i≤d and (B i ) i≤d is O(M(d)) operations modulo p α ; this dominates the cost of recovering h. The remaining cost is that of computing (C i ) i≤d , which is reported in Subsection 5.4 in terms of numbers of operations modulo p α . The previous estimate on N(p, d) concludes the proof. 2 Remark that a more general question is to compute α,β (T -q(α, β)), where q is any polynomial in k[X, Y ], and where the product is taken as before over all roots α of f and β of g. Here, we dealt with q = X + Y ; the case q = XY is actually simpler, and is discussed in [START_REF] Dvornicich | Newton symmetric functions and the arithmetic of algebraically closed fields[END_REF]5,6]. However, in the general case, even in characteristic zero, no quasi-linear time algorithm is known as of now.

Experimental results

We implemented the composed sum algorithm over F 2 (i.e., p = 2 here). We used the NTL C++ package as a basis [START_REF] Shoup | NTL: A library for doing number theory[END_REF]. Since NTL does not implement bivariate resultants, we also used Magma [4] for comparison with the resultant method. All timings are obtained on an AMD Athlon 64 with 5GB of RAM.

Figure 11 gives detailed timings for our algorithm; each colored area gives the time of one of the main tasks. The less costly step is the first, the conversion from the original polynomials to their Newton sums. Then, we give the time needed to compute all the power series roots needed for our multiplication algorithm, followed by the evaluationinterpolation process itself; finally, we give the time necessary to recover h from its power sums. Altogether, the practical behavior of our algorithm matches the quasi-linear complexity estimates. The steps we observe correspond to the increase in the number of variables in our multivariate polynomials, and are the analogues of the steps observed in classical FFT.

Figure 12 gives timings obtained in Magma, using the built-in resultant function, on the same set of problems as above. As predicted by the complexity analysis, the results are significantly slower (about two orders of magnitude for the larger problems).

Conclusion

Several questions remain open after this work. First, it should be stressed that for triangular sets in few variables (say up to 4), the approach discussed here is not expected to be competitive with the general purpose algorithms: our advantages appear in structured situations. Indeed, the most challenging open problem remains how to unconditionally get rid of all exponential factors in multiplication algorithms for triangular sets. More immediate questions may be the following: at the fine tuning level, adapting the idea of the Truncated Fourier Transform [START_REF] Van Der Hoeven | The Truncated Fourier Transform and applications[END_REF] should enable us to reduce the step effect in the timings of the previous section. Besides, it will be worthwhile to investigate what other applications can be dealt with using the "homotopy multiplication" model, such as the product of matrices with entries defined modulo a triangular set, or further tasks such as modular inversion or modular composition.

Fig. 1 .

 1 Fig. 1. Lifting all roots of a bivariate polynomial.

  with d = max i≤n d i . In view of the estimates C(d) ∈ O(M(d) log(d)) and M(d) ∈ O(d lg(d) lg lg(d)), we also have the upper bound L(d) ∈ O(lg(δ d ) 3 ), which shows that L(d) is of polylogarithmic growth in δ d .

Fig. 2 .

 2 Fig. 2. Evaluation algorithm.

Proposition 3 .

 3 Given an equiprojectable set V of multi-degree d that supports interpolation, one can construct the associated triangular set T in time O(δ d L(d)).

Fig. 4 .

 4 Fig. 4. Associated triangular set of V .

Proposition 4 .

 4 Let V ⊂ R n be an equiprojectable set of multi-degree d = (d 1 , . . . , d n ) that supports interpolation, and let T be the associated triangular set. Then one can perform multiplication modulo T in time O(δ d L(d)).

Fig. 5 .

 5 Fig. 5. Multiplication algorithm.

Fig. 6 .

 6 Fig. 6. Lifting the roots of S.

Proposition 7 .

 7 Given A, B, T and V , one can compute AB mod T in time O(δ d L(d)M(r)) ⊂ O˜(δ d r), with r = H(2d 1 -2, . . . , 2d n -2).

Fig. 7 .

 7 Fig. 7. Multiplication algorithm.

Fig. 9 .

 9 Fig. 9. A view of the reduction tree in the coordinates (f1, f2).

  and we use the upper bound L(d) ≤ nC(d)/d of Equation (3). 2

  1 and (x 1 , . . . , x p ) = (0, . . . , p -1). Note as well that we can take n ∈ O(log p (d)), and that δ d = p n ≤ pd. By Proposition 9, we deduce that the cost of computing C, and thus all c 0 , . . . , c d , is O(d log(d) M(p) M(p log p (d))). If p is fixed, we obtain the estimate O(d log(d)M(log(d))). This is not as good as the estimate O(M(d)) = O(d log(d) log log(d)) we obtained in characteristic zero, but quite close.

Proposition 11 .

 11 Given f and g, one can compute h usingO (M(d) + d log(d) M(p) M(p log d (p))) N(p, d) bit operations, with N(p, d) = O(M Z (log(p)) log(log(p)) + M Z (log(d))).After simplification, this cost is seen to be O˜(dp 2 ) bit operations. Also, if we consider p fixed, the cost becomesO (M(d) + d log(d)M(log(d))) M Z (log(d)) ,that is, quasi-linear.

Fig. 10 .

 10 Fig. 10. Composed sum in small characteristic.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Detailed timings for our algorithm.

  i in R[X].We discuss an application of this question in the next section.The naive algorithm has cost O(d 2 ). If 1, . . . , d are units in R, the former equation takes the form

	i≤d

This research was supported by ANR Gecko, the joint Inria-Microsoft Research Lab, MITACS, NSERC and the Canada Research Chair program.