§1. Discrete Fourier Transform  Lecture I  Page 1

Lecture I BASIC ARITHMETIC

This lecture considers the arithmetic operations (addition, subtraction, multiplication and division) in three basic algebraic structures: polynomials, integers, matrices. These operations are the basic building blocks for other algebraic operations, and hence are absolutely fundamental in algorithmic algebra. Strictly speaking, division is only defined in a field. But there are natural substitutes in general rings: it could be always be replaced by the divisibility predicate. In a domain, we can define exact division (the exact division of u by v is carried out iff the v divides u). In case of Euclidean rings (Lecture II), it could be replaced by the quotient and remainder functions. The ring of polynomials over a field and the integers are Euclidean,

Complexity of Multiplication. In most algebraic structures of interest, the obvious algorithms for addition and subtraction take linear time and are easily seen to be optimal. Hence we have no more to say about them. As for the division-substitutes, they turn out to be reducible to multiplication. Hence the term “complexity of multiplication” can be regarded a generic term to cover such operations as well. After these considerations, what remains to be addressed is multiplication itself. The pervading influence of Schönhage and Strassen in all these results cannot be overstated.

We use some other algebraic structures in addition to the ones introduced in Lecture 0, §1:

GF (p^{m}) = Galois field of order p^{m}, p prime,
n  = integers modulo n ≥ 1,
M_{m,n}(R) = m by n matrices over a ring R,
M_{n}(R)  = M_{n,n}(R).

Finite structures such as GF (p^{m}) and _{n} have independent interest, but they also turn out to be critical for algorithms in infinite structures such as .

§1. The Discrete Fourier Transform

The key to fast multiplication of integers and polynomials is the discrete Fourier transform.

Roots of unity. In this section, we work with complex numbers. A complex number α ∈ is an nth root of unity if α^{n} = 1. It is a primitive nth root of unity if, in addition, α^{m} ≠ 1 for all m = 1,…, n − 1. In particular,

e^{\frac{2π}{n} ı} = cos \frac{2π}{n} + ı sin \frac{2π}{n}

(ı = \sqrt{−1}) is a primitive nth root of unity. There are exactly φ(n) primitive nth roots of unity where φ(n) is the number of positive integers less than or equal to n that are relatively prime to n (note that φ(n) is Euler’s phi-function).

Example: A primitive 8th root of unity is ω = e^{\frac{2π}{8} ı} = \frac{1}{\sqrt{2}} + ı\frac{1}{\sqrt{2}}. It is easy to check the only other primitive roots are ω^{3}, ω^{5} and ω^{7} (so φ(8) = 4). [These roots are easily visualized in the complex plane.]

Let ω denote any primitive nth root of unity. We note a basic identity.

Lemma 1 (Cancellation Property) If s ≢ 0 mod n then ∑_{j=0}^{n−1} ω^{js} = 0.

To see this, consider the identity x^{n} − 1 = (x − 1)(∑_{j=0}^{n−1} x^{j}). Substituting x = ω^{s} makes the left-hand side equal to zero. The right-hand side becomes (ω^{s} − 1)(∑_{j=0}^{n−1} ω^{js}). Since ω^{s} ≢ 1 for s ≢ 0 mod n, the result follows.

Let F (ω) = F_{n}(ω) denote the matrix

( \begin{matrix}1 & 1 & 1 & ⋯ & 1 \\ 1 & ω & ω^{2} & ⋯ & ω^{n−1} \\ 1 & ω^{2} & ω^{4} & ⋯ & ω^{2(n−1)} \\ ⋮ & & & & ⋮ \\ 1 & ω^{n−1} & ω^{2(n−1)} & ⋯ & ω^{(n−1)^{2}} \end{matrix} ) .

Definition 1 (The DFT and its inverse) Let a = (a_{0},…, a_{n−1})^{T} ∈ ^{n}. The discrete Fourier transform (abbr. DFT) of a is DFT_{n}(a) := A = (A_{0},…, A_{n−1})^{T} where A_{i} = ∑_{j=0}^{n−1} a_{j}ω^{ij}, for i = 0,…, n − 1. Thus,

DFT_{n}(a) = F (ω) · ( a_{0} ; a_{1} ; ⋮ ; a_{n−1} ) = ( A_{0} ; A_{1} ; ⋮ ; A_{n−1} ) .

