More On Squaring and Multiplying Large Integers

ABSTRACT

This is an expansion of the paper On Squaring and Multiplying Large Integers presented at the IEEE Symposium on Computer Arithmetic (ARITH-11) in July 1993. This paper presents more of the ideas behind these methods; some new results; and some speculations for areas of future work.

Much of the work presented here is based on methods that appear in Don Knuth's classic work The Art of Computer Programming as well as the work of A. L. Toom.

In actual experiments, some of the simpler methods of squaring and multiplying are shown to be best in a surprisingly large number of cases. Some new methods are also discussed which are useful out to quite large numbers. In light of these results it seems that methods such as the Schönhage and Strassen FFT multiply, while of theoretical interest, may never be best for any reasonably sized numbers.

The problem to be discussed here is how to find the best (fastest) ways to square and multiply large numbers in software. The methods presented here were implemented in C and assembly language on an HP-9000/720 but the general observations and conclusions should be true in wider areas of application.

The approach used was to write a collection of routines for squaring and multiplying w-word numbers producing 2w-word results for various specific values of w. These were general purpose routines in the sense that w was a parameter.

In the course of writing these routines, it became clear that some common operations on large integers stored as arrays of unsigned 32-bit words would be needed. These operations were written in assembly language and performed the functions detailed in Table 1.

A collection of routines for squaring and multiplying "small" integers were also written in assembly. These provided a basis upon which the larger routines could be built.

All other routines were written in C and compiler optimized.

Running time was measured for w in the range 1 to 1,175,553 words (37,617,696 bits) by counting instruction cycles on a 50MHz HP-9000/720.

We will define T_Si(w) as the time required for method i to square a number of length w and T_Mi(w) as the time required for method i to multiply a number of length w. (We will use the expression T_i(w) when discussing either method or method i in general.)

Method i is considered best for some length w if

