Brice Boyer
email: boyer@ncsu.edu.

Jean-Guillaume Dumas
email: guillaume.dumas@imag.fr.

Matrix multiplication over word-size prime fields using Bini's approximate formula

Bini's approximate formula (or border rank) for matrix multiplication achieves a better complexity than Strassen's matrix multiplication formula. In this paper, we show a novel way to use the approximate formula in the special case where the ring is Z/pZ. Besides, we show an implementation à la FFLAS-FFPACK, where p is a word-size prime number, that improves on state-of-the-art Z/pZ matrix multiplication implementations.

Introduction

A fast reliable matrix multiplication implementation over Z/pZ is crucial in exact linear algebra. Indeed, many algorithm rely on fast matrix multiplication as a building block (computation of factorised forms, characteristic polynomial, black box methods,. . . ). But also matrix multiplication on other rings such as the integer ring Z, polynomials Z/pZ[X] or Galois fields F q reduce to matrix multiplication over Z/pZ.

In [START_REF] Dumas | Dense linear algebra over word-size prime fields: the FFLAS and FFPACK packages[END_REF], matrix product over Z/pZ is computed efficiently with Strassen-Winograd's subcubic matrix multiplication algorithm ( [START_REF] Volker Strassen | Gaussian elimination is not optimal[END_REF][START_REF] Winograd | On multiplication of 2x2 matrices[END_REF]). A cascade algorithm is used: first a few recursive steps of Strassen's algorithm are performed, second, below a threshold, numerical BLAS matrix multiplication routines are called. The numerical routines are used in an exact manner and reductions modulo p are delayed as much as possible.

We use similar techniques in this article. Our main contribution is the use of an approximate algorithm to compute the exact matrix multiplication. We first recall Bini's algorithm ( [START_REF] Bini | O(n 2.7799 ) complexity for matrix multiplication[END_REF]) in Section 2. Then in Section 3, we show two fashions for the application of Bini's algorithm to get an exact matrix multiplication: first via rounding, second via p-adic expansion. We finally provide 1 in Section 4 a new schedule for the algorithm minimizing memory usage and compare in Section 5 with existing exact matrix multiplication implementations.

Bini's formula

2.1 Algorithm.

We first recall Bini's approximate formula. Let A ∈ K m×k , B ∈ K k×n , and let C = A × B be the matrix product. We consider the special case (m, k, n) = (3, 2, 2) and use a parameter ε. Bini's approximate formula computes a matrix

C ε = A × B + εD(ε), (1)
where D is a polynomial on K 3×2 (ε), with 10 multiplications only. We write a bilinear formula for the multiplication C = A × B:

C = µ r=1 A, X r B, Y r Z r (2)
where •, • is the Frobenius (term-wise) product and X, Y and Z are given 1 in the Table 1. We now write Bini's (3, 2, 2)-algorithm (Table 2) using 2 Table 1 and Equation (2). We use the following notations troughout the rest of the paper: we divide the matrix A of size m × k into six submatrices Aij (or A ij ) of equal size m/3 × k/2, where i and j correspond to row and column indices.

