Jean-Guillaume Dumas

Bruno Grenet

In-place accumulation of fast multiplication formulae

This paper deals with simultaneously fast and in-place algorithms for formulae where the result has to be linearly accumulated: some of the output variables are also input variables, linked by a linear dependency. Fundamental examples include the in-place accumulated multiplication of polynomials or matrices, C += AB. The difficulty is to combine in-place computations with fast algorithms: those usually come at the expense of (potentially large) extra temporary space, but with accumulation the output variables are not even available to store intermediate values. We first propose a novel automatic design of fast and in-place accumulating algorithms for any bilinear formulae (and thus for polynomial and matrix multiplication) and then extend it to any linear accumulation of a collection of functions. For this, we relax the in-place model to any algorithm allowed to modify its inputs, provided that those are restored to their initial state afterwards. This allows us, in fine, to derive unprecedented in-place accumulating algorithms for fast polynomial multiplications and for Strassen-like matrix multiplications.

Introduction

Multiplication is one of the most fundamental arithmetic operations in computer science and in particular in computer algebra and symbolic computation. In terms of arithmetic operations, for instance, 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][START_REF] Volker Strassen | Gaussian elimination is not optimal[END_REF], many sub-quadratic (resp. sub-cubic) algorithms were developed for polynomial (resp. matrix) multiplication. But these fast algorithms usually come at the expense of (potentially large) extra temporary space to perform the computation. On the contrary, classical, quadratic (resp. cubic) algorithms, when computed sequentially, quite often require very few (constant) extra registers. Further work then proposed simultaneously "fast" and "in-place" algorithms, for matrix or polynomial operations [START_REF] Boyer | Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm[END_REF][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] 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 extend the latter line of work for accumulating algorithms. Actually, one of the main ingredient of the latter line of work is to use the (free) space of the output as intermediate storage. But when the result has to be accumulated, i.e., if the output is also part of the input, this free space does not even exist. To be able to design accumulating in-place algorithms we thus relax the in-place model to allow algorithms to also modify their input, therefore to use them as intermediate storage, provided that they are restored to their initial state after completion of the procedure. This is in fact a natural possibility in many programming environments. Furthermore, this restoration allows for recursive combinations of such procedures, as the (non concurrent) recursive calls will not mess-up the state of their callers. We thus propose a generic technique transforming any bilinear algorithm into an in-place algorithm under this model. This then directly applies to accumulating polynomial and matrix multiplication algorithms, including fast ones. Further, the technique actually generalizes to any linear accumulation, i.e. not only bilinear formulae, provided that the input of the accumulation can be itself reversibly computed in-place (therefore also potentially in-place of some of its own input if needed).

Next, we give our model for in-place computations and recall classical in-place algorithms in Section 2. We then detail in Section 3 our novel technique for in-place accumulation. Finally we apply this technique and further optimizations in order to derive new fast and in-place algorithms for the accumulating multiplication of matrices, Section 4, and of polynomials, Section 5.

In-place model with mutable input

We here call in-place, an algorithm using only the space of its inputs, its outputs, and at most O(1) extra space; the algorithm can also modify its (mutable) inputs, provided that these inputs are restored to their initial state afterwards. This is a less powerful model than when the input is purely readonly, but it turns out to be crucial in our case, especially when we have accumulation operations. Indeed, the algorithms we describe are in-place with accumulation, and the archetypical example is a multiply-accumulate operation c += a × b. For such an algorithm, the condition is that a and b are restored to their initial states at the end of the computation, while c (which is also part of the input) is replaced by c + ab. Also, as a variant, by over-place, we mean an algorithm where the output replaces (parts of) its input (e.g., like a ← b• a). In the following we signal by a "Read-only:" tag the parts of the input that the algorithm is not allowed to modify (the other parts are mutable as long as they are restored).

For recursive algorithms, some space may be required to store the recursive call stack (this space is bounded by the recursion depth of the algorithms, and managed in practice by the compiler). In our complexity summaries (Tables 2 and3 in Appendix C), the space complexity includes the stack. Nonetheless, we call in-place a recursive algorithm whose only extra space is the call stack. For more details on these models, we refer to [START_REF] Buhrman | Computing with a full memory: catalytic space[END_REF][START_REF] Daniel | Space-and time-efficient polynomial multiplication[END_REF][START_REF] Giorgi | Generic reductions for in-place polynomial multiplication[END_REF]. We end this section by recalling in Alg. 1, that classical algorithms for matrix or polynomial operations can be performed strictly in-place.

Algorithm 1 Quadratic/cubic in-place accumulating polynomial/matrix multiplication Input: A, B, C polynomials of degrees m, n, m+n. Read-only: A, B. Output: C(X) += A(X)B(X)

1: for i, j do 2:

C[i+j] += A[i]B[j];

Input: A, B, C matrices of dimensions m×ℓ, ℓ×n, m×n. Read-only: A, B. Output: C += AB 1: for i, j, k do

C ij += A ik B kj ;

3 In-place linear accumulation Karatsuba polynomial multiplication [START_REF] Karatsuba | Multiplication of multidigit numbers on automata[END_REF] and Strassen matrix multiplication [START_REF] Volker Strassen | Gaussian elimination is not optimal[END_REF] are famous optimizations of bilinear formulae on their inputs: results are linear combinations of products of bilinear combinations of the inputs. To compute recursively such a formula in-place, we perform each product one at a time. For each product, both factors are then linearly combined in-place into one of the entry beforehand and restored afterwards. The product of both entries is at that point accumulated in one part of the output and then distributed to the other parts. The difficulty is to perform this distribution in-place, without recomputing the product. Our idea is to presubtract one output from the other, then accumulate the product to one output, and finally re-add the newly accumulated output to the other one: overall both outputs just have accumulated the product, in-place. Potential constant factors can also be dealt with pre-divisions and postmultiplications. Basically we need two kind of in-place operations, and their combinations. First, as shown in Eq. ( 1), an in-place accumulation of a quantity multiplied by a (known in advance) invertible constant:

{c /= µ; c += m; c * = µ; } computes in-place c ← c + µ • m. (1)
Second, as shown in Eq. ( 2), an in-place distribution of the same quantity, without recomputation, to several outputs:

{d -= c; c += m; d += c; } computes in-place c ← c + m; d ← d + m. (2)
Example 1 shows how to combine several of these operations, while also linearly combining parts of the input.

Example 1. Suppose, that for some inputs/outputs a, b, c, d, r, s one wants to compute an intermediate product p = (a + 3b) * (c + d) only once and then distribute and accumulate that product to two of its outputs (or results), such that we have both r ← r + 5p and s ← s + 2p. To perform this in-place, first accumulate a += 3b and c += d, then pre-divide r by 5, as in Eq. (1). Now we directly have p = ac and it can be computed once, and then accumulated to r, and to s, if the latter is prepared: divide it by 2, and pre-subtract r. This is s /= 2, s -= r followed by r += ac. After this, we can reciprocate (or unroll) the precomputations: this distributes the product to the other result and restores the read-only inputs to their initial state (for s, another possibility is to directly pre-subtract 2r and to post-add 2r). This is summarized as:

a += 3b; c += d; r /= 5 ; s -= 2r; r += ac; s += 2r; r * = 5 ; c -= d; a -= 3b; computes in-place r ← r + 5(a + 3b)(c + d); s ← s + 2(a + 3b)(c + d).
Algorithm 2 shows how to implement this in general, taking into account the constant (or read-only) multiplicative coefficients of all the linear combinations. We suppose that inputs are in three distinct sets: left-hand sides, a, right-hand sides, b, and those accumulated to the results, c. We denote by ⊙ the point-wise multiplications of left-hand sides by right-hand sides. Then Alg. 2 computes c += µ m, for m = (α a) ⊙ (β b), with α, β and µ matrices of constants.

Algorithm 2 In place bilinear formula

Input: a ∈ F m , b ∈ F n , c ∈ F s ; α ∈ F t×m , β ∈ F t×n , µ ∈ F s×t , without zero-rows in α, β, µ. Read-only: α, β, µ. Output: c += µ m, for m = (α a) ⊙ (β b). 1: for ℓ = 1 to t do 2:
Let i s.t. α ℓ,i = 0; ai * = α ℓ,i ;

3:

for λ = 1 to m, λ = i, α ℓ,λ = 0 do ai += α ℓ,λ a λ end for 4:
Let j s.t. β ℓ,j = 0; bj * = β ℓ,j ;

5:

for λ = 1 to n, λ = j, β ℓ,λ = 0 do bj += β ℓ,λ b λ end for 6: Let k s.t. µ k,ℓ = 0; c k /= µ k,ℓ ; 7:
for λ = 1 to s, λ = k, µ λ,ℓ = 0 do c λ -= µ λ,ℓ c k end for for λ = 1 to s, λ = k, µ λ,ℓ = 0 do c λ += µ λ,ℓ c k end for {undo 7}

10:

c k * = µ k,ℓ ; {undo 6}
11:

for λ = 1 to n, λ = j, β ℓ,λ = 0 do bj -= β ℓ,λ b λ end for {undo 5}

12:

bj /= β ℓ,j ; {undo 4}

13:

for λ = 1 to m, λ = i, α ℓ,λ = 0 do ai -= α ℓ,λ a λ end for {undo 3}

14:

ai /= α ℓ,i ; {undo 2}

15: return c.

Remark 2. Lines 2 to 7 and 9 to 14 of Alg. 2 are acting on independent parts of the input, a and b, and of the output c. If needed they could therefore be computed in parallel or in different orders, and even potentially grouped or factorized across the main loop (on ℓ).

To simplify the counting of operations, we denote by ADD both the addition or subtraction of elements, += or -= ; by MUL the (tensor) product of elements, ⊙; and by SCA the scaling by constants, * = or /= . We also denote by #x (resp. ♯x) the number of non-zero (resp. ∈ {0, 1, -1}) elements in a matrix x. Theorem 3. Algorithm 2 is correct, in-place, and requires t MUL, 2(#α + #β + #µ) -5t ADD and 2(♯α + ♯β + ♯µ) SCA operations.

Proof. First, as the only used operations ( += , -= , * = , /= ) are in-place ones, the algorithm is inplace. Second, the algorithm is correct both for the input and the output: the input is well restored, as (α ℓ,i a i + α ℓ,λ a λ -α ℓ,λ a λ )/α ℓ,i = a i and (β ℓ,j b

j + β ℓ,λ b λ -β ℓ,λ b λ )/β ℓ,j = b j ; the output is correct as c λ -µ λ,ℓ c k /µ k,ℓ + µ λ,ℓ (c k /µ k,ℓ + a i b j ) = c λ + µ λ,ℓ a i b j and (c k /µ k,ℓ + a i b j )µ k,ℓ = c k + µ k,ℓ a i b j . Third,
for the number of operations, Lines 2 and 3 require one multiplication by a constant for each non-zero element a λ in the row and one less addition. But multiplications and divisions by 1 are no-op, and by -1 can be dealt with subtraction. This is #α -t additions and ♯α constant multiplications. Lines 4 and 5 (resp. Lines 6 and 7) are similar for each non-zero element in b λ (resp. in µ). Finally Line 8 performs t multiplications of elements and t additions. The remaining lines double the number of ADD and SCA.