T_i(w) <= T_j(w') for all j and for all w' >= w.

That is, a method is best for a given length if there is no faster way of squaring or multiplying numbers at least as large as this one.

To start the ball rolling, basis routines for operating on small integers were written using the obvious method.

Roughly speaking, if the algorithm for a multiply is

FOR i = w-1 to 0 :

FOR j = w-1 to 0 :

(r[i+j,i+j+1] and CarryOut) += a[i] × b[i],

then squaring can be made strictly faster than multiply by accumulating the upper half off-diagonal elements separately and combining them with the diagonals according to the formula

Result = 2 × Offdiagonals + Diagonals.

Obviously, the running time, T₁(w), increases quadratically. T₁(w), in cycles, is well approximated on the HP-9000/720 by the least squares formulas

T_S1(w) = 4.53w² + 10w + 30.71 and
T_M1(w) = 9.82w² - 9.06w + 90.26.

These routines turned out to be best for all w up to 33 for square (23 for multiply).

It seems that the first discussion of a method of multiplication (squaring, actually) that was better than w² in the width of the operand was a short paper written by Karatsuba and Ofman that appeared in the Soviet journal Doklady in 1963. The method involves splitting the operand into 2 parts: the high order bits and the low order bits, based on the polynomial

(A₁x + A₀)² = A₁²(x² - x) + (A₁ + A₀)²x + A₀²(1 - x)

where x = 2ⁿ.

Knuth presented a minor improvement on this based on the slightly different polynomial

(A₁x + A₀)² = A₁²(x² + x) - (A₁ - A₀)²x + A₀²(x + 1).

Both of these may be regarded as special cases of solving for the C_i in the polynomial

(A₁x + A₀)² = C₂x² + C₁x + C₀.

Karatsuba and Ofman solve this equation by letting x take on the values {Infinity, 1, 0} where x = Infinity corresponds to

lim x --> Infinity

(A1x + A0)2/x2 =

lim x --> Infinity

(C2x2 + C1x + C0)/x2

in this context.

Rather than express this system in the usual polynomial form, I will express it as a linear transformation of the operand in the following way

V2 V1 V0

=

A1 A1 + A0 A0

=

1 0 1 1 0 1

A1 A0

(Please forgive the fairly primitive look of these matrix expressions. It is the best I have been able to do under the limitations of html.)

After squaring each element of the vector, the solution to the resulting system

1 0 0 1 1 1 0 0 1

C2 C1 C0

=

S2 S1 S0

=

A12 (A1 + A0)2 A02

may be expressed as

C2 C1 C0

=

1 0 0 -1 1 -1 0 0 1

S2 S1 S0

This method has a small problem in that the term A₁ + A₀ requires one more bit than the other two. Karatsuba and Ofman address this problem by showing how you can square an n+1-bit number using an n-bit algorithm together with a couple of shifts and adds.

Knuth solves the equation by letting x take on the values {Infinity, -1, 0}. A similar transformation

V2 V1 V0

=

A1 A1 - A0 A0

=

1 0 1 -1 0 1

A1 A0

results, after squaring, in a similar system

1 0 0 1 -1 1 0 0 1

C2 C1 C0

=

S2 S1 S0

=

A12 (A1 - A0)2 A02

which has the solution

C2 C1 C0

=

1 0 0 1 -1 1 0 0 1

S2 S1 S0

Knuth presented this method as Karatsuba and Ofman's but Knuth's variation leads to a real improvement. By rearranging the algebra a bit, he was able to replace the A₁ + A₀ term with A₁ - A₀. Since the only use of the term is to square it, its sign is unimportant. Therefore, one can always subtract the lesser from the greater and save that extra bit.

The test and branch involved results in much less overhead than Karatsuba and Ofman's method.

This method can be turned into a multiply method by replacing the S_i with the P_i in the formula

1 0 0 1 -1 1 0 0 1

C2 C1 C0

=

P2 P1 P0

=

A1B1 (A1 - A0)(B1 - B0) A0B0

and doing the implied extra work.

These methods are diagrammed in Figure 1.

Assuming that the time for an add (or subtract) of length w is O(w) , the time to square a number (or multiply two numbers) of length 2w would be

T₂(2w) = 3T_best(w) + O(w).

The overhead involved, one add before the squaring step (two adds before the multiplying step) and two after, is strictly larger than the overhead for a w² method. Thus, although a w² squaring (or multiply) of length 2w would take

T₁(2w) = 4T_best(w) + O(w),

there is a definite crossover point as w increases where T₁(w) is roughly the same as T₂(w). This point occurs where the time involved in the extra overhead of Knuth's method (which I will call the 2-way method from this point on) is roughly equivalent to the time for the 4^th multiply in the w² method.

Below this point the w² method wins because of its simplicity and above this point the 2-way method wins because of its better asymptotic behavior.

Figure 2 illustrates this by plotting the ratio T₂(w) / T₁(w) for some 2-way routines. (T₁(w) is actually the least squares approximation.)

Ratio of T₂(w), the running time of Knuth's 2-way multiply to T₁(w), the running time of the w² basis.

In the limit, the running time of the 2-way method triples each time the size of the problem doubles. Therefore, its time increases as w^{log(3) / log(2)} = w^1.585.

The running time of the 2-way routines (in cycles) is approximated by the least squares formulas

T_S2(w) = 26.64w^{log(3) / log(2)} - 54.63w + 227 and
T_M2(w) = 46.2w^{log(3) / log(2)} - 84.11w + 979.

Building on a basis set of w² routines up to 33 words for square (23 words for multiply), the 2-way method takes over at 34 words (24 words) and generally dominates in these experiments up to 768 words (384 words).

Mark Shand has pointed out that this crossover point is a strong function of the ratio of the multiply time to add time on a given machine. The HP-9000/720 performs an unsigned 32 × 32 --> 64 multiply in 2 cycles (5 to 7 cycles if you include load/store overhead). On a machine with a relatively slower multiply the crossover point can be much lower.

We now have enough background to show the general idea behind these squaring (and multiplying) methods.

In each squaring case, the operand is split into k parts A_j. These parts are transformed via a Q_ij (derived from the original polynomial) into m parts V_i = Q_ijA_j. Then, each of the V_i is squared with the best method available, producing S_i = V_i².

In each multiply case, the operands are split into k parts A_j and B_j. These parts are transformed via Q_ij into m parts V_i = Q_ijA_j and W_i = Q_ijB_j. Then, each of the V_i is multiplied by W_i with the best method available, producing P_i = V_iW_i.

Now, the solution to the system M_ijC_j = S_i ( M_ijC_j = P_i) (derived from the square of the original polynomial) can be expressed as C_i = D_ik^-1R_kjS_j ( C_i = D_ik^-1R_kjP_j) with all the coefficients in the integers. (Up to now, D_ik^-1, which expresses the integer divides, has not been needed as it was the identity matrix.) Finally, the C_i are shifted and added up to produce the result.

Toom showed that circuits could be constructed for squaring integers where the size of the circuit was bounded by O(nc^Sqrt(log(n))) and the delay was bounded by O(c^Sqrt(log(n))).

Cook showed that an algorithm for squaring integers could be realized which had a running time of O(n2^{5Sqrt(log(n))}).

This algorithm, which is actually a collection of algorithms, has come to be known as the Toom-Cook method.

The first method in this collection is equivalent to Karatsuba's method.

The second method divides a number into three parts and involves solving for the C_i in the polynomial

(A₂x² + A₁x + A₀)² = C₄x⁴ + C₃x³ + C₂x² + C₁x + C₀.

Cook solves this equation by letting x take on the values {4, 3, 2, 1, 0} in the polynomial, leading to the transformation

V4 V3 V2 V1 V0

=

16 4 1 9 3 1 4 2 1 1 1 1 0 0 1

A2 A1 A0

which, after squaring, yields the system

256 64 16 4 1 81 27 9 3 1 16 8 4 2 1 1 1 1 1 1 0 0 0 0 1

C4 C3 C2 C1 C0

=

S4 S3 S2 S1 S0

=

(16A2 + 4A1 + A0)2 (9A2 + 3A1 + A0)2 (4A2 + 2A1 + A0)2 (A2 + A1 + A0)2 A02

the solution of which I will express as

C4 C3 C2 C1 C0

=

-1 24 0 0 0 0 0 12 0 0 0 0 0 24 0 0 0 0 0 12 0 0 0 0 0 1

1 -4 6 -4 1 -3 14 -24 18 -5 11 -56 114 -104 35 -3 16 -36 48 -25 0 0 0 0 1

S4 S3 S2 S1 S0

The diagonal terms in the first matrix are the denominators in divides of intermediate results. In practice, there turn out to be four divides by 3 since the rest amounts to shifts of 2 or 3 bits.

Knuth solves this with x taken from {2, 1, 0, -1, -2} leading to the transformation

V4 V3 V2 V1 V0

=

4 2 1 1 1 1 0 0 1 1 -1 1 4 -2 1

A2 A1 A0

and the system

16 8 4 2 1 1 1 1 1 1 0 0 0 0 1 1 -1 1 -1 1 16 -8 4 -2 1

C4 C3 C2 C1 C0

=

S4 S3 S2 S1 S0

=

(4A2 + 2A1 + A0)2 (A2 + A1 + A0)2 A02 (A2 - A1 + A0)2 (4A2 - 2A1 + A0)2

which has the somewhat more symmetrical solution

C4 C3 C2 C1 C0

=

-1 24 0 0 0 0 0 12 0 0 0 0 0 24 0 0 0 0 0 12 0 0 0 0 0 1

1 -4 6 -4 1 1 -2 0 2 -1 -1 16 -30 16 -1 -1 8 0 -8 1 0 0 1 0 0

S4 S3 S2 S1 S0

Although both of these solutions look similar, Knuth's is somewhat better than the previous one. The sign symmetry simplifies the system and greatly reduces the overhead (in the form of add, subtract, shift, and the like).

If, however, one uses reciprocal symmetry as well as sign symmetry in the choices of x, as in the set {Infinity, 2, 1, 1/2, 0}, then the greater symmetry in the system results in even greater simplification and greater reduction in overhead. (In general, x = p/q corresponds to q²ⁿ(A_n(p/q)ⁿ + ... + A₀)² = q²ⁿ(C_2n(p/q)²ⁿ + ... + C₀).)

This system

V4 V3 V2 V1 V0

=

1 0 0 4 2 1 1 1 1 1 2 4 0 0 1

A2 A1 A0

and

1 0 0 0 0 16 8 4 2 1 1 1 1 1 1 1 2 4 8 16 0 0 0 0 1

C4 C3 C2 C1 C0

=

S4 S3 S2 S1 S0

=

A22 (4A2 + 2A1 + A0)2 (A2 + A1 + A0)2 (A2 + 2A1 + 4A0)2 A02

has the solution

C4 C3 C2 C1 C0

=

-1 1 0 0 0 0 0 6 0 0 0 0 0 2 0 0 0 0 0 6 0 0 0 0 0 1

1 0 0 0 0 -21 2 -12 1 -6 7 -1 10 -1 7 -6 1 -12 2 -21 0 0 0 0 1

S4 S3 S2 S1 S0

(Oddly enough, it appears that reciprocal symmetry is more important than sign symmetry in the sense that, if one must break one or the other, the system seems simpler if sign symmetry is broken rather than reciprocal symmetry.)

This solution has only 2 divides by 3 and 3 shifts as well as much less overhead in the computation of the intermediate results.

The computation of S₂ requires 2 extra bits and S₁ and S₃ each require 3.

The squaring method implied by this solution is shown in Figure 3.

Notice that it comes out looking like a complicated version of the 2-way method in that there is an initial transformation stage, a squaring stage, and a final transformation stage.

All of the operations shown can be made O(w). Shifts can be done separately or as part of the combined adds and subtracts.

The time to square (multiply) a number of length 3w becomes

T₃(3w) = 5T_best(w) + O(w).

This suggests that the asymptotic running time of the 3-way method increases as w^{log(5) / log(3)} = w^1.465.

The improvement is not nearly so dramatic as the improvement in the 2-way method over the w² method. The advantages of more complex routines will be even less dramatic for much more work.

The running time of the 3-way method is approximated by the formulas

T_S3(w) = 64.59w^{log(5) / log(3)} - 243.51w + 33913 and
T_M3(w) = 90.15w^{log(5) / log(3)} + 75.91w - 44016.

As before, to turn this method into a multiply, one duplicates those operations before the squaring stage and multiplies the corresponding terms instead of squaring them.

The running time of a Toom-Cook, Knuth, and k-way methods will all have three components: m = 2k - 1 squarings of elements of size w/k; some O(w) overhead; and some fixed overhead.

Assuming that one starts with a basis routine of length w₀ that takes t₀ time, we have

T_k(w₀) = t₀
T_k(k^r+1w₀) = mT_k(k^rw₀) + c₁k^rw₀ + c₀
T_k(kw₀) = mt₀ + c₁w₀ + c₀
T_k(k²w₀) = m²t₀ + c₁w₀(m + k) + c₀(m + 1)
T_k(k³w₀) = m³t₀ + c₁w₀(m² + mk + k²) + c₀(m² + m + 1)
...
T_k(k^rw₀) = m^rt₀ + c₁w₀(m^r - k^r) / (m - k) + c₀(m^r - 1) / (m - 1).

Rearranging terms in this last equation gives

T_k(k^rw₀) = (t₀ + c₁w₀ / (m - k) + c₀ / (m - 1))m^r - c₁w₀k^r / (m - k) - c₀ / (m - 1)

or, in terms of w = k^rw₀ becomes

T_k(w) = (t₀ + c₁w₀ / (m - k) + c₀ / (m - 1))(w/w₀)^{log(m) / log(k)} - c₁w / (m - k) - c₀ / (m - 1).

This formula has the expected form: a dominant (w/w₀)^{log(m) / log(k)} term, an O(w) term, and a constant term, both of which may be neglected for large w.

But, and this is important, the dominant term contains terms proportional to the overhead. Thus, while the overhead terms may be neglected for large w, the overhead itself may not as it directly contributes to the magnitude of the dominant term.

Extension to a k-way method involves solving for the C_i in the polynomial

(A_kx^k + ... + A₀)² = C_2kx^2k + ... + C₀.

Cook would solve this polynomial by letting x take on the values {0, ..., 2k}.

Knuth would solve this polynomial by letting x take on the values {-k, ..., k}.

But if one chooses x from the reciprocally symmetric set {1, Infinity, 0, 2, 1/2, -2, -1/2, 3, 1/3, -3, -1/3, 3/2, 2/3, ...}, that is, the set of small rational numbers, +/- p/q, such that the GCD(p,q) = 1, the resulting system of equations is very symmetrical and the algorithm is simpler.

I will illustrate the various methods for k = 4.

Cook:

V6 V5 V4 V3 V2 V1 V0

=

216 36 6 1 125 25 5 1 64 16 4 1 27 9 3 1 8 4 2 1 1 1 1 1 0 0 0 1

A3 A2 A1 A0

46656 7776 1296 216 36 6 1 15625 3125 625 125 25 5 1 4096 1024 256 64 16 4 1 729 243 81 27 9 3 1 64 32 16 8 4 2 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1

C6 C5 C4 C3 C2 C1 C0

=

S6 S5 S4 S3 S2 S1 S0

=

(256A3 + 36A2 + 6A1 + A0)2 (125A3 + 25A2 + 5A1 + A0)2 (64A3 + 16A2 + 4A1 + A0)2 (27A3 + 9A2 + 3A1 + A0)2 (8A3 + 4A2 + 2A1 + A0)2 (A3 + A2 + A1 + A0)2 A02

C6 C5 C4 C3 C2 C1 C0

=

-1 720 0 0 0 0 0 0 0 240 0 0 0 0 0 0 0 144 0 0 0 0 0 0 0 48 0 0 0 0 0 0 0 360 0 0 0 0 0 0 0 60 0 0 0 0 0 0 0 1

1 -6 15 -20 15 -6 1 -5 32 -85 120 -95 40 -7 17 -114 321 -484 411 -186 35 -15 104 -307 496 -461 232 -49 137 -972 2970 -5080 5265 -3132 812 -10 72 -225 400 -450 360 -147 0 0 0 0 0 0 1

S6 S5 S4 S3 S2 S1 S0

Knuth:

V6 V5 V4 V3 V2 V1 V0

=

27 9 3 1 8 4 2 1 1 1 1 1 0 0 0 1 -1 1 -1 1 -8 4 -2 1 -27 9 -3 1

A3 A2 A1 A0

729 243 81 27 9 3 1 64 32 16 8 4 2 1 1 1 1 1 1 1 1 0 0 0 0 0 0 1 1 -1 1 -1 1 -1 1 64 -32 16 -8 4 -2 1 729 -243 81 -27 9 -3 1

C6 C5 C4 C3 C2 C1 C0

=

S6 S5 S4 S3 S2 S1 S0

=

(27A3 + 9A2 + 3A1 + A0)2 (8A3 + 4A2 + 2A1 + A0)2 (A3 + A2 + A1 + A0)2 A02 (-A3 + A2 - A1 + A0)2 (-8A3 + 4A2 - 2A1 + A0)2 (-27A3 + 9A2 - 3A1 + A0)2

C6 C5 C4 C3 C2 C1 C0

=

-1 720 0 0 0 0 0 0 0 240 0 0 0 0 0 0 0 144 0 0 0 0 0 0 0 48 0 0 0 0 0 0 0 360 0 0 0 0 0 0 0 60 0 0 0 0 0 0 0 1

1 -6 15 -20 15 -6 1 1 -4 5 0 -5 4 -1 -1 12 -39 56 -39 12 -1 -1 8 -13 0 13 -8 1 2 -27 270 -490 270 -27 2 1 -9 45 0 -45 9 -1 0 0 0 1 0 0 0

S6 S5 S4 S3 S2 S1 S0

4-way:

V6 V5 V4 V3 V2 V1 V0

=

1 0 0 0 8 4 2 1 -8 4 -2 1 1 1 1 1 1 -2 4 -8 1 2 4 8 0 0 0 1

A3 A2 A1 A0

1 0 0 0 0 0 0 64 32 16 8 4 2 1 64 -32 16 -8 4 -2 1 1 1 1 1 1 1 1 1 -2 4 -8 16 -32 64 1 2 4 8 16 32 64 0 0 0 0 0 0 1

C6 C5 C4 C3 C2 C1 C0

=

S6 S5 S4 S3 S2 S1 S0

=

A32 (8A3 + 4A2 + 2A1 + A0)2 (-8A3 + 4A2 - 2A1 + A0)2 (A3 + A2 + A1 + A0)2 (A3 - 2A2 + 4A1 - 8A0)2 (A3 + 2A2 + 4A1 + 8A0)2 A02

This last solution was implemented as shown in the rather busy Figure 4.

Except for its complexity, it is much the same as the other methods.

The computation of the S₃ term requires 2 extra bits and S₁, S₂, S₄, and S₅ each require 4.

The 4-way method is asymptotically w^{log(7) / log(4)} = w^1.404 and is approximated on the HP-9000/720 by the formula

T_S4(w) = 120.79w^{log(7) / log(4)} - 858.7w + 976623.

Figure 5 shows the equations necessary to implement a 5-way method. This method would be asymptotically w^{log(9) / log(5)} = w^1.365 but has not yet been implemented. Out of a sense of mercy to the reader, I won't even attempt a graph of the solution.

The best known method of squaring or multiplying integers is the FFT multiply discovered by Schönhage and Strassen. That is, this is the method with the best known asymptotic behavior of O(wlog(w)log(log(w))).

I shall not present this method here in any great detail but I recommend the excellent description in Chapter 7 of Aho, Hopcroft, and Ullman.

In brief, the FFT multiply of order k splits a number into 2^k-1 parts for increasingly large values of k, performs an order k FFT on it, squares (or multiplies) those 2^k elements, and performs an inverse FFT on the result. All arithmetic is performed in a ring with a solution to the equation z^2k-1 = -1.

If the ring chosen is the integers mod some base, b, the equation becomes z^2k-1 = -1 mod b and b must be larger than the result coefficients.

It is common to choose b = 2^m2k-1 + 1 with m2^k-1 > 2w(log(w) / log(2)) so that z = 2^m, although other choices are possible. (This transform is also known as a Number Theoretic Transform or NTT.)

The FFT multiply has the advantages that: it is easily extensible to any value of k; both the initial and final transformations can be made to take place in stages that have identical connectivity; and, since all the divides are by powers of 2, they can be done as shifts.

But, it has a major disadvantage in that, since the ring must be large enough to represent the result coefficients, all of the basic operations are twice as large as in a corresponding Toom-Cook style method.

Therefore, the overhead involved can be large compared to the k-way methods.

For example, with z = 2^m and b = 2^2m + 1, we have z² = -1 mod b and the transformation for an FFT-2 would look like

V3 V2 V1 V0

=

A0 + A1 A0 - A1 A0 + A1z A0 - A1z

=

1 0 1 -1 1 z 1 -z

A1 A0

and

z -1 -z 1 -z -1 z 1 -1 1 -1 1 1 1 1 1

C3 C2 C1 C0

=

S3 S2 S1 S0

=

V32 V22 V12 V02