We stress that the formula (2) automatically yields (by duality, [START_REF] Hopcroft | Duality applied to the complexity of matrix multiplication[END_REF]) a (2, 2, 3) formula by simply considering (Y ⊤ r , X ⊤ r , Z ⊤ r ), and a (2, 3, 2) formula using (Z ⊤ r , X r , Y ⊤ r ). Actually, we use (εZ ⊤ r , X r , 1 /εY ⊤ r ) in our algorithm so that the pre-addition phase contain only scalings by ε, and divisions by ε only occur in the post-addition phase, just like in Table 2. Combining those three formulae yields a (12, 12, 12) square matrix algorithm via the (12, 12, 12) = (4

•3, 6•2, 6•2), then (4, 6, 6) = (2 • 2, 2 • 3, 3 • 2) and
(2, 2, 3) variants. This algorithm has an log 12 (1000) ≈ 2.780 exponent (smaller than Strassen's log 2 (7) ≈ 2.807).

We represent in Figure 1 the dependency graph for Bini's (3, 2, 2) algorithm. In this figure, we show scalar division by ε with dotted lines while we draw dashed lines for scalar multiplication by ε. We notice that similarly to Bodrato's algorithm for the Strassen multiplication ([Bod10]), Bini's algorithm has symmetries. These symmetries can be used in the scheduling of the algorithm, taking for instance advantage of independent pre-additions for parallelization, or for possibly saving operations for the squaring in the (12, 12, 12) algorithm.

1 This table corrects some typos in the matrix w (s) r (ε) in [Bin80, eq. (5.2)] (namely signs in Z 4 , Z 7 and Z 9 ).

2 We have actually permuted the indices r = 0, . . . , 9 in Zr for the products Pr.

X0 = 1 0 0 Y0 = ε 0 0 1 Z0 = 1/ε 0 0 1 0 0 X1 = 0 0 0 Y1 = 0 0 -1 -1 Z1 = 1/ε 0 0 0 0 0 X2 = 1 0 0 Y2 = 0 0 0 1 Z2 = -1/ε -1/ε 0 0 0 0 X3 = 0 ε 0 0 Y3 = -ε 0 1 0 Z3 = 1/ε 0 1 0 0 0 X4 = 1 ε 0 0 Y4 = 0 ε 0 1 Z4 = 0 1/ε 0 -1 0 0 X5 = 0 1 0 Y5 = 1 0 0 ε Z5 = 0 0 1 0 0 1/ε X6 = 0 1 0 Y6 = -1 -1 0 0 Z6 = 0 0 0 0 0 1/ε X7 = 0 0 0 Y7 = 1 0 0 0 Z7 = 0 0 0 0 -1/ε -1/ε X8 = 0 1 ε Y8 = 0 1 0 -ε Z8 = 0 0 0 1 0 1/ε X9 = 0 0 ε Y9 = 1 0 ε 0 Z9 = 0 0 -1 0 1/ε 0
S1 ← A11 + A22 T1 ← B22 + ε • B11 T2 ← B21 + B22 S3 ← A32 + ε • A31 T3 ← B11 + ε • B21 S4 ← A22 + ε • A12 T4 ← B21 -ε • B11 S5 ← A11 + ε • A12 T5 ← B22 + ε • B12 S6 ← A21 + A32 T6 ← B11 + ε • B22 T7 ← B11 + B12 S9 ← A21 + ε • A31 T9 ← B12 -ε • B22 B P0 ← A11 × B22 P1 ← S1 × T1 P2 ← A22 × T2 P3 ← S3 × T3 P4 ← S4 × T4 P5 ← S5 × T5 P6 ← S6 × T6 P7 ← A21 × T7 P8 ← A32 × B11 P9 ← S9 × T9 B C11 ← (P1 -P2 + P4 -P0)/ε C12 ← (P5 -P0)/ε C21 ← P4 -P3 + P6 C22 ← P1 -P5 + P9 C31 ← (P3 -P8)/ε C32 ← (P6 -P7 + P9 -P8)/ε 2.2 Complexity.
We give an idea of the number of operations in Strassen-Winograd's and Bini's algorithm for one level of recursion only. On the one hand, we have W (m, k, n) = 7W ( m /2, k /2, n /2) + 4 m /2 k /2 + 4 k /2 n /2 + 7 m /2 n /2, that is 7 multiplications, 4 preadditions on each side, and 7 post-additions, see [START_REF] Winograd | On multiplication of 2x2 matrices[END_REF] and [BDPZ09, Alg. 1]). On the other hand, one has B(m, k, n) = 10W ( m /3, k /2, n /2) + 10 m /3 k /2 + 14 k /2 n /2 + 16 m /3 n /2, that is 10 multiplications, 10 additions or scalings on the left side, 14 additions or scalings on the right, and 16 final additions or scalings by ε.

We write the relative change B(6m,6k,6n)-W (6m,6k,6n) W (6m,6k,6n)

≈ -1 21 . For 6m = 6k = 6n = 3 000, and one level of recursion, one would get ≈ -4.6%. We will not consider two levels of recursions of Bini's algorithm, for reasons that will appear later, so the second level will be a call to Strassen-Winograd's algorithm. In that case, the relative change will be ≈ -4.5%.

Application to exact matrix multiplication

We first recall the transformation in [START_REF] Bini | Relations between exact and approximate bilinear algorithms[END_REF] from an approximate formula (Equation (1)) to and exact algorithm C = A × B: it consists in finding a special set of d + 1 scalars α i , where d = deg ε (εD (ε)), and d + 1 pair-wise distinct scalars ε i . Then it is possible (invertible Vandermonde matrix in ε j i ) to make sure that d+1 i=1 α i = 1 and for j = 1, . . . , d that

d+1 i=1 α i ε j i = 0. Then d+1 i=1 α i C εi breaks into d+1 i=1 α i A×B + d+1 i=1 α i ε i D(ε i ).
The first part reduces to A × B, while the second one gives 0 (after switching the summation operand between

d+1 i=1 and ε i D(ε i ) = d j=1 d j ε j i ).
It is not practicable, at least for the size of matrices we consider, to perform three ε-approximate multiplication while competing against Strassen-Winograd's algorithm implementation, as exemplified in our experiments of Table 6: we get only 5-10% speed-up for one call to the approximate algorithm compared to one call of Strassen-Winograd's. However, we apply a different method requiring only one call to the approximate multiplication, for the special case where the field is Z/pZ.

We have used this approximate algorithm in the following two cases, both for a storage on double floating point machine words. First, we consider the case ε = 2 -27 (Section 3.1), where the idea is to store two exact integers in one double as x + εy. Then any term in ε 2 will be neglected, as ε 2 approaches the machine precision, and a final rounding will remove the ε-approximations. Then in Section 3.2, we take ε = p within the numerical routines and the p-approximations are removed by a final reduction modulo p.

Case

ε = 2 -27 .
We recover the correct result by rounding to the nearest integer, but we need to make sure that the error is not too large and that no exact digit is lost during the approximate computations. We consider one level of recursion.

We need to loose no digit in the computation of the exact part C. In the upper and lower band of of matrix C (i.e. submatrices C1* and C3*) in Table 2, we need that the ε part of the P i 's product is correct. In the central band, we need the integral part to be correct. More precisely, we need to make sure that:

• No digit is lost in the ε part of the computation of P 3 = (A 32 + ε • A 31 ) × (B 11 + ε • B 21 ) (from
examination of all the products one can see that P 3 is the worst case);

• No digit is lost during the post addition phase

• The residual errors vanish when rounded to the nearest.

We represent elements modulo p in the standard representation {0, . . . , p -1}.

Let M = ⌊ k /2⌋(p -1) 2 be the largest element that can occur in a product of A ij by B kl . This is a dot product between ⌊ k /2⌋-long rows and columns of a submatrix of A and B whose entries are bounded by p -1, in the case where k is not even, we would take the floor of k /2 and apply pealing techniques. We then consider the number λ = M + 2M ε + M ε 2 that is the largest possible entry, occurring for instance in P 3 , the number λ 0 = M + 2M ε which is of interest and the error λ 1 = M ε 2 . If M fills 27 bits, then λ 0 will use all the bits in its double representation. Indeed, for M = 2 27 -1, one has λ 0 = 2 27 1 + 2 -27 -2 -53 , which fills all the bits in the mantissa of a double). Then, λ 1 cannot change any bit (the exponent in λ 0 is 27 so 2 -54 is below precision). For any smaller value of M , the error term λ 1 will of course not affect λ 0 . Now, all the digits in P 3 that need to be correct are correct, and so is the middle band of C. The same reasoning applies for the other entries. We now require that ⌊ k /2⌋(p -1) 2 = M < 2 27 . Now, we consider the approximation error. For one recursion level, the components of the error polynomial is either of the form

D 1,1 (ε) = -A 12 × B 11 , or the form D 2,1 (ε) = -εA 12 × B 11 -εA 31 × B 21 + A 12 × B 21 + A 21 × B 22 - A 22 × B 11 -A 31 × B 11 -A 32 × B 21 + A 32 × B 22 .
In fact, the term in ε 2 disappears as explained previously. We examine carefully the coefficient of degree ε in the εD 2,1 expression: the terms A 12 , A 21 and B 22 are only involved additively and the term B 11 only appears with a -. If we want to bound by above, we may assume the first ones are p -1 and the second one is 0. There remains -A 32 B 21 + (A 32 + B 21 )(p -1) + ε(p -1) 2 . This expression is maximum when

A 32 = B 21 = (p -1) /2. So the polynomial D(ε), for ε > 0, is bounded: D(ε) 7 4 ⌊ k /2⌋(p -1) 2 .
The same arguments give a similar bound from below: D(ε) -3⌊ k /2⌋(p -1) 2 For our rounding purposes, we require that εD(ε) ∞ < 1 2 , so we require ⌊ k /2⌋(p -1) 2 = M < 1 6ε = 1 6 2 27 . We have thus proved the following: Proposition 1. For ε = 2 -27 and an (m, k, n) matrix multiplication over a ring Z/pZ, the rounding to the nearest integer of the output C ε in Equation (1) of one call to Bini's (3, 2, 2)-approximate algorithm with double floating point arithmetic, gives the exact result C, provided that:

