Fast Fourier Transform¶

From Computability, Complexity & Algorithms — Charles Brubaker and Lance Fortnow

Introduction: Convolution and Polynomial Multiplication¶

The Fast Fourier Transform (FFT) enables efficient convolution of sequences, reducing the naive \(O(nm)\) cost to \(O(N \log N)\) where \(N = n + m - 1\) is the output length.

Convolution arises in signal processing, image filtering, and polynomial multiplication. The \(k\)-th coefficient of the convolution of sequences \(a\) and \(b\) is:

\[c_k = \sum_{j} a_j \cdot b_{k-j}\]

ConvolutionStart — Initial alignment of the two sequences for convolution at \(k = 0\).¶

ConvolutionEnd — Final state after sliding the reversed sequence through the convolution.¶

QuadTerm — Alignment demonstrating calculation of the \(x^2\) term in polynomial multiplication.¶

Convolution of \(a\) and \(b\) is equivalent to multiplying the corresponding polynomials \(A(x)\) and \(B(x)\).

Polynomial Representations¶

A degree-\(n\) polynomial \(A(x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_{n-1} x^{n-1}\) can be represented in two ways:

Coefficient form: the vector \((a_0, a_1, \ldots, a_{n-1})\).
Value form: the values \(A(x_0), A(x_1), \ldots, A(x_{N-1})\) at \(N\) distinct points (requires \(N \geq n\)).

Conversion from coefficients to values at \(N\) points uses the Vandermonde matrix:

\[\begin{split}\begin{pmatrix} 1 & x_0 & x_0^2 & \cdots \\ 1 & x_1 & x_1^2 & \cdots \\ \vdots & & & \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ \vdots \end{pmatrix} = \begin{pmatrix} A(x_0) \\ A(x_1) \\ \vdots \end{pmatrix}\end{split}\]

VandermondeDef — Vandermonde matrix converting polynomial coefficients to point values.¶

Key observation: In value form, pointwise multiplication of two degree-\(N\) polynomials costs only \(O(N)\):

CoefficientVsValue — Multiplication is \(O(N^2)\) in coefficient form but \(O(N)\) in value form.¶

The Three-Step FFT Process¶

TheProcess — **Evaluate → Multiply → Interpolate**: convert to values, multiply pointwise, convert back to coefficients.¶

Evaluate \(A\) and \(B\) at \(N\) chosen points → value representations.
Multiply pointwise: \(C(x_i) = A(x_i) \cdot B(x_i)\) for each \(i\) — cost \(O(N)\).
Interpolate the product values back to coefficients of \(C\).

MultQuiz — *Quiz:* Multiply two polynomials using value representation.¶

TheProcessBigO — Time complexity of each stage — the bottleneck is evaluation and interpolation.¶

The naive evaluation at \(N\) arbitrary points costs \(O(N^2)\) (one Horner evaluation per point). The FFT achieves \(O(N \log N)\) by choosing the evaluation points cleverly.

Even-Odd Decomposition¶

Split \(A(x)\) into even- and odd-indexed coefficients:

\[A_e(x) = a_0 + a_2 x + a_4 x^2 + \cdots \qquad A_o(x) = a_1 + a_3 x + a_5 x^2 + \cdots\]

Then:

\[\begin{split}A(x) &= A_e(x^2) + x \cdot A_o(x^2) \\ A(-x) &= A_e(x^2) - x \cdot A_o(x^2)\end{split}\]

EvenOddTerms — Evaluating at \(+x\) and \(-x\) shares the computation of \(A_e(x^2)\) and \(A_o(x^2)\).¶

Saving: One recursive call at \(x^2\) gives both \(A(x)\) and \(A(-x)\) — halving the work at each level.

Roots of Unity¶

To apply the even-odd trick recursively, choose evaluation points as the \(N\)-th roots of unity:

\[\omega_N = e^{2\pi i / N}, \qquad \omega_N^0, \omega_N^1, \ldots, \omega_N^{N-1}\]

These satisfy the key cancellation property:

\[\omega_N^{j + N/2} = -\omega_N^j\]

so the positive-negative pairing needed by the even-odd split is automatically present.

RootsOfUnity8 — The 8th roots of unity arranged on the unit circle in the complex plane.¶

RootsOfUnity4 — Squaring the 8th roots of unity yields the 4th roots of unity.¶

Recursive FFT Structure¶

FFTFirstLayer — First decomposition level: evaluate at 8th roots → two sub-problems at 4th roots.¶

FFTSecondLayer — Second decomposition level: 4th roots → two sub-problems at 2nd roots (base case).¶

ApplyFFTQuiz — *Quiz:* Apply FFT to evaluate a polynomial at the four 4th roots of unity.¶

FFT Algorithm¶

FFT(a, ω):          # a = coefficient vector, ω = primitive N-th root
    if N == 1: return a[0]
    a_e = [a[0], a[2], ..., a[N-2]]
    a_o = [a[1], a[3], ..., a[N-1]]
    y_e = FFT(a_e, ω²)
    y_o = FFT(a_o, ω²)
    for j = 0 to N/2 - 1:
        y[j]       = y_e[j] + ω^j · y_o[j]
        y[j + N/2] = y_e[j] - ω^j · y_o[j]
    return y

Recurrence:

\[T(N) = 2T(N/2) + \Theta(N) \;\Longrightarrow\; T(N) = \Theta(N \log N)\]

by the Master Theorem.

Butterfly Network¶

The FFT computation can be visualised as a butterfly network — a circuit of \(\log N\) layers each performing \(N/2\) butterfly operations.

ButterflyRightLeft — Right-to-left paths encode the **bit-reversal** permutation of input indices.¶

ButterflyLeftRight — Left-to-right paths encode the **twiddle-factor** \(\omega^j\) multiplications.¶

UpdatedProcess — Updated running time: evaluation and interpolation both cost \(O(N \log N)\), giving overall \(O(N \log N)\) polynomial multiplication.¶

Inverse FFT (Interpolation)¶

Given the values \(A(\omega^0), \ldots, A(\omega^{N-1})\), recover the coefficients using the inverse DFT. The Vandermonde matrix at roots of unity has a convenient inverse:

\[a_j = \frac{1}{N} \sum_{k=0}^{N-1} A(\omega^k) \cdot \omega^{-jk}\]

This is simply FFT evaluated at the conjugate roots \(\bar{\omega} = \omega^{-1}\), followed by scaling by \(1/N\).

VandermondeAtRootsOfUnity — Vandermonde matrix structure at the \(N\)-th roots of unity and its inverse.¶

InverseFFTAlg — Inverse FFT algorithm: run FFT with \(\omega^{-1}\) then divide by \(N\).¶

ApplyInverseFFTQuiz — *Quiz:* Recover polynomial coefficients from their values at roots of unity.¶

Complete Process¶

FinalProcess — Full polynomial multiplication pipeline achieving \(O(N \log N)\) overall.¶

Summary of costs:

Stage	Naive cost	FFT cost
Evaluate \(A\), \(B\) at \(N\) points	\(O(N^2)\)	\(O(N \log N)\)
Pointwise multiply	\(O(N)\)	\(O(N)\)
Interpolate \(C\) from values	\(O(N^2)\)	\(O(N \log N)\)
Total	\(O(N^2)\)	\(O(N \log N)\)

Key Formulas¶

Even-odd split:

\[A(x) = A_e(x^2) + x \cdot A_o(x^2), \qquad A(-x) = A_e(x^2) - x \cdot A_o(x^2)\]

Primitive \(N\)-th root of unity:

\[\omega_N = e^{2\pi i / N}\]

Cancellation property: