Jean-Guillaume Dumas

Bruno Grenet

In-place fast polynomial modular remainder

We consider the simultaneously fast and in-place computation of the Euclidean polynomial modular remainder R(X) ≡ A(X) mod B(X) with A and B of respective degrees n and m ≤ n. But fast algorithms for this usually come at the expense of (potentially large) extra temporary space. To remain in-place a further issue is to avoid the storage of the whole quotient Q(X) such that A = BQ + R. If the multiplication of two polynomials of degree k can be performed with M(k) operations and O(k) extra space, and if it is allowed to use the input space of A or B for intermediate computations, but putting A and B back to their initial states after the completion of the remainder computation, we here propose an in-place algorithm (that is with its extra required space reduced to O(1) only) using at most O n m M(m) log(m) arithmetic operations, if M(m) is quasi-linear, or O n m M(m) otherwise. We also propose variants that compute -still in-place and with the same kind of complexity bounds -the over-place remainder A(X) ≡ A(X) mod B(X), the accumulated remainder R(X) += A(X) mod B(X) and the accumulated modular multiplication R(X) += A(X)C(X) mod B(X). To achieve this, we develop techniques for Toeplitz matrix operations which output is also part of the input. Fast and in-place accumulating versions are obtained for the latter, and thus for convolutions, and then used for polynomial remaindering. This is realized via further reductions to accumulated polynomial multiplication, for which fast in-place algorithms have recently been developed.

Introduction

Modular methods with dense univariate polynomials over a finite ring are of central importance in computer algebra and symbolic computation. For instance, they are largely used with lifting or Chinese remaindering techniques as the bases to compute at a larger precision. There, the quotient of the Euclidean division is not needed, but is very often computed anyway along the algorithm.

In terms of arithmetic operations, from the work of [START_REF] Moenck | Fast modular transforms via division[END_REF][START_REF] Kung | On computing reciprocals of power series[END_REF] to the more recent results of [4], many sub-quadratic algorithms were developed for this task, based on fast polynomial multiplication [7,17,[START_REF] Von | Modern Computer Algebra[END_REF]. But these fast algorithms come at the expense of (potentially large) extra temporary space to perform the computation. On the contrary, classical, quadratic algorithms, when computed sequentially, quite often require very few (constant) extra registers [START_REF] Monagan | In-place arithmetic for polynomials over zn[END_REF]. Further work then proposed simultaneously "fast" and "in-place" algorithms, for both matrix and polynomial operations [5,[START_REF] Giorgi | Generic reductions for in-place polynomial multiplication[END_REF][START_REF] Giorgi | Fast in-place algorithms for polynomial operations: division, evaluation, interpolation[END_REF].

We here propose fast algorithms to extend the latter line of work. In the polynomial setting, we compute in-place the remainder only of the Euclidean division. This means that, e.g., with respect to [START_REF] Giorgi | Fast in-place algorithms for polynomial operations: division, evaluation, interpolation[END_REF]Alg. 3], we obtain the remainder without needing any space for the quotient. As polynomials and Toeplitz matrices are indeed different representations of the same objects, see, e.g., [START_REF] Bini | Fast parallel polynomial division via reduction to triangular Toeplitz matrix inversion and to polynomial inversion modulo a power[END_REF][START_REF] Bini | Polynomial and matrix computations, 1st Edition[END_REF][START_REF] Giorgi | Generic reductions for in-place polynomial multiplication[END_REF], we also develop fast methods for Toeplitz matrix operations, in-place with accumulation or over-place: those are difficult cases where the output is also part of the input. The difficulty when the result overwrites (parts) of the input, is that in-place methods start with absolutely no margin for extra space. Thus, for instance, the generic recursive techniques

In place and over-place quadratic algorithms

We recall here that classical, quadratic, algorithms for polynomial remaindering and triangular matrix operations can be performed in-place. For any field F we have for instance the following over-place algorithms for triangular matrix operations, given in Alg. 1.

Algorithm 1 Over-place quadratic triangular matrix operations (left: matrix-vector multiplication; right: triangular system solve)

Input: U ∈ F m×m upper triangular and v ∈ F m Read-only:

U Output: v ← U • v 1: for i = 1 to m do 2:
for j = 1 to i -1 do 3:

v j += U j,i v i 4:
end for 5:

v i ← U i,i v i ; 6: end for Output: v ← U -1 • v 1: for i = m down-to 1 do 2:
for j = m down-to i + 1 do 3:

v i -= U i,j v j 4:
end for 5:

v i ← U -1 i,i v i ; {If U i,i ∈ F * } 6: end for
To any vector v = [v 1 , . . . , v m ], we associate the polynomial of degree ≤ m -1, with these coefficients: V (X) = m-1 i=0 v i+1 X i . We then give a quadratic in-place remaindering in Alg. 2, as the classical long division. By abuse of notation, we will often use vectors for polynomials and vice versa. For n coefficients, we will then speak of a polynomial or a vector of size n and degree n -1.

q ← r M-1 • b -1 M ; {Local quotient is r M-1 • b -1 M } 5:
r ← [a i , r 1 , . . . , r M-2 ];

6:

r -= q • b; 7: end for 8: return R(X) = M-2 i=0 r i+1 X i . Algorithm 2 can be made over-place of its input a by considering that r is just a "pointer" to some position in a (with q in r M-1 ): this is then close to the in-place quadratic version given in [START_REF] Monagan | In-place arithmetic for polynomials over zn[END_REF].

Fast algorithms

If M(M ) is a sub-multiplicative complexity bound on an algorithm computing the multiplication of polynomials of degree M , we suppose in the following that there exists (not in-place) algorithms, see, e.g., [START_REF] Giorgi | Fast in-place algorithms for polynomial operations: division, evaluation, interpolation[END_REF]17,[START_REF] Von | Modern Computer Algebra[END_REF], such that M ×M Toeplitz matrix-vector multiplication as well as triangular Toeplitz system solve require less than O(M(M )) operations and O(M ) extra space. Further, for sparse inputs, we suppose that M ×M sparse Toeplitz (constant number of non-zero bands) matrix-vector multiplication as well as triangular sparse Toeplitz system solve requires less than O(M ) operations and O(M ) extra space.

In-place algorithms with accumulation

In an in-place setting, another difficulty arises when performing accumulating operations: the result is not even available to store intermediate computations. Recently, in [START_REF] Dumas | In-place accumulation of fast multiplication formulae[END_REF], a bilinear technique was developed for the fast accumulated in-place polynomial multiplication, preserving an O(M(M )) complexity bound. Thus we suppose in the following that there exist fast routines computing C(X) += A(X)B(X) in-place. In the following, we reduce other in-place accumulating routines (Toeplitz, circulant, convolutions operations) to this fundamental building block. This allows us, in fine, to obtain an in-place fast polynomial modular remainder.

Fast in-place convolution with accumulation

From algorithms for polynomial multiplications at a lower level, one can devise an algorithm for the generalized accumulated convolution C += AB mod (X nf ). To compute such convolutions, one can reduce the computations to full in-place products. Those accumulated in-place polynomial multiplications (and their cost), will also be denoted by M. Then, the initial idea is to unroll a first recursive iteration and then to call any in-place accumulated polynomial multiplication on halves.

Algorithm 3 shows this, first with a Karatsuba-like iteration, using any in-place accumulated polynomial multiplications (e.g. classical, Karatsuba, Toom-k, DFT, etc.). Apart from the full polynomial multiplications, only a linear number of operations will be needed. Therefore, the overall complexity bound will be that of the accumulated polynomial multiplications.

Suppose indeed first that n is even. Let t = n/2, and

A(X) = a 0 (X)+X t a 1 (X), B(X) = b 0 (X)+X t b 1 (X), C(X) = c 0 (X) + X t c 1 (X), all of degree n -1 = 2t -1 (so that all of a 0 , a 1 , b 0 , b 1 , c 0 and c 1 are of degree at most t -1). Then let τ 0 + X t τ 1 = a 0 b 1 + a 1 b 0 . Since X 2t = X n ≡ f mod X n -f , we have that: C + AB mod (X n -f ) = C + a 0 b 0 + X t τ 0 + f • τ 1 + f • a 1 b 1 .
This can be computed with 4 full accumulated sequential product, each of degree no more than 2t -2 = n -2, and by exchanging the lower and upper parts when accumulating a 0 b 1 X t and a 0 b 1 X t : this comes from the fact that (τ

0 +X t τ 1 )X t ≡ X t τ 0 + f • τ 1 mod X n -f .
À la Karatsuba, a more efficient version could use only 3 full polynomial products: let instead m 0 = m 00 + m 01 X t = a 0 b 0 , m 1 = m 10 + m 11 X t = (a 0 + a 1 )(b 0 + b 1 ) and m 2 = m 20 + m 21 X t = a 1 b 1 be these 3 full products, to be computed and accumulated in-place. As

X 2t = X n ≡ f mod X n , then C + AB mod (X n -f ) is : c 0 + m 00 + f m 20 + f (m 11 -m 01 -m 21 ) +(c 1 + m 01 + f m 21 + (m 10 -m 00 -m 20 ))X t (1)
Equation ( 1) is an in-place accumulating linear computation, in the sense of [START_REF] Dumas | In-place accumulation of fast multiplication formulae[END_REF]: [START_REF] Dumas | In-place accumulation of fast multiplication formulae[END_REF]Alg. 3] derives an accumulating in-place algorithm, using 2 × 2 invertible blocks, that computes