2⌊ k /2⌋(p -1) 2 < 1 3 2 27 .
Remark 2. These bounds are tight. Indeed, the proof gives explicit matrices A and B, for instance A ⊤ = * * p-1 0 p-1 p-1 and B = p-1 * p-1 0 , or in the balanced representation, if we denote p m = p-1 2 , A ⊤ = * pm pm -pm pm pm and B = pm * pm -pm . Remark 3. If we use a balanced representation (i.e. where an element modulo an odd p is represented between 1-p 2 and p-1 2 ), we can prove that the bound is reduced to D(ε) ∞ 3 2 ⌊ k /2⌋(p -1) 2 , thus gaining in Proposition 1:

⌊ k /2⌋(p -1) 2 < 1 3 2 27 .
3.2 Case ε = p.

We consider again one level of recursion. The elements S i for all i and T i for i = 4, 9 are 0 and 2(p -1) or (p -1) + ε(p -1) = (p -1)(p + 1). Elements T 4 and T 9 are -p(p -1) and p -1. A careful examination of the coefficients in C, before any division by ε, shows that the elements are bounded by ±⌊ k /2⌋(p -1) 2 (p + 1) 2 .

Finally, our algorithm clearly computes an exact multiplication C ε = A × B since the remainder is pD(p) ≡ 0 (mod p). We have then proved: Proposition 4. For ε = p and an (m, k, n) matrix multiplication over a ring Z/pZ, the reduction modulo p of the output C ε in Equation (1) of one call to Bini's (3, 2, 2)-approximate algorithm with double floating point arithmetic, gives the exact result C, provided that:

⌊ k /2⌋(p -1) 2 (p + 1) 2 < 2 53 .
Remark 5. With a balanced representation and an odd prime, a careful examination of the bounds then gives:

1 /2⌊ k /2⌋(p -1) 2 p(p + 1) < 2 53 .

Memory usage and scheduling

In this section, we provide schedules for Bini's approximate multiplication, in a similar fashion to [START_REF] Boyer | Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm[END_REF]. These schedules are then implemented. We try to require as little extra memory (temporaries) as possible.

4.1 Scheduling for the (3, 2, 2) multiplication.

We created the next schedule with only two temporaries (Table 3). In that table, we write e for the parameter ε. A star (*) represents a matrix multiplication while a dot (.) represents a scalar multiplication. The first column represents the operation number, the third column the algorithmic variable (cf. Table 2).

The extra memory used by X is m /3 × k /2 and the extra memory occupied by Y is k /2 × n /2 except for the last but one (step 34) where actually m /3 × n /2 is needed instead. But, at step 34, X is also usable so A 32 B 11 can be stored in a combination of X and Y if k /2 × ( m /3 + n /2) m /3 × n /2. We have proven:

Lemma 1. The extra memory used for one level of recursion in Table 3

is max ( m /3 × n /2, ( m /3 + n /2) × k /2).
This requirement is smaller than Strassen-Winograd's (cf. [START_REF] Boyer | Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm[END_REF]) where X is of size m /2 × max ( k /2, n /2) and Y has size k /2 × n /2. For instance, for m = n = k, we have 5 12 m 2 for Bini's and 1 2 m 2 for Strassen-Winograd's. Lemma 2. Two temporaries are a minimum in Table 3.

algorithm (see [BDPZ09, HLJJ + 96]

). This can also be proven "by hand", considering a finished pebble game on Figure 1 and trying to use only one extra pebble only for the previous moves until no 'previous' move is possible. Remark 6. There is a priori no clear reason whether we should put more effort into reducing the number of operations, the number of allocations or a balance between the two. For instance, in the (2, 3, 2) algorithm, elements A12 and A22, are only involved with a product with ε, and three times each; the same happens for B. One could save operations by creating a temporary eA12 := e . A12 and reuse it three times. Our implementations show that recomputing does not affect negatively the timings.

Schedules for other shapes and other properties