The inverse discrete Fourier transform of A = (A_{0},…, A_{n−1})^{T} is DFT_{n}^{−1}(A) = \frac{1}{n}F (ω^{−1}) · A. Thus,

DFT_{n}^{−1}(A) := \frac{1}{n} ( \begin{matrix}1 & 1 & 1 & ⋯ & 1 \\ 1 & ω^{−1} & ω^{−2} & ⋯ & ω^{−n+1} \\ ⋮ & & & & ⋮ \\ 1 & ω^{−n+1} & ω^{−2(n−1)} & ⋯ & ω^{−(n−1)^{2}} \end{matrix} ) · ( A_{0} ; A_{1} ; ⋮ ; A_{n−1} ) .

Note that ω^{−1} = ω^{n−1}. We will omit the subscript ‘n’ in DFT_{n} when convenient. The following shows that the two transforms are indeed inverses of each other:

Lemma 2 We have F (ω^{−1}) · F (ω) = F (ω) · F (ω^{−1}) = nI_{n} where I_{n} is the identity matrix.

Proof. Let F (ω^{−1}) · F ω) = [c_{j,k}]_{j,k=0}^{n−1} where

c_{j,k} = ∑_{i=0}^{n−1} ω^{−ji}ω^{ik} = ∑_{i=0}^{n−1} ω^{i(j−k)}.

If j = k, then c_{j,k} = ∑_{i=0}^{n−1} ω^{0} = n. Otherwise, −n < j − k < j − k < n and j − k ≠ 0 implies c_{j,k} = 0, using lemma 1. Similarly, F (ω) · F (ω^{−1}) = nI_{n}. Q.E.D.

Connection to polynomial evaluation and interpolation. Let a be the coefficient vector of the polynomial P (X) = ∑_{i=0}^{n−1} a_{i}X^{i}. Then computing DFT(a) amounts to evaluating the polynomial P (X) at all the nth roots of unity, at

X = 1, X = ω, X = ω^{2},…, X = ω^{n−1}.

Similarly, computing DFT^{−1}(A) amounts to recovering the polynomial P (X) from its values (A_{0},…, A_{n−1}) at the same n points. In other words, the inverse discrete Fourier transform interpolates, or reconstructs, the polynomial P (X) from its values at all the n roots of unity. Here we use the fact (Lecture IV.1) that the interpolation of a degree n − 1 polynomial from its values at n distinct points is unique. (Of course, we could also have viewed DFT as interpolation and DFT^{−1} as evaluation.)

The Fast Fourier Transform. A naive algorithm to compute DFT and DFT^{−1} would take Θ(n^{2}) complex arithmetic operations. In 1965, Cooley and Tukey discovered a method that takes O(n log n) operations. This has come to be known^{1} as the fast Fourier transform (FFT). This algorithm is widely used. Actually, the basic ideas of the FFT were known in some form prior to 1965 (e.g., Runge and König, 1924).

Let us now present the FFT algorithm to compute DFT(a) where a = (a_{0},…, a_{n−1}). In fact, it is a fairly straightforward divide-and-conquer algorithm. To simplify discussion, let n be a power of 2. It is convenient to be able to interchangeably talk of the polynomial P (X) of degree at most n − 1 whose vector of coefficients is a. As noted, computing DFT(a) amounts to computing the n values

P (1), P (ω), P (ω^{2}),…, P (ω^{n−1}).  (1)

First, let us express P (X) as the sum of its odd part and its even part:

P (X) = P_{e}(X^{2}) + XP_{o}(X^{2})

where P_{e}(Y), P_{o}(Y) are polynomials of degrees at most \frac{n}{2} and \frac{n−1}{2}, respectively. E.g., for P (X) = 3X^{6} − X^{4} + 2X^{3} + 5X − 1, we have P_{e}(Y) = 3Y^{3} − Y^{2} − 1, P_{o}(Y) = 2Y + 5. Thus we have reduced the problem of computing the values in (1) to the following:

FFT Algorithm:

Input: a polynomial P (X) with coefficients given by an n-vector a, and ω, a primitive nth root of unity.

Output: DFT_{n}(a).

1. Evaluate P_{e}(X^{2}) and P_{o}(X^{2}) at X^{2} = 1, ω^{2},…, ω^{n−2}, ω^{n}, ω^{n+2},…, ω^{2n−2}.

2. Multiply P_{o}(ω^{2j}) by ω^{j} for j = 0,…, n − 1.

3. Add P_{e}(ω^{2j}) to ω^{j}P_{o}(ω^{2j}), for j = 0,…, n − 1.

1This is, of course, a misnomer since the FFT is not a transform but an algorithm.

§2. Polynomial Multiplication  Lecture I  Page 4

Analysis. Note that in step 1, we have 1 = ω^{n}, ω^{2} = ω^{n+2},…, ω^{n−2} = ω^{2n−2}. So it suffices to evaluate P_{e} and P_{o} at only n/2 values, X = 1, ω^{2},…, ω^{n−2}, i.e., at all the (n/2)th roots of unity. But this is equivalent to the problem of computing DFT_{n/2}(P_{e}) and DFT_{n/2}(P_{o}). We accomplish this recursively. Steps 2 and 3 takes n multiplications and n additions respectively. Overall, if T (n) is the number of complexity additions and multiplications, we have

T (n) = 2T (n/2) + 2n

which has the solution, T (n) = 2n log n.

Since the same method can be applied to the inverse discrete Fourier transform, we have shown:

Theorem 3 (Complexity of FFT) Assuming the availability of a primitive nth root of unity, the discrete Fourier transform DFT_{n} and its inverse can be computed in O(n log n) complex arithmetic operations.

Exercise 1.1: Show that the number of multiplications in step 2 can be reduced to n/2. HINT: Then half of the additions in step 3 becomes subtractions.

Remark: There are several closely related fast transform methods which has the same framework. See “A unified treatment of discrete fast unitary transforms”, B.J Fino and V.R. Algazi (SIAM J.Comp.,6:4(1977)700–717).

§2. Polynomial Multiplication

We consider the multiplication of complex polynomials. To exploit the FFT algorithm, we make a fundamental connection.

Convolution and polynomial multiplication. Assume n ≥ 2 is even. The convolution of two n-vectors a = (a_{0},…, a_{n−1})^{T} and b = (b_{0},…, b_{n−1})^{T} is the n-vector

c = a ∗ b :=(c_{0},…, c_{n−1})^{T}

where c_{i} = ∑_{j=0}^{i} a_{j}b_{i−j}. Let P (X) and Q(X) be polynomials of degrees less than n/2. Then R(X) := P (X)Q(X) is a polynomial of degree less than n−1. Let a and b denote the coefficient vectors of P and Q (padded out with initial zeroes to make vectors of length n). Then it is not hard to see that a ∗ b gives the coefficient vector of R(X). Thus convolution is essentially polynomial multiplication.

Theorem 4 (Convolution Theorem) Let a, b be n-vectors (n even) whose initial n/2 entries are zeroes. Then

DFT^{−1}(DFT(a) · DFT(b)) = a ∗ b.

§3. Modular FFT  Lecture I  Page 5

Proof. Suppose DFT(a) = (A_{0},…, A_{n−1})^{T} and DFT(b) = (B_{0},…, B_{n−1})^{T}. Let C = (C_{0},…, C_{n−1})^{T} where C_{i} = A_{i}B_{i}. From the evaluation interpretation of DFT, it follows that C_{i} is the value of the polynomial R(X) = P (X)Q(X) at X = ω^{i}. Hence DFT^{−1}(C) is the interpolation of polynomial P (X)Q(X) from its values at all the nth roots of unity. This reconstructs R(X) correctly (Lecture IV§1). So C is the coefficient vector of R(X). But the coefficient vector of R(X) is just the convolution of a and b. Q.E.D.

This theorem reduces the problem of convolution (equivalently, polynomial multiplication) to two DFT and one DFT^{−1} computations. We immediately conclude from the FFT result (Theorem 3):

Theorem 5 (Algebraic complexity of polynomial multiplication) Assuming the availability of a primitive nth root of unity, we can compute the product P Q of two polynomials P, Q ∈ [X] of degrees less than n in O(n log n) complex operations.

Remark: If the coefficients of our polynomials are not complex numbers but some other ring, then a similar result holds provided the ring contains an analogue to the roots of unity. Such a situation arises in our next section.

Exercise 2.1: Show that the quotient P div Q and remainder P mod Q can be computed in O(n log n) complex operations.

Exercise 2.2: Let q = p^{m} where p ∈ is prime, m ≥ 1. Show how to multiply in GF (q) in O(mL(m)) operations of _{p} and compute inverses in O(mL^{2}(m)) operations of _{p}. HINT: use the fact that GF (q) is isomorphic to GF (p)[X]/(F (X)) where F (X) is any polynomial of degree m that is irreducible over GF (p).

Exercise 2.3: Let q = p^{m} as above. Show how to multiply two degree n polynomials over GF (q) in O(nL^{2}(n)) operations of GF (q). and compute the GCD of two such polynomials in O(nL^{2}(n)) operations of GF (q).

§3. Modular FFT

To extend the FFT technique to integer multiplication, a major problem to overcome is how one replaces the complex roots of unity with some discrete analogue. One possibility is to carry out the complex arithmetic to a suitable degree of accuracy. This was done by Strassen in 1968, achieving a time bound that satisfies the recurrence T (n) = O(nT (log n)). For instance, this implies T (n) = O(n log n(log log n)^{1+ϵ}) for any ϵ > 0. In 1971, Schönhage and Strassen managed to improved this to T (n) = O(n log n log log n). While the complexity improvement can be said to be strictly of theoretical interest, their use of modular arithmetic to avoid approximate arithmetic has great interest. They discovered that the discrete Fourier transform can be defined, and the FFT efficiently implemented, in _{M} where

M = 2^{L} + 1,  (2)

for suitable values of L. This section describes these elegant results.

First, we make some general remarks about _{M} for an arbitrary modulus M > 1. An element x ∈ _{M} is a zero-divisor if there exists y ≠ 0 such that x · y = 0; a (multiplicative) inverse of x is y such that xy = 1. For example, in _{4}, the element 2 has no inverse and 2 · 2 = 0. We claim:

An element x ∈ _{M} has a multiplicative inverse (denoted x^{−1}) if and only if x is not a zero-divisor.

To see this claim, suppose x^{−1} exists and x · y = 0. Then y = 1 · y = x^{−1}x · y = 0. Conversely, if x is not a zero-divisor then the elements in the set {x · y : y ∈ _{M}} are all distinct because if x · y = x · y′ then x(y − y′) = 0 and y − y′ ≠ 0, contradiction. Hence, by pigeon-hole principle, 1 occurs in the set. This proves our claim. We have two basic consequences: (i) If x has an inverse, the inverse is unique. [In proof, if x · y = 1 = x · y′ then x(y − y′) = 0 and so y = y′.] (ii) _{M} is a field iff M is prime. [In proof, if M has the proper factorization xy then x is a zero-divisor. Conversely, if M is prime then every x ∈ _{M} has an inverse because the extended Euclidean algorithm implies there exist s, t ∈ _{M} such that sx+tM = 1, i.e., s = x^{−1}(mod M).]

In the rest of this section and also the next one, we assume M has the form in Equation (2). Then 2^{L} ≡ −1(mod M) and 2^{2L} = (M − 1)^{2} ≡ 1(mod M). We also use the fact that every element of the form 2^{i} (i ≥ 0) has an inverse in _{M}, viz., 2^{2L−i}.

Representation and basic operations modulo M. We clarify how numbers in _{M} are represented. Let 2^{L} ≡ −1(mod M) be denoted with the special symbol \overline{1}. We represent each element of _{M} − \{\overline{1}\} in the expected way, as a binary string (b_{L−1},…, b_{0}) of length L; the element \overline{1} is given a special representation. It is relatively easy to add and subtract in _{M} under this represention using a linear number of bit operations, i.e., O(L) time. Of course, special considerations apply to \overline{1}.

Exercise 3.1: Show that addition and subtraction take O(L) bit operations.

We will also need to multiply by powers of 2 in linear time. Intuitively, multiplying a number X by 2^{j} amounts to left-shifting the string X by j positions; a slight complication arises when we get a carry to the left of the most significant bit.

Example: Consider multiplying 13 = (1101) by 2 = (0010) in _{17}. Left-shifting (1101) by 1 position gives (1010), with a carry. This carry represents 16 ≡ −1 = \overline{1}. So to get the final result, we must add \overline{1} (equivalently, subtract 1) from (1010), yielding (1001). [Check: 13 × 2 ≡ 9(mod 17) and 9 = (1001).]

In general, if the number represented by the string (b_{L−1},…, b_{0}) is multiplied by 2^{j} (0 < j < L), the result is given as a difference:

(b_{L−j−1}, b_{L−j−2},…, b_{0}, 0,…, 0) − (0,…, 0, b_{L−1}, b_{L−2},…, b_{L−j}).

But we said that subtraction can be done in linear time. So we conclude: in _{M}, multiplication by 2^{j} takes O(L) bit operations.

Primitive roots of unity modulo M. Let K = 2^{k} and K divides L. We define

ω:= 2^{L/K}.

For instance, in _{17} we have L = 4. With K = 2, we get ω^{i} = 4, 16, 13, 1. So ω is a primitive 4th root of unity.

Lemma 6 In _{M}, ω is a primitive (2K)th root of unity.

Proof. Note that ω^{K} = 2^{L} ≡ −1(mod M). Thus ω^{2K} ≡ 1(mod M), i.e., it is a (2K)th root of unity. To show that it is in fact a primitive root, we must show ω^{j} ≢ 1 for j = 1,…, (2K − 1). If j ≤ K then ω^{j} = 2^{Lj/K} ≤ 2^{L} < M so clearly ω^{j} ≢ 1. If j > K then ω^{j} = −ω^{j−K} where j − K = 1,…, K − 1. Again, ω^{j−K} < 2^{L} ≡ −1 and so −ω^{j−K} ≢ 1. Q.E.D.

We next need the equivalent of the cancellation property (Lemma 1). The original proof is invalid since _{M} is not necessarily an integral domain.

Lemma 7 The cancellation property holds: for any integer s not divisible by 2K,

∑_{j=0}^{2K−1} ω^{js} ≡ 0(mod M).

Proof. Let (s mod 2K) = 2^{p}q where q is odd, 0 < 2^{p} < 2K and let r = 2K · 2^{−p} > 1. Then

∑_{j=0}^{2K−1} ω^{js} = ∑_{j=0}^{r−1} ω^{js} + ∑_{j=r}^{2r−1} ω^{js} + ⋯ + ∑_{j=2K−r}^{2K−1} ω^{js}

= ∑_{j=0}^{r−1} ω^{js} + ω^{rs} ∑_{j=0}^{r−1} ω^{js} + ⋯ + ω^{rs(2^{p}−1)} ∑_{j=0}^{r−1} ω^{js}

= 2^{p} ∑_{j=0}^{r−1} ω^{js},

since ω^{rs} ≡ 1. Note that ω^{rs/2} = ω^{Kq} ≡ (−1)^{q} = −1. The lemma follows since

∑_{j=0}^{r−1} ω^{js} = ∑_{j=0}^{\frac{r}{2} −1} (ω^{sj} + ω^{s(j+ \frac{r}{2} )}) = ∑_{j=0}^{\frac{r}{2} −1} (ω^{sj} − ω^{sj}) = 0.

Q.E.D.

Using ω, we define the discrete Fourier transform and its inverse in _{M} as usual: DFT_{2K}(a) := F (ω) · a and DFT_{2K}^{−1}(A) := \frac{1}{2K} F (ω^{−1}) · A. To see that the inverse transform is well-defined, we should recall that \frac{1}{2K} and ω^{−1} both exist. Our proof that DFT and DFT^{−1} are inverses (Lemma 2) goes through. We obtain the analogue of Theorem 3:

Theorem 8 The transforms DFT_{2K}(a) and DFT_{2K}^{−1}(A) for (2K)-vectors a, A ∈ ( _{M})^{2K} can be computed using the Fast Fourier Transform method, taking O(KL log K) bit operations.

§4. Integer Multiplication  Lecture I  Page 8

Proof. We use the FFT method as before (refer to the three steps placed in a box in §1). View a as the coefficient vector of the polynomial P (X). Observe that ω^{2} is a primitive Kth root of unity in _{M}. Of course, we can compute ω^{2} easily from ω in O(L/K) bit operations. This allows us to implement step 1 recursively, by calling DFT_{K} twice, once on the even part P_{e}(Y) and again on the odd part P_{o}(Y). In step 2, we need to compute ω^{j} (which is easy) and multiply it to P_{o}(ω^{2j}) (also easy), for j = 0,…, 2K − 1. Step 2 takes O(KL) bit operations. Finally, we need to add ω^{j}P_{o}(ω^{2j}) to P_{e}(ω^{2j}) in step 3. This also takes O(KL) bit operations. Thus the overall number of bit operations T (2K) satisfies the recurrence

