Some Toom-Cook Methods For Even Long Integers

Alberto Zanoni

Centro “Vito Volterra” – Università di Roma “Tor Vergata” Via Columbia 2 – 00133 Roma, Italy

                                         zanoni@volterra.uniroma2.it

Abstract. We present a new approach to evaluation and interpolation phases of some (balanced and unbalanced) Toom-Cook multiplication methods for long integers when at least one factor is even. Some other optimization are also indicated.

AMS Subject Classification: 11A05, 11A25, 11K65, 11Y70

Keywords and phrases: Polynomial, long integers, multiplication, Karatsuba, Toom, interpolation

1 Introduction

We describe an ad hoc approach to evaluation and interpolation phases of some Toom-Cook multiplication methods for long integers, when at least one factor is even. In particular, the classical Toom-3 method (with balanced and unbalanced factors) and the Toom-3.5 method (with slightly and very unbalanced factors) are presented. General possible optimizations of Toom-4.5 and Toom-5 are also shown.

2 Toom-3 classical method

Long integer multiplications is standardly reduced to polynomial multiplications by considering as coefficients a_{i}, b_{i} the digits of a certain base B expansion of the factors (for computer applications, typically B = 2^{32k} for some k). The classical Toom-Cook method [5], [3] – Toom-3 for short – applies when B is such that the obtained polynomials have both degree 2 (balanced case) or 3 and 1, respectively (unbalanced case). Coefficients multiplication is similarly treated by recursion, up to a certain threshold when Karatsuba or high school multiplication methods are more effective.

2.1 Complexity issues

We note that if a_{i}, b_{i} length (number of bits in base 2 expansion) in evaluation phase is n, in interpolation phase the coefficients length is about 2n. In order to analyze complexity, we consider the possible availability of an ad hoc function: consider, for e ∈ ℤ, the following two equivalent processes (I) and (II)

(I) \begin{matrix}T ← 2^{e}X; \\ Z ← Y ± T; \end{matrix}  ;  (II) Z ← Y ± 2^{e}X;

for which one could write a shift-add function sa(X, Y, e, op) ↦ (Y + sign(op)2^{e}X) (with op ∈ {−1, 1}) performing (II) process, which reads X, Y just once and uses no temporary variable, taking benefit of code locality. As memory access is much more time consuming than computing, the less we read/write data, the better it is.

A reasonable computational model to analyze the complexity of interpolation was proposed by Bodrato and Zanoni in [2] – where unbalanced Toom–(n+1/2) methods are also introduced – considering operations costs represented by the constants reported

Operation  Time  Operation  Time
Sign change unitary − ∼ 0  Shift  ≪ , ≫  S
Addition  +  A  Division by k  /  D_{(k)}
Subtraction  −  A  Shift-add  sa A + 1 2

in the aside table, referring to the execution time of additions/subtractions, shifts (multiplications/divisions by power of 2), exact divisions by small constants and sa, whose real values depend on the particular chosen architecture. The indicated costs are relative to the operands (maximum) length: as all operations are linear, costs are proportional to it. This means for example that there is a multiplicative factor 2 distinguishing evaluation and interpolation costs, as operands in interpolation phase have more or less double length with respect to the ones in the evaluation phase.

Obviously we have 1 2 ⩽ S. If sa is not available, process (I) must be used, so that 1 2 = S and one more temporary variable is possibly needed. As the shift operation must read only one operand in memory instead of the two ones needed by additions/subtractions, we reasonably also suppose that S ⩽ A.

2.2 Method description

The balanced version of Toom-3 method considers two quadratic polynomials:

a(x) = a_{2}x^{2} + a_{1}x + a_{0}  ;  b(x) = b_{2}x^{2} + b_{1}x + b_{0}

For the unbalanced version (when a(x) and b(x) have different degrees) we instead have

a(x) = a_{3}x^{3} + a_{2}x^{2} + a_{1}x + a_{0}  ;  b(x) = b_{1}x + b_{0}

To find the coefficients of their product c(x) = a(x)b(x) by using the EMI scheme, we first consider the values w_{i} obtained by evaluating a(x)b(x) in the five interpolation points {∞, 2, −1, 1, 0}, as shown below:

Balanced  Unbalanced
M_{1}  a_{2}b_{2}  a_{3}b_{1}  = w_{4} =  c_{4}
M_{2} (4a_{2} + 2a_{1} + a_{0})(4b_{2} + 2b_{1} + b_{0}) (8a_{3} + 4a_{2} + 2a_{1} + a_{0})(2b_{1} + b_{0}) = w_{3} = 16c_{4} + 8c_{3} + 4c_{2} + 2c_{1} + c_{0}
M_{3} (a_{2} − a_{1} + a_{0})(b_{2} − b_{1} + b_{0})  (a_{3} − a_{2} + a_{1} − a_{0})(b_{1} − b_{0}) = w_{2} = c_{4} − c_{3} + c_{2} − c_{1} + c_{0}
M_{4} (a_{2} + a_{1} + a_{0})(b_{2} + b_{1} + b_{0})  (a_{3} + a_{2} + a_{1} + a_{0})(b_{1} + b_{0}) = w_{1} = c_{4} + c_{3} + c_{2} + c_{1} + c_{0}
M_{5}  a_{0}b_{0}  a_{0}b_{0}  = w_{0} =  c_{0}

In both cases we have

w = Mc  with  M = (\begin{matrix}1 & 0 & 0 & 0 & 0 \\ 16 & 8 & 4 & 2 & 1 \\ 1 & −1 & 1 & −1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \end{matrix} ) ; w = (\begin{matrix}w_{4} \\ w_{3} \\ w_{2} \\ w_{1} \\ w_{0} \end{matrix}) ; c = (\begin{matrix}c_{4} \\ c_{3} \\ c_{2} \\ c_{1} \\ c_{0} \end{matrix})

so that c = M^{−1}w (M_{i} corresponds to the i^{th} line of the matrix).

In [1] Bodrato found the aside shown evaluation procedure for balanced factors in the three more “delicate” points {1, −1, 2} (we show the results u_{i} for a; values v_{i} for factor b are similarly obtained – we have w_{i} = u_{i}v_{i}): the total evaluation cost for both factors is 2(5A + 1 2) = 10A + 2( 1 2).

\begin{aligned}& & & & a_{2} & a_{1} & a_{0} \\ 1) & u_{2} & = a_{2} + a_{0} & [ & 1 & 0 & 1 ] & A \\ 2) & u_{1} & = u_{2} + a_{1} & [ & 1 & 1 & 1 ] & A \\ 3) & u_{2} & = u_{2} − a_{1} & [ & 1 & −1 & 1 ] & A \\ 4) & u_{3} & = u_{1} + a_{2} & [ & 2 & 1 & 1 ] & A \\ 5) & u_{3} & = (u_{3} ≪1) − a_{0} & [ & 4 & 2 & 1 ] & A + 1 2 \end{aligned} (1)

The unbalanced version is instead realized as follows, with a slightly bigger cost: [7A+3( 1 2)]+3A = 10A+3( 1 2)

\begin{aligned}& & & & a_{3} & a_{2} & a_{1} & a_{0} \\ 1) & u_{3} & = a_{2} + a_{0} & [ & 0 & 1 & 0 & 1 & ] A \\ 2) & u_{2} & = a_{3} + a_{1} & [ & 1 & 0 & 1 & 0 & ] A \\ 3) & u_{1} & = u_{2} + u_{3} & [ & 1 & 1 & 1 & 1 & ] A \\ 4) & u_{2} & = u_{2} − u_{3} & [ & 1 & −1 & 1 & −1 & ] A \\ 5) & u_{3} & = a_{2} + (a_{3} ≪1) & [ & 2 & 1 & 0 & 0 & ] A + 1 2 \\ 6) & u_{3} & = a_{1} + (u_{3} ≪1) & [ & 4 & 2 & 1 & 0 & ] A + 1 2 \\ 7) & u_{3} & = a_{0} + (u_{3} ≪1) & [ & 8 & 4 & 2 & 1 & ] A + 1 2 \end{aligned}  \begin{aligned}& & & & b_{1} & b_{0} \\ 8) & v_{1} & = b_{1} + b_{0} & [ & 1 & 1 & ] A \\ 9) & v_{2} & = b_{1} − b_{0} & [ & 1 & −1 & ] A \\ 10) & v_{3} & = v_{1} + b_{1} & [ & 2 & 1 & ] A \end{aligned}  (2)

