Notes on fast multiplication of large numbers

Bill Casselman Department of Mathematics University of British Columbia

cass@math.ubc.ca

If one multiplies two numbers of N digits by the elementary method taught in grammar school it takes time roughly proportional to N^{2}. If N is relatively small this is reasonably efficient, but for large numbers there are ways to speed this calculation up enormously. For very large N, the best techniques use the fast Fourier transformation and its inverse as intermediate steps, reducing the amount of time roughly to O(N log N). These techniques are not difficult to implement in practical programs.

Using the fast Fourier transform for this purpose was apparently first suggested by Schönhage and Strassen. The basic reference for this now is the second edition of volume II in Knuth’s series on the art of computer programming. The technique which I will explain most carefully has the remarkable virtue that it doesn’t use floating point calculations. This elegant modification of the original idea comes from Schönhage, who did it in the context of arithmetic modulo 2^{n} +1. The basic ideas of this procedure became somewhat clearer when Nussbaumer used them to do convolution of polynomials. The first explicit description of Schönhage-Nussbaumer that I am aware of is the solution to exercise 4.6.4.59 in Knuth’s book (although it apparently gets the order of FFT transformations wrong). For me, none of these accounts of Schönhage-Nussbaumer are terribly clear, and it was not until Joe Buhler explained to me an exposition of Daniel Bernstein’s that I understood what was going on.

In fact, I have benefitted from conversations with Buhler throughout writing this note and the accompanying programs.

Polynomials and multi-precision arithmetic

Suppose we are doing arithmetic in base B. Every non-negative number is a sum

p = ∑_{0}^{N−1} p_{i}B^{i}

with 0 ≤ p_{i} ≤ B − 1. If P is the polynomial we obtain by replacing B by an indeterminate z then the map z ↦ B takes P (z) to p = P (B). If P maps to p, Q to q then R = P Q maps to r = pq. Replacing multiplication of integers by multiplication of polynomials is the first step towards fast multiplication of large numbers.

At first sight, the transition from the polynomial product to the integer product is a bit complicated. This is because the coefficients of the polynomial product will not usually be in the range [0, B), so some mild reduction must be done. This replaces carrying in the usual multiplication technique.

For a practical program, we shall need some estimate of the sizes of the product coefficients. They are given by the convolution sums

r_{k} = ∑_{0}^{k} p_{i}q_{k−i} (k = 0, …, 2N − 1)

and we can use this to estimate that the coefficient r_{k} of degree k satisfies

0 ≤ r_{k} ≤ C (B − 1)^{2}, C = min(k + 1, 2N − k)

and that all coefficients are at most N(B − 1)^{2}. They will not generally be in the range [0, B), but they can be expressed as numbers in base B with a fixed precision depending on N, and the process of reducing the product again to base B notation will not be difficult. Essentially it will amount to addition of multi-precision numbers of a fixed length.

For example, to multiply the decimal numbers 97 and 86 we multiply the polynomials 9z + 7 and 8z + 4 to get

(9z + 7)(8z + 6) = 72z^{2} + (56 + 54)z + 42

which we reduce according to the tableau

4 2 1 1 0 7 2 8 3 4 2

For computers the most efficient base to work with is a power of two. If B = 2^{n} then any positive integer no larger than N(B − 1)^{2} can be represented by less than 2n + log_{2} N bits. If n is the number of bits in a machine word then for practical purposes this will be at most 3n since the word size in a machine controls both the size of integer representation as well as the size of addresses of data. (We must be able to address the array holding the coefficients in n bits also.) If we know the product of the two polynomials, we can therefore reconstruct the product pq by arithmetic addition (including carrying) in numbers of 3n bits, or 3 machine words.

As I have mentioned, most if not all fast multi-precision arithmetic uses the translation from numbers into polynomials and back again, so the first component in a complete package is a routine to reduce a polynomial with limited precision coefficients back to a number in base B.

General fast methods of multiplication

The one basic trick in nearly all methods of performing fast multi-precision multiplication is the combination of evaluation and interpolation of polynomials in some guise or another. If P and Q are two polynomials of degree k − 1 and ℓ − 1 then P Q has degree k + ℓ − 2. To calculate P Q we evaluate P and Q at k + ℓ − 1 points z_{i}, and interpolate the product P Q from its values P (z_{i})Q(z_{i}).

One of the simplest techniques doing this is Karatsuba’s. Suppose we are given two polynomials P (z) = a_{0} +a_{1}z and Q(z) = b_{0} +b_{1}z. Then the product R(z) = P (z)Q(z) has degree at most two and can be written