it computes c += µ m where m = (α a) ⊙ (β b), for α = β = 1 0 1 1 0 1 ∈ F 3×2 and µ = 1 -f -1 1 0 f 1 0 f -f -1 f ∈ F 2×6 . If f / ∈ {0, 1} and µ = [M 0 |M 1 |M 2 ], M -1 0 = (1 -f ) -1 1 f 1 1 , M -1 1 = 0 1 f -1 0 , and M -1 2 = (f 2 -f ) -1 f f 1 f . From this,
[ ci cj ] += M [ ρ0 ρ1 ], via the equivalent algorithm: [ ci cj ] * = M -1 ; [ ci cj ] += [ ρ0 ρ1 ]; [ ci cj ] * = M .
After some simplifications described below, we obtain Algorithm 3.

Algorithm 3 In-place even degree accumulating f -convolution

Input: A(X), B(X), C(X) polynomials of odd degree n -1; f ∈ F \ {0, 1}. Output: C += AB mod (X n -f ) 1: if n ≤ Threshold then {Constant-time if Threshold ∈ O(1)} 2:
return the quadratic in-place polynomial multiplication. Let t = n/2; {n is even}, 5:

Let A = a 0 + X t a 1 ; B = b 0 + X t b 1 ; C = c 0 + X t c 1 ; 6:
c 1 += c 0 ;

7:

c 1 /= (1 -f ); {Simplification (i)} 8: c 0 += f • c 1 ; 9:
c 1 c 0 -= a 1 • b 1 ; {m 2 via M with a 1 b 1 = m 20 m 21 } 12: c 0 -= c 1 ; {this is c 0 /f + m 00 /f -m 01 + m 20 -m 21 } 13: c 1 * = (1 -f ); {Simplification (iv)} 14: c 1 -= f • c 0 ; {this is c 1 -m 00 + m 01 -m 20 + f m 21 } 15:
a 0 += a 1 ;

16:

b 0 += b 1 ; 17: c 1 c 0 += a 0 • b 0 ; {m 1 via M with (a 0 + a 1 )(b 0 + b 1 ) = m 10 m 11 } 18: b 0 -= b 1 ; 19:
a 0 -= a 1 ; 20:

c 0 * = f ; {Simplification (v)} 21: end if
The algorithm is obtained from the output of [START_REF] Dumas | In-place accumulation of fast multiplication formulae[END_REF]Alg. 3] after applying the following simplifications:

(i) As 1 + f (1 -f ) -1 = (1 -f ) -1 , then M -1 0 [ c0 c1 ] = (1 -f ) -1 1 f 1 1 [ c0 c1 ] can be sequentially computed via c 1 += c 0 ; c 1 /= (1 -f ); c 0 += f c 1 . (ii) As M -1 2 • M 0 = 0 -1 -f -1 0 = -1 0 0 f -1 [ 0 1
1 0 ], then computing M 2 right after M 0 allows to simplify the intermediate in-place operations into a variable swap, negations and multiplication by f -1 .

(iii) That negation can be delayed, since (-M

• c) += ρ is equivalent to -((M • c) -= ρ). (iv) Similarly, -M -1 1 • M 2 = 1 -f -1 1
, thus computing M 1 after M 2 again simplifies to the sequential computation:

c 0 -= c 1 ; c 1 -= (1 -f ); c 1 -= f c 0 . (v) Finally, M 1 = f 0 0 1 [ 0 1
1 0 ] and thus this combination is again just a swap of variables and a multiplication by f . Overall, we obtain Alg. 3 that works for even n with f / ∈ {0, 1}. We then need to design variant algorithms for the other cases. Algorithm 4 deals with the case where f = 1 (as M 0 and M 2 are not invertible when f = 1). We present here only the version with 4 full products, but a more efficient version with 3 products only like Alg. 3, could also be derived.

Algorithm 4

In-place even degree accumulating 1-convolution Input: A(X), B(X), C(X) polynomials of odd degree n -1;