T (2K) = 2T (K) + O(KL)

which has solution T (2K) = O(KL log K), as claimed. Q.E.D.

Remarks: It is not hard to show (exercise below) that if M is prime then L is a power of 2. Generally, a number of the form 2^{2^{n}} + 1 is called Fermat number. The first 4 Fermat numbers are prime which led to Fermat to the rather unfortunate conjecture that they all are. No other primes have been discovered so far and many are known to be composite (Euler discovered in 1732 that the 5th Fermat number 2^{2^{5}} + 1 is divisible by 641). Fermat numbers are closely related to a more fortunate conjecture of Mersenne, that all numbers of the form 2^{p} − 1 are prime (where p is prime): although the conjecture is false, at least there is more hope that there are infinitely many such primes.

Exercise 3.2: (a) If a^{L} + 1 is prime where a ≥ 2, then a is even and L is a power of two. (b) If a^{L} − 1 is prime where L > 1, then a = 2 and L is prime.

Exercise 3.3: Show that Strassen’s recurrence T (n) = n · T (log n) satisfies

T (n) = O ((∏_{i=0}^{k−1} log^{(i)} n) (log^{(k)} n)^{1+ϵ})  (3)

for any k < log^{∗}(n).

Exercise 3.4: The “original” fast algorithm for integer multiplication uses the fact that if U = 2^{L}U_{0} + U_{1} and V = 2^{L}V_{0} + V_{1} where U_{i}, V_{i} are L-bit numbers, then W = UV = 2^{2L}U_{0}V_{0} + 2^{L}(U_{0}V_{1} + U_{1}V_{0}) + U_{1}V_{1}, which we can rewrite as 2^{2L}W_{0} + 2^{L}W_{1} + W_{2}. But if we compute (U_{0} + U_{1})(V_{0} + V_{1}), W_{0}, W_{2}, we also obtain W_{1}. Show that this leads to a

time bound of T (n) = O(n^{lg 3}).

§4. Fast Integer Multiplication

The following result of Schönhage and Strassen [3] is perhaps the fundamental result of the algorithmic algebra.

Theorem 9 (Complexity of integer multiplication) Given two integers u, v of sizes at most n bits, we can form their product uv in O(n log n log log n) bit-operations.

For simplicity, we prove a slightly weaker version of this result, obtaining a bound of O(n log^{2.6} n) instead.

A simplified Schönhage-Strassen algorithm. Our goal is to compute the product W of the positive integers U, V. Assume U, V are N-bit binary numbers where N = 2^{n}. Choose K = 2^{k}, L = 3 · 2^{ℓ} where

k:= ⌊\frac{n}{2}⌋ , ℓ:= ⌈\frac{n}{2}⌉ .

We breakup U, V into K pieces, each of size at most 2^{ℓ} bits; this gives us the (2K)-vectors:

\overline{U} = (0,…, 0, U_{K−1},…, U_{0}), \overline{V} = (0,…, 0, V_{K−1},…, V_{0}),

where U_{j}, V_{j} are 2^{ℓ}-bit strings. Now regard \overline{U}, \overline{V} as the coefficient vectors of the polynomials P (X) = ∑_{j=0}^{K−1} U_{j}X^{j} and Q(X) = ∑_{j=0}^{K−1} V_{j}X^{j}. Let

\overline{W} = (W_{2K−1},…, W_{0})

be the convolution of \overline{U} and \overline{V}. Note that each W_{i} in \overline{W} satisfies the inequality

0 ≤ W_{i} ≤ K · 2^{2·2^{ℓ}}  (4)

since it is the sum of at most K products of the form U_{j}V_{i−j}. Hence

0 ≤ W_{i} < 2^{3·2^{ℓ}} < M

for sufficiently large input size N. In this case, \overline{W} is computed exactly in _{M}. It follows that \overline{W} is the coefficient vector of the product R(X) = P (X)Q(X). Since P (2^{2^{ℓ}}) = U and Q(2^{2^{ℓ}}) = V, it follows that R(2^{2^{ℓ}}) = UV = W. Hence