c_{0} + c_{1}z + c_{2}z^{2} .

It is determined by its values at three distinct points. This suggests that we can form the polynomial product by evaluating the product polynomial at three points by three multiplications and then recombining. Note that the straightforward formulas for the product coefficients

\begin{aligned}c_{0} & = a_{0}b_{0} \\ c_{0} & = a_{0}b_{1} + a_{1}b_{0} \\ c_{2} & = a_{1}b_{1} \end{aligned}

require in contrast four multiplications.

This idea can be extended to deal with homogeneous linear polynomials a_{0}x + a_{1}y and b_{0}x + b_{1}y, and this gives the most useful formulas. Explicitly, we have the formulas

(a_{0}x + a_{1}y)(b_{0}x + b_{1}y) = c_{0}x^{2} + c_{1}xy + c_{2}y^{2} .

We can write

\begin{aligned}P (0) & = a_{0} \\ P (1) & = a_{0} + a_{1} \\ P (−1) & = a_{0} − a_{1} \\ P (∞) & = a_{1} \end{aligned}

and hence derive, effectively interpolating from the values at 0, ±1, and ∞

c_{0} = P (0)Q(0) = a_{0}b_{0}
c_{1} = P (1)Q(1) − P (0)Q(0) − P (∞)Q(∞)
= (a_{0} + a_{1})(b_{0} + b_{1}) − a_{0}b_{0} − a_{1}b_{1}
c_{2} = P (∞)Q(∞) = a_{1}b_{1}

or also in slightly different form

c_{0} = P (0)Q(0) = a_{0}b_{0}
c_{1} = P (0)Q(0) + P (∞)Q(∞) − P (−1)Q(−1)
= a_{0}b_{0} + a_{1}b_{1} − (a_{0} − a_{1})(b_{0} − b_{1})
c_{2} = P (∞)Q(∞) = a_{1}b_{1}

Both these formulas require only three multiplications. Both of these can be used recursively on numbers of 2^{n} words to calculate products using 3^{n} rather than 4^{n} multiplications. This is very practical, and in practice timings are indeed cut by a factor (3/4)^{n}.

By the way, having at your disposal a method of multiplication which is less than quadratic in the size of multiplicands changes drastically your choice of algorithms for many operations. For example, to calculate large factorials in this scheme you would incorporate an algorithm to calculate the product of successive integers i … j − 1 for arbitrary positive i and j, which calls itself by recursion after dividing the range in two halves.

I find that in practice an efficient implementation of the standard algorithm works best for numbers of size up to about eight words, but that above that Karatsuba works better. If calculating products of words of length 2n an amount of temporary storage equal to 2n is also needed, and since Karatsuba is recursive it requires in input a total amount of storage equal to 2n + n + n/2 + … which for all practical purposes equals 4n.

So the second component in a package for fast multiplication is a routine to implement Karatsuba’s technique. In the following code, the product of arrays p and q, is to be stored in r, using tmp for scratch work. The arrays p and q have length 2n, r and tmp have length 4n. The array p_{0} is the lower half of p, p_{1} its upper half.

karatsuba(r, p, q, tmp, n)

if n < limit then do ordinary multiplication

otherwise

call karatsuba to store the product p_{0}q_{0} in the lower end of r, using the lower half of tmp as scratch, and passing the upper half as scratch storage in the call

store |p_{1} − p_{0}| in the third quarter r_{2} of r and |q_{1} − q_{0}| in the fourth quarter r_{3}, taking note of signs for later use

call karatsuba to store the product |p_{1} − p_{0}||q_{1} − q_{0}| in the lower half of tmp, using the upper half as storage in the call now call karatsuba to store p_{1}q_{1} in the top half of tmp, using the top half of r as storage in the call

add (or subtract) the pieces into r, carefully

In writing my own programs to carry this process out, I have incorporated a few tricks of memory allocation from the low-level package libI distributed with LiDiA.

The fast Fourier transform

Karatsuba’s method amounts to evaluation and interpolation involving three points. For numbers of medium size it probably can’t be bettered. But for large calculations it is more efficient to use N-th roots of unity instead, where one can use the Fourier transform to do both evaluation and interpolation.

Suppose P to be a polynomial of degree less than N, or

P (z) = p_{0} + p_{1}z + p_{2}z^{2} + ⋯ + p_{N−1}z^{N−1} .

The evaluation of P at powers of ω amounts to finding the inverse Fourier transform of the sequence p_{i} modulo N since

P (ω^{k}) = p_{0} + p_{1}ω^{k} + p_{2}ω^{2k} + ⋯ + p_{N−1}ω^{(N−1)k} .