We have implemented schedules for the (2, 2, 3) shape (essentially the same as Table 3 up to exchanging Aij with Bji, and Cij with Cji) and for the (2, 3, 2) shape. In the later case, we created 4 temporaries in our implementation. Just as in [START_REF] Boyer | Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm[END_REF], we can allow the overwriting of (part of) A and/or B. The following Table 4 shows how to perform the multiplication with no extra memory, overwriting a third of matrix A (namely A11 and A12). A similar in-place schedule overwriting B could be written up. The schedule in Table 4 has strong requirements on the sizes of the inputs, for instance, the size of a submatrix of C needs to be equal to one of A (operation # 2 and # 34), larger than one of B (# 3), and the size of a submatrix of A needs to be larger that that of one of B (# 14). In order to relax these requirements, one can write an in-place schedule that overwrites parts of A and B. This phenomenon is common with [START_REF] Boyer | Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm[END_REF].

Finally, we could also imagine schedules for the product with accumulation C ← αA × B + βC that require less memory than the naive approach, and may save operations.

16-36

Same as Table 3, replacing X:= with A11:= and Y:= with A12:=.

Implementation

We implemented in the FFLAS-FFPACK 3 (cf. [START_REF] Dumas | Dense linear algebra over word-size prime fields: the FFLAS and FFPACK packages[END_REF][START_REF] Boyer | Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm[END_REF]) library the algorithm corresponding to the schedule in Table 3. Only one level of recursion is performed, which means that the matrix multiplication within the algorithm are calls to the implementation of Strassen-Winograd's algorithm, that is the fgemm routine of FFLAS.

Achievable size of p.

First of all, according to Propositions 1 and 4, we can take moduli as large as those reported in the Table 5: Then, we want to compare against the implementation of Winograd's algorithm in the FFLAS-FFPACK library. Our implementation uses the double type only, but we need to compete against the float and double implementation of Winograd's algorithm in FFLAS. Indeed, when changed form double to float, our bounds in Propositions 1 and 4 would employ 2 12 instead of 2 27 and (resp.) 2 24 instead of 2 53 , but the bounds linking k and p would be too restrictive. Hence the implementation on double type only. However, the library FFLAS provides a very efficient matrix multiplication implementation over Z/pZ for small p taking advantage of a float representation, since BLAS sgemm routine is close to twice faster than dgemm. Therefore, we compare to both the float and double implementations of FFLAS.

One can find precise bounds linking the prime p, dimension k, the number of recursive calls l, the floating point representation, and the field representation in [START_REF] Dumas | Dense linear algebra over word-size prime fields: the FFLAS and FFPACK packages[END_REF]. For instance, with k = 3 000 and one level of recursion on a float type, one can perform no intermediate modular reduction for p < 39 (full delayed reduction), which is rather small. When p is larger and modular reductions cannot be delayed, FFLAS on float then uses smaller, more numerous, BLAS blocks and performs more reductions. This can still prove to be faster than FFLAS on double, because we get the benefit from faster BLAS on float type. So there is a trade-off between more reductions and smaller blocking with faster BLAS on float and fewer reductions and larger blocking with slower BLAS of double.

Timings.

The timings presented in Table 6 were performed on a x86_64 Gentoo Linux laptop with a 2,3GHz Intel ® Core ™ i7 and 6Gb of RAM, using atlas-3.10.1 BLAS, averaging on four runs (in seconds). The compilation was done with g++-4.8.2 enabling -Ofast and AVX extensions (support introduced in the latest SVN-1.8.0 version of FFLAS-FFPACK). The first column represents the (m, k, n) dimensions of matrices, in thousands; or simply (n) when all dimensions are equal. The entries are created randomly in the ring Z/pZ. The second column is the modulo used for the multiplication (not necessarily prime). The third and fourth column contain our matrix multiplication algorithm from Section 3.1 and (resp.) Section 3.2. The symbol † signifies that the (2, 3, 2) variant was benchmarked. Then, the following two columns are the reference implementation of Winograd algorithm in FFLAS, using the double and float types. Finally, the relative change (in percent) is computed as the ratio (t b -t w )/t w • 100 where t b is this paper implementation's timings and t w if the best of FFLAS reference implementations. The self optimised FFLAS threshold for double was 1 000, while it was 1 640 for float. This threshold essentially determines when Strassen-Winograd's algorithm takes advantage over the naïve (BLAS) algorithm on a given representation (float or double). We did exactly one recursive level for Bini's algorithm (hence no threshold there). The symbol "n/a" (not available) corresponds entries in the table that cannot be filled, i.e. the corresponding algorithm could not be implemented because the modulo p is too large. The columns titled 'M' correspond to a standard modular representation of a field, while 'MB' refers to a balanced representation. The boldface for timings emphasizes the best time for a modular (M) and balanced (MB) representation among all tested algorithms for a given row. matrix multiplication needs to consider a threshold (via automatic benchmarking, cf. [BDG + 14]) below which switch to float representation is preferable (if memory allows). This is not done yet in FFLAS. The balanced representation allows to gain an ≈ 10% speed-up on size 3 900 where the standard representation could not be used. The best speed-up of ≈ 15% around sizes 2 700 to 3300 could be explained by optimal size BLAS block calls.

The ability to adapt to the sizes for the recursion (one may choose to divide m or n or k by 3 instead of 2) makes it possible to always perform better than fgemm; that would not be the case had we not implemented the (2, 3, 2) shape. We notice that timings are almost identical for a constant product mnk, which is not the case for fgemm. Therefore, our "triple cascading" (Bini+Winograd+BLAS) improves on the standard cascading (Wino-grad+BLAS) in a generic fashion: we automatically get better speed-up by pluging in a new algorithm that performs better on some domain of its parameters, and that can plug-in (recursively) the best available routines.

We may also gain speed-ups for larger k and p than those described in Table 5 by first doing a few recursive levels of Winograd's algorithm and finishing by the algorithm discussed in this paper, when the dimension k becomes small enough in a recursive step. This would be a (Winograd+Bini+(Winograd+BLAS)) cascade.

Figure 1 :

 1 Figure 1: Bini's algorithm dependencies

  = e . B11 + B22 T1 11 C21 := X * Y P1 12 C22 := C21 -C22 C22 13 C11 := C21 -C11 C11 14 A12 := e . A12 + A22 S4 15 C32 := B21 -e . B11 T4

Table 1 :

 1 Bini's approximate (3, 2, 2) formula.

Table 2 :

 2 Bini's algorithm

Table 3 :

 3 Schedule for Bini's algorithm using only two temporaries

Table 4 :

 4 Schedule for Bini's algorithm, inplace, overwriting a third of A.

Table 5 :

 5 Largest moduli for delayed reduction in Bini's algorithm

		ε = 2 -27		ε = p
	size Modular Balanced Modular Balanced
	1 000	212	599	2 060	2 450
	2 000	150	424	1 732	2 060
	3 000	123	346	1 565	1 861
	4 000	106	300	1 456	1 732

3 See http://linalg.org/projects/fflas-ffpack/.

Table 6 :

 6 Timings (s) of Bini's algorithm vs.Winograd's (fgemm) on Z/pZ using Modular and ModularBalanced representation. /a n/a 0.31 0.31 0.32 0.32 1.50 1.52 -2.76 -2.98 (2.1) 1 001 n/a n/a 0.81 0.82 0.85 0.85 1.07 1.10 -5.64 -4.60 (2.1) 1 501 n/a n/a 0.80 0.82 0.86 0.86 4.33 4.37 -5.77 -4.62 (2.7) 1 001 n/a n/a 1.70 1.72 2.03 2.02 2.42 2.43 -16.3 -15.0 (2.7) 1 501 n/a n/a 1.71 1.71 2.03 2.03 9.25 9.25 -15.9 -16.0 (3.3) 1 001 n/a n/a 3.08 3.08 3.47 3.47 4.48 4.47 -11.1 -11.0 (3.3) 1 501 n/a n/a 3.09 3.09 3.46 3.45 17.0 16.8 -10.7 -10.6 (3.9) 1 001 n/a n/a 4.94 4.97 5.39 5.54 6.89 6.96 -8.34 -10.5 (3.9) 1 501 n/a n/a n/a 4.94 5.39 5.51 27.8 28.1 n/a -10.4

			Sec. 3.1	Sec. 3.2	FFLAS	FFLAS	rel. change
	dimensions (×1 000)	p	ε = 2 -27 M MB	ε = p M MB	double M MB	float M MB	M	(%)	MB
	(1.5)	141 0.31 0.31 0.31 0.31 0.32 0.32 0.25 0.25 +21.4 +22.2
	(1.5)	451	n/a 0.31 0.31 0.31 0.32 0.32 0.32 0.32 -2.99 -3.28
	(1.5)	1 001 n/a n/a 0.31 0.31 0.32 0.32 0.41 0.43 -2.95 -2.92
	(1.5)	1 501 n

DiscussionThe timings show that our implementation is competitive with Winograd's algorithm implementation, usually providing an ≈ 5% speed-up, and it is always faster than FFLAS using a double representation. As expected, for small primes, the float representation performs better, but this phenomenon is only relevant for small moduli (≈ 400 and less). A general fast efficient implementation of

This material is based on work supported in part by the National Science Foundation under Grant CCF-1115772 (Kaltofen). This material is based on work supported in part by the Agence Nationale pour la Recherche under Grant ANR-11-BS02-013 HPAC (Dumas).