W = ∑_{i=0}^{2K−1} 2^{2^{ℓ}j}W_{j}.

We can easily obtain each summand in this sum from \overline{W} by multiplying each W_{j} with 2^{2^{ℓ}j}. As each W_{j} has k + 2 · 2^{ℓ} < L non-zero bits, we illustrate this summation as follows:

Figure 1. Illustrating forming the product W = UV.

From this figure we see that each bit of W is obtained by summing at most 3 bits plus at most 2 carry bits. Since W has at most 2N bits, we conclude:

Lemma 10 The product W can be obtained from \overline{W} in O(N) bit operations.

It remains to show how to compute \overline{W}. By the convolution theorem,

\overline{W} = DFT^{−1}(DFT(\overline{U}) · DFT(\overline{V})).

The three transforms takes O(KL log K) = O(N log N) bit operations (Theorem 8). The scalar product DFT(\overline{U})·DFT(\overline{V}) requires 2K multiplications of L-bit numbers, which is accomplished recursively. Thus, if T (N) is the bit-complexity of this algorithm, we obtain the recurrence

T (N) = O(N log N) + 2K · T (L).  (5)

Write t(n) := T (N)/N where N = 2^{n}. The recurrence becomes

t(n) = O(n) + 2\frac{K}{N} T (L)
= O(n) + 2 · \frac{3}{L}T (L)
= O(n) + 6 · t(\frac{n}{2} + c),

for some constant c. We shift the domain of t(n) by defining s(n) := t(n + 2c). Then s(n) = (n + 2c) + 6t((n/2) + 2c) = O(n) + 6s(n/2). This gives the solution s(n) = O(n^{lg 6}). Backsubstituting, we obtain

T (N) = O(N log^{α} N),  α = lg 6 < 2.5848.  (6)

Refinements. Our choice of L = 3 · 2^{ℓ} is suboptimal. A slight improvement (attributed to Karp in his lectures) is to compute each W_{i} (i = 0,…, 2K − 1) in two parts: let M′:= 2^{2·2^{ℓ}} + 1 and M′′:= K. Since M′, M′′ are relatively prime and W_{i} < M′M′′, it follows that if we have computed W′_{i}:= W_{i} mod M′ and W′′_{i}:= W_{i} mod M′′, then we can recover W_{i} using the Chinese remainder theorem (Lecture IV). It turns out that computing all the W′′_{i}’s and the reconstruction of W_{i} from W′_{i}, W′′_{i} can be accomplished in linear time. The computation of the W′_{i}’s proceeds exactly as the above derivation. The new recurrence we have to solve is

t(n) = n + 4t(n/2)

which has the solution t(n) = O(n^{2}) or T (N) = O(N log^{2} N). To obtain the ultimate result, we have to improve the recurrence to t(n) = n + 2t(n/2). In addition to the above ideas (Chinese remainder, etc), we must use a variant convolution called “negative wrapped convolution” and DF T_{K} instead of DF T_{2K}. Then W_{i}’s can be uniquely recovered.

Exercise 4.1: Carry out the outline proposed by Karp.

Integer multiplication in other models of computation. Recall that M_{B}(n) denotes the Turing complexity of multiplying two n-bit integers. It is tedious but possible to carry out the Schönhage-Strassen algorithm on a Turing machine, in the same time complexity. Thus we conclude

M_{B}(n) = O(n log n log log n) = nL(n).

§5. Matrix Multiplication  Lecture I  Page 11

This bound on M_{B}(n) can be improved for more powerful models of computation. Schönhage [2] has shown that linear time is sufficient on pointer machines. Using general simulation results, this translates to O(n log n) time on logarithmic-cost successor RAMs (Lecture I). In parallel models, O(log n) time suffices.

Extending the notation of M_{B}(n), let

M_{B}(m, n)

denote the Turing complexity of multiplying two integers of sizes (respectively) at most m and n bits. Thus, M_{B}(n) = M_{B}(n, n). It is straightforward to extend the bound on M_{B}(n) to M_{B}(m, n).

Exercise 4.2: Show that M_{B}(m, n) = max{m, n}L(min{m, n}).