Output: C += AB mod (X n -1) 1: if n ≤ Threshold then {Constant-time if Threshold ∈ O(1)} 2:
return the quadratic in-place polynomial multiplication. Let t = n/2; {n is even, let Y = X t , so that Y 2 ≡ 1}, 5:

Let A = a 0 + X t a 1 ; B = b 0 + X t b 1 ; C = c 0 + c 1 X t ; 6: c 0 c 1 += a 0 • b 0 {Via M, of degree 2t -2 ≤ n -1} 7: c 0 c 1 += a 1 • b 1 {Via M, since Y 2 ≡ 1} 8: c 1 c 0 += a 0 • b 1 {Via M, since (u + vY )Y ≡ v + uY } 9: c 1 c 0 += a 1 • b 0 {Via M, since (u + vY )Y ≡ v + uY } 10: end if
Algorithm 5 deals with the odd-n case where f is invertible. For the sake of simplicity, we also only present the version with 4 full products. The additional difficulty here is that the degrees of the lower and upper parts are different. Therefore, Alg. 5 ensures that there is always the correct space to accumulate.

Finally, we also have to deal with the case f = 0, that is with the in-place short product C += AB mod X n . If the field cardinality is not 2, then a simple way is to compute two f -convolutions, with two different values for f , as shown in Remark 1.

Remark 1. If AB = QX n + R, of overall degree less than 2n -1, then AB = Q(X n -1) + (R + Q) = Q(X n -g) + (R + gQ)
, for any g, and thus the respective remainders, by X n -1, and respectively by X ng, share the same quotient Q. Now for any λ ∈ {0, 1}, taking g = λ(λ -1) -1 , we obtain that

R = λR + (1 -λ)R together with 0 = λ•1 + (1 -λ)g. Therefore C += AB mod X n can be computed by C += λAB mod (X n -1) followed by C += (1 -λ)AB mod (X n -g).

This is two calls among Algs. 3 to 5, with an overall complexity bound still of O(M(n)).

There remains to compute a short product in the field with 2 elements. For this, we split the inputs and output in 3 blocks of size close to n/3 to obtain the 3 × 14 matrix µ of Eq. ( 2): this is the 2 × 2 expansion of [8, Eq. ( 24)] modulo 2, recalled in Eq. (3), again letting

m i = m i0 + Y m i1 . 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 1 0 1 0 0 0 0 0 0 0 1 1 0 1 1 1 1 1 1 0 1 0 •   m00 m01 . . . m61   (2)
Algorithm 5 In-place odd degree accumulating f -convolution

Input: A(X), B(X), C(X) polynomials of even degree n -1; f ∈ F * . Output: C += AB mod (X n -f ) 1: if n ≤ Threshold then {Constant-time if Threshold ∈ O(1)} 2:
return the quadratic in-place polynomial multiplication.

3: end if 4: Let t = (n + 1)/2;
{n is odd, so that X 2t ≡ f X},

5: Let A = a 0 + X t a 1 ; B = b 0 + X t b 1 ; C = c 0 + c 1 X t ; {a 1 ,b 1 ,c 1 of degree n -1 -t = t -2} 6: c 0 c 1 += a 0 • b 0 {Via M, of degree 2t -2 ≤ n -1} 7: c 1..(t-1) c 1 += a 1 • b 1 {Via M, of degree (2t -4) + 1 ≤ n -1} 8: c 0..(t-2) /= f ; 9: c t..(2t-2) c 0..(t-2) += a 0 • b 1 {Via M, since 2t-3<n-1 (u+vX t-1 )X t ≡f v+uX t } 10: c t..(2t-2) c 0..(t-2) += a 1 • b 0 {Via M, since 2t-3<n-1 (u+vX t-1 )X t ≡f v+uX t } 11: c 0..(t-2) * = f ; A(Y ) = Y 2 a 2 + Y a 1 + a 0 ; B(Y ) = Y 2 b 2 + Y b 1 + b 0 ; C(Y ) = Y 2 c 2 + Y c 1 + c 0 ; m 0 = a 0 • b 0 ; m 1 = (a 0 + a 1 + a 2 ) • (b 0 + b 1 + b 2 ); m 2 = a 2 • b 2 ; m 3 = (a 0 + a 2 ) • (b 0 + b 2 ); m 4 = (a 1 + a 2 ) • (b 1 + b 2 ); m 5 = (a 1 + a 2 ) • (b 0 + b 1 ); m 6 = (a 0 + a 2 ) • (b 1 + b 2 ); t 0 = c 0 + m 0 ; t 1 = c 1 + m 1 + m 2 + m 3 + m 4 ; t 2 = m 1 + m 3 + m 4 + m 5 + m 6 ; then C + AB mod Y 3 ≡ Y 2 t 2 + Y t 1 + t 0 mod 2 (3)
The five first pairs of columns of Eq. ( 2) can be dealt with the technique of [START_REF] Dumas | In-place accumulation of fast multiplication formulae[END_REF]Alg. 3]. The last two pairs are not full rank and thus cannot be handled like this (in other words, there is no space to put the high order terms of these two products, even though they are needed for the final result). They, however, just represent (m 50 + m 60 ) mod X n/3 and thus can be performed via 2 recursive calls of degree n/3. There remains to handle the cases where n is not a multiple of 3. For this, let t = ⌊n/3⌋ and Y = X t and start by computing C += AB mod X 3t via the algorithm sketched above. Now the remaining one or two coefficients, of degree 3t and 3t + 1, when applicable, are just accumulated directly, scalar multiplications by scalar multiplications. The whole process is shown in Alg. 6.

Overall, Alg. 6 uses 5 polynomial multiplications with polynomials of degree n/3, then 2 recursive calls of short products of degree n/3 and O(n) extra computations. The overall complexity is thus still O(M(n)), as shown in Theorem 2. Finally, Algs. 3 to 6 all together provide a complete solution for the fast in-place convolution with accumulation, as given in Alg. 7.

Theorem 2. Using an in-place polynomial multiplication with complexity bounded by M(n), Alg. 7 is correct, in-place and has complexity bounded by O(M(n)).

Proof. Correctness of Alg. 3 comes from that of Eq. (1). Correctness of Algs. 4 and 5 is direct from the degrees of the sub-polynomials (as given in the comments, line by line). Correctness of Alg. 6 comes from that of Eq. [START_REF] Bini | Polynomial and matrix computations, 1st Edition[END_REF].

Also, all four algorithms are in-place as they use only in-place atomic operations or in-place polynomial multiplication.

Algorithm 6 In-place accumulating short product (0-convolution)