From the equations

p̂_{0} = P (1)
= p_{0} + p_{1} + p_{2} + ⋯ p_{N−1}
p̂_{1} = P (ω)
= p_{0} + p_{1}ω + p_{2}ω^{2} + ⋯ p_{N−1}ω^{N−1}
…
p̂_{k} = P (ω^{k})
= p_{0} + p_{1}ω^{k} + p_{2}ω^{2k} + ⋯ p_{N−1}ω^{(N−1)k}
…

We can recover the original coefficients upon multiplying by various powers of ω and adding up:

p_{0} = N^{−1}(p̂_{0} + p̂_{1} + ⋯ + p̂_{N−1})
p_{1} = N^{−1}(p̂_{0} + p̂_{1}ω^{−1} + ⋯ + p̂_{N−1}ω^{−(N−1)})
…
p_{k} = N^{−1}(p̂_{0} + p̂_{1}ω^{−k} + ⋯ + p̂_{N−1}ω^{−(N−1)k})
…

We can think of this process as determining a new polynomial

P̂ = ∑ p̂_{i}x^{i}

from the original P. Let E be the evaluation map P ↦ P̂. What the formulas for inversion mean is that applying E twice takes P into a polynomial N P^{γ} = N \overline{P} where P^{γ}(z) = P (z^{−1}) for all N-th roots of unity z.

Now suppose N to be a power of 2, and let ω be a primitive N-th root of unity. Since N is a power of 2 this is equivalent to the condition that ω^{N/2} = −1.

We begin the fast version of the Fourier transform by rewriting P as

P (z) = P_{0}(z^{2}) + P_{1}(z^{2})z

where

\begin{aligned}P_{0}(z) & = p_{0} + p_{2}z + ⋯ + p_{N−2}z^{M−1} \\ P_{1}(z) & = p_{1} + p_{3}z + ⋯ + p_{N−1}z^{M−1} \end{aligned}

and M = N/2 is also a power of two. Since ω^{M} = −1 we have

\begin{matrix}P (ω^{k}) = P_{0}(ω^{2k}) + P_{1}(ω^{2k})ω^{k} \\ P (ω^{k+M}) = P_{0}(ω^{2k}) − P_{1}(ω^{2k})ω^{k} \end{matrix}

where we may assume 0 ≤ k < M. This can be interpreted as saying that the array of N items P (ω^{k}) (k = 0, …, N − 1) can be calculated from the array of items P_{0}((ω^{2})^{k}), P_{1}((ω^{2})^{k}) (k = 0, …, N/2 − 1) by means of certain simple linear combinations. The point now is that we can repeat this procedure with the polynomials P_{0} and P_{1} of degree one-half that of P, and with ω^{2}, which is a primitive M-th root of unity, instead of ω.

This means we write

\begin{aligned}P_{0}(z) & = P_{00}(z^{2}) + P_{01}(z^{2})z \\ P_{1}(z) & = P_{10}(z^{2}) + P_{11}(z^{2})z \end{aligned}

where

\begin{aligned}P_{00}(z) & = p_{0} + p_{4}z + ⋯ \\ P_{01}(z) & = p_{2} + p_{6}z + ⋯ \\ P_{10}(z) & = p_{1} + p_{5}z + ⋯ \\ P_{11}(z) & = p_{3} + p_{7}z + ⋯ \end{aligned}

We arrive finally at polynomials of degree 0 which are coefficients of P. For example if

P (z) = p_{0} + p_{1}z + p_{2}z^{2} + p_{3}z^{3} + p_{4}z^{4} + p_{5}z^{5} + p_{6}z^{6} + p_{7}z^{7}

then we get the polynomials

\begin{aligned}P_{0} & = p_{0} + p_{2}z + p_{4}z^{2} + p_{6}z^{3} \\ P_{1} & = p_{1} + p_{3}z + p_{5}z^{2} + p_{7}z^{3} \\ P_{00} & = p_{0} + p_{4}z \\ P_{01} & = p_{2} + p_{6}z \\ P_{10} & = p_{1} + p_{5}z \\ P_{11} & = p_{3} + p_{7}z \\ P_{000} & = p_{0} \\ P_{001} & = p_{4} \\ P_{010} & = p_{2} \\ P_{011} & = p_{6} \\ P_{100} & = p_{1} \\ P_{101} & = p_{5} \\ P_{110} & = p_{3} \\ P_{111} & = p_{7} \end{aligned}

The rule for assigning the final set of constants is bit inversion: if N = 2^{n} then we get constants P_{s} where s = s_{0}s_{1} … s_{n−1} is a sequence of bits and P_{s} = p_{t} where t is the integer in the range [0, 2^{n}) whose binary expression in bits is the reverse of the sequence s. From the sequence of these constants we get a sequence of polynomials of degree one by linear combinations with coefficients 1 and −1. From these in turn we get polynomials of degree three by linear combinations with coefficients 1, i, and so on.

In the above example we start with these constants and eventually fill in in this tableau:

\begin{aligned}P_{000} & = p_{0} & P_{00}(1) & = P_{000} + P_{001} & P_{0}(1) & = P_{00}(1) + P_{01}(1) & P (1) & = P_{0}(1) + P_{1}(1) \\ P_{001} & = p_{4} & P_{00}(−1) & = P_{000} − P_{001} & P_{0}(i) & = P_{00}(−1) + P_{01}(−1)i & P (ω) & = P_{0}(i) + P_{1}(i)ω \\ P_{010} & = p_{2} & P_{01}(1) & = P_{010} + P_{011} & P_{0}(−1) & = P_{00}(1) − P_{01}(1) & P (ω^{2}) & = P_{0}(−1) + P_{1}(−1)i \\ P_{011} & = p_{6} & P_{01}(−1) & = P_{010} − P_{011} & P_{0}(−i) & = P_{00}(−1) − P_{01}(−1)i & P (ω^{3}) & = P_{0}(−i) + P_{1}(−i)ω^{3} \\ P_{100} & = p_{1} & P_{10}(1) & = P_{100} + P_{101} & P_{1}(1) & = P_{10}(1) + P_{11}(1) & P (−1) & = P_{0}(1) − P_{1}(1) \\ P_{101} & = p_{5} & P_{10}(−1) & = P_{100} − P_{101} & P_{1}(i) & = P_{10}(−1) + P_{11}(−1)i & P (−ω) & = P_{0}(i) − P_{1}(i)ω \\ P_{110} & = p_{3} & P_{11}(1) & = P_{110} + P_{111} & P_{1}(−1) & = P_{10}(1) − P_{11}(1) & P (−ω^{2}) & = P_{0}(−1) − P_{1}(−1)i \\ P_{111} & = p_{7} & P_{11}(−1) & = P_{110} − P_{111} & P_{1}(−i) & = P_{10}(−1) − P_{11}(−1)i & P (−ω^{3}) & = P_{0}(−i) − P_{1}(−i)ω^{3} \end{aligned}

We can do this in two steps, first applying the permutation from the array p_{i} to p_{σ(i)}, where σ is bit inversion, and then applying the linear transformations directly to the array y_{i} = p_{σ(i)}. In the first step we transform the polynomial P to the polynomial P^{σ} with permuted coefficients, and then apply a map I call F to this permuted (inverted? perverted?) polynomial.

We can summarize what we know so far:

E = Fσ, EE = γ

where

\begin{aligned}EP_{i} & = P (ω^{i}) \\ P^{γ}(ω^{i}) & = P (ω^{−i}) \\ P_{i}^{σ} & = P_{σ(i)} \end{aligned}

and F is the map we obtain from the process exhibited in the tableau above. It is best described in terms of the algorithm producing it (y is an array of length N):

F (y, N)

k := 1

while k < N (being in turn 1, 2, …, N/2) w = ω_{2k}

a := 1
for j = 0 to k − 1
(a = w^{j})
i := j
while i < N
y_{i}, y_{i+k} := y_{i} + ay_{i+k}, y_{i} − ay_{i+k}
i := i + 2k
a := aw
k := 2k

The inverse to this map can also be constructed very simply by using the identity

[ \begin{matrix}1 & ζ \\ 1 & −ζ \end{matrix} ]^{−1} = [ \begin{matrix}1 & 1 \\ ζ^{−1} & −ζ^{−1} \end{matrix} ] .

Here is explicit code:

F^{∗}(y, N)

k = N/2

while k > 0
a := 1
w = ω_{k}^{−1}
for j = 0 to k − 1
(a = w^{j})
i = j
while i < N
y_{i}, y_{i+k} = (1/2)[y_{i} + y_{i+k}], (1/2)[a(y_{i} + y_{i+k})]
i = i + 2k
a := aw
k = k/2

The map F can only be described explicitly in terms of the effect of applying it after permutation of coefficients. The map F^{∗} however has a simple direct interpretation which plays a fundamental role in fast multiplication. If we apply F^{∗} to the coefficients of a polynomial P we get the sequence of evaluations P (ω^{σ(i)}). This follows from the formulas above relating the composition E^{2} of evaluation with itself, bit inversion σ, and conjugation γ.

