On Multiplication: I · II · III · IV · V

The Fourier Transform: From Signals to Arithmetic

Why We Need a New Idea

At the end of the Toom-Cook story in Vol. I we noted a frustrating ceiling: as we increase the splitting parameter $k$, the exponent $\log_k(2k-1)$ drifts toward 1, but it never reaches it. Worse, the constant factor in the $O(n)$ additive work (evaluation and interpolation matrices of size $k \times k$, coefficient blowup at large evaluation points) grows so fast that no finite $k$ delivers practical gains beyond Toom-4 or so.

To break through this wall we need to abandon the strategy of evaluating at a handful of ad-hoc points and instead evaluate at a structured infinite family of points whose algebraic symmetry enables a radically faster algorithm. That family is the complex roots of unity, and the algorithm that exploits their symmetry is the Fast Fourier Transform.

Before we can state Schönhage and Strassen’s multiplication algorithm, we must build the FFT from the ground up. This requires four conceptual layers: (1) the observation that integer multiplication is polynomial convolution, (2) the Discrete Fourier Transform as an evaluation map, (3) the Cooley-Tukey decomposition that makes the DFT fast, and (4) the inverse transform that recovers coefficients from point values.

Multiplication Is Convolution

Polynomials and digit vectors

We have already seen that an $n$-digit integer in base $B$ can be written as a polynomial evaluated at $B$. Let us now make the multiplication side precise. Given two integers with digit vectors $\mathbf{a} = (a_0, a_1, \ldots, a_{n-1})$ and $\mathbf{b} = (b_0, b_1, \ldots, b_{n-1})$, their product is a new integer whose digit vector $\mathbf{c} = (c_0, c_1, \ldots, c_{2n-2})$ satisfies:

(with the convention that $a_j = 0$ for $j \geq n$ and similarly for $b_j$). This is exactly the definition of the linear convolution of the sequences $\mathbf{a}$ and $\mathbf{b}$, or equivalently, the coefficient vector of the product polynomial $A(x) \cdot B(x)$ where:

Computing this convolution directly requires computing each of the $2n - 1$ output coefficients, and for each we sum up to $n$ products. The total cost is $O(n^2)$ – schoolbook multiplication in disguise.

The Convolution Theorem

The single most important theorem in this entire story is:

Convolution Theorem. Let $\mathcal{F}$ denote the Discrete Fourier Transform (defined below). If $\mathbf{c} = \mathbf{a} \ast \mathbf{b}$ is the convolution of two sequences, then:

\[\mathcal{F}(\mathbf{a} \ast \mathbf{b}) = \mathcal{F}(\mathbf{a}) \cdot \mathcal{F}(\mathbf{b})\]

where $\cdot$ denotes pointwise (component-by-component) multiplication.

In words: convolution in the “coefficient domain” becomes pointwise multiplication in the “frequency domain.” Since pointwise multiplication of two length-$N$ vectors costs only $O(N)$, the entire cost of multiplication reduces to the cost of two forward transforms and one inverse transform. If each transform costs $O(N \log N)$, the total cost is $O(N \log N)$ – an exponential improvement over $O(N^2)$.

This is not a hand-wave. Let us now build the DFT rigorously and prove why it can be computed in $O(N \log N)$.

The Discrete Fourier Transform

Complex exponentials and Euler’s formula

Before defining the DFT, we need the language of complex numbers. Recall Euler’s formula:

This elegant identity tells us that the complex exponential $e^{i\theta}$ traces out the unit circle in the complex plane as $\theta$ varies from $0$ to $2\pi$. Every point on the unit circle can be written as $e^{i\theta}$ for some angle $\theta$.

The $N$-th roots of unity

Fix a positive integer $N$. The $N$-th roots of unity are the $N$ complex numbers that satisfy $z^N = 1$. They are:

These are $N$ points spaced equally around the unit circle, like the vertices of a regular $N$-gon inscribed in the circle $\lvert z \rvert = 1$. The primitive $N$-th root of unity is:

so that $\omega_N^k = (\omega_N)^k$. We will usually drop the subscript $N$ when the context is clear.

The roots of unity possess remarkable algebraic properties that are the engine of the FFT:

1. Periodicity: $\omega^{k+N} = \omega^k$ for all $k$. The roots cycle with period $N$.

2. Cancellation (Half-turn symmetry): $\omega^{k + N/2} = -\omega^k$ when $N$ is even. Geometrically, the point diametrically opposite $\omega^k$ on the unit circle is $-\omega^k$. This is the single most important property for the FFT.

3. Summation: $\displaystyle\sum_{k=0}^{N-1} \omega^{jk} = \begin{cases} N & \text{if } N \mid j \cr 0 & \text{otherwise} \end{cases}$

   This orthogonality relation is what makes the inverse DFT work.

4. Squaring (Halving): If $N$ is even, then $\lbrace (\omega_N^k)^2 : k = 0, \ldots, N-1 \rbrace$ gives exactly the $N/2$-th roots of unity, each appearing twice. That is, $(\omega_N)^2 = \omega_{N/2}$.

Definition of the DFT

Given a vector $\mathbf{x} = (x_0, x_1, \ldots, x_{N-1})$, its Discrete Fourier Transform is the vector $\mathbf{X} = (X_0, X_1, \ldots, X_{N-1})$ defined by:

Equivalently, $X_k = P(\omega^k)$ where $P(z) = \sum_{j} x_j z^j$ is the polynomial whose coefficients are the entries of $\mathbf{x}$. The DFT is nothing more than evaluating the polynomial at all $N$-th roots of unity simultaneously.

The DFT matrix

We can express the DFT as a matrix-vector product $\mathbf{X} = F_N \, \mathbf{x}$ where the DFT matrix $F_N$ has entries:

Explicitly, for $N = 4$ with $\omega = \omega_4 = e^{2\pi i/4} = i$:

Computing $\mathbf{X} = F_N \mathbf{x}$ by brute-force matrix-vector multiplication costs $O(N^2)$. The entire point of the FFT is to exploit the structure of $F_N$ to compute this product in $O(N \log N)$.

The Inverse DFT

The orthogonality of the roots of unity (property 3 above) guarantees that the DFT is invertible. The inverse is:

or in matrix form, $\mathbf{x} = \frac{1}{N} F_N^{-1} \mathbf{X}$ where $F_N^{-1}$ is the matrix with entries $\omega_N^{-jk}$ – the same as $F_N$ but with $\omega$ replaced by $\omega^{-1} = \overline{\omega}$ (the complex conjugate). In other words, the inverse DFT is computed by the same algorithm as the forward DFT, just with the twiddle factors conjugated and the output scaled by $1/N$. Any algorithm that computes the forward DFT efficiently also computes the inverse efficiently, for free.

Proof of inversion. We need to verify that $\frac{1}{N}F_N^{-1} F_N = I_N$, i.e., that:

When $j = \ell$, every term in the sum is $\omega^0 = 1$, so the sum is $N$ and we get $N/N = 1$. When $j \neq \ell$, let $m = \ell - j \not\equiv 0 \pmod{N}$. Then $\sum_{k=0}^{N-1} \omega^{mk}$ is a geometric series with ratio $r = \omega^m \neq 1$:

since $\omega^N = 1$ by definition. $\blacksquare$

The Fast Fourier Transform (Cooley-Tukey, Radix-2)

The core idea: divide and conquer on even and odd indices

The breakthrough of Cooley and Tukey (1965) – though the idea traces back to Gauss (1805) – is to observe that when $N$ is even, the DFT of size $N$ can be decomposed into two DFTs of size $N/2$ plus $O(N)$ additional work.

Write $N = 2M$. Split the input sequence $\mathbf{x}$ into its even-indexed and odd-indexed elements:

Now consider the DFT sum for an arbitrary output index $k$:

Separate the even-indexed and odd-indexed terms:

Using the squaring property $\omega_N^2 = \omega_M$ (where $M = N/2$), each of these inner sums is itself a DFT of size $M$:

So the full DFT decomposes as:

But we need $X_k$ for $k = 0, 1, \ldots, N-1$, while $E_k$ and $D_k$ are periodic with period $M = N/2$ (they are DFTs of size $M$). For the “upper half” $k + M$, the cancellation property $\omega_N^{k+M} = -\omega_N^k$ gives:

These two equations together constitute the butterfly operation: from the pair $(E_k, D_k)$ and the twiddle factor $\omega_N^k$, we compute both $X_k$ and $X_{k+M}$ using one complex multiplication and two complex additions.

The butterfly diagram

Each butterfly takes two inputs, multiplies one by a twiddle factor, and produces two outputs via addition and subtraction:

  E_k ──────────┬──── (+) ──── X_k
                 │
          ω^k    ×
                 │
  D_k ──────────┴──── (−) ──── X_{k+M}

For $M = N/2$ values of $k$ (namely $k = 0, 1, \ldots, M-1$), we perform $M$ butterflies, each costing $O(1)$. The total additional work at this level of recursion is $O(N)$.

The recurrence and complexity

Let $T(N)$ denote the cost of computing a DFT of size $N$ (assumed a power of 2). The Cooley-Tukey decomposition gives:

By the Master Theorem (case 2, with $a = 2$, $b = 2$, $f(N) = \Theta(N)$):

compared to $O(N^2)$ for the naive DFT. For $N = 2^{20} \approx 10^6$, this is the difference between $\sim 10^{12}$ operations and $\sim 2 \times 10^7$ – a speedup of 50,000$\times$.

The full recursion tree

When $N = 2^s$, the recursion unfolds $s = \log_2 N$ levels deep. At each level, we partition the data, recurse on two halves, and combine with $N$ butterfly operations. The total work is:

Each level performs exactly $N/2$ butterflies (each costing one complex multiplication and two additions), for a total of $\frac{N}{2}\log_2 N$ complex multiplications.

A worked example: 8-point FFT

Let $N = 8$ and $\omega = \omega_8 = e^{2\pi i/8} = e^{i\pi/4} = \frac{1+i}{\sqrt{2}}$. Consider the input $\mathbf{x} = (1, 1, 1, 1, 0, 0, 0, 0)$ (a rectangular pulse).

Level 0 (split): Separate into even and odd:

Level 1: Each half recursively splits again. For $\mathbf{e}$:

These are 2-point DFTs: $\text{DFT}_2(a, b) = (a + b, \; a - b)$. So:

Level 1 butterfly (for $\mathbf{e}$): Combine with twiddle factors $\omega_4^0 = 1$ and $\omega_4^1 = i$:

So $\text{DFT}_4(\mathbf{e}) = (2, \; 1+i, \; 0, \; 1-i)$.

An identical calculation gives $\text{DFT}_4(\mathbf{d}) = (2, \; 1+i, \; 0, \; 1-i)$.

Level 0 butterfly: Combine with twiddle factors $\omega_8^k$ for $k = 0, 1, 2, 3$:

The key observation is not the specific numerical values but the structure: $3$ levels of $4$ butterflies each $= 12$ complex multiplications, versus $8^2 = 64$ for the naive DFT. The ratio $12/64 \approx 19\%$ – and this ratio improves as $N$ grows.

The visualization below runs this exact example through the 8-point butterfly network, stage by stage. The input goes in bit-reversed, the butterflies span 1, then 2, then 4, and every value shown is computed live. Hover or tap any value to see which two numbers and which twiddle factor produced it.

Putting It Together: FFT-Based Polynomial Multiplication

We now have all the pieces. To multiply two polynomials $A(x)$ and $B(x)$ of degree $n-1$ (and thereby two $n$-digit integers):

Step 1. Pad. The product polynomial has degree $2n - 2$, so we need at least $2n - 1$ evaluation points. Choose $N = 2^{\lceil \log_2(2n) \rceil}$ (the next power of 2 at or above $2n$). Pad both coefficient vectors with zeros to length $N$.

Step 2. Forward FFT. Compute $\hat{\mathbf{a}} = \text{FFT}_N(\mathbf{a})$ and $\hat{\mathbf{b}} = \text{FFT}_N(\mathbf{b})$. Cost: $O(N \log N)$ each.

Step 3. Pointwise multiply. Compute $\hat{c}_k = \hat{a}_k \cdot \hat{b}_k$ for $k = 0, \ldots, N-1$. Cost: $O(N)$.

Step 4. Inverse FFT. Compute $\mathbf{c} = \text{IFFT}_N(\hat{\mathbf{c}})$. Cost: $O(N \log N)$.

Step 5. Carry propagation. The entries of $\mathbf{c}$ are the exact convolution coefficients (no rounding – yet). Each $c_k$ may exceed the base $B$, so we perform a single left-to-right carry pass: set $c_k \leftarrow c_k \bmod B$ and add $\lfloor c_k / B \rfloor$ to $c_{k+1}$. Cost: $O(N)$.

Total cost: $O(N \log N) = O(n \log n)$.

At this point you might ask: if the FFT already gives us $O(n \log n)$ multiplication, why do we need Schönhage-Strassen? The answer lies in a subtle but critical issue: numerical precision.

Next: Vol. III: Schönhage-Strassen and Harvey-van der Hoeven

On Multiplication: Vol. I · Vol. II · Vol. III · Vol. IV · Vol. V