Input: A(X), B(X), C(X) polynomials of degree < n in F[X]; Output: C += AB mod X n . 1: if n ≤ Threshold then {Constant-time if Threshold ∈ O(1)} 2:
return the quadratic in-place polynomial multiplication.

3: end if 4: if F \ {0, 1} = ∅ then {Use Remark 1} 5:
Choose λ ∈ {0, 1} and let g = λ(λ -1) -1 ;

6:

C += λAB mod (X n -1); {Alg. 4 or Alg. 5}

7:

return C += (1 -λ)AB mod (X n -g). {Alg. 3 or Alg. 5} 8: end if {The rest is thus just when F ≃ F 2 } 9: Let t = ⌊n/3⌋; 10: Let A mod X 3t = a 0 +X t a 1 +X 2t a 2 ; B mod X 3t = b 0 +X t b 1 +X 2t b 2 ; C mod X 3t = c 0 +c 1 X t +X 2t c 2 ; 11: Let m i = [ mi0 mi1 ] = m i0 + m i1 X t , i = 0..6, as in Eq. (3); 12:   c 0 c 1 c 2   += 1 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 1 0 1 0 0 0 1 1 0 1 1 1 1 1 • [ m0 m1 m2 m3 m4 ] ⊺ {[11, Alg. 3]} 13: a 1 += a 2 ; b 0 += b 1 ; 14: c 2 += a 1 • b 0 mod X t {Recursive call: c 2 += m 5 mod X t } 15: b 0 -= b 1 ; a 1 -= a 2 ; 16: a 0 += a 2 ; b 1 += b 2 ; 17: c 2 += a 0 • b 1 mod X t {Recursive call: c 2 += m 6 mod X t } 18: b 1 -= b 2 ; a 0 -= a 2 ; 19: if n ≥ 3t + 1 then {Scalar accumulations} 20:
for i = 0 to 3t do c 3t += A i

for i = 0 to 3t + 1 do c 3t+1 += A i • B 3t+1-i end for 24: end if Algorithm 7 In-place convolution with accumulation Input: A(X), B(X), C(X) polynomials of degree < n; f ∈ F. Output: C += AB mod (X n -f ) 1: if f = 0 then return Alg. 6 end if 2: if n is odd then return Alg. 5 end if {f = 0} 3: if f = 1 then return Alg. 4 end if {n is even} 4: return Alg. 3. {f ∈ {0, 1}}
Then, Algs. 3 to 5 use only a linear number of operations plus up to four calls to polynomial multiplications of half-degrees. Their complexity bound thus stems from the fact that M(n)/n is non-decreasing and that the number of operations is less than O(n) + 4M(n/2).

Finally, Alg. 6 uses a linear number of atomic operations, five degree-n/3 polynomial multiplications and only two recursive calls with degree-n/3 polynomials. Thus its complexity bound satisfies T (n) ≤ 2T (n/3) + 5M(n/3) + O(n). Again, the non-decreasingness of M(n)/n gives that T (n) = O(M(n)).

Circulant and Toeplitz matrix operations with accumulation

Toeplitz matrix-vector multiplication can be reduced to circulant matrix-vector multiplication, via an embedding into a double-size circulant matrix. But this is not immediately in-place, since doubling the size requires a double space. We see in the following how we can instead double the operations while keeping the same dimension. We start by the usual definitions, also extending circulant matrices to f-circulant matrices, following, e.g., [21, Theorem 2.6.4].

Accumulating & in-place f -circulant

Definition 3. For a ∈ F m , C( a) is the circulant matrix represented by a, that is, the m × m matrix (C ij ), such that C 1j = a j and the (i + 1)-th row is the cyclic right shift by 1 of the i-th row. Definition 4. For f ∈ F and a ∈ F m , the (lower) f-circulant matrix represented by a, C f ( a), is the m × m matrix (Γ ij ), such that:

for C = C( a), Γ ij = C ij if i ≤ j, Γ ij = f • C ij otherwise.
The algebra of f -Circulant matrices is in fact isomorphic to the algebra of polynomials modulo X n -f [21, Theorem 2.6.1]. This means that the product of an f -circulant matrix by a vector is obtained by the convolution of this vector by the vector representing the f -circulant matrix, and thus via Alg. 7.

Further, it is also well known that circulant matrices are diagonalized by a discrete Fourier transform, and hence can be manipulated via fast Fourier transforms as C 1 ( a) = F -1 m diag(F m a)F m , for F m a DFT-matrix, see, e.g., [15, § 4.7.7]. This gives an alternative way to compute circulant matrix-vector multiplication inplace (using and restoring afterwards both the matrix and the vector). Indeed, the (even truncated) Fourier transform and its inverse can be computed in-place [START_REF] Daniel | Space-and time-efficient polynomial multiplication[END_REF][START_REF] Harvey | An in-place truncated Fourier transform and applications to polynomial multiplication[END_REF][START_REF] Coxon | An in-place truncated Fourier transform[END_REF]. This gives us an in-place algorithm to compute the accumulation c += C 1 ( a) • b as:

1. a ← F m a, b ← F m b and c ← F m c; 2. c += diag( a) • b; 3. c ← F -1 m c; b ← F -1 m b and a ← F -1 m a.
Now, this diagonalization allows to compute fast in-place accumulated circulant matrix-vector only when primitive roots exists, and only via this DFT algorithm. On the contrary, Alg. 7 has no restriction on f and is a reduction to any accumulated in-place polynomial multiplication (including DFT ones, as, e.g., [START_REF] Dumas | In-place accumulation of fast multiplication formulae[END_REF]Alg. 7]).

Accumulating & in-place Toeplitz

Definition 5. For a ∈ F 2m-1 , T( a) is the (square) Toeplitz matrix represented by a, that is, the m × m matrix (T ij ), such that T ij = a m+j-i . Similarly, for a ∈ F m+n-1 , we denote by T m,n ( a) the m×n rectangular Toeplitz matrix defined by its first column, a 1..m bottom to top, and its first row, a m..(m+n-1) , left to right.

Again, there is an isomorphism between rectangular Toeplitz matrices and the middle product of polynomials, see, e.g., [13, § 3.1]. We also immediately see that with these notations, we have for instance

C 1 ( a) = C( a), C 0 ( a) = 1 2 (C 1 ( a) + C -1 ( a)), or also C f ( a) = T([f • a 2..m , a])
, where [ u, v] denotes the vector obtained by concatenation of u and v.

Fast algorithms for f-circulant matrices then allows us to build algorithms, by reduction, for accumulation with triangular and square Toeplitz matrices first, as sums of f-circulant in Algs. 8 and 9, and then for any Toeplitz matrix, again as sums of triangular Toeplitz matrices in Alg. 10.

Algorithm 8

In-place accumulating Upp. Triang. Toeplitz m-v. mult.