Zimmermann, in GMP library [4] version 4.2.1, proposed the following sequence of operations (inversion sequence, or IS for short) to invert M.

\begin{matrix}M = & (\begin{matrix}1 & 0 & 0 & 0 & 0 \\ 16 & 8 & 4 & 2 & 1 \\ 1 & −1 & 1 & −1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹^{M_{2}+(2)M_{3}} (\begin{matrix}1 & 0 & 0 & 0 & 0 \\ 18 & 6 & 6 & 0 & 3 \\ 1 & −1 & 1 & −1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹_{M_{3}+M_{4}}^{M_{2}/(3)} ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 6 & 2 & 2 & 0 & 1 \\ 2 & 0 & 2 & 0 & 2 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \end{matrix} ) \\ ⟹_{M_{3}≫(1)}^{M_{2}+M_{5}} & ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 6 & 2 & 2 & 0 & 2 \\ 1 & 0 & 1 & 0 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹_{M_{2}−(2)M_{1}}^{M_{2}≫(1)} ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹^{M_{4}−M_{2}} (\begin{matrix}1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{matrix}) ⟹^{M_{2}−M_{3}} (\begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{matrix}) ⟹_{M_{3}−M_{5}}^{M_{3}−M_{1}} I \end{matrix}

with computational cost

cost_{GMP} = 8A + D_{(3)} + 2S + 2( 1 2)

More efficient inversion sequences have already been found by exhaustive search in [2] with smaller cost:

cost_{BZ} = 8A + D_{(3)} + 2S + ( 1 2)

3 Toom-3 with (at least) an even factor

We consider here the case in which at least one of a_{0}, b_{0} is even, so that c_{0} = a_{0}b_{0} is even, too. This happens in 75 % of all possible cases, and its detection is quite fast (just test the least meaningful bit of a_{0} and possibly b_{0}). We will distinguish the balanced and unbalanced versions.

3.1 Evaluation phase: the balanced case

Factors have here the same degree. To fix ideas, without loss of generality let’s suppose a_{0} is even. Then, from

M_{2} : (4a_{2} + 2a_{1} + a_{0})(4b_{2} + 2b_{1} + b_{0}) = 16c_{4} + 8c_{3} + 4c_{2} + 2c_{1} + c_{0}

dividing by 2 both sides we still obtain integer numbers

M′_{2} : (2a_{2} + a_{1} + \frac{a_{0}}{2} )(4b_{2} + 2b_{1} + b_{0}) = 8c_{4} + 4c_{3} + 2c_{2} + c_{1} + \frac{c_{0}}{2}

The new evaluation sequence (ES, for short) for factor a is shown aside. It has exactly the same computational complexity of (1): note infact that we do not have to explicitly divide by 2. This will help us in the interpolation phase, permitting to save one shift.

\begin{aligned}& & & & a_{2} & a_{1} & a_{0} \\ 1) & u_{2} & = a_{2} + a_{0} & [ & 1 & 0 & 1 & ] A \\ 2) & u_{1} & = u_{2} + a_{1} & [ & 1 & 1 & 1 & ] A \\ 3) & u_{2} & = u_{2} − a_{1} & [ & 1 & −1 & 1 & ] A \\ 4) & u_{3} & = u_{2} + a_{2} & [ & 2 & 1 & 1 & ] A \\ 5) & u_{3} & = u_{3} − (a_{0} ≫1) & [ & 2 & 1 & \frac{1}{2} & ] A + 1 2 \end{aligned}

3.2 Evaluation phase: the unbalanced case

In the unbalanced case we have deg(a) = 3 and deg(b) = 1. We have two asymmetrical subcases:

– When a_{0} is even we have

M′_{2} : (4a_{3} + 2a_{2} + a_{1} + \frac{a_{0}}{2} )(2b_{1} + b_{0}) = 8c_{4} + 4c_{3} + 2c_{2} + c_{1} + \frac{c_{0}}{2}

and the new evaluation of the first factor is slightly different but not worse than before: the whole evaluation cost does therefore not change with respect to the balanced case.

a_{3} a_{2} a_{1} a_{0}
1′) u_{3} = a_{2} + a_{0}  [ 0 1 0 1 ] A
2′) u_{2} = a_{3} + a_{1}  [ 1 0 1 0 ] A
3′) u_{1} = u_{2} + u_{3}  [ 1 1 1 1 ] A
4′) u_{2} = u_{2} − u_{3}  [ 1 −1 1 −1 ] A
5′) u_{3} = a_{2} + (a_{3} ≪1) [ 2 1 0 0 ] A + 1 2
6′) u_{3} = a_{1} + (u_{3} ≪1) [ 4 2 1 0 ] A + 1 2
7′) u_{3} = u_{3} + (a_{0} ≫1) [ 4 2 1 \frac{1}{2} ] A + 1 2

– When b_{0} is even, the situation gets unfortunately a bit worse:

M_{2} : (8a_{3} + 4a_{2} + 2a_{1} + a_{0})(b_{1} + \frac{b_{0}}{2} ) = 8c_{4} + 4c_{3} + 2c_{2} + c_{1} + \frac{c_{0}}{2}

The ES cost for the second factor grows then to 3A + ( 1 2), so that the whole evaluation cost is 10A + 4( 1 2).

b_{1} b_{0}
9′) v_{1} = b_{1} + b_{0}  [ 1 1 ] A
10′) v_{2} = b_{1} − b_{0}  [ 1 −1 ] A
11′) v_{3} = b_{1} + (b_{0} ≫1) [ 1 \frac{1}{2} ] A + 1 2

Obviously, if both a_{0} and b_{0} are even, it is therefore preferrable to aproach the first subcase.

3.3 Interpolation phase

We can therefore consider a different matrix M′ to be inverted, with the second line divided by 2.

M′ = ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 8 & 4 & 2 & 1 & \frac{1}{2} \\ 1 & −1 & 1 & −1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \end{matrix} )

We propose for it the following inversion sequence, someway inspired by Zimmermann’s:

M′ ⟹^{M′_{2}+M′_{3}} ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 9 & 3 & 3 & 0 & \frac{3}{2} \\ 1 & −1 & 1 & −1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹_{M′_{3}+M′_{4}}^{M′_{2}/(3)} (\begin{matrix}1 & 0 & 0 & 0 & 0 \\ 3 & 1 & 1 & 0 & \frac{1}{2} \\ 2 & 0 & 2 & 0 & 2 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \end{matrix}) ⟹_{M′_{2}−(2)M′_{1}}^{M′_{2}+( \frac{1}{2} )M′_{5}} (\begin{matrix}1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 1 \\ 2 & 0 & 2 & 0 & 2 \\ 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \end{matrix}) ⟹_{M′_{3}≫(1)}^{M′_{4}−M′_{2}} (\begin{matrix}1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{matrix}) ⟹^{M′_{2}−M′_{3}}(\begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{matrix}) ⟹_{M′_{3}−M′_{5}}^{M′_{3}−M′_{1}} I

Its computational cost is 8A + D_{(3)} + S + 2( 1 2), smaller than cost_{GMP}. We could here avoid a shift because we implicitly did it in the evaluation phase, but with no (small, when only b_{0} is even in the unbalanced case) extra cost.

We report below an implementation in gp-pari: the three possible ES’s and the common IS.

\\ Evaluation: Unbalanced case (a even).               \\ Evaluation: Unbalanced case (b even).
a = a3*x^3 + a2*x^2 + a1*x + a0; b = b1*x + b0;        a = a3*x^3 + a2*x^2 + a1*x + a0; b = b1*x + b0;
w0 = a2 + a0;            w4 = b1 - b0;                 w0 = a2 + a0;            w4 = b1 - b0;
w1 = a3 + a1;                                          w1 = a3 + a1;
w3 = w1 - w0;                                          w3 = w1 - w0;
w2 = w4*w3; \\ Evaluation in (-1)                      w2 = w4*w3; \\ Evaluation in (-1)
w3 = w1 + w0;            w4 = b1 + b0;                 w3 = w1 + w0;            w4 = b1 + b0;
w1 = w3*w4; \\ Evaluation in (1)                       w1 = w3*w4; \\ Evaluation in (1)
w0 = a2 + (a3<<1);       w4 = w4 + b1;                 w0 = a2 + (a3<<1);       w4 = b1 + (b0>>1);
w0 = a1 + (w0<<1);                                     w0 = a1 + (w0<<1);
w0 = w0 + (a0>>1);                                     w0 = a0 + (w0<<1);
w3 = w0*w4; \\ Evaluation in (2) divided by 2.         w3 = w0*w4; \\ Evaluation in (2) divided by 2.
w0 = a0*b0; \\ Evaluation in (0)                       w0 = a0*b0; \\ Evaluation in (0)
w4 = a3*b1; \\ Evaluation in (1/0)                     w4 = a3*b1; \\ Evaluation in (1/0)
\\ Evaluation: Balanced case (a even)                  \\ Interpolation
a = a2*x^2 + a1*x + a0; b = b2*x^2 + b1*x + b0;        w3 = w3 + w2;       \\ A        (9 3 3 0 3/2)
                                                       w3 = w3 / 3;        \\ D        (3 1 1 0 1/2)
w0 = a2 + a0;            w4 = b2 + b0;                 w2 = w2 + w1;       \\ A        (2 0 2 0 2)
w1 = w0 - a1;            w3 = w4 - b1;                 w3 = w3 + (w0>>1);  \\ A + _1_2 (3 1 1 0 1)
w2 = w1*w3; \\ Evaluation in (-1)                      w3 = w3 - (w4<<1);  \\ A + _1_2 (1 1 1 0 1)
w0 = w0 + a1;            w4 = w4 + b1;                 w1 = w1 - w3;       \\ A        (0 0 0 1 0)
w1 = w0*w4; \\ Evaluation in (1)                       w2 = w2 >> 1;       \\ S        (1 0 1 0 1)
w0 = w0 + a2;            w4 = w4 + b2;                 w3 = w3 - w2;       \\ A        (0 1 0 0 0)
w0 = w0 - (a0>>1);       w4 = (w4<<1) - b0;            w2 = w2 - w0;       \\ A        (1 0 1 0 0)
w3 = w0*w4; \\ Evaluation in (2) divided by 2.         w2 = w2 - w4;       \\ A        (0 0 1 0 0)
w0 = a0*b0; \\ Evaluation in (0)
w4 = a2*b2; \\ Evaluation in (1/0)                     c = w4*x^4 + w3*x^3 + w2*x^2 + w1*x + w0;

We point out that one could equivalently use −2 as interpolation value instead of 2. This reduces the probability of carry in M_{2} computation, but one has then to cope with more negative values.

4 Toom-3.5 with (at least) an even factor

When deg(c) = deg(a) + deg(b) = 5 we may apply the Toom-3.5 method, which is intrinsically unbalanced. The used interpolation values are {∞, −2, 2, 1, −1, 0}. There are two versions: the slightly and the very unbalanced one, and for each of them we have to consider two cases, depending on which among a_{0}, b_{0} is even.

4.1 Evaluation phase: the slightly unbalanced case

When a(x) = a_{3}x^{3} + a_{2}x^{2} + a_{1}x + a_{0} and b(x) = b_{2}x^{2} + b_{1}x + b_{0}, if sa is available, the ES is reported below: the cost for a is 8A + 4( 1 2), while the cost for b is 6A + 2( 1 2). The total cost is 14A + 6( 1 2).

\begin{aligned}& & & & a_{3} & a_{2} & a_{1} & a_{0} \\ 1) & u_{2} & = a_{3} + a_{1} & [ & 1 & 0 & 1 & 0 & ] A \\ 2) & u_{3} & = a_{2} + a_{0} & [ & 0 & 1 & 0 & 1 & ] A \\ 3) & u_{4} & = u_{2} − u_{3} & [ & −1 & 1 & −1 & 1 & ] A \\ 4) & u_{5} & = u_{2} + u_{3} & [ & 1 & 1 & 1 & 1 & ] A \\ 5) & u_{1} & = a_{0} + (a_{2} ≪2) & [ & 0 & 4 & 0 & 1 & ] A + 1 2 \\ 6) & u_{3} & = a_{1} + (a_{3} ≪2) & [ & 4 & 0 & 1 & 0 & ] A + 1 2 \\ 7) & u_{2} & = u_{1} − (u_{3} ≪1) & [ & −8 & 4 & −2 & 1 & ] A + 1 2 \\ 8) & u_{3} & = u_{1} + (u_{3} ≪1) & [ & 8 & 4 & 2 & 1 & ] A + 1 2 \end{aligned}  \begin{aligned}& & & & b_{2} & b_{1} & b_{0} \\ 9) & v_{5} & = b_{2} + b_{0} & [ & 1 & 0 & 1 & ] A \\ 10) & v_{4} & = v_{5} − b_{1} & [ & 1 & −1 & 1 & ] A \\ 11) & v_{5} & = v_{5} + b_{1} & [ & 1 & 1 & 1 & ] A \\ 12) & v_{2} & = v_{5} + b_{2} & [ & 2 & 1 & 1 & ] A \\ 13) & v_{3} & = (v_{2} ≪1) − b_{0} & [ & 4 & 2 & 1 & ] A + 1 2 \\ 14) & v_{2} & = v_{3} − (b_{1} ≪2) & [ & 4 & −2 & 1 & ] A + 1 2 \end{aligned} (3)

If sa is not available, the ES is reported below: the cost for a is 8A+3S, while the cost for b is 6A+2S (the proposed ES for b is different from the one that could straightforwardly be obtained from the above one, in order to reduce carry presence probability). The total cost is 14A + 5S.

\begin{aligned}& & a_{3} & a_{2} & a_{1} & a_{0} \\ 1) - 4) as in eq. (3 ) \\ 5) & u_{1} & = a_{2} ≪2 & [ & 0 & 4 & 0 & 0 & ] S \\ 6) & u_{1} & = u_{1} + a_{0} & [ & 0 & 4 & 0 & 1 & ] A \\ 7) & u_{3} & = a_{3} ≪2 & [ & 4 & 0 & 0 & 0 & ] S \\ 8) & u_{3} & = u_{3} + a_{1} & [ & 4 & 0 & 1 & 0 & ] A \\ 9) & u_{3} & = u_{3} ≪1 & [ & 8 & 0 & 2 & 0 & ] S \\ 10) & u_{2} & = u_{1} − u_{3} & [ & −8 & 4 & −2 & 1 & ] A \\ 11) & u_{3} & = u_{1} + u_{3} & [ & 8 & 4 & 2 & 1 & ] A \end{aligned}  \begin{aligned}& & & & b_{2} & b_{1} & b_{0} \\ 12) & v_{5} & = b_{2} + b_{0} & [ & 1 & 0 & 1 & ] A \\ 13) & v_{4} & = v_{5} − b_{1} & [ & 1 & −1 & 1 & ] A \\ 14) & v_{5} & = v_{5} + b_{1} & [ & 1 & 1 & 1 & ] A \\ 15) & v_{2} & = b_{2} ≪2 & [ & 4 & 0 & 0 & ] S \\ 16) & v_{2} & = v_{2} + b_{0} & [ & 4 & 0 & 1 & ] A \\ 17) & v_{1} & = b_{1} ≪1 & [ & 0 & 2 & 0 & ] S \\ 18) & v_{3} & = v_{2} + v_{1} & [ & 4 & 2 & 1 & ] A \\ 19) & v_{2} & = v_{2} − v_{1} & [ & 4 & −2 & 1 & ] A \end{aligned}  (4)