Exercise 4.3: Show that we can take remainders u mod v and form quotients u div v of integers in the same bit complexity as multiplication.

Exercise 4.4: Show how to multiply in _{p} (p ∈ a prime) in bit complexity O(log pL(log p)), and form inverses in _{p} in bit complexity O(log pL^{2}(log p)).

§5. Matrix Multiplication

For arithmetic on matrices over a ring R, it natural that our computational model is algebraic programs over the base comprising the ring operations of R. Here the fundamental discovery by Strassen (1968) [4] that the standard algorithms for matrix multiplication is suboptimal started off intense research for over a decade in the subject. Although the final word is not yet in, rather substantial progress had been made. These results are rather deep and we only report the current record, due to Coppersmith and Winograd (1987) [1].

Proposition 11 (Algebraic complexity of matrix multiplication) The product of two matrices in M_{n}(R) can be computed of O(n^{α}) operations in the ring R, where α = 2.376. In other words,

MM(n) = O(n^{α}).

It is useful to extend this result to non-square matrix multiplication. Let MM(m, n, p) denote the number of ring operations necessary to compute the product of an m × n matrix by a n × p matrix. So MM(n) = MM(n, n, n).

Theorem 12 Let MM(n) = O(n^{α}) for some α ≥ 2. Then

MM(m, n, p) = O(mnp · k^{α−3})

where k = min{m, n, p}.

Proof. We extend the bound on multiplying square matrices to non-square matrices. Suppose A is a m × n matrix, B a n × p matrix. First assume m = p but n is arbitrary. Then the bound in our theorem amounts to:

MM(m, n, m) = \{ \begin{matrix}O(nm^{α−1}) & if m ≤ n \\ O(m^{2}n^{α−2}) & if n ≤ m. \end{matrix}

We prove this in two cases. Case: m ≤ n. We partition A into r = ⌈n/m⌉ matrices, A = [A_{1}|A_{2}| ⋯ |A_{r}] where each A_{i} is an m-square matrix except possibly for A_{r}. Similarly partition B into r square matrices, B^{T} = [B_{1}^{T} |B_{2}^{T} | ⋯ |B_{r}^{T}]. Then

AB = A_{1}B_{1} + A_{2}B_{2} + ⋯ + A_{r}B_{r}.

We can regard A_{r}B_{r} as a product of two m-square matrices, simply by padding out A_{r} and B_{r} with zeroes. Thus each A_{i}B_{i} can be computed in m^{α} operations. To add the products A_{1}B_{1},…, A_{r}B_{r} together, we use O(rm^{2}) = O(rm^{α}) addition operations. Hence the overall complexity of computing AB is O(rm^{α}) = O(nm^{α−1}) as desired. Case: n ≤ m. We similarly breakup the product AB into r^{2} products of the form A_{i}B_{j}, i, j = 1,…, r, r = ⌈m/n⌉. This has complexity O(r^{2}n^{α}) = O(m^{2}n^{α−2}). This completes the proof for the case m = p.

Next, since the roles of m and p are symmetric, we may assume m < p. Let r = ⌈p/m⌉. We have two cases: (1) If m ≤ n then MM(m, n, p) ≤ rMM(m, n, m) = O(pnm^{α−2}). (2) If n < m, then MM(m, n, p) ≤ rMM(m, n, m) = O(rm^{2}n^{α−2}) = O(pmn^{α−2}). Q.E.D.

Notice that this result is independent of any internal details of the O(n^{α}) matrix multiplication algorithm. Webb Miller (“Computational complexity and numerical stability”, SIAM J.Comp., 4:2(1975)97–107) has shown that under sufficient conditions for numerical stability, any algorithm for matrix multiplication over a ring requires n^{3} multiplications.

References

[1] Don Coppersmith and Shmuel Winograd. Matrix multiplication via arithmetic progressions. Journal of Symbolic Computation, 9:251–280, 1990. Extended Abstract: Symposium on Theory of Computing, Vol.19, 1987, pp.1-6.

[2] A. Schönhage. Storage modification machines. SIAM Journal Computing, 9:490–508, 1980.

[3] A. Schönhage and V. Strassen. Schnelle Multiplikation grosser Zahlen. Computing, 7:281– 292, 1971.

[4] Volker Strassen. Gaussian elimination is not optimal. Numerische Mathematik, 14:354–356, 1969.