Input: a, b, c ∈ F m . Output: c += T([ 0, a]) • b. 1: return c += C 0 ( a) • b. {Alg. 7}
Remark 6. We here present f-circulant matrices where the coefficient acts on the lower left part of the matrix (excluding the diagonal). In the same manner, one can design an in-place fast algorithm for the other type of f-circulant matrices, where the coefficient would act on the upper right part of the matrix (excluding the diagonal). One could thus derive the transposed version of Alg. 8, with a single call to (the transposed version of ) Alg. 7. Now it is also possible to use the same algorithm but on the transposed matrix and on reversed vectors: letting b be the reversed vector of b, one can see by inspection that:

C f ( a) ⊺ • b is the reverse of C f ( a) • b (4)
Algorithm 9 In-place accumulating square Toeplitz m-v. mult. Proof. The complexity bound comes from that of Alg. 7. Correctness is obtained directly looking at the values of the matrices. First, for Alg. 8, we have that T([ 0, a]) = C 0 ( a). Second, for Alg. 9, Remark 6 and Eq. ( 4) show that it is possible to compute

Input: a 1 ∈ F m , a 2 , b, c ∈ F m+1 , Output: c += T([ a 1 , a 2 ]) • b. 1: c += C 0 ( a 2 ) • b; {Alg. 7} 2
C 0 ( a 1 ) ⊺ • b 1 via the reverse of C 0 ( a 1 ) • b 1 . T([ a 1 , a 2 ]) = C 0 ( a 2 ) + 0 ⊺ 0 C0( a1) ⊺ 0 .
From this, we give in Alg. 10 an in-place rectangular Toeplitz matrix-vector multiplication.

Proposition 8. Algorithm 10 is correct and requires less than O max{m,n} min{m,n} M(min{m, n}) operations.

Proof. If m>n, then there are first ⌊m/n⌋ square n × n calls to Alg. 9. This is bounded by O((m/n)M(n)) ≤ O(M(m)) operations. The remaining recursive call is then negligible. This is similar when m<n.

These in-place accumulated Toeplitz matrix-vector multiplications, in turns, allow us to obtain both over-place triangular Toeplitz multiplication or system solve, given in Algs. 11 and 12.

Algorithm 10 In-place accumulating rectangular Toeplitz matrix-vector multiplication

Input: a ∈ F m+n-1 , b ∈ F n , c ∈ F m , Output: c += T m,n ( a) • b. 1: if m = n then return c += T( a) • b. end if {Algorithm 9} 2: if m > n
Input: a, b ∈ F m , s.t. a 1 ∈ F * . Output: b ← T([ a, 0]) • b. 1: if m ≤ Threshold then {Constant-time if Threshold ∈ O(1)} 2:
return the quadratic in-place triang. m-v. mult. Proof. First, let T = T([ 0, a]) and consider it as blocks

{Alg. 1} 3: end if 4: Let k = ⌈m/2⌉, b 1 = b 1..k and b 2 = b (k+1)..m ; 5: b 2 ← T([ a (k+1)..m , 0]) • b 2 ; {Recursive call} 6: b 2 += T m-k,k ([a 1 , . . . , a m-1 ]) • b 1 ; {Algorithm 7: b 1 ← T([a (m-k+1)..m , 0]) • b 1 ; {Recursive call}
T 1 = T([ 0, a 1..k ]), T 2 = T([ 0, a (k+1)..m ]) and G = T k,m-k ( a 2..m ). Then T = T1 G 0 T2 . Now define H, s.t. T -1 = T -1 1 H 0 T -1 2 . Then H satisfies T -1 1 G + HT 2 = 0. Also, we have T -1 b = [ T -1 1 b1+Hb2 T -1 2 b2 ] ⊺ . Let b2 = T -1 2 b 2 and b1 = T -1 1 b 1 + Hb 2 . Then b1 = T -1 1 b 1 + HT 2 b2 = T -1 1 b 1 -T -1 1 G b2 = T -1 1 b 1 -G b2
and this shows that the algorithm is correct. Now for the complexity bound, the cost function is We now have over-place Toeplitz methods. Next, we reduce the extra space for polynomial remaindering. Eventually, we combine these two techniques to obtain in-place polynomial remaindering.

T (m) ≤ 2T (m/2) + O(M(m)). Again, this is ≤ O(M(m) log(m)) if M(m)
Algorithm 12 Over-place triang. Toeplitz system solve

Input: a, b ∈ F m , s.t. a 1 ∈ F * . Output: b ← T([ 0, a]) -1 • b. 1: if m ≤ Threshold then {Constant-time if Threshold ∈ O(1)} 2:
return the quadratic in-place triang. syst. solve. {Alg. 1} 3: end if

4: Let k = ⌈m/2⌉, b 1 = b 1..k and b 2 = b (k+1)..m ; 5: b 2 ← T([ 0, a (k+1)..m ]) -1 • b 2 ; {Recursive call} 6: b 1 -= T k,m-k ( a 2..m ) • b 2 ; {Algorithm 10} 7: b 1 ← T([ 0, a 1..k ]) -1 • b 1 ; {Recursive call} 8: return b.

In-place modular remainder

We consider the fast in-place (resp. over-place) computation of the Euclidean polynomial modular remainder R(X) ≡ A(X) mod B(X) (resp. A(X) ≡ A(X) mod B(X)) with A and B of respective degrees n and m ≤ n. Standard algorithms for the remainder require O n m M(m) arithmetic operations and, apart from that of A and B, at least O(nm) extra memory [START_REF] Giorgi | Fast in-place algorithms for polynomial operations: division, evaluation, interpolation[END_REF]. This extra space is notably usually used to store the whole quotient Q(X)

such that A = BQ + R with deg R < deg B.
We first show how to avoid the storage of the whole of this quotient, and propose an algorithm still using O n m M(m) arithmetic operations but only O(m) extra space (when the divisor B is sparse with a constant number of non-zero terms, the arithmetic complexity bound reduces to O(n)).

Second, we combine this with the techniques of Sections 2 and 3 and use the input space of A or B for intermediate computations in order to propose in-place and over-place algorithms for the modular remainder using at most O(M(m) log(m)) arithmetic operations if M(m) is quasi-linear, and O(M(m)) otherwise.

Successively over-writing the quotients

With two polynomials A and B of respective degrees N and M , the computation of the Euclidean division remainder R of degree strictly less than M such that A = BQ + R with quotient Q, can be rewritten as R ≡ A -BQ mod X M . This is therefore enough to compute the quotient only up to the degree M -1: let

A M ≡ A mod X M and Q M ≡ Q mod X M , then R ≡ A M -BQ M mod X M .
This observation is the ingredient that allows to compute the remainder using an extra space only of the order of the degree of the divisor B. One can also see this as the long division algorithm applied to blocks of dimension M .

Let us write the Euclidean equation A = BQ + R in a Toeplitz matrix form. In Eqs. ( 5) and ( 6), we view the polynomials A, Q and R as vectors [a 0 , . . . , a N ], [q 0 , . . . , q N -M ], [r 0 , . . . , r M-1 , 0, . . . , 0] and then B as a Toeplitz matrix B = T([0, . . . , 0, b M , . . . , b 0 , 0, . . . , 0]). Then, focusing on the last N -M + 1 rows of Eq. ( 5) we obtain directly the upper triangular (N -M + 1) × (N -M + 1) Toeplitz system of equations whose solution is only the quotient, as shown in Eq. ( 6).

      a0 . . . . . . . . . aN       =       b0 0 . . . . . . bM . . . b0 . . . . . . 0 bM       •    q0 . . . . . . qN-M    +      r0 . . . rM-1 0      . ( 5
)    aM . . . . . . aN    =       bM ... b0 0 . . . . . . . . . . . . . . . b0 0 . . . . . . bM       •    q0 . . . . . . qN-M    . (6)
We now let n = N -M + 1 and suppose that B is really of degree M , that is, its leading coefficients b M is invertible in the coefficient domain. For the sake of simplicity we also assume that n is a multiple of M (otherwise, for now, just complete the polynomial A with virtual leading zero coefficients up to the next multiple of M ) and let µ = n M . We then denote the M ×M blocks of the Toeplitz matrix in Eq. ( 6) by:

T = T([ 0 M-1 , b M..1 ]) and G = T([b M-1..0 , 0 M-1 ]), that is: T = bM ... b1 . . . . . . bM and G = b0 . . . . . . bM-1 ... b0 . (7)
This in turns gives a way to access only the first coefficients of Q in an upper triangular Toeplitz system:

   q0 . . . . . . qM-1    = I M 0    T G 0 . . . . . . . . . G 0 T    -1    aM . . . . . . aN    (8)
with a 2-block band structure where I M 0 is the concatenation of the M ×M identity matrix and the (N -M )×(N -M ) zero matrix. Finally, recovering the remainder from the first M rows of Eq. ( 5) equations is just like multiplying the quotient by G and thus R = A -BQ mod X M can be written as:

R =    a0 . . . . . . aM-1    -G 0    T G 0 . . . . . . . . . G 0 T    -1    aM . . . . . . aN    (9)
In fact, as we need only the first M coefficient of the Toeplitz system we just need the first block-row of the (upper triangular) inverse of the upper triangular 2-band Toeplitz matrix. Now this first block-row, of the inverse, of an upper triangular Toeplitz matrix U , is given by a direct formula, obtained from either of the equations U • U -1 = I or U -1 • U = I (see, e.g., [9, Eq.( 1)] for the scalar case). If we denote by H i the blocks of that row, we have Eq. ( 10):

     H 1 = T -1 H i-1 G + H i T = 0, i = 2..µ T H i + GH i-1 = 0, i = 2..µ (10)
which solves as:

H i = T -1 (-GT -1 ) i = (-GT -1 ) i T -1 .
We have shown:

Lemma 11.    T G 0 . . . . . . . . . G 0 T    -1 = T -1 .     I -GT -1 ... (-GT -1 ) µ-1 . . . . . . . . . . . . -GT -1 0 I    
Now denote by [ a 0 , a 1 , . . . , a µ ] the decomposition into blocks of dimension M of [a 0 , . . . , a N ]. Combining Eq. ( 9) and Lemma 11, we obtain now that:

R =    a0 . . . . . . aM-1    -GT -1 µ i=1 -GT -1 i-1 a i = µ i=0 -GT -1 i a i (11)
From Eq. ( 11) we thus can immediately deduce the following Alg. 13 that uses only O(M ) extra memory space in a Horner-like fashion of the polynomial in (-GT -1 ) of Eq. [START_REF] Dumas | In-place accumulation of fast multiplication formulae[END_REF]. Note that this algorithm does not modify its input along its course: both A(X) and B(X) are for now read-only (in particular virtually padding a with virtual zeroes, Line 3, is therefore not an issue).

Algorithm 13 Overwritten-quotient Euclidean remainder Input: A(X), B(X) in F[X] of respective degrees N and M . Read-only: A(X), B(X). Output: R(X) ≡ A(X) mod B(X) of degree at most M -1.

1: if M > N then return A. end if 2: Let n = N -M + 1, µ = n M ; 3: Let [ a 0 , . . . , a µ ] = [a 0 , . . . , a N , 0]; {Blocks of dimension M } 4: Let T = T([ 0 M-1 , b M , . . . , b 1 ]), G = T([b M-1 , . . . , b 0 , 0 M-1 ]); 5: r = a µ ; { r in-place of the result} 6: for i = µ -1 down-to 0 do 7: t = T -1 • r;
{Triang. Toeplitz solve}

8: r = (-G) • t; {Triang. Toeplitz m-v. mult.} 9:
r += a i ; 10: end for

11: return R = M-1 i=0 r i X i .
Theorem 12. Algorithm 13 is correct and requires O N M M(M ) arithmetic operations and O(M ) extra memory space. If the polynomial B is sparse with a constant number of non-zero coefficients, the arithmetic complexity is reduced to O(N ).

Proof. Correctness is given by Lemma 11 and Eq. [START_REF] Dumas | In-place accumulation of fast multiplication formulae[END_REF]. For the complexity bounds, we use Section 1. • Line 8, r = (-G) • t, corresponds to recovering the lower part of the multiplication of that current quotient by B;

• Line 9, r += a i , updates the next M coefficients of the current dividend (the leading ones being zero by construction of the current quotient).

An in-place atomic (and thus quadratic) long division is given in [START_REF] Monagan | In-place arithmetic for polynomials over zn[END_REF]. In-place long division by blocks is also sketched for instance in the proof of [START_REF] Giorgi | Fast in-place algorithms for polynomial operations: division, evaluation, interpolation[END_REF]Lemma 2.1]. The latter Lemma gives about 3λM(N ) operations and (2 + s)M extra space. In fact a refined analysis should also give the same (better) complexity bound as that of Theorem 12, that is less than 2λ N M M(M ) operations and (1 + s)M extra space.

In-place remainders via Toeplitz techniques

We now derive algorithms that use only O(1) extra memory space in the in-place model of Section 1.1: modifying the inputs is possible if and only if all inputs are restored to their initial state after the completion of the algorithm. This allows us to store some intermediate results, over-writing the input, provided that we can afterwards recompute the initial inputs in their entirety. Further, this enables recursive calls, as intermediate values are used but restored along the recursive descent. The general idea is then to combine Sections 2, 3 and 4.1. in fact, in the following, we present an in-place polynomial remaindering, IPER (and two variants OPER, APER in Appendix A):

• IPER: for R(X) ≡ A(X) mod B(X), creates a small remainder from a large polynomial using only the output space and that of the modulus B (i.e. A is read-only and B is restored to its initial state after completion);

• OPER: for A(X) ← A(X)/B(X); A(X) mod B(X) , updates a large polynomial in-place by its quotient and small remainder (with B restored after completion);

• APER: for R(X) += A(X) mod B(X), accumulates the remainder of a large polynomial into a small one (with both A and B restored after completion).

We can now present in Alg. 14, a fully in-place remainder where only B(X) is modified but restored: this variant replaces only Lines 7 to 9 of Alg. 13 by their over-place variants, Algs. 11 and 12, that modify and restore some parts of B.

Algorithm 14 IPER(R, A, B): In-place Polynomial Euclidean Remainder Input: A(X), B(X) in F[X] of respective degrees N and M . Read-only: A(X).

Output: R(X) ≡ A(X) mod B(X) of degree at most M -1.

1: if M > N then return A. r += a i ; 8: end for 9: return R = M-1 i=0 r i X i .

From Alg. 14 one can also obtain, Alg. 15 (given in Appendix A) which directly updates the dividend over-place, while also remaining fully in-place. In Alg. 14, A is read-only and can thus be virtually padded with zeroes, when N -M + 1 is not a multiple of M . This is not the case anymore in Alg. 15. This is why this variant takes some special care of the last s = (N + 1) mod M coefficients.

Finally, one can also obtain the variant Alg. [START_REF] Harvey | An in-place truncated Fourier transform and applications to polynomial multiplication[END_REF], where the remainder is accumulated in-place (also given in Appendix A). As in Alg. 15 the dividend is updated over-place, but in fact enough information is preserved so as to put it back into its initial state after the accumulation of the actual remainder. This suffices to prove the Theorem's claims for Alg. 14. These costs also dominates the costs of Algs. 9 and 10 by Lemma 7 and Proposition 8, and therefore that of Alg. 15: its Lines 5 to 8 are the over-place first iteration of Alg. 14, for s < M . Finally, Alg. 16 is also similar, using the same sub-algorithms, and only roughly doubling the overall cost (by symmetry before and after the accumulation).

Conclusion

We have presented novel reductions of accumulated in-place f-circulant and Toeplitz matrix-vector multiplications to in-place polynomial multiplication. This allows us to derive novel algorithms for accumulated or over-place Toeplitz multiplication or Toeplitz system solving. We also present algorithms that reduce the extra storage required to compute the remainder only when dividing polynomials. Eventually, we combine these techniques to propose the first in-place, over-place and accumulating algorithms computing only the remainder of the polynomial Euclidean division. Appendix B also shows a direct application of these techniques for the multiplication in a polynomial extension.

A Over-place & accumulated remainders

We here give the variants of Alg. 14 where, in Alg. 15, the polynomial is overwritten by its quotient and remainder, see Remark 13; and where, in Alg. 16, the remainder is accumulated.

Algorithm 15 OPER(A, B): Over-place Polynomial Euclidean Quotient and Remainder Input: A(X), B(X) in F[X] of respective degrees N and M . Output: A is replaced by both A(X)/B(X) and A(X) mod B(X) of degrees N -M and at most M -1.

1: if M > N then return A.