– When a_{0} is even, dividing by 2 for the interpolation value x = 2 (and similarly for x = −2) we have

M′_{2} : (4a_{3} + 2a_{2} + a_{1} + \frac{a_{0}}{2} )(4b_{2} + 2b_{1} + b_{0}) = 16c_{5} + 8c_{4} + 4c_{3} + 2c_{2} + c_{1} + \frac{c_{0}}{2}

The ES for a changes, but the cost does not. We show it in both cases, when sa is and is not available, respectively.

\begin{aligned}& & a_{3} & a_{2} & a_{1} & a_{0} \\ 1′) - 4′) as 1) - 4) in eq. (3 ) \\ 5′) & u_{1} & = a_{0} + (a_{2} ≪2) & [ & 0 & 4 & 0 & 1 & ] & A + 1 2 \\ 6′) & u_{3} & = a_{1} + (a_{3} ≪2) & [ & 4 & 0 & 1 & 0 & ] & A + 1 2 \\ 7′) & u_{2} & = (u_{1} ≫1) − u_{3} & [ & −4 & 2 & −1 & \frac{1}{2} & ] & A + 1 2 \\ 8′) & u_{3} & = (u_{1} ≫1) + u_{3} & [ & 4 & 2 & 1 & \frac{1}{2} & ] & A + 1 2 \end{aligned}  \begin{aligned}a_{3} & a_{2} & a_{1} & a_{0} \\ 1′) - 4′) as 1) - 4) in eq. (4) \\ 5′) & u_{1} & = a_{2} ≪1 & [ & 0 & 2 & 0 & 0 & ] & S \\ 6′) & u_{2} & = a_{0} ≫1 & [ & 0 & 0 & 0 & \frac{1}{2} & ] & S \\ 7′) & u_{1} & = u_{1} + u_{2} & [ & 0 & 2 & 0 & \frac{1}{2} & ] & A \\ 8′) & u_{3} & = a_{3} ≪2 & [ & 4 & 0 & 0 & 0 & ] & S \\ 9′) & u_{3} & = u_{3} + a_{1} & [ & 4 & 0 & 1 & 0 & ] & A \\ 10′) & u_{2} & = u_{1} − u_{3} & [ & −4 & 2 & −1 & \frac{1}{2} & ] & A \\ 11′) & u_{3} & = u_{1} + u_{3} & [ & 4 & 2 & 1 & \frac{1}{2} & ] & A \end{aligned}  (5)

– When b_{0} is even, dividing by 2 for the interpolation value x = 2 (and similarly for x = −2) we instead have

M′_{2} : (8a_{3} + 4a_{2} + 2a_{1} + a_{0})(2b_{2} + b_{1} + \frac{b_{0}}{2} ) = 16c_{5} + 8c_{4} + 4c_{3} + 2c_{2} + c_{1} + \frac{c_{0}}{2}

The ES for b changes, but the cost does not. We show both cases, when sa is and is not available, respectively.

\begin{aligned}& & & & b_{2} & b_{1} & b_{0} \\ 9′) & v_{5} & = b_{2} + b_{0} & [ & 1 & 0 & 1 & ] A \\ 10′) & v_{4} & = v_{2} − b_{1} & [ & 1 & −1 & 1 & ] A \\ 11′) & v_{5} & = v_{2} + b_{1} & [ & 1 & 1 & 1 & ] A \\ 12′) & v_{1} & = v_{5} + b_{2} & [ & 2 & 1 & 1 & ] A \\ 13′) & v_{1} & = v_{5} − (b_{0} ≫1) & [ & 2 & 1 & \frac{1}{2} & ] A + 1 2 \\ 14′) & v_{2} & = v_{1} − (b_{1} ≪1) & [ & 2 & −1 & \frac{1}{2} & ] A + 1 2 \end{aligned}  \begin{aligned}& & & & b_{2} & b_{1} & b_{0} \\ 12′) & v_{5} & = b_{2} + b_{0} & [ & 1 & 0 & 1 & ] & A \\ 13′) & v_{4} & = v_{5} − b_{1} & [ & 1 & −1 & 1 & ] & A \\ 14′) & v_{5} & = v_{5} + b_{1} & [ & 1 & 1 & 1 & ] & A \\ 15′) & v_{1} & = b_{2} ≪1 & [ & 2 & 0 & 0 & ] & S \\ 16′) & v_{2} & = b_{0} ≫1 & [ & 0 & 0 & \frac{1}{2} & ] & S \\ 17′) & v_{1} & = v_{1} + v_{2} & [ & 2 & 0 & \frac{1}{2} & ] & A \\ 18′) & v_{3} & = v_{2} + b_{1} & [ & 2 & 1 & \frac{1}{2} & ] & A \\ 19′) & v_{2} & = v_{2} − b_{1} & [ & 2 & −1 & \frac{1}{2} & ] & A \end{aligned}  (6)

4.2 Evaluation phase: the very unbalanced case

When a(x) = a_{4}x^{4} + a_{3}x^{3} + a_{2}x^{2} + a_{1}x + a_{0}, b(x) = b_{1}x + b_{0}, if sa is available, the ES is reported below: the cost for a is 10A + 5( 1 2), for b is 4A. The total cost is 14A + 5( 1 2).

\begin{aligned}& & & & a_{4} & a_{3} & a_{2} & a_{1} & a_{0} \\ 1) & u_{2} & = a_{3} + a_{1} & [ & 0 & 1 & 0 & 1 & 0 & ] A \\ 2) & u_{1} & = a_{4} + a_{0} & [ & 1 & 0 & 0 & 0 & 1 & ] A \\ 3) & u_{1} & = u_{1} + a_{2} & [ & 1 & 0 & 1 & 0 & 1 & ] A \\ 4) & u_{4} & = u_{1} − u_{2} & [ & 1 & −1 & 1 & −1 & 1 & ] A \\ 5) & u_{5} & = u_{1} + u_{2} & [ & 1 & 1 & 1 & 1 & 1 & ] A \\ 6) & u_{1} & = a_{1} + (a_{3} ≪2) & [ & 0 & 4 & 0 & 1 & 0 & ] A + 1 2 \\ 7) & u_{2} & = a_{2} + (a_{4} ≪2) & [ & 4 & 0 & 1 & 0 & 0 & ] A + 1 2 \\ 8) & u_{3} & = a_{0} + (u_{2} ≪2) & [ & 16 & 0 & 4 & 0 & 1 & ] A + 1 2 \\ 9) & u_{2} & = u_{3} − (u_{1} ≪1) & [ & 16 & −8 & 4 & −2 & 1 & ] A + 1 2 \\ 10) & u_{3} & = u_{3} + (u_{1} ≪1) & [ & 16 & 8 & 4 & 2 & 1 & ] A + 1 2 \end{aligned}  \begin{aligned}& & & & b_{1} & b_{0} \\ 11) & v_{4} & = b_{0} − b_{1} & [ & −1 & 1 & ] A \\ 12) & v_{5} & = b_{0} + b_{1} & [ & 1 & 1 & ] A \\ 13) & v_{2} & = v_{4} − b_{1} & [ & −2 & 1 & ] A \\ 14) & v_{3} & = v_{5} + b_{1} & [ & 2 & 1 & ] A \end{aligned}  (7)

If sa is not available, the ES for a is reported below (the ES of b does not change): its cost is 10A + 4S. The total cost is 14A + 4S.

\begin{aligned}& & a_{4} & a_{3} & a_{2} & a_{1} & a_{0} \\ 1) - 5) as in eq. (7 ) \\ 6) & u_{1} & = a_{3} ≪2 & [ & 0 & 4 & 0 & 0 & 0 & ] S \\ 7) & u_{1} & = u_{1} + a_{1} & [ & 0 & 4 & 0 & 1 & 0 & ] A \\ 8) & u_{1} & = u_{1} ≪1 & [ & 0 & 8 & 0 & 2 & 0 & ] S \\ 9) & u_{2} & = a_{4} ≪2 & [ & 4 & 0 & 0 & 0 & 0 & ] S \end{aligned}  \begin{aligned}& & & & a_{4} & a_{3} & a_{2} & a_{1} & a_{0} \\ 10) & u_{2} & = u_{2} + a_{2} & [ & 4 & 0 & 1 & 0 & 0 & ] A \\ 11) & u_{3} & = u_{2} ≪2 & [ & 16 & 0 & 4 & 0 & 0 & ] S \\ 12) & u_{3} & = u_{3} + a_{0} & [ & 16 & 0 & 4 & 0 & 1 & ] A \\ 13) & u_{2} & = u_{3} − u_{1} & [ & 16 & −8 & 4 & −2 & 1 & ] A \\ 14) & u_{3} & = u_{3} + u_{1} & [ & 16 & 8 & 4 & 2 & 1 & ] A \end{aligned}  (8)

– When a_{0} is even, dividing by 2 for the interpolation value x = 2 (and similarly for x = −2) we have

M′_{2} : (8a_{4} + 4a_{3} + 2a_{2} + a_{1} + \frac{a_{0}}{2} )(2b_{1} + b_{0}) = 16c_{5} + 8c_{4} + 4c_{3} + 2c_{2} + c_{1} + \frac{c_{0}}{2}

The ES for a changes, but the cost does not. We show both cases, when sa is and is not available, respectively. The first 5 steps are as in equation (7).

\begin{aligned}& & & & a_{4} & a_{3} & a_{2} & a_{1} & a_{0} \\ 6′) & u_{1} & = a_{1} + (a_{3} ≪2) & [ & 0 & 4 & 0 & 1 & 0 & ] A + 1 2 \\ 7′) & u_{2} & = a_{2} + (a_{4} ≪2) & [ & 4 & 0 & 1 & 0 & 0 & ] A + 1 2 \\ 8′) & u_{3} & = a_{0} + (u_{2} ≪2) & [ & 16 & 0 & 4 & 0 & 1 & ] A + 1 2 \\ 9′) & u_{2} & = (u_{3} ≫1) − u_{1} & [ & 8 & −4 & 2 & −1 & \frac{1}{2} & ] A + 1 2 \\ 10′) & u_{3} & = (u_{3} ≫1) + u_{1} & [ & 8 & 4 & 2 & 1 & \frac{1}{2} & ] A + 1 2 \end{aligned}  \begin{aligned}& & & & a_{4} & a_{3} & a_{2} & a_{1} & a_{0} \\ 6′) & u_{1} & = a_{3} ≪2 & [ & 0 & 4 & 0 & 0 & 0 & ] & S \\ 7′) & u_{1} & = u_{1} + a_{1} & [ & 0 & 4 & 0 & 1 & 0 & ] & A \\ 8′) & u_{2} & = a_{4} ≪3 & [ & 8 & 0 & 0 & 0 & 0 & ] & S \\ 9′) & u_{3} & = a_{2} ≪1 & [ & 0 & 0 & 2 & 0 & 0 & ] & S \\ 10′) & u_{2} & = u_{2} + u_{3} & [ & 8 & 0 & 2 & 0 & 0 & ] & A \\ 11′) & u_{3} & = a_{0} ≫1 & [ & 0 & 0 & 0 & 0 & \frac{1}{2} & ] & S \\ 12′) & u_{3} & = u_{2} + u_{3} & [ & 8 & 0 & 2 & 0 & \frac{1}{2} & ] & A \\ 13′) & u_{2} & = u_{3} − u_{1} & [ & 8 & −4 & 2 & −1 & \frac{1}{2} & ] & A \\ 14′) & u_{3} & = u_{3} + u_{1} & [ & 8 & 4 & 2 & 1 & \frac{1}{2} & ] & A \end{aligned}  (9)

– When b_{0} is even, dividing by 2 for the interpolation value x = 2 (and similarly for x = −2) we instead have

M′_{2} : (16a_{4} + 8a_{3} + 4a_{2} + 2a_{1} + a_{0})(b_{1} + \frac{b_{0}}{2} ) = 16c_{5} + 8c_{4} + 4c_{3} + 2c_{2} + c_{1} + \frac{c_{0}}{2}

In this case the situation gets unfortunately a bit worse: The new (partial) ES cost for the second factor grows to 4A + ( 1 2) is sa is available, or 4A + S if it is not, respectively.

\begin{aligned}& & & & b_{1} & b_{0} \\ 12′) & v_{4} & = b_{0} − b_{1} & [ & −1 & 1 & ] A \\ 13′) & v_{5} & = b_{0} + b_{1} & [ & 1 & 1 & ] A \\ 14′) & v_{2} & = (b_{0} ≫1) − b_{1} & [ & −1 & \frac{1}{2} & ] A + 1 2 \\ 15′) & v_{3} & = b_{0} − v_{2} & [ & 1 & \frac{1}{2} & ] A \end{aligned}  \begin{aligned}& & & & b_{1} & b_{0} \\ 12′) & v_{4} & = b_{0} ≫1 & [ & 0 & \frac{1}{2} & ] & S \\ 13′) & v_{2} & = v_{4} − b_{1} & [ & −1 & \frac{1}{2} & ] & A \\ 14′) & v_{3} & = v_{4} + b_{1} & [ & 1 & \frac{1}{2} & ] & A \\ 15′) & v_{4} & = b_{0} − b_{1} & [ & −1 & 1 & ] & A \\ 16′) & v_{5} & = b_{0} + b_{1} & [ & 1 & 1 & ] & A \end{aligned}  (10)

The total ES cost amounts then to 14A + 6( 1 2) and 14A + 5S, respectively.

4.3 Interpolation phase

The IS proposed in [1] for the general case has a cost of 12A + 2S + D_{(6)} + D_{(12)} + 2( 1 2). Strictly following the EMI scheme, we can manage the IS so that the cost becomes one of the following – when sa is available or not, respectively. Note that at least one division is now by a different constant.^{1}

cost′ = 12A + S + D_{(3)} + D_{(12)} + 3( 1 2)  ;  cost′′ = 12A + 4S + 2D_{(3)}

The new M matrix resulting from the evaluation values {∞, −2, 2, −1, 1, 0} with second and third line divided by 2 is

M = (\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 \\ −16 & 8 & −4 & 2 & −1 & \frac{1}{2} \\ 16 & 8 & 4 & 2 & 1 & \frac{1}{2} \\ −1 & 1 & −1 & 1 & −1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} )

If sa is available, the proposed IS is:

\begin{matrix}M ⟹_{M_{4}+M_{5}}^{M_{2}+M_{3}} & (\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 16 & 0 & 4 & 0 & 1 \\ 16 & 8 & 4 & 2 & 1 & \frac{1}{2} \\ 0 & 2 & 0 & 2 & 0 & 2 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹_{M_{4}≫1}^{M_{3}−( \frac{1}{2} )M_{2}}(\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 16 & 0 & 4 & 0 & 1 \\ 16 & 0 & 4 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 & 0 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹_{M_{5}−M_{4}}^{M_{2}−M_{6}}(\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 16 & 0 & 4 & 0 & 0 \\ 16 & 0 & 4 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹_{M_{4}−M_{6}}^{M_{3}−M_{5}}(\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 16 & 0 & 4 & 0 & 0 \\ 15 & 0 & 3 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} ) \\ ⟹_{M_{3}/(3)}^{M_{2}−(4)M_{4}} & ( \begin{matrix}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 12 & 0 & 0 & 0 & 0 \\ 5 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹_{M_{3}−(4)M_{1}}^{M_{2}/(12)}( \begin{matrix}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹_{M_{5}−M_{3}}^{M_{4}−M_{2}}( \begin{matrix}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹^{M_{3}−M_{1}} I \end{matrix}

If sa is not available, the IS is slightly different, and a small “trick” has to be used in order not to use any extra temporary variable: evaluation and interpolation phases have to be interlaced someway. We don’t report explicitly this second IS: it can be easily deduced from the code reported in appendix A or B.

If we instead mix a bit IS and ES we can also work as follows: rewrite ES so that the interpolation matrix is the shown aside M: that is, simply reorganize the ES in order to put the evaluation in −2 (divided by 2) in M_{6}, where c_{0} = a_{0}b_{0} – the evaluation in 0 – should be, and do not compute c_{0}. Note that M_{2} line remains for now undefined. The IS (when sa is available) becomes then the following one, otherwise the above considerations apply:

M = (\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 \\ ⟨not computed value⟩ \\ 16 & 8 & 4 & 2 & 1 & \frac{1}{2} \\ −1 & 1 & −1 & 1 & −1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ −16 & 8 & −4 & 2 & −1 & \frac{1}{2} \end{matrix} )

M ⟹_{M_{5}−M_{4}}^{M_{2}=M_{3}+M_{6}} (\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 16 & 0 & 4 & 0 & 1 \\ 16 & 8 & 4 & 2 & 1 & \frac{1}{2} \\ −1 & 1 & −1 & 1 & −1 & 1 \\ 2 & 0 & 2 & 0 & 2 & 0 \\ −16 & 8 & −4 & 2 & −1 & \frac{1}{2} \end{matrix} ) ⟹_{M_{6}=a_{0}b_{0}}^{M_{3}−M_{6}}(\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 16 & 0 & 4 & 0 & 1 \\ 32 & 0 & 8 & 0 & 2 & 0 \\ −1 & 1 & −1 & 1 & −1 & 1 \\ 2 & 0 & 2 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} ) ⟹ ⋯ ⟨as in [2]⟩ ⋯ ⟹ I

The cost becomes one of the following – when sa is available or not, respectively. Note that at least one division is now by a different constant.

cost′ = 12A + S + D_{(6)} + D_{(12)} + 2( 1 2)  ;  cost′′ = 12A + 3S + D_{(3)} + D_{(6)}

We don’t provide the code in this case, as it can be easily deduced from the one relative to the precedent case.

1 As GMP has an optimized function to divide a long integer by 3, it can be directly used, gaining efficiency.

5 Some savings in Toom-4.5 and Toom-5

Good IS’s for Toom-4.5 and Toom-5 were introduced in [2]. They were not proven to be optimal (in the considered model): because of their too big dimension it was impossible to find the optimal IS by exhaustive search. By slightly changing the model (in particular, “interlacing” ES and IS), it is possible to avoid one 1 2 (or a shift) for Toom-4.5 and two for Toom-5, for whatever a, b.

5.1 Saving in Toom-4.5

It uses \{∞, −1, −2, \frac{1}{2}, 1, 2, −\frac{1}{2}, 0\} as interpolation points (with lines corresponding to ±\frac{1}{2} opportunely multiplied by 2^{7}). The matrix is

M =(\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ −1 & 1 & −1 & 1 & −1 & 1 & −1 & 1 \\ −128 & 64 & −32 & 16 & −8 & 4 & −2 & 1 \\ 1 & 2 & 4 & 8 & 16 & 32 & 64 & 128 \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 128 & 64 & 32 & 16 & 8 & 4 & 2 & 1 \\ 1 & −2 & 4 & −8 & 16 & −32 & 64 & −128 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} )

The proposed IS contains the “decoupling” instructions for lines 3 and 6, related to values −2 and 2: the partial cost of these two operations is 2A + 1 2.

\begin{aligned}M_{3} & = M_{3} − M_{6} & ; & ( & −256 & 0 & −64 & 0 & −16 & 0 & −4 & 0 & ) \\ M_{6} & = (M_{6} ≪ 1) − M_{3} & ; & ( & 0 & 128 & 0 & 32 & 0 & 8 & 0 & 2 & ) \end{aligned}

By slightly mixing evaluation and interpolation phases it is possible to reduce this cost to 2A, as follows:

– First perform the ES for Toom-4.5 giving M, but modified such in a way that the obtained interpolation matrix is the aside shown M′: that is, simply reorganize the ES in order to put the evaluation in −2 in M_{8}, where c_{0} = a_{0}b_{0} – the evaluation in 0 – should be, and do not compute c_{0}. Note that M_{3} line remains for now undefined.

M′ =(\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ −1 & 1 & −1 & 1 & −1 & 1 & −1 & 1 \\ & ⟨not computed value⟩ \\ 1 & 2 & 4 & 8 & 16 & 32 & 64 & 128 \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 128 & 64 & 32 & 16 & 8 & 4 & 2 & 1 \\ 1 & −2 & 4 & −8 & 16 & −32 & 64 & −128 \\ −128 & 64 & −32 & 16 & −8 & 4 & −2 & 1 \end{matrix} )

– Then decouple and complete the ES

\begin{aligned}M_{3} & = M_{6} − M_{8} & ; & ( & 256 & 0 & 64 & 0 & 16 & 0 & 4 & 0 & ) \\ M_{6} & = M_{6} + M_{8} & ; & ( & 0 & 128 & 0 & 32 & 0 & 8 & 0 & 2 & ) \\ M_{8} & = a_{0}b_{0} & ; & ( & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & ) \end{aligned}

– Finally complete the IS following [2].

5.2 Saving in Toom-5

The corresponding matrix M is shown below: the IS proposed in [2] contained two decouplings: for values ±2 (lines 2 and 5) and ±\frac{1}{2} (lines 3 and 8). In this case mixing ES and IS results in a saving of 2( 1 2).

M =(\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 256 & −128 & 64 & −32 & 16 & −8 & 4 & −2 & 1 \\ 1 & 2 & 4 & 8 & 16 & 32 & 64 & 128 & 256 \\ 4^{8} & 4^{7} & 4^{6} & 4^{5} & 256 & 64 & 16 & 4 & 1 \\ 256 & 128 & 64 & 32 & 16 & 8 & 4 & 2 & 1 \\ 1 & −1 & 1 & −1 & 1 & −1 & 1 & −1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & −2 & 4 & −8 & 16 & −32 & 64 & −128 & 256 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} )

Alberto Zanoni  9

– First perform the ES for Toom-5 giving M, but modified such in a way that the obtained interpolation matrix is the aside shown M′: that is, simply reorganize the ES in order to put the evaluation in −2 in M_{1}, the evaluation in −\frac{1}{2} (multiplied by 2^{8}) in M_{9}, and do not compute c_{0} and c_{9}. Note that M_{2} and M_{8} lines remain for now undefined.

– Then decouple and adapt some IS steps to the new situation, as it must be slightly modified. The IS is a bit involved, we report gp-pari code for it:

M′ =( \begin{matrix}256 & −128 & 64 & −32 & 16 & −8 & 4 & −2 & 1 \\ & ⟨not computed value⟩ \\ 1 & 2 & 4 & 8 & 16 & 32 & 64 & 128 & 256 \\ 4^{8} & 4^{7} & 4^{6} & 4^{5} & 256 & 64 & 16 & 4 & 1 \\ 256 & 128 & 64 & 32 & 16 & 8 & 4 & 2 & 1 \\ 1 & −1 & 1 & −1 & 1 & −1 & 1 & −1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ & ⟨not computed value⟩ \\ 1 & −2 & 4 & −8 & 16 & −32 & 64 & −128 & 256 \end{matrix} )

\\ Evaluation of W0,...,W8 such that                         W1 = W1 + W7;
W0 <-- Evaluation in (-1/2) times 2^8                        W4 = W4 + W6;
W1 <-- Still not inizialized                                 W1 = W1 - (80*W3);
W2 <-- Evaluation in (1)                                     W6 = W6 - (510*W0);
W3 <-- Evaluation in (-1)                                    W5 = W5 - W7;
W4 <-- Evaluation in (2)                                     W6 = 3*W6 + W4;
W5 <-- Evaluation in (1/2) times 2^8                         W1 = W1/180;
W7 <-- Still not inizialized                                 W4 = W4 + (378*W2);
W8 <-- Evaluation in (-2)                                    W7 = W7>>2;
                                                             W3 = W3 - W7;
\\ Interpolation (and evaluation completion)                 W4 = W4 /(-72);
W7 = W4 - W8;                                                W6 = W6 /(-360);
W4 = W4 + W8;                                                W7 = W7 - W1;
W1 = W6 - W0;                                                W2 = W2 - W6;
W6 = W6 + W0;                                                W5 = W5 - (W4<<8);
W8 = Evaluation in (1/0);                                    W6 = W6 - W4;
W0 = Evaluation in (0);                                      W5 = W5 - (W6<<12);
W2 = W2 + W3;                                                W5 = W5 - (W2<<4);
W6 = W6 - W2;                                                W5 = W5 + (W3<<8);
W2 = W2 >> 1;                                                W3 = W3 + W7;
W3 = W2 - W3;                                                W7 = (W7*180) + W5;
W5 = W5 - W0;                                                W7 = W7/11340;
W5 = W5 - (W8<<16);                                          W5 = W5 + (720*W3);
W4 = W4 - (W2<<9);                                           W5 = W5/(-2160);
W2 = W2 - W8;                                                W3 = W3 - W5;
W2 = W2 - W0;                                                W1 = W1 - W7;

6 Conclusions

We showed how some Toom-Cook methods can be modified when one factor is even, permitting some savings. Lowerlevel implementations should be realized in order to understand the real gain over the classical methods. Possible optimizations of the interpolation phase of Toom-4.5 and Toom-5 were also suggested.

7 Acknowledgements

The author wishes to thank Marco Bodrato for many helpful discussions and fruitful suggestions and considerations. This work is dedicated to the memory of the author’s beloved grandfather, Lino Broffoni.

References

1. Marco Bodrato. Towards optimal Toom-Cook multiplication for univariate and multivariate polynomials in characteristic 2 and 0. In Claude Carlet and Berk Sunar, editors, WAIFI ’07 proceedings, volume 4547 of Lecture Notes in Computer Science. Springer, June 2007.

2. Marco Bodrato and Alberto Zanoni. Integer and polynomial multiplication: Towards optimal Toom-Cook matrices. In Christopher W. Brown, editor, Proceedings of the ISSAC 2007 Conference. ACM press, July 2007. URL:

http://bodrato.it/papers/#ISSAC2007.

3. Stephen A. Cook. On the minimum computation time of functions. PhD thesis, Department of Mathematics, Harvard University, 1966.

4. The GNU multiple precision (GMP) library documentation. URL: http://gmplib.org/#DOC.

5. Andrei L. Toom. The complexity of a scheme of functional elements realizing the multiplication of integers. Soviet Mathematics

Doklady, 3:714–716, 1963. URL: http://www.de.ufpe.br/~toom/articles/engmat/MULT-E.PDF.

A Toom-3.5 code when a_{0} is even

1. Slightly unbalanced case : We propose two complete procedures: when sa is or is not available. To avoid the use of any extra temporaries, the version without sa is sensibly different: evaluation and interpolation sequence are mixed. The values for c_{0} is computed during the IS, as soon it is needed, so that the variable that will contain it can meanwhile be used as a temporary one. In the code presented when sa is available we suppose that 2( 1 2) ⩽ S. If

  it is not the case, one should have  W0 = W0>>1; W4 = W0 + W3; W0 = W0 - W3;  instead of  {W4 = (W0>>1) + W3; \atop W0 = (W0>>1) - W3;}

With sa  Without sa

a = a3*x^3 + a2*x^2 + a1*x + a0;                          a = a3*x^3 + a2*x^2 + a1*x + a0;
b =          b2*x^2 + b1*x + b0;                          b =          b2*x^2 + b1*x + b0;
\\\\\\\\\\\\\\ Evaluation                                 \\\\\\\\\\\\\\ Evaluation
W0 = a0 + a2;        W1 = b0 + b2;                        W0 = a0 + a2;        W1 = b0 + b2;
W4 = a1 + a3;        W5 = W1 - b1;                        W4 = a1 + a3;        W5 = W1 - b1;
W3 = W0 - W4;                                             W3 = W0 - W4;
W2 = W3*W5;  \\ Evaluation in (-1)                        W2 = W3*W5;  \\ Evaluation in (-1)
W3 = W0 + W4;        W5 = W1 + b1;                        W3 = W0 + W4;        W5 = W1 + b1;
W1 = W3*W5;  \\ Evaluation in (1)                         W1 = W3*W5;  \\ Evaluation in (1)
W0 = a0 + (a2<<2);   W5 = W5 + b2;                        W0 = a2<<1;          W5 = W5 + b2;
W3 = a1 + (a3<<2);   W5 = (W5<<1) - b0;                   W3 = a0>>1;          W5 = W5<<1;
W4 = (W0>>1) + W3;                                        W0 = W0 + W3;        W5 = W5 - b0;
W0 = (W0>>1) - W3;                                        W3 = a3<<2;
W3 = W4*W5;  \\ Evaluation in (2) divided by 2.           W3 = W3 + a1;
                     W5 = W5 - (b1<<2);                   W4 = W0 + W3;
W4 = W0*W5;  \\ Evaluation in (-2) divided by 2.          W0 = W0 - W3;
W0 = a0*b0;  \\ Evaluation in (0)                         W3 = W4*W5;  \\ Evaluation in (2) divided by 2.
W5 = a3*b2;  \\ Evaluation in (1/0)                                            W4 = b1<<2;
                                                                               W5 = W5 - W4;
\\\\\\\\\\\\ Interpolation                                W4 = W0*W5;  \\ Evaluation in (-2) divided by 2.
W4 = W4 + W3;                                             W5 = a3*b2;  \\ Evaluation in (1/0)
W2 = W2 + W1;
W3 = W3 - (W4>>1);                                        \\\\\\\\\\\\ Interpolation
W2 = W2>>1;                                               W4 = W4 + W3;
W4 = W4 - W0;                                             W2 = W2 + W1;
W1 = W1 - W2;                                             W0 = W4>>1;
W3 = W3 - W1;                                             W3 = W3 - W0;
W2 = W2 - W0;                                             W2 = W2>>1;
W4 = W4 - (W2<<2);                                        W1 = W1 - W2;
W3 = W3/3;                                                W3 = W3 - W1;
W4 = W4/12;                                               W3 = W3/3;
W3 = W3 - (W5<<2);                                        W0 = W5<<2;
W1 = W1 - W3;                                             W3 = W3 - W0;
W2 = W2 - W4;                                             W1 = W1 - W3;
W3 = W3 - W5;                                             W3 = W3 - W5;
                                                          W0 = a0*b0;  \\ Evaluation in 0.
c = W5*x^5 + W4*x^4 + W3*x^3 + W2*x^2 + W1*x + W0;        W4 = W4 - W2;
                                                          W4 = W4/3;
                                                          W2 = W2 - W0;
                                                          W4 = W4 - W2;
                                                          W4 = W4>>2;
                                                          W2 = W2 - W4;
                                                          c = W5*x^5 + W4*x^4 + W3*x^3 + W2*x^2 + W1*x + W0;

Alberto Zanoni  11

2. Very unbalanced case : We propose just the evaluation phases when sa is (complete ES) or is not (not complete ES) available. The interpolation (pure or mixed) and the final reconstruction are the same as above, respectively.

With sa  Without sa

a = a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0;                 a = a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0;
b =                            b1*x + b0;                 b =                            b1*x + b0;
\\\\\\\\\\\\\\ Evaluation                                 \\\\\\\\\\\\\\ Evaluation
W0 = a0 + a4;                                             W0 = a0 + a4;
W0 = W0 + a2;                                             W0 = W0 + a2;
W4 = a1 + a3;                                             W4 = a1 + a3;
W3 = W0 - W4;        W5 = b0 - b1;                        W3 = W0 - W4;        W5 = b0 - b1;
W2 = W3*W5;  \\ Evaluation in (-1)                        W2 = W3*W5;  \\ Evaluation in (-1)
W3 = W0 + W4;                                             W3 = W0 + W4;
W1 = a2 + (a4<<2);                                        W1 = a4<<2;
W1 = a0 + (W1<<2);                                        W1 = W1 + a2;
W4 = a1 + (a3<<2);                                        W1 = W1<<2;
W0 = (W1>>1) + W4;                                        W1 = W1 + a0;
W1 = (W1>>1) - W4;   W5 = W5 - b1;                        W4 = a3<<2;
W4 = W1*W5;  \\ Evaluation in (-2) divided by 2.          W4 = W4 + a1;
                     W5 = b0 + b1;                        W1 = W1>>1;
W1 = W3*W5;  \\ Evaluation in (1)                         W0 = W1 + W4;
                     W5 = W5 + b1;                        W1 = W1 - W4;        W5 = W5 - b1;
W3 = W0*W5;  \\ Evaluation in (2) divided by 2.           W4 = W1*W5;  \\ Evaluation in (-2) divided by 2.
W0 = a0*b0;  \\ Evaluation in (0)                                              W5 = b0 + b1;
W5 = a4*b1;  \\ Evaluation in (1/0)                       W1 = W3*W5;  \\ Evaluation in (1)
                                                                               W5 = W5 + b1;
                                                          W3 = W0*W5;  \\ Evaluation in (2) divided by 2.
                                                          W5 = a4*b1;  \\ Evaluation in (1/0)

B Toom-3.5 code when b_{0} is even

1. Slightly unbalanced case : We propose just the evaluation phases when sa is (complete ES) or is not (not complete ES) available. The interpolation (pure or mixed) and the final reconstruction are the same as in the case when a_{0} is even, respectively.

With sa  Without sa

a = a3*x^3 + a2*x^2 + a1*x + a0;                          a = a3*x^3 + a2*x^2 + a1*x + a0;
b =          b2*x^2 + b1*x + b0;                          b =          b2*x^2 + b1*x + b0;
\\\\\\\\\\\\\\ Evaluation                                 \\\\\\\\\\\\\\ Evaluation
W0 = a0 + a2;        W1 = b0 + b2;                        W0 = a0 + a2;        W1 = b0 + b2;
W4 = a1 + a3;        W5 = W1 - b1;                        W4 = a1 + a3;        W5 = W1 - b1;
W3 = W0 - W4;                                             W3 = W0 - W4;
W2 = W3*W5;  \\ Evaluation in (-1)                        W2 = W3*W5;  \\ Evaluation in (-1)
W3 = W0 + W4;        W5 = W1 + b1;                        W3 = W0 + W4;        W5 = W1 + b1;
W1 = W3*W5;  \\ Evaluation in (1)                         W1 = W3*W5;  \\ Evaluation in (1)
W0 = a0 + (a2<<2);   W5 = W5 + b2;                        W0 = a2<<2;          W5 = W5 + b2;
W3 = a1 + (a3<<2);   W5 = W5 - (b0>>1);                   W0 = W0 + a0;        W4 = b0>>1;
W4 = W0 + (W3<<1);                                        W3 = a3<<2;          W5 = W5 - W4;
W0 = W0 - (W3<<1);                                        W3 = W3 + a1;
W3 = W4*W5;  \\ Evaluation in (2) divided by 2.           W3 = W3<<1;
                     W5 = W5 - (b1<<1);                   W4 = W0 + W3;
W4 = W0*W5;  \\ Evaluation in (-2) divided by 2           W0 = W0 - W3;
W0 = a0*b0;  \\ Evaluation in (0)                         W3 = W4*W5;  \\ Evaluation in (2) divided by 2.
W5 = a3*b2;  \\ Evaluation in (1/0)                                            W4 = b1<<1;
                                                                               W5 = W5 - W4;
                                                          W4 = W0*W5;  \\ Evaluation in (-2) divided by 2.
                                                          W5 = a3*b2;  \\ Evaluation in (1/0)

2. Very unbalanced case : We propose just the evaluation phases when sa is (complete ES) or is not (not complete ES) available. The interpolation (pure or mixed) and the final reconstruction are the same as in the case when a_{0} is, respectively.

With sa  Without sa

a = a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0;                 a = a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0;
b =                            b1*x + b0;                 b =                            b1*x + b0;
\\\\\\\\\\\\\\ Evaluation                                 \\\\\\\\\\\\\\ Evaluation
W0 = a0 + a4;                                             W0 = a0 + a4;
W0 = W0 + a2;                                             W0 = W0 + a2;
W5 = a1 + a3;                                             W4 = a1 + a3;
W3 = W0 - W5;        W4 = b0 - b1;                        W3 = W0 - W4;        W5 = b0 - b1;
W2 = W3*W4;  \\ Evaluation in (-1)                        W2 = W3*W5;  \\ Evaluation in (-1)
W3 = W0 + W5;        W4 = b0 + b1;                        W3 = W0 + W4;        W5 = b0 + b1;
W1 = W3*W4;  \\ Evaluation in (1)                         W1 = W3*W5;  \\ Evaluation in (1)
W0 = a2 + (a4<<2);                                        W0 = a4<<2;
W0 = a0 + (W0<<2);                                        W0 = W0 + a2;
W3 = a1 + (a3<<2);                                        W0 = W0<<2;
                                                          W0 = W0 + a0;
W4 = W0 + (W3<<1);                                        W3 = a3<<2;
W0 = W0 - (W3<<1);   W5 = (b0>>1) + b1;                   W3 = W3 + a1;
W3 = W4*W5;  \\ Evaluation in (2) divided by 2.           W3 = W3<<1;
                     W5 = b0 - W5;                        W4 = W0 + W3;        W5 = b0>>1;
W4 = W0*W5;  \\ Evaluation in (-2) divided by 2.          W0 = W0 - W3;        W5 = W5 + b1;
                                                          W3 = W4*W5;  \\ Evaluation in (2) divided by 2.
W0 = a0*b0;  \\ Evaluation in (0)                                              W5 = b0 - W5;
W5 = a4*b1;  \\ Evaluation in (1/0)                       W4 = W0*W5;  \\ Evaluation in (-2) divided by 2.
                                                          W5 = a4*b1;  \\ Evaluation in (1/O).