The Fourier Transform

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 with $2k-1$ rows, coefficient blowup at large evaluation points) grows so fast that implementations stop at small $k$: GMP’s largest variant is Toom-8½, and above that size it switches to an FFT.

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 family of points whose algebraic symmetry enables a much 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$. 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 the integer $\sum_k c_k B^k$, where the coefficient 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 of length $N$ (defined below). If $\mathbf{c} = \mathbf{a} \ast \mathbf{b}$ is the cyclic convolution of two length-$N$ sequences, $c_k = \sum_{j=0}^{N-1} a_j b_{(k-j) \bmod N}$, 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.” Cyclic convolution wraps indices around modulo $N$, but if $\mathbf{a}$ and $\mathbf{b}$ are zero-padded to length $N \geq 2n - 1$, nothing wraps and it equals the linear convolution above. 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)$, a vast improvement over $O(N^2)$.

This is not a hand-wave. We now build the DFT, prove the theorem, and show why the DFT 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 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 standard primitive $N$-th root of unity (one whose powers give all $N$ roots) is:

so that $\omega_N^k = (\omega_N)^k$. We will usually drop the subscript $N$ when the context is clear. The positive exponent follows CLRS; signal-processing texts usually write $e^{-2\pi i/N}$, and nothing below depends on the choice.

Four algebraic properties of the roots of unity drive 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, a half turn around the circle takes $\omega^k$ to the diametrically opposite point $-\omega^k$. The FFT relies on this property more than on any other.

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. This proves the Convolution Theorem: cyclic convolution computes the coefficients of $A(z)B(z)$ with $z^N$ replaced by $1$, and since $(\omega^k)^N = 1$, that replacement does not change the value at $\omega^k$. So the $k$-th entry of $\mathcal{F}(\mathbf{a} \ast \mathbf{b})$ is $A(\omega^k)B(\omega^k)$, the product of the $k$-th entries of $\mathcal{F}(\mathbf{a})$ and $\mathcal{F}(\mathbf{b})$.

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} \overline{F_N} \mathbf{X}$ where $\overline{F_N}$ 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). So $F_N^{-1} = \frac{1}{N} \overline{F_N}$. In other words, the inverse DFT is computed by the same algorithm as the forward DFT, just with every power of $\omega$ conjugated and the output scaled by $1/N$. Any fast algorithm for the forward DFT is therefore also a fast algorithm for the inverse.

Proof of inversion. We need to verify that $\frac{1}{N}\overline{F_N} F_N = I_N$, i.e., that for all $j$ and $\ell$:

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$, so $E_{k+M} = E_k$ and $D_{k+M} = D_k$). 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 about 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 the halves with $N/2$ butterflies, which produce $N$ outputs. Counting one unit of work per output, 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 is $12/64 \approx 19\%$, and it 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. In exact arithmetic, the entries of $\mathbf{c}$ are the convolution coefficients themselves. Each $c_k$ may exceed the base $B$, so we perform a single carry pass from $c_0$ upward: add $\lfloor c_k / B \rfloor$ to $c_{k+1}$, then set $c_k \leftarrow c_k \bmod B$. Cost: $O(N)$.

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

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 is numerical precision.