What this means in practice is that we can forget about bit-inversion in order to do fast multiplication. In order to calculate products of polynomials P and Q, say of degree N − 1 with N a power of 2, we (1) calculate the array of evaluations P (ω^{−σ(i)}), Q(ω^{−σ(i)}) at the 2N-th roots of unity (in permuted order) by applying F^{∗}; (2) multiply the evaluations; (3) apply the inverse interpolation F of degree 2N; and then (4) multiply through by 1/2N.

The procedure of evaluation works in any ring at all containing primitive 2N-th roots of unity. Inversion is possible whenever we can divide by two in the ring.

Nussbaumer’s algorithm

Nussbaumer made two valuable observations (implicitly made earlier by Schönhage):

Proposition 1. If M is a power of 2 then in the polynomial ring

A_{M} = ℤ[z]/(z^{M} + 1)

the image of z is a primitive 2M-th root of unity.

This suggests that we can create roots of unity out of thin air!

This is simple, since z^{M} = −1 in this ring.

Proposition 2. If N = MR then the ring A_{N} is isomorphic to A_{R}[Z]/(Z^{M} − z).

This is because any polynomial in A_{N}

P = ∑_{0}^{MR−1} p_{i}z^{i}

can be rewritten as

∑_{i=0}^{M−1} P_{i}(z^{M})z^{i}, P_{i}(y) = ∑_{j=0}^{R−1} p_{jM+i}z^{Mj} = ∑_{j=0}^{R−1} p_{jM+i}y^{j}

where y = z^{M}. We now introduce a change of variables. The isomorphism we are looking for maps P (z) to the polynomial ∑ P_{i}(z)Z^{i} in Z with coefficients P_{i}(z) in

ℤ[z]/(z^{R} + 1) = A_{R} .

As a consequence of the first result, we know an efficient way to compute products in the ring A_{N}[z] in terms of products in A_{N}, by means of evaluation at 2N-th roots of unity, or any 2M-th root of unity for M|N. This will be especially efficient since multiplication by a root of unity just amounts to a shift.

The second result lets us proceed by recursion. Let N be a power of two. We set M = ⌊N/2⌋, R = ⌈N/2⌉. Then N = MR and M|R. In fact R/M = 1 or 2. Suppose we want to multiply two polynomials P and Q in A_{N}. We identify them with polynomials \overline{P} and \overline{Q} in A_{R}[Z]/(Z^{M} − z), where A_{R} = ℤ[z]/(z^{R} + 1). We can find the product in A_{R}[Z] by means of the Fourier evaluation & interpolation process, the one of degree 2M since the element z^{R/M} will be a 2M-th root of unity in A_{R}, applying recursion to A_{R}. We can then turn the product back into a polynomial in z by substituting z^{M} for z and z for Z.

Implementation

There are are several tricky points in any practical implementation. One of the most interesting is how to handle overflow in the arithmetic. We know that if we multiply polynomials whose coefficients are nonnegative and have n bits then in practical circumstances the product will have coefficients whicgh are also non-negative and 3n bits in size. But the various Fourier transforms in between may involve much larger integers which can even be negative. This is a relatively mild inflation, but requires some care in handling efficiently. The natural thing to do is perform calculations modulo 2^{3n} − 1, say, where size inflation cannot occur. An apparently better one is to use three separate moduli of one-word size. The advantage is that multiplying two three-word integers involves nine machine multiplications, whereas multiplying three single words for the smaller moduli involves only three. Using three moduli will involve some extra arithmetic at the beginning and end of the process, but this will be small if the moduli are chosen conveniently. In particlar conversion out of three-moduli arithmetic can be done easily using Garner’s method (p. 274 of Knuth, ACP II). Modulus arithmetic becomes particularly efficient if one uses moduli for which machine calculations can be done quickly. In practice this means moduli of the form 2^{m} − 1.

In my own programs n = 32, and the moduli I choose are 2^{32} − 1, 2^{31} − 1, and 2^{29} − 1.

For small moduli the Fourier transform involves what seems to be an inordinate amount of calculation. Therefore for multiplying small polynomails (of degree up to about seven) I use a kind of Karatsuba convolution.

Here are some sample times for my Sparc IPC to multiply two numbers of 2^{n} words of 32 bits (on a good day):

n  time in seconds
8  0.2
9  0.4
10  0.8
11  2.0
12  4.6
13  10.6
14  23.4
15  51.0