a µ ← T -1 1 • a µ ; {Algorithm 7: Let G 1 = T M,s ([b M-1 , . . . , b 0 , 0 s-1 ]); {Left s columns of G} 8: a µ-1 -= G 1 • a µ ; {Algorithm 9: end if 10: for i = µ-1 down-to 1 do 11: a i ← T -1 • a i ; {Algorithm 12:
a i-1 -= G • a i ; {Algorithm 13: end for 14: return The Quotient [ a 1 , . . . , a µ ] and the Remainder a 0 .

A nice property of OPER is that it is reversible: obviously, by computing in the reverse order a j-1 += G• a j ; a j ← T • a j one reverses the algorithm and recovers A. We denote by OPER -1 this recovery.

B In-place modular multiplication

As a direct application of this possibility to compute modular remainder in-place, we can obtain for instance that modular multiplications, of the form R(X) += A(X)C(X) mod B(X), can be computed in-place. Therefore one can use Alg. 16, but computing some of the coefficients of AC only when they are needed. These latter computations can be performed in-place of C, provided that they can be reverted in the end. This is actually always possible if A is smaller than both B and C, as explained, thereafter, in Alg. 17 and Proposition 15. Otherwise, if C is smaller than both A and B, then of course just compute R += CA mod B. Finally in the remaining cases, we show in Alg. [START_REF] Kung | On computing reciprocals of power series[END_REF] and Proposition 16 that one can always start by reducing the smaller of A or C by B, in-place using Alg. 15, as the latter reduction is revertible.

With the notations of Alg. 17, Eq. ( 12) becomes: A0 Aµ • c = [ G T ] • q + r which gives both A µ • c 2 = T • q and A 0 • c = G• q+ r. Then, as M ≥ L implies N ≥ q, both A µ and T are invertible upper triangular matrices (with Algorithm 16 APER(R, A, B): Accumulated in-place Polynomial Euclidean Remainder Input: R(X), A(X), B(X) in F[X] of resp. degrees M -1, N and M . Output: R(X) += A(X) mod B(X) of degree at most M -1.

1: if M > N then return R += A. the leading polynomial coefficient in the diagonal). The residue to be added is thus

A 0 • c -G • T -1 • A µ • c 2 ,
as in Alg. 17. Now for the complexity, A µ , G and T are matrices with dimensions smaller that M ×M . Toeplitz operations with those thus require less than O(M(M ) log(M )) operations if M(M ) is quasi-linear, and O(M(M )) otherwise by Propositions 8 to 10. Finally, A 0 is an M ×N Toeplitz so that the cost of applying it to a vector is bounded by O( N M M(M )), by Proposition 8. With Alg. 17, we can obtain for instance that multiplications in a polynomial extension of a finite field, can be computed in-place. Indeed, this is a special case where L ≤ N < M , as operations are performed modulo the irreducible polynomial of the extension and elements of the extension field are polynomials of degree upper bounded by that of this irreducible polynomial.

Finally, Algorithm 18 shows that one can also compute R += AC mod B, without any constraint on the degrees. For this we combine Algs. 15 and 17: if A is larger than C, compute R += CA mod B instead; if A is larger than B, start by reducing it in-place via OPER (Alg. 15), then compute the multiplication by Alg. 17 and eventually restore A via OPER -1 . Proof. This is the combination of Theorem 14 and Proposition 15.

Algorithm 2 1 i=0 a i+1 X i end if 2 :

 212 In-place quadratic polynomial remainderInput: a ∈ F N , b ∈ F M with b M ∈ F * Read-only: a, b Output: R(X) = A(X) mod B(X) 1: if N < M then return R(X) = A(X) = N -Let b = [b 1 , . . . , b M-1 ], n = N -M + 1 andr ← [a n+1 , . . . , a N ]; 3: for i = n down-to 1 do 4:

c 0 c 1 = m 00 m 01 } 10 :

 110 += a 0 • b 0 ; {m 0 via M with a 0 b 0 c 0 /= f ; {Simplifications (ii) and (iii)} 11:

: 3 : 5 :Lemma 7 .

 357 Let b 1 = b 1..m and c 2 = c 2..m+1 ; In-place reverse a 1 (resp. b 1 , c 2 ) into a 1 (resp. b 1 , c 2 ) 4: c 2 += C 0 ( a 1 ) • b 1 ; {Eq. (4) and Alg. 7} In-place reverse a 1 (resp. b 1 , c 2 ) into a 1 (resp. b 1 , c 2 ) 6: return c. Algorithms 8 and 9 are correct and have complexity bounded by O(M(n)).

  then

3 : 4 :c 1 5 :c 2 else 7 : 9 :c

 3415279 Let c 1 = c 1..n and c 2 = c (n+1)..m ; += T( a (m-n+1)..(m+n-1) ) • b; {Algorithm 9} += T m-n,n ( a 1..(m-1) ) • b; {Recursive call} 6: Let b 1 = b 1..m and b 2 = b (m+1)..n ; 8: c += T( a 1..(2m-1) ) • b 1 ; {Algorithm 9} += T m,n-m ( a (m+1)..(m+n-1) ) • b 2 ; {Recursive call} 10: end if 11: return c. Algorithm 11 Over-place triang. Toeplitz m-v. mult.

8: return b. Proposition 9 .

 9 Algorithm 11 is correct and requires less than O(M(m) log(m)) operations if M(m) is quasi-linear, and O(M(m)) otherwise. Proof. For the correctness, let T = T([ a, 0]) and consider it as blocks T 1 = T([a (m-k+1)..m , 0]), T 2 = T([ a (k+1)..m , 0]) and G = T m-k,k ([a 1 , . . . , a m-1 ]). Then T = T1 0 G T2 . Thus T b = T1b1 Gb1+T2b2 . Let b1 = T 1 b 1 , b2 = T 2 b 2 and b2 = Gb 1 + T 2 b 2 . Then b2 = b2 + Gb 1 and the algorithm is correct. Now for the complexity bound, the cost function is T (m) ≤ 2T (m/2) + O(M(m)). If M(m) is quasi-linear this is ≤ O(M(m) log(m)), and ≤ O(M(m)) otherwise. Proposition 10. Algorithm 12 is correct and requires less than O(M(m) log(m)) operations if M(m) is quasi-linear, and O(M(m)) otherwise.

  is quasi-linear, and O(M(m)) otherwise.

  3: for each block, the triangular Toeplitz system solve and the Toeplitz m-v. mult. require respectively λ σ M(m) and λ τ M(m) operations and, sequentially, max{σ; τ }M extra space. Apart from this space, we only need one extra vector, t, to store intermediate results. Overall we thus perform µ((λ σ + λ τ )M(M ) + M ) operations. With µ = n/M and n = N -M + 1, this is N -M+1 M ((λ σ + λ τ )M(M ) + M ) = O( N M M(M )) operations, using (1 + max{σ; τ })M extra space. Now if B is sparse with a constant number of non-zero elements, each triangular Toeplitz system solve and Toeplitz matrix-vector multiplication can be performed with only O(M ) operations with the same extra memory space. Thus the overall arithmetic bound becomes O(µM ) = O(N ). Remark 13. Algorithm 13 is in fact just the long division polynomial algorithm applied to sub-blocks of the polynomial of size M : • Line 7, t = T -1 • r, corresponds to computing the quotient of the current leading coefficients, of the dividend, by B;

2 :

 2 Let n, µ, [ a 0 , . . . , a µ ], T , G as in Alg. 13; 3: r = a µ ; { r in-place of the result} 4: for i = µ -1 down-to 0 do 5: r ← T -1 • r; {Algorithm 6: r ← (-G) • r; {Algorithm 7:

Theorem 14 .

 14 Algorithms 14, 15 and 16 are correct, in-place and require less than O N M M(M ) log(M ) operations if M(M ) is quasi-linear, and O N M M(M ) otherwise. Proof. Algorithm 14 calls µ= N -M+1 M = O N M times Algs. 11 and 12, each call requiring less than O(M(M ) log(M )) operations by Propositions 9 and 10 if M(M ) is quasi-linear, and O(M(M )) otherwise.

2 : 3 : 5 :

 235 Let n, µ, T , G as in Alg. 14 and s = (N +1) mod M ; Let [ a 0 , . . . , a µ-1 ] = [a 0 , . . . , a N -s ] and a µ = a (N -s+1)..N ; 4: if s = 0 thenLet T 1 = T([ 0 s-1 , b M , . . . , b M-s+1 ]);{s×s upper left of T } 6:

Proposition 15 .

 15 Algorithm 17 is correct and requires less thanO log(M ) + N M M(M ) operations if M(M ) is quasi-linear,and O NM M(M ) otherwise. Proof. Replacing the left-hand side of Eq. (5) by AC we see that AC = BQ + R can be written as Eq. (12):

Proposition 16 .

 16 Algorithm 18 is correct and requires less than O log(M ) + L+N M M(M ) operations if M(M ) is quasi-linear, and O L+N M M(M ) otherwise.

2 :

 2 Let n, µ, s, T , T 1 , G, G 1 , [ a 0 , . . . , a µ ] as in Alg. 15, and r = [r 0 , . . . , r M-1 ];3: if s = 0 then -= G 1 • a µ ; {Alg. 6: end if 7: for i = µ -1 down-to 1 do -= G • a i ;{Alg. 9} 10: end for 11: r += a 0 ; 12: for i = 1 to µ -1 do {Undo Lines 8 and 9 via Algs. 9 and 11}13:a i-1 += G • a i ; a i ← T •a i ; 14: end for 15: if s = 0 then {Undo Lines 4 and 5 via Algs. 10 and 11} 16: a µ-1 += G 1 • a µ ; a µ ← T 1 • a µ ; 17: end if 18: return r.

	4:	a µ ← T -1 1 • a µ ;	{Alg. 12}
	5: a µ-1 8: a i ← T -1 • a i ;	{Alg. 12}
	9:	a i-1