This is t + 2(#α + #β + #µ -3t) = 2(#α + #β + #µ) -5t ADD.
Remark 4. Similarly, slightly more generic accumulation operations of the form c ← γ ⊙ c + µ m, for a vector γ ∈ F s , can also be computed in-place: precompute first c ← γ ⊙ c, then call Alg. 2.

For instance, to use Alg. 2 with matrices or polynomials, each product m ℓ is in fact computed recursively. Further, in an actual implementation of a fixed formula, one can combine more efficiently the pre-and postcomputations over the main loop on ℓ, as in Remark 2. See Sections 4 and 5 for examples of recursive calls, together with sequential optimizations and combinations.

In fact the method for accumulation, computing each bilinear multiplication once is generalizable. With the notations of Alg. 2, any algorithm of the form c += µ m can benefit from this technique, provided that each m j can be obtained from a function that can be computed in-place. Let F j : Ω → F be such a function on some inputs from a space Ω, for which an in-place algorithm exists. Then we can accumulate it in-place, if it satisfies the following constraint, that it is not using its output space as an available intermediary memory location. Further, this function can be in-place in different models: it can follow our model of Section 2, if there is a way to put its input back into their initial states, or some other model, again provided that it follows the above constraint.

Then, the idea is just to keep from Alg. 2 the Lines 6 to 10, replacing Line 8 by the in-place call to F j , potentially surrounding that call by manipulations on the inputs of F j (just like the one performed on a and b in Alg. 2). We have shown: Theorem 5. Let c ∈ F s and µ ∈ F s×t , without zero-rows. Let F = (F j : Ω → F) j=1..t be a collection of functions and ω ∈ Ω. If all these functions are computable in-place, without using their output space as an intermediary memory location, then there exists an in-place algorithm computing c += µ F (ω) in-place, requiring a single call to each F j , together with (2#µ -t) ADD and 2♯µ SCA ops. This is close to catalytic machines and transparent space [START_REF] Buhrman | Computing with a full memory: catalytic space[END_REF]; but here we allow only the input and output as catalytic space, and we do preserve the (not in-place) time complexity bound, up to a (quasi)linear overhead.

We give examples of the application of the generalized method of Theorem 5 to non bilinear formulae in Appendix B.

In-place Strassen matrix multiplication with accumulation

Considered as 2 × 2 matrices, the matrix product with accumulation C += A • B could be computed using Strassen-Winograd (S.-W.) algorithm by performing the computations of the following Eq. (3).

ρ 1 ← a 11 b 11 , ρ 3 ← (-a 11 -a 12 + a 21 + a 22 )b 22 , ρ 2 ← a 12 b 21 , ρ 4 ← a 22 (-b 11 + b 12 + b 21 -b 22 ), ρ 5 ← (a 21 + a 22 )(-b 11 + b 12 ), ρ 6 ← (-a 11 + a 21 )(b 12 -b 22 ), ρ 7 ← (-a 11 + a 21 + a 22 )(-b 11 + b 12 -b 22 ) [ c11 c12 c21 c22 ] += ρ1+ρ2 ρ1-ρ3+ρ5-ρ7 ρ1+ρ4+ρ6-ρ7 ρ1+ρ5+ρ6-ρ7 (3)
This algorithm uses 7 multiplications of half-size matrices and 24 + 4 additions (that can be factored into only 15 + 4, see [START_REF] Winograd | La complexité des calculs numériques[END_REF]: 4 involving A, 4 involving B and 7 involving the products, plus 4 for the accumulation). This can be used recursively on matrix blocks, halved at each iteration, to obtain a sub-cubic algorithm. To save on operations, it is of course interesting to compute the products only once, that is store them in extra memory chunks. To date, up to our knowledge, the best versions that reduced this extra memory space (also overwriting the input matrices but not putting them back in place) were proposed in [START_REF] Boyer | Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm[END_REF]: their best sub-cubic accumulating product used 2 temporary blocks per recursive level, thus a total of extra memory required to be 2 3 n 2 . With Alg. 2 we instead obtain a sub-cubic algorithm for accumulating matrix multiplication with O(1) extra space requirement. From Eq. ( 3) indeed (see also the representation in [START_REF] Nader | On the additive complexity of 2 × 2 matrix multiplication[END_REF]), we can extract the matrices µ, α and β to be used in Alg. 2 as given in Eq. ( 4).

µ = 1 1 0 0 0 0 0 1 0 -1 0 1 0 -1 1 0 0 1 0 1 -1 1 0 0 0 1 1 -1 α =    1 0 0 0 0 1 0 0 -1 -1 1 1 0 0 0 1 0 0 1 1 -1 0 1 0 -1 0 1 1    β =    1 0 0 0 0 0 1 0 0 0 0 1 -1 1 1 -1 -1 1 0 0 0 1 0 -1 -1 1 0 -1    (4)
All coefficients being 1 or -1 the resulting in-place algorithm can of course compute the accumulation C += AB without constant multiplications. It thus requires 7 recursive calls and, from Theorem 3, at most 2(#α + #β + #µ -3t) = 2(14 + 14 + 14 -3 * 7) = 42 block additions. Just like the 24 additions of Eq. ( 3) can be factored into 15, one can optimize also the in-place algorithm. For instance, looking at α we see that performing the products in the order ρ 6 , ρ 7 , ρ 3 , ρ 5 and accumulating in a 21 allows to perform all additions/subtractions in A with only 6 operations (this is in fact optimal, see the following Prop. 8). This is similar for β if the order ρ 6 , ρ 7 , ρ 4 , ρ 5 is used and accumulation is in b 12 . Thus ordering for instance ρ 6 , ρ 7 , ρ 4 , ρ 3 , ρ 5 will reduce the number of block additions to 26. Now looking at µ (more precisely at its transpose, see [START_REF] Kaminski | Addition requirements for matrix and transposed matrix products[END_REF]), a similar reduction can be obtained, e.g., if one of the orders (ρ 6 , ρ 7 , ρ 1 , ρ 5 ) or (ρ 5 , ρ 7 , ρ 1 , ρ 6 ) is used and accumulation is in c 22 .

Therefore, using the ordering ρ 6 , ρ 7 , ρ 1 , ρ 4 , ρ 3 , ρ 5 , ρ 2 requires only 18 additions, as shown with Alg. 7, in Appendix A. Thus, without thresholds and for powers of two, the dominant term of the overall arithmetic cost is 8n log 2 (7) , for the in-place version, roughly a third more operations than the 6n log 2 (7) dominant term of the cost for the version using extra temporaries.

Any bilinear algorithm for matrix multiplication (see, e.g., https://fmm.univ-lille.fr/) can be dealt with similarly. Further, even the accumulating version of the non bilinear algorithm of [START_REF] Dumas | Some fast algorithms multiplying a matrix by its adjoint[END_REF] can benefit from our techniques of in-place accumulation (see Appendix B).

We now prove that 18 additions is the minimal number of additions required by an in-place algorithm resulting from any bilinear algorithm for matrix multiplication using only 7 multiplications. For this we need to consider elementary operations on variables (similar to elementary linear algebra operators): variableswitching (swapping variable i and variable j); variable-multiplying (multiplying a variable by a constant); variable-addition (adding one variable, potentially multiplied by a constant, to another variable). An elementary program is a program making use of only these three kind of operations. Then we consider the in-place implementation of a linear function on its input: for α ∈ F t×m and a ∈ F m , we want to compute at least once each of the t coefficients of α a, but using only elementary operations and using only the variables of a as temporary variables. We start by proving in Lemma 6 that in any bilinear algorithm for matrix multiplication using only 7 multiplications, the columns of the associated matrices α, β, µ (as in Eq. ( 4)) cannot contain too many zeroes.

Lemma 6. If (α, β, µ) ∈ F 7×4 × F 7×4 × F 4×7 is
a bilinear algorithm for matrix multiplication, then none of α, β, µ ⊺ contain a zero column vector, nor a multiple of a vector of the canonical basis.

Proof. The dimensions of the matrices indicate that the multiplicative complexity of the algorithm is 7.

From [START_REF] De Groote | On varieties of optimal algorithms for the computation of bilinear mappings II[END_REF] we know that all such bilinear algorithms can be obtained from one another. Following [START_REF] Nader | On the additive complexity of 2 × 2 matrix multiplication[END_REF]Lemma 6], then any associated α, β, µ ⊺ matrix is some row or column permutation, or the multiplication by some G ⊗ H (the Kronecker product of two invertible 2×2 matrices), of the matrices of Eq. ( 4). By duality [START_REF] Hopcroft | Duality applied to the complexity of matrix multiplication and other bilinear forms[END_REF], see also [START_REF] Nader | On the additive complexity of 2 × 2 matrix multiplication[END_REF]Eq. (3)], it is also sufficient to consider any one of the 3 matrices. We thus let

G = g1,1 g1,2 g2,1 g2,2 , H = h1,1 h1,2
h2,1 h2,2 and K = G ⊗ H be their Kronecker product. Then any column of K is of the form ux uy vx vy

⊺ , for u ∈ {g 1,1 , g 1,2 }, v ∈ {g 2,1 , g 2,2
}, with u and v both in the same column of G and x ∈ {h 1,1 , h 1,2 }, y ∈ {h 2,1 , h 2,2 }, with x and y both in the same column of H. Further as G is invertible, u and v cannot be both zero simultaneously and, similarly, x and y cannot be both zero simultaneously. Now consider for instance α • K, with α of Eq. ( 4). Then any column θ of α • K is of the form: ux uy -ux -uy + vx + vy vy vx + vy -ux + vx -ux + vx + vy ⊺ . For such a column to be a multiple of a vector of the canonical basis or the zero vector, at least 6 of its 7 coefficients must be zero. For instance, this means that at least two out of rows 1, 2 and 4 must be zero: or that at least two of ux, uy or vy must be zero. This limits us to three cases: (1) u = 0, (2) y = 0 or (3) x = v = 0. If u = 0, then θ = v 0 0 x + y y x + y x x + y ⊺ ; at least one of rows 4 or 6 has to be zero, thus, w.l.o.g. suppose x = 0, we obtain that θ = vy 0 0 1 1 1 0 1 ⊺ with none of v nor y being zero (otherwise G or H is not invertible); such a column cannot be a multiple of a vector of the canonical basis nor the zero vector. Similarly, if y = 0, then θ = x u 0 -u + v 0 v -u + v -u + v ⊺ ; at least one of rows 1 or 5 has to be zero, thus, w.l.o.g. suppose v = 0, we obtain that θ = ux 1 0 -1 0 0 -1 -1 ⊺ ; such a column cannot be a multiple of a vector of the canonical basis nor the zero vector. Finally, if x = v = 0, then θ = uy 0 1 -1 0 0 0 0 ⊺ ; again that column cannot be a multiple of a vector of the canonical basis nor the zero vector. Now we show that any in-place elementary algorithm requires at least 1 extra operation to put back the input in its initial state.

Lemma 7. Let a ∈ F m and α ∈ F t×m with at least one row which is neither the zero row, nor a vector of the canonical basis. Now suppose that, without any constraints in terms of temporary registers, k is the minimal number of elementary operations required to compute α a. Then any algorithm computing all the t values of α a, in-place of a, requires at least k + 1 elementary operations.

Proof. Consider an in-place algorithm realizing α a in f operations. Any zero or canonical basis vector row can be realized without any operations on a. Now take this algorithm at the moment where the last of the other rows of α are realized (at that point all the t values are realized). Then this last realization (a non trivial linear combination of the initial values of a) has to have been stored in one variable of a, say a i . Therefore, at this point, the in-place algorithm has to perform at least one more operation to put back a i to its initial state. Therefore, by replacing all the in-place computations by operations on extra registers and omitting the operation(s) that restore this a i , we obtain an algorithm with less than f -1 elementary operations that realizes α a and thus: (f -1) ≥ k. Proposition 8. For the in-place realization of each of the two linear operators α and β, of any bilinear matrix multiplication algorithm using only 7 multiplications, and the restoration of the initial states of their input, at least 6 operations are needed.

Proof. A bilinear matrix multiplication algorithm has to compute α a, with a the entries of the left input of the matrix multiplication, while β deals with the right input. These α and β matrices cannot contain a (4-dimensional) zero row: otherwise there would exist an algorithm using less than 6 multiplications, but 7 is minimal [START_REF] Winograd | On multiplication of 2×2 matrices[END_REF]. If α or β contain at least 5 rows that are not vectors of the canonical basis, then they require at least 5 non-trivial operations to be computed, and therefore at least 6 elementary operations with an in-place algorithm, by Lemma 7. The matrices also cannot contain more than 3 multiple of vectors of the canonical basis, by [7, Lemma 8] (and thus require at least 4 operations to be computed).

There thus remains now only to consider matrices with exactly 3 rows that are multiple of vectors of the canonical basis. Let M be the 4×4 sub-matrix obtained from α (or β) by removing those 3 vectors of the canonical basis. By Lemma 6, no column of M can be the zero column: otherwise a 7-dimensional column of α (or β) would be either a multiple of a canonical basis vector, or the zero vector. This means that every variable of a has to be used at least once to realize the 4 operations of M a. Now suppose that there exists an in-place algorithm realizing M a in 5 elementary operations. Any operations among these 5 that, as its results, puts back a variable into its initial state, does not realize any row of M a (because putting back a variable to its initial state is the trivial identity on this initial variable, and this would be represented by a 4-dimensional vector of the canonical basis, which M do not contain, by construction). Therefore, at most one among these 5 operations puts back a variable of a into its initial state (otherwise M a, and therefore α a or β a, would be realizable in strictly less than 4 operations). This means that at most one variable of a can be modified during the algorithm (otherwise the algorithm would not be able to put back all its input variables into their initial state).

W.l.o.g suppose this only modified variable is a 1 . Finally, as all the other 3 variables must be used in at least one of the 5 elementary operations, at least 3 operations are of the form a 1 += λ i a i for i = 2, 3, 4 and some constants λ i . After those, to put back a 1 into its initial state, each one of these 3 independent variables, a 2 , a 3 and a 4 , must be "removed" from a 1 at some point of the elementary program. But, with a total of 5 operations, there remains only 2 other possible elementary operations, each one of those modifying only a 1 . Therefore not all 3 variables can be removed and thus no in-place algorithm can use only 5 operations.

Finally, there remains to consider the linear combinations of the 7 multiplications to conclude that Alg. 7 realizes the minimal number of operations for any in-place algorithm with 7 multiplications. Theorem 9. At least 25 additions are required to compute in-place any bilinear matrix multiplication algorithm using only 7 multiplications and to restore its two input matrices to their initial states afterwards.

Proof. Proposition 8 shows that at least 6 operations are required to realize α (or β). For µ, we in fact compute c += µ ρ, so we need to consider the matrix P = [ I4 µ ] ∈ F 4×11 and the vector ξ = c ρ . Consider now an elementary program that realizes P ξ, in-place of c only. This implies for instance that if ρ is zero, c should be put back to its initial state at the end of the program. Finally, thus consider the transposed program P ⊺ c: it must be in-place of c, while putting back c to its initial state afterwards. Now, by Prop. 8, µ ⊺ , and thus P ⊺ ∈ F 11×4 , requires at least 6 elementary operations to be performed. By Tellegen's transposition principle, see also [START_REF] Kaminski | Addition requirements for matrix and transposed matrix products[END_REF]Theorem 7], computing the transposed program thus requires at least 6 + (11 -4) = 13 operations. This gives a total of at least 6 + 6 + 13 = 25 additions.

Theorem 9 thus shows that our Alg. 7 with 18 elementary additions and 7 recursive calls (thus 7 more, and a total of 18 + 7 = 25 additions) is an optimal in-place bilinear matrix multiplication algorithm using only 7 multiplications.

In-place polynomial multiplication with accumulation

Algorithm 2 can also be used for polynomial multiplication. An additional difficulty is that this does not completely fits the setting, as multiplication of two size-n inputs will in general span a (double) size-2n output. This is not an issue until one has to distribute separately the two halves of this 2n values (or more generally to different parts of different outputs). In the following we show that this can anyway always be done for polynomial multiplications.

In-place accumulating Karatsuba

For instance, we immediately obtain an in-place Karatsuba polynomial multiplication since it writes as in Eq. ( 5), from which we can extract the associated µ, α, β matrices shown in Eq. ( 6).

(Y a 1 + a 0 )(Y b 1 + b 0 ) = a 0 b 0 + Y 2 (a 1 b 1 ) + Y ((a 1 + a 0 )(b 1 + b 0 ) -a 0 b 0 -a 1 b 1 ) (5)
µ = 1 0 0 -1 1 -1 0 0 1 α = 1 0 1 1 0 1 β = 1 0 1 1 0 1 (6)
Then, with Y = X t and a i , b i , c i polynomials in X (and a 0 , b 0 , c 0 of degree less than t), this is detailed, with accumulation, in Eq. ( 7):

A(Y ) = Y a 1 + a 0 ; B(Y ) = Y b 1 + b 0 ; C(Y ) = Y 3 c 11 + Y 2 c 10 + Y c 01 + c 00 ; m 0 = a 0 • b 0 = m 01 Y + m 00 ; m 2 = a 1 • b 1 = m 21 Y + m 20 ; m 1 = (a 0 + a 1 ) • (b 0 + b 1 ) = m 11 Y +
To deal with the distributions of each half of the products of Eq. ( 7), each coefficient in µ in Eq. ( 6) can be expanded into 2×2 identity blocks, and the middle rows combined two by two, as each tensor product actually spans two sub-parts of the result; we obtain Eq. ( 8):

µ (2) = I2 02 02 02 02 02 + 0 0 0 -I2 I2 -I2 0 0 0 + 02 02 02 02 02 I2 = 1 0 0 0 0 0 -1 1 1 0 -1 0 0 -1 0 1 1 -1 0 0 0 0 0 1 (8)
Finally, Eq. ( 7) then translates into an in-place algorithm thanks to Alg. 2 and Eqs. ( 6) and [START_REF] Buhrman | Computing with a full memory: catalytic space[END_REF]. The first point is that products double the degree: this corresponds to a constraint that the two blocks have to remain together when distributed. In other words, this means that the µ (2) matrix needs to be considered two consecutive columns by two consecutive columns. This is always possible if the two columns are of full rank 2. Indeed, consider a 2 × 2 invertible submatrix M = [ v w

x y ] of these two columns. Then computing

[ ci cj ] += M [ ρ0 ρ1
] is equivalent to computing a 2 × 2 version of Eq. ( 1):

[ ci cj ] * = M -1 ; [ ci cj ] += [ ρ0 ρ1 ] ; [ ci cj ] * = M . (9)
The other rows of these two columns can be dealt with as before by pre-and post-multiplying/dividing by a constant and pre-and post-adding/subtracting the adequate c i and c j . Now to apply a matrix M = a b c d to a vector of results u v , it is sufficient that one of its coefficient is invertible. W.l.o.g suppose that its upper left element, a, is invertible. Thus a b c d =

1 0 ca -1 1 a b 0 d-ca -1 .
Then the in-place evaluation of Eq. ( 10) performs this application, using the two (known in advance) constants x = ca -1 and y = d -ca -1 b:

u * = a u += b • v v * = y v += x • u      computes in-place: u v ← a b c d ⊙ u v = a u+b v c u+d v for x = ca -1 and y = d -xb (10)
Remark 10. In practice for 2 × 2 blocks, if a is not invertible, permuting the rows is sufficient since c has to be invertible for the matrix to be invertible: for

J = [ 0 1 1 0 ], if M = c d 0 b = J•M , then M = J• M and M -1 = M -1 •J so that Eq. (9) just becomes: [ ci cj ] * = J; [ ci cj ] * = M -1 ; [ ci cj ] += [ ρ0 ρ1 ] ; [ ci cj ] * = M ; [ ci cj ] * = J.
We now have all the tools to produce in-place polynomial algorithms. We start, in Alg. 3, with a version of Alg. 2 that regroups the intermediate computations into consecutive blocks. Theorem 11. Algorithm 3 is correct, in-place, and requires t MUL-2D, 2(#α + #β + #µ -t) ADD and 2(♯α + ♯β + ♯µ + 2t) SCA operations.

Proof. Thanks to Eqs. ( 9) and [START_REF] Coxon | An in-place truncated Fourier transform[END_REF] and Remark 10, correctness is similar to that of Alg. 2 in Theorem 3. Then, Eq. ( 10) requires 4 SCA and 2 ADD operations and is called 2t times. The rest is similar to Alg. 2 and amounts to 2t + 2(#α -t + #β -t + #µ -2t) + (2t)2 ADD and 2(♯α + ♯β + ♯µ -2t) + (2t)4 SCA operations. Let i s.t. α ℓ,i = 0; ai * = α ℓ,i ;

Algorithm 3 In place bilinear 2 by 2 formula

Input: a ∈ F m , b ∈ F n , c ∈ F s ; α ∈ F t×m , β ∈ F t×n , µ ∈ F s×(2t) = [ M 1 ... Mt ],
3:

for λ = 1 to m, λ = i, α ℓ,λ = 0 do ai += α ℓ,λ a λ end for 4:

Let j s.t. β ℓ,j = 0; bj * = β ℓ,j ;

5:

for λ = 1 to n, λ = j, β ℓ,λ = 0 do bj += β ℓ,λ b λ end for

6: Let k, f s.t. M = µ k,2ℓ µ k,2ℓ+1 µ f,2ℓ µ f,2ℓ+1
is invertible;

7:

c k c f ← M -1 c k c f
{Via Eq. ( 10) and Remark 10}

8:

for λ = 1 to s, λ ∈ {f, k}, µ λ,2ℓ = 0 do c λ -= µ λ,2ℓ c k end for 9:

for λ = 1 to s, λ ∈ {f, k}, µ λ,2ℓ+1 = 0 do c λ -= µ λ,2ℓ+1 c f end for 10:

c k c f += ai • bj {This is the accumulation of the product m k m f } 11:
for λ = 1 to s, λ ∈ {f, k}, µ λ,2ℓ+1 = 0 do c λ += µ λ,2ℓ+1 c f end for 12:

for λ = 1 to s,λ ∈ {f, k}, µ λ,2ℓ = 0 do c λ += µ λ,2ℓ c k end for 13:

c k c f ← M c k c f
{Via Eq. ( 10) and Remark 10, undo 7}

14:

for λ = 1 to n, λ = j, β ℓ,λ = 0 do bj -= β ℓ,λ b λ end for ; bj /= β ℓ,j ;

15:

for λ = 1 to m, λ = i, α ℓ,λ = 0 do ai -= α ℓ,λ a λ end for ; ai /= α ℓ,i ;

16: return c.

There remains to use a double expansion of the output µ ∈ F s×t to simulate the double size of the intermediate products (MUL-2D), producing µ (2) ∈ F s×(2t) , as in Eq. ( 8), that is used as an input in Alg. 3. This double expansion matrix is obtained by the following duplication: µ (2) (i, 2j) = µ(i, j) and µ (2) (i + 1, 2j + 1) = µ(i, j) for i = 1..s and j = 1..t. We prove, in Lemma 12, that in fact any such double expansion of a representative matrix is suitable for the in-place computation of Alg. 3. Lemma 12. If µ does not contain any zero column, then each pair of columns of µ (2) , resulting from the expansion of a single column in µ, contains an invertible lower triangular 2×2 submatrix.

Proof. Consider the top most non-zero element of a column. It is expanded as a 2×2 identity matrix whose second row is merged with the first row of the next identity matrix: in picture, [ a b ] is expanded to

a 0 b a * b .
For instance with m 00 + Y m 01 = a 0 b 0 = ρ 0 + Y ρ 1 , consider the upper left 2 × 2 block of µ (2) in Eq. ( 8), that is

M = 1 0 -1 1 , whose inverse is M -1 = [ 1 0 1 1 ]. One has first to precompute M -1 [ c00 c01 ]
, that is nothing on c 00 and c 01 += c 00 for the second coefficient. Then, afterwards, the third row, for c 10 , will just be -m 01 : for this just pre-add c 10 += c 01 , and post-subtract c 10 -= c 01 after the product actual computation. This example is exactly lines 14 to 16 of Alg. 4 thereafter. One could also consider instead the first and last rows, but in this particular case without any advantage in term of computations. To complete Eq. ( 7), the computation of m 2 is dealt with in the same manner, while that of m 1 is direct in the results (and requiring pre and post additions of its inputs). This gives then the whole of Alg. 4. The second point is to deal with unbalanced dimensions and degrees for Y = X t and recursive calls: first separate the largest polynomial in two parts, so that two sub-products are performed: a large balanced one, and, recursively, a smaller unbalanced one. Second, for the balanced case, the idea is to ensure that three out of four parts of the result, t 00 , t 01 and t 10 , have the same size and that the last one t 11 is smaller. This ensures that all accumulations can be performed in-place. Again the details can be found in Alg. 4. Proposition 13. Alg. 4 is correct and requires O(mn log 2 (3)-1 ) operations.

Proof. With the above analysis, correctness comes from that of Alg. 3 applied to Eq. ( 6 return the quadratic in-place multiplication. {Alg. 1} 3: else if m > n then 4:

Let A(X) = A 0 (X) + X n+1 A 1 (X)

Let A = a 0 + X t a 1 ; B = b 0 + X t b 1 ;
13:

Let C = c 00 + c 01 X t + c 10 X 2t + c 11 X 3t ; {d • c 11 = 2n -3t} 14:
c 01 += c 00 ; c 10 += c 01 ;

15:

[ c00 c01 ] += a 0 • b 0 {Recursive call for m 0 }
a 0 += a 1 ; b 0 += b 1 ; {d • a 0 = t ≥ n -t = d • a 1 } 18: [ c01 c10 ] += a 0 • b 0 {Recursive call for m 1 } 19: b 0 -= b 1 ; a 0 -= a 1 ; 20:
c 10 += c 11 ; c 01 += c 10 ;

21:

[ c10 c11 ] += a 1 • b 1 {Recursive call for m 2 }
) ≤ O( κ-1 i=1 u i u log 2 (3)-1 i+1
) operations. But, u i+1 ≤ u 2 = n and we let s i = u i + u i+1 so that u i = s i -u i+1 ≤ s i . Now, from [16, Corollary 2.6], we have that s i ≤ s 1 (2/3) i-1 . Thus the number of operations is bounded by O(

κ-1 i=1 s i n log 2 (3)-1 ) ≤ O(n log 2 (3)-1 s 1 ( 1 1-(2/3) -1)) = O(n log 2 (3)-1 (m + n)) = O(mn log 2 (3)-1 ).
In order to compare the complexity when m = n, note that all coefficients of α, β and µ (2) being 1 or -1, Alg. 4 does compute the accumulation C += AB without constant multiplications. Also, the de-duplication of each recursive output has enabled some natural reuse, so in fact there is a cost of 2(#α -t + #β -t) = 2(4 -3 + 4 -3) with t = 3, and 2(2(#µ -t) = 4(5 -3) = 2(#µ (2) -2t), for a total of 3 recursive accumulating calls and at most 12 block additions. Thus, without thresholds, for powers of two and with m = n, the dominant term of the overall cost only goes from 10n log 2 (3) , for the simple application of Eq. ( 7), to 14n log 2 (3) , for the fully in-place version in Alg. 4.

Further bilinear polynomial multiplications

We have shown that any bilinear algorithm can be transformed into an in-place version. This approach thus also works for any Toom-k algorithm using 2k -1 interpolations points instead of the three points of Karatsuba (Toom-2).

For instance Toom-3 uses interpolations at 0, 1, -1, 2, ∞. Therefore, α and β are the Vandermonde matrices of these points for the 3 parts of the input polynomials and µ is the inverse of the Vandermonde matrix of these points for the 5 parts of the result, as shown in Eq. ( 11) thereafter.

µ = 1 0 0 0 0 1 1 1 1 1 1 -1 1 -1 1 1 -2 4 -8 16 0 0 0 0 1 -1 =     1 0 0 0 0 1 2 1 3 -1 1 6 -2 -1 1 2 1 2 0 -1 - 1 2 1 6 1 2 - 1 6 2 0 0 0 0 1     ; α = β = 1 0 0 1 1 1 1 -1 1 1 -2 4 0 0 1 (11)
With the same kind of duplication as in Eq. ( 8), apart from the recursive calls, the initially obtained operation count is 2(11+11-2 * 5)+2(2(16-5)) = 68 additions and 2(2+2+2( 11)) = 52 scalar multiplications. Following the optimization of [START_REF] Bodrato | Towards optimal Toom-Cook multiplication for univariate and multivariate polynomials in characteristic 2 and 0[END_REF], we see in α and β that the evaluations at 1 and -1 (second and third rows) share one addition. As they are successive in our main loop, subtracting one at the end of the second iteration, then followed by re-adding it at the third iteration can be optimized out. This is 2 less operations. Together with shared coefficients in the rows of µ, some further optimizations of [START_REF] Bodrato | Towards optimal Toom-Cook multiplication for univariate and multivariate polynomials in characteristic 2 and 0[END_REF] can probably also be applied, where the same multiplicative constants appear at successive places.

Fast bilinear polynomial multiplication

When sufficiently large roots of unity exist, polynomial multiplications can be computed fast in our in-place model via a discrete Fourier transform and its inverse, as shown in Alg. 5, for power of two dimensions, and in Alg. 6, for general dimensions.

Let F ∈ D[X] of degree < n and ω be a principal n-th root of unity, where n = 2 p . The discrete Fourier transform of F at ω is defined as

DFT n (F, ω) = (F (ω 0 ), F (ω 1 ), . . . , F (ω n-1 )). The map is in- vertible, of inverse DFT -1 n (•, ω) = 1 n DFT n (•, ω -1
). Further, the DFT can be computed over-place, replacing the input by the output [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF]. Actually, for over-place algorithms and their extensions to the truncated Fourier transform, it is more natural to work with the bit-reversed DFT defined by brDFT

n (F, ω) = (F (ω [0]p ), F (ω [1]p ), . . . , F (ω [n-1]p )) where [i] p = p-1 j=0 d j 2 p-j is the length-p bit reversal of i = p-1 j=0 d j 2 j , d j ∈ {0, 1}. If π denotes the bit-reversal permutation (that is π(i) = [i] p ), we have brDFT n (•, ω) = π • DFT n (•, ω). Its inverse becomes brDFT -1 n (•, ω) = 1 n DFT n (•, ω -1 ) • π = 1 n DFT n (π(•), ω -1
) since π is an involution. Remark 14. The Fast Fourier Transform (FFT) algorithm has two main variants: decimation in time (DIT) and decimation in frequency (DIF). Both algorithms can be performed over-place, replacing the input by the output. Without applying any permutation to the entries of the input/output vector, the over-place DIF-FFT algorithm naturally computes brDFT n (•, ω), while the over-place DIT-FFT algorithm on ω -1 computes n • brDFT -1 n (•, ω).

Algorithm 5 In-place power of two accumulating multiplication

Input: a, b and c of length 2 L , 2 L and 2 L+1 , containing the coefficients of A, B, C ∈ D[X] respectively; ω ∈ D primitive 2 L+1 -th root of unity.

Output: c contains the coefficients of C + A • B. 1: Let n = 2 L ; 2: c ← brDFT 2n ( c, ω); {over-place} 3: a ← brDFT n ( a, ω 2 ); b ← brDFT n ( b, ω 2 ) {over-place} 4: for i = 0 to n -1 do c i += a i × b i end for 5: a ← brDFT -1 n ( a, ω 2 ); b ← brDFT -1 n ( b, ω 2 ) {Undo 3} 6: for i = 0 to n -1 do a i * = ω i ; b i * = ω i end for 7: a ← brDFT n ( a, ω 2 ); b ← brDFT n ( b, ω 2 ) {over-place} 8: for i = 0 to n -1 do c i+n += a i × b i end for 9: a ← brDFT -1 n ( a, ω 2 ); b ← brDFT -1 n ( b, ω 2 ) {Undo 7} 10: for i = 0 to n -1 do a i /= ω i ; b i /= ω i end for {Undo 6} 11: return c ← brDFT -1 2n ( c, ω)
Theorem 15. Using an over-place brDFT algorithm with complexity bounded by O(n log n), Alg. 5 is correct, in-place and has complexity bounded by O(n log n).

Proof. Alg. 5 follows the pattern of the standard FFT-based multiplication algorithm. Our goal is to compute brDFT 2n (A, ω), brDFT 2n (B, ω) and brDFT 2n (C, ω), then obtain brDFT 2n (C + AB, ω) and finally C + AB using an inverse brDFT. Computations on C and then C + AB are performed over-place using any standard over-place brDFT algorithm. The difficulty happens for A and B that are stored in length-n arrays. We use the following property of the bit reversed order: for

k < n/2, [k] p = 2[k] p-1 , and for k ≥ n/2, [k] p = 2[k -n/2] p-1 + 1. Therefore, the first n coefficients of brDFT 2n (A, ω) are (A(ω 2[0]p-1 ), . . . A(ω 2[ n 2 -1]p-1 )) = brDFT n (A, ω 2 ).
Similarly, the next n coefficients are brDFT n (A(ωX), ω 2 ). Therefore, one can compute brDFT(A, ω 2 ) and brDFT(B, ω 2 ) in a and b respectively, and update the first n entries of c. Next we restore a and b using brDFT -1 n (•, ω 2 ). We compute A(ωX) and B(ωX) and again brDFT(A(ωX), ω 2 ) and brDFT(B(ωX), ω 2 ) to update the last n entries of c. Finally, we restore a and b and perform the inverse brDFT on c. The cost is dominated by the ten brDFT ±1 computations.

Since there are 2 size-2n and 8 size-n brDFT ±1 computations, the dominant term of the cost is 18n log n ring operations, twice as large as the dominant term for the standard (not in-place) algorithm [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF].

The case where n is not a power of two is loosely similar, using as a routine a truncated Fourier transform (TFT) rather than a DFT [START_REF] Van Der Hoeven | The Truncated Fourier Transform and Applications[END_REF]: let ω be an N -th root of unity for some N = 2 p . The length-n (bit-reversed) TFT of a polynomial

F ∈ D[X], n < N , is brTFT n (F, ω) = (F (ω [0]p ), . . . , F (ω [n-1] p
)), that is the n first coefficients of brDFT N (F, ω). As for the DFT, the (bit-reversed) TFT and its inverse can be computed over-place [START_REF] Harvey | An in-place truncated Fourier transform and applications to polynomial multiplication[END_REF][START_REF] Daniel | Efficient Computation with Sparse and Dense Polynomials[END_REF][START_REF] Arnold | A new truncated Fourier transform algorithm[END_REF][START_REF] Coxon | An in-place truncated Fourier transform[END_REF].

Given inputs A and B ∈ D[X] of respective lengths m and n and an output C ∈ D[X] of length m + n -1, we aim to replace C by C + AB. The idea is first to replace C by brTFT m+n-1 (C, ω) where ω is a 2 p -th principal root of unity, 2 p ≥ m + n -1. That is, the vector c now contains as its i-th entry the value C(ω [i]p ). The goal is then to replace C(ω

[i]p ) by C(ω [i]p ) + A(ω [i]p )B(ω [i]p
), for i = 0 to m + n -2. We cannot compute the length m + n -1 brTFT's of A and B since this takes too much space. Instead, we will progressively compute some parts of these brTFT's by means of (standard) brDFT's, and update c accordingly. The starting point of this strategy is the following lemma. Lemma 16 ([18,[START_REF] Daniel | Efficient Computation with Sparse and Dense Polynomials[END_REF]). Let F ∈ D[X], ℓ and s be two integers such that 2 ℓ divides s and ω be a 2 p -th principal root of unity. Define F s,ℓ (X) = F (ω [s]p X) mod X 2 ℓ -1 . Then brDFT 2 ℓ (F s,ℓ , ω 2 p-ℓ ) = (F (ω [s]p ), . . . , F (ω [s+2 ℓ -1]p )).

Proof. Let ω ℓ = ω 2 p-ℓ . This is a principal 2 ℓ -th root of unity since ω is a principal 2 p -th root of unity. In particular, for any i < 2 ℓ , F s,ℓ (ω

[i] ℓ ℓ ) = F (ω [s]p ω [i] ℓ ℓ ). Now, ω [i] ℓ ℓ = ω [i]p since 2 p-ℓ [i] ℓ = [i] p . Furthermore, [s] p + [i] p = [s + i] p since i < 2 ℓ and 2 ℓ divides s. Finally, F s,ℓ (ω [i] ℓ ℓ ) = F (ω [s+i]p ).
Corollary 17. Let F ∈ D[X] stored in an array f of length n. Let also ℓ and k be two integers and ω be a 2 p -th principal root of unity, with 2 ℓ ≤ n and (k+1)2 ℓ ≤ 2 p . There exists an algorithm, partTFT k,ℓ ( f , ω), that replaces the first 2 ℓ entries of f by Proof. Algorithm partTFT k,ℓ ( f , ω) is the following:

F (ω [k•2 ℓ ]p ), . . . , F (ω [(k+1)•2 ℓ -1]p ),
1: for i = 0 to n -1 do f i * = ω i[k•2 ℓ ]p end for 2: for i = 2 ℓ to n -1 do f i-2 ℓ += f i end for 3: f 0..2 ℓ -1 ← brDFT 2 ℓ ( f 0..2 ℓ -1 , ω 2 p-ℓ )
Its correctness is ensured by Lemma 16. Its inverse algorithm partTFT -1 k,ℓ ( f , ω) does the converse:

1: f 0..2 ℓ -1 ← brDFT -1 2 ℓ ( f 0..2 ℓ -1 , ω 2 p-ℓ ) 2: for i = 2 ℓ to n -1 do f i-2 ℓ -= f i end for 3: for i = 0 to n -1 do f i /= ω i[k•2 ℓ ]p end for
In both algorithms, the call to brDFT ±1 has cost O(ℓ • 2 ℓ ), and the two other steps have cost O(n).

To implement the previously sketched strategy, we assume that m ≤ n for simplicity. We let ℓ, t be such that 2 ℓ ≤ m < 2 ℓ+1 and 2 ℓ+t ≤ n < 2 ℓ+t+1 . Using partTFT ±1 , we are able to compute (A(ω [k•2 ℓ ]p ), . . . , A(ω [(k+1)•2 ℓ -1]p )) for any k, and restore A in its initial state afterwards. Similarly, it is possible to compute (B(ω [k•2 ℓ+t ]p ), . . . , B(ω [(k+1)•2 ℓ+t -1]p )) and restore B. ℓ ← ⌊log 2 min{r, m}⌋; t ← ⌊log 2 min{r, n}⌋ -ℓ;

5: k ← m + n -1 -r {over-place: B(ω [k•2 ℓ+t ]p] ), . . . , B(ω [(k+1)•2 ℓ+t -1]p )} 6: b ← partTFT k,ℓ+t ( b, ω) 7:
for s = 0 to 2 t -1 do {over-place:

A(ω [(k•2 t +s)2 ℓ ]p] ), . . . , A(ω [(k•2 t +s+1)2 ℓ -1]p )} 8: a ← partTFT s+k•2 t ,ℓ ( a, ω) 9:
for Proof. The fact that the algorithm is in-place comes from Corollary 17. The only slight difficulty is to produce, fast and in-place, the relevant roots of unity. This is actually dealt with in the original over-place TFT algorithm [START_REF] Harvey | An in-place truncated Fourier transform and applications to polynomial multiplication[END_REF] and can be done the same way here.

i = 0 to 2 ℓ -1 do c i+(k•2 t +s)2 ℓ += a i b i+s•2 ℓ end for 10: a ← partTFT -1 s+k•2 t ,ℓ ( a,
To assess its correctness, first note that the values of Lines 4 and 5 are computed so that 2 ℓ ≤ r, m and 2 ℓ+t ≤ r, n. One iteration of the while loop updates the entries c k to c k+2 ℓ+t -1 where k = m + n -1 -r. To this end, we first compute B(ω [k•2 ℓ+t ]p] ) to B(ω [(k+1)•2 ℓ+t -1]p ) in b using partTFT. Then, since a may be too small to store 2 ℓ+t values, we compute the corresponding evaluations of A by groups of 2 ℓ , using a smaller partTFT. After each computation in a, we update the corresponding entries in c and restore a. Finally, at the end of the iteration, entries k to k + 2 ℓ+t -1 of c have been updated and b can be restored. This proves the correctness of the algorithm.

We now bound its complexity. Since m ≤ n, our aim is to bound it by O(n log n). Let us first bound the number of iterations of the while loop. We identify two phases, first iterations where r ≥ n and then iterations with r < n. During the first phase, 2 ℓ+t > n 2 entries of c are updated at each iteration, hence the first phase has at most 3 iterations. In the second phase, 2 ℓ+t > r 2 entries are updated per iteration. The second phase starts with r < n and each iteration decreases r by half, hence the second phase has at most log 2 n iterations.

The cost of an iteration is dominated by the calls to partTFT ±1 . The cost of a call to partTFT ±1 k,ℓ with a size-m input is the sum of a linear term O(m) and a non-linear term O(ℓ • 2 ℓ ). At each iteration, there are two calls to partTFT ±1 on b and 2 t+1 calls to partTFT ±1 on a. The linear terms sum to O

(n + m • 2 t ) = O(n) since m • 2 t < 2 ℓ+1+t ≤ 2n.
Over the log 2 n iterations, the global cost due to these linear terms is O(n log n).

The cost due to the non-linear terms in one iteration is O((ℓ + t) • 2 ℓ+t ). In the first iterations, 2 ℓ+t ≤ n and these costs sum to O(n log n). In the next iterations, 2 ℓ+t ≤ r < n. Since r is halved at each iteration, the non-linear costs in these iterations sum to O i n

2 i log n 2 i = O(n log n).

Conclusion

We here provide a generic technique mapping any bilinear formula (and more generally any linear accumulation) into an in-place algorithm. This allows us for instance to provide the first accumulating in-place Strassen-like matrix multiplication algorithm. This also allows use to provide fast in-place accumulating polynomial multiplications algorithms. Many further accumulating algorithm can then be reduced to these fundamental building blocks, see for instance Toeplitz, circulant, convolutions or remaindering operations in [START_REF] Dumas | In-place fast polynomial modular remainder[END_REF]. For the in-place realization of each of linear operators arising from the sparsification of any bilinear matrix multiplication algorithm using only 7 multiplications, and for the restoration of the initial states of their input, at least 6 operations are needed.

Proof. The alternate basis method of [START_REF] Karstadt | Matrix multiplication, a little faster[END_REF], consist in sparsifying the matrices of Eq. ( 4), via right multiplication by 4×4 invertible matrices. The sparsest obtained matrices are given in

   0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 1 -1 0 -1 1 0 0 0 -1 0 1    ,    0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 -1 0 1 0 1 -1 0 -1 1 0 0    , 0 0 0 1 1 0 0 0 0 1 0 1 -1 1 0 1 0 0 0 0 1 1 0 0 0 0 -1 0 (12)
We then follow the same line of reasoning as in Prop. 8, where we mostly need to adapt Lemma 6. W.l.o.g, we consider a 4×4 transformation M of the left matrix in Eq. [START_REF] Dumas | Some fast algorithms multiplying a matrix by its adjoint[END_REF].

If the resulting product matrix has only 3 rows that are multiple of a canonical vector, then 6 multiplications are minimal by Prop. 8. The only other possibility is thus, as in Eq. ( 12), that it contains exactly 4 rows that are multiples of a canonical vector. Now, each one of the four columns of the resulting product has the form:

a b c d c -b c -d a -c ⊺ .
Such a column can never be a zero-column, nor a multiple of a canonical vector. Therefore, as in Prop. 8 all 4 variables must appear in the 3 rows that are not canonical and at most one variable can be modified. So again at least 3 operations are of the form a 1 += λ i a i for i = 2, 3, 4 and some constants λ i . But then again at least 3 operations are required to put back a 1 in its initial state and no in-place algorithm can use strictly less than 6 operations.

B Fast In-place Square & Symmetric Rank-k Update

Thanks to Alg. 7 and with some care on transposes, the same technique can be adapted to, e.g., [START_REF] Dumas | Some fast algorithms multiplying a matrix by its adjoint[END_REF]Alg. 12], which performs the multiplication of a matrix by its transpose. With an accumulation, this is a classical Symmetric Rank-k Update (or SYRK): C ← αAA ⊺ + βC.

Following the notations of the latter algorithm, which is not a bilinear algorithm on its single input matrix, the in-place accumulating version is shown in Alg. 8, for α = β = 1, using any (fast to apply) skewunitary Y ∈ F n×n . It has been obtained automatically by the method of Theorem 5, and it thus preserves the need of only 5 multiplications P 1 to P 5 . It has then been scheduled so as to reduce the number of extra operations.

Algorithm 8 requires 3 recursive calls, 2 multiplications of two independent half matrices, 4 multiplications by a skew-unitary half matrix, 8 additions (of half inputs), 12 semi-additions (of half triangular outputs). Provided that the multiplication by the skew-unitary matrix can be performed in-place in negligible time, this gives a dominant term of the complexity bound for Alg. 8 of a fraction 2 2 ω -3 of the cost of the full in-place algorithm. This is a factor 1 2 , when Alg. 7 is used for the two block multiplications of independent matrices (P 4 and P 5). Now, the skew-unitary matrices used in [START_REF] Dumas | Some fast algorithms multiplying a matrix by its adjoint[END_REF], are either a multiple of the identify matrix, or the Kronecker product of a b -b a by the identity matrix, for a 2 + b 2 = -1 and a = 0. The former is easily performed in-place in time O(n 2 ). For the latter, it is sufficient to use Eq. ( 10): the multiplication a b -b a u can be realized in place by the algorithm: u 1 * = a; u 1 += b•u 2 ; u 2 * = (a + b 2 a -1 ); u 2 += -ba -1 •u 1 .

The same technique can be used on the symmetric algorithm for the square of matrices given in [START_REF] Bodrato | A Strassen-like matrix multiplication suited for squaring and higher power computation[END_REF]. The resulting in-place algorithm is given in Alg. 9.

C Reduced memory algorithms for polynomial multiplication

We compare in Table 2 the procedure given in Alg. 4 (obtained via the automatic application of Alg. 3) with previous Karatsuba-like algorithms for polynomial multiplications, designed to reduce their memory footprint (see also [START_REF] Giorgi | Efficient algorithms and implementation in exact linear algebra[END_REF]Table 2.2]).

  ). When m = n, with 3 recursive calls and O(n) extra operations, the algorithm thus requires overall O(n log 2 (3) ) operations. Otherwise, it requires m n equal degree calls, then a recursive call with n and (m mod n). Now, let u 1 = Algorithm 4 In-place Karatsuba polynomial multiplication with accumulation Input: A, B, C polynomials of degrees m, n, m + n with m ≥ n. Output: C += AB 1: if n ≤ Threshold then {Constant-time if Threshold ∈ O(1)} 2:

16 :

 16 c 10 -= c 01 ; c 01 -= c 00 ; {c 10 -m 01 and c 01 + m 01 -m 00 } 17:

22 :c

 22 01 -= c 10 ; c 10 -= c 11 ; {t 01 and t 10 in Eq. (7)} 23: return C. m and u 2 = n and if the Euclidean algorithm on them requires k steps, let u i for i = 1..k denote the successive residues within this Euclidean algorithm (and u k ∈ {0, 1}). Let κ = k -1 if u k = 0 and κ = k otherwise. Now, Alg. 4 requires less than O( κ-1 i=1 ⌊ ui ui+1 ⌋u log 2 (3) i+1

  and an inverse algorithm partTFT -1 k,ℓ that restores f to its initial state. Both algorithms use O(1) extra space and have complexity bounded by O(n + ℓ • 2 ℓ ).

Algorithm 6

 6 In-place fast accumulating polynomial multiplication Input: a, b and c of length m, n and m + n -1, m ≤ n, containing the coefficients of A, B, C ∈ D[X] respectively; ω ∈ D principal 2 p -th root of unity with 2 p-1 < m + n -1 < 2 p Output: c contains the coefficients of C + A • B. 1: c ← brTFT m+n-1 ( c, ω); {over-place} 2: r ← m + n -1 3: while r ≥ 0 do 4:

13 :Theorem 18 .

 1318 return c ← brTFT -1 m+n-1 ( c, ω) Algorithm 6 is correct and in-place. If the algorithm brDFT used inside partTFT has complexity O(n log n), then the running time of Alg. 6 is O((m + n) log(m + n)).

Algorithm 8

 8 In-place accumulating multiplication by its transposeInput: A = [ a11 a12 a21 a22 ] ∈ F m×2n ; symmetric C = c11 c21 ⊺ c21 c22 ∈ F m×m . Output: Low (C) += Low (A • A ⊺ ).{update bottom left triangle} Low(C 22 ) := Low(C 22 ) -Low(C 11 ); Low(C 21 ) := Low(C 21 ) -Low(C 11 ); Up(C 21 ) := Up(C 21 ) -Low(C 11 ) ⊺ ; Low(C 11 ) := Low(C 11 ) + Low(A 11 * A 11 ⊺ ); # P 1 Rec. Up(C 21 ) := Up(C 21 ) + Low(C 11 ) ⊺ ; Low(C 21 ) := Low(C 21 ) + Low(C 11 ); Low(C 22 ) := Low(C 22 ) + Low(C 11 ); Low(C 11 ) := Low(C 11 ) + Low(A 12 * A 12 ⊺ ); # P 2 Rec. A 11 := A 11 * Y; A 21 := A 21 * Y; A 11 := A 11 -A 21 ; A 21 := A 21 -A 22 ; Low(C 22 ) := Low(C 22 ) -Low(C 21 ); Low(C 22 ) := Low(C 22 ) -Low(C 21 ⊺ ); C 21 := C 21 + A 11 * A 21 ⊺ ; # P 4 (e.g. Alg. 7) Low(C 22 ) := Low(C 22 ) + Low(C 21 ⊺ ); A 21 := A 21 -A 11 ; Up(C 21 ) := Up(C 21 ) -Low(C 21 ) ⊺ ; Low(C 21 ) := Low(C 21 ) + Low(A 21 * A 21 ⊺ ); # P 5 Rec. Up(C 21 ) := Up(C 21 ) + Low(C 21 ) ⊺ ; Low(C 22 ) := Low(C 22 ) + Low(C 21 ); A 21 := A 21 + A 12 ; C 21 := C 21 + A 22 * A 21 ⊺ ; # P 3 (e.g. Alg. 7) A 21 := A 21 -A 12 ; A 21 := A 21 + A 11 ; A 21 := A 21 + A 22 ; A 11 := A 11 + A 21 ; A 21 := A 21 * Y -1 ; A 11 := A 11 * Y -1 ;

Algorithm 9

 9 In-place accumulating S.-W. matrix-squareInput: A = [ a11 a12 a21 a22 ], C = [ c11 c12 c21 c22 ]. Output: C += A 2 .A 22 := A 22 -A 21 ; C 12 := C 12 + C 22 ; C 22 := C 22 + A 22 * A 22 ; A 22 := A 22 + A 12 ; A 22 := A 22 -A 11 ; C 12 := C 12 -A 22 * A 12 ; C 21 := C 21 -A 21 * A 22 ; C 21 := C 21 -C 22 ; A 22 := A 22 + A 11 ; C 22 := C 22 -A 22 * A 22 ; C 11 := C 11 + C 22 ; C 22 := C 22 -A 12 * A 21 ; A 22 := A 22 + A 21 ; C 12 := C 12 -C 22 ; C 11 := C 11 -C 22 ; C 22 := C 22 + A 22 * A 22 ; A 22 := A 22 -A 12 ; C 21 := C 21 + C 22 ; C 11 := C 11 + A 11 * A 11 ;

  m 10 ; t 00 = c 00 + m 00 ; t 01 = c 01 + m 01 + m 10 -m 00 -m 20 ; t 10 = c 10 + m 11 + m 20 -m 01 -m 21 ; t 11 = c 11 + m 21 ; then C + AB = Y 3 t 11 + Y 2 t 10 + Y t 01 + t 00

  with no zero-rows in α, β, µ, s.t. (ai • bj ) fits two result variables c k , c l and s.t. Mi ∈ F s×2 is of full-rank 2 for i = 1..t. Read-only: α, β, µ. Output: c += µ m, for m = (α a) ⊙ (β b)

1: for ℓ = 1 to t do 2:

Table 1 :

 1 Algorithm 7 In-place accumulating S.-W. matrix-multiplication A 21 := A 21 -A 11 ; B 12 := B 12 -B 22 ; C 21 := C 21 -C 22 ; C 22 := C 22 + A 21 * B 12 ; A 21 := A 21 + A 22 ; B 12 := B 12 -B 11 ; C 12 := C 12 -C 22 ; C 22 := C 22 -A 21 * B 12 ; C 11 := C 11 -C 22 ; C 22 := C 22 + A 11 * B 11 ; C 11 := C 11 + C 22 ; B 12 := B 12 + B 21 ; C 21 := C 21 + C 22 ; C 21 := C 21 + A 22 * B 12 ; B 12 := B 12 + B 22 ; B 12 := B 12 -B 21 ; A 21 := A 21 -A 12 ; C 12 := C 12 -A 21 * B 22 ; A 21 := A 21 + A 12 ; A 21 := A 21 + A 11 ; C 22 := C 22 + A 21 * B 12 ; C 12 := C 12 + C 22 ; B 12 := B 12 + B 11 ; A 21 := A 21 -A 22 ; C 11 := C 11 + A 12 * B 21 ; Reduced-memory accumulating S.-W. multiplication Alg.

	Input: A = [ a11 a12 a21 a22 ], B = b11 b12 b21 b22 , C = [ c11 c12 c21 c22 ]. Output: C += AB.

Table 2 :

 2 Reduced-memory algorithms for Karatsuba polynomial multiplication Alg.Then, Alg. 6 is compared with previous FFT-based algorithms for polynomial multiplications designed to reduce their memory footprint in Table3(see also[START_REF] Giorgi | Efficient algorithms and implementation in exact linear algebra[END_REF] Table 2.2]).

		Memory	inputs	accumulation
	[30]	n + 5 log n read-only	✗
	[27, 28]	5 log n	read-only	✗
	[14]	O(1)	read-only	✗
	Alg. 4	5 log n	mutable	D

Table 3 :

 3 Reduced-memory algorithms for FFT polynomial multiplication Alg.

		Memory	inputs	accumulation
	[9]	2n	read-only	✗
	[27]	O(2 ⌈log 2 n⌉ -n) read-only	✗
	[18]	O(1)	read-only	✗
	Alg. 6	O(1)	mutable	D

We here give an in-place version of Strassen-Winograd algorithm for matrix multiplication. We first directly apply our Alg. 2 to the classical, not in-place Strassen-Winograd algorithm, following the specific scheduling strategy of Section 4. This strategy allows to reduce the number of additions obtained when calling Alg. 2, from 42 + 7 to 18 + 7: mostly remove successive additions/subtractions that are reciprocal on either submatrices. This optimized version is given in Alg. 7 and reaches the minimal possible number of extra additions/subtractions, as shown in Theorem 9.

The memory footprint of Alg. 7 is compared to that of previously known algorithms in Table 1.

To go beyond our minimality result for operations, one could try an alternate basis of [24]. But an argument similar to that of Prop. 8 shows that alternate basis does not help for the in-place case.