Multiplying polynomials in Θ(n lg n) by changing what a polynomial is — and one of the most consequential algorithms ever written.
Multiplying two degree-n polynomials by the schoolbook method costs Θ(n²). Chapter 30 gets it down to Θ(n lg n), and the mechanism is a change of representation rather than a cleverer multiplication. In coefficient form multiplication is a convolution and is slow; in point-value form it is pointwise and is trivially linear. The fast Fourier transform converts between the two representations in Θ(n lg n), and the whole algorithm is: convert, multiply pointwise, convert back.
| Coefficient form | Point-value form | |
|---|---|---|
| Stores | a = (a₀, a₁, …, an-1) | {(x₀,y₀), …, (xn-1, yn-1)} with yₖ = A(xₖ) |
| Evaluate at a point | Θ(n) by Horner’s rule | Θ(n²) — must interpolate first |
| Add | Θ(n) | Θ(n) |
| Multiply | Θ(n²) — convolution | Θ(n) — pointwise! |
In coefficient form the product coefficients are a convolution:
cⱼ = ∑ from k=0 to j of aₖ bj-k // every pair contributes: Θ(n²)A(xₖ) = yₖ and B(xₖ) = zₖ, then the product satisfies C(xₖ) = yₖ zₖ. One multiplication per point, Θ(n) in total. The whole chapter exists to make the conversion between representations cheap enough to exploit that.n polynomials has degree bound 2n - 1, so it needs 2n point-value pairs to be determined. Both inputs must therefore be evaluated at 2n points, padded with zero coefficients up to length 2n. Forgetting this is the standard implementation bug.coefficient form ──── evaluate at 2n points ────▶ point-value form
a, b Θ(n lg n) via FFT y, z
│
│ pointwise
│ multiply, Θ(n)
▼
coefficient form ◀─── interpolate ─────────────── point-value form
c Θ(n lg n) via inverse FFT y·zThree steps, Θ(n lg n) + Θ(n) + Θ(n lg n) = Θ(n lg n). Everything now depends on doing evaluation and interpolation fast, and in general both are Θ(n²).
0, 1, 2, …. Choose the n complex nth roots of unity and the evaluation acquires a recursive structure that divide-and-conquer can exploit.The n complex nth roots of unity are the solutions of zⁿ = 1, evenly spaced around the unit circle:
ωₙk = e2πik/n for k = 0, 1, …, n-1
where ωₙ = e2πi/n is the principal nth root of unityn nth roots produces the n/2 (n/2)th roots, each hit exactly twice.n > 0 is even, then the squares of the n complex nth roots of unity are the n/2 complex (n/2)th roots of unity. This is what makes the recursion work — the subproblem is the same problem at half the size, with the same kind of evaluation points.Two supporting facts: the cancellation lemma ωdkdn = ωkn, and the summation lemma that ∑ₖ (ωₙk)ʲ = 0 for any j not divisible by n — which is what makes the inverse transform work.
n roots of unity is the discrete Fourier transform:yₖ = A(ωₙk) = ∑ from j=0 to n-1 of aⱼ ωₙkj
y = DFTₙ(a)Computed directly this is n evaluations at Θ(n) each, so Θ(n²). The FFT computes the same thing in Θ(n lg n).
Split the coefficients by parity into two polynomials of half the degree:
A[0](x) = a₀ + a₂x + a₄x² + … // even-indexed coefficients
A[1](x) = a₁ + a₃x + a₅x² + … // odd-indexed coefficients
A(x) = A[0](x²) + x · A[1](x²)A at the n nth roots reduces to evaluating two half-size polynomials at the squares of those roots — which by the halving lemma are just the n/2 (n/2)th roots. So the subproblem is the same problem at half the size.RECURSIVE-FFT(a)
1 n = a.length // n is a power of 2
2 if n == 1
3 return a
4 ωₙ = e2πi/n; ω = 1
5 a[0] = (a₀, a₂, …, an-2)
6 a[1] = (a₁, a₃, …, an-1)
7 y[0] = RECURSIVE-FFT(a[0])
8 y[1] = RECURSIVE-FFT(a[1])
9 for k = 0 to n/2 - 1
10 yₖ = yₖ[0] + ω · yₖ[1] // the butterfly
11 yk+n/2 = yₖ[0] - ω · yₖ[1]
12 ω = ω · ωₙ
13 return yT(n) = 2T(n/2) + Θ(n) = Θ(n lg n) — the merge-sort recurrence, solved by master case 2.
ω · yₖ[1], using the fact that ωₙk+n/2 = -ωₙk. Half the work is free by symmetry. Iterative implementations arrange the input in bit-reversed order and then run lg n levels of butterflies in place, with no recursion and no extra array.Interpolation — going from point-value back to coefficients — looks like a different problem, and is not.
The DFT is multiplication by the Vandermonde matrix Vₙ with entries vⱼₖ = ωₙkj. Its inverse has a startlingly simple closed form:
(Vₙ⁻¹)ⱼₖ = ωₙ-kj / n
aⱼ = (1/n) · ∑ from k=0 to n-1 of yₖ ωₙ-kjωₙ for ωₙ⁻¹, divide the result by n, and RECURSIVE-FFT computes the inverse DFT. One piece of code serves both directions, which is why FFT libraries expose a single routine with a direction flag.| Application | How the FFT is used |
|---|---|
| Signal processing | Convert a waveform between time and frequency domains. This is the original motivation and still the dominant use. |
| Audio and image compression | MP3, JPEG, and their descendants transform to the frequency domain and discard coefficients the eye or ear cannot detect. |
| Large-integer multiplication | Digits are polynomial coefficients, so multiplying big integers is polynomial multiplication. The Schönhage-Strassen and Harvey-van der Hoeven algorithms build on this. |
| Convolution generally | Any convolution becomes pointwise multiplication in the transformed domain — used for filtering, correlation, and pattern matching. |
| Polynomial arithmetic | Multiplication, division, and evaluation in computer algebra systems. |
n is a power of 2; real libraries handle arbitrary n by mixed-radix decomposition or by padding.Θ(n²); point-value form makes it Θ(n), pointwise. The chapter is about converting cheaply between them.2n points, multiply pointwise, interpolate back. The product needs 2n points, so pad both inputs.nth roots of unity.n nth roots are the n/2 (n/2)th roots. That self-similarity is what makes divide-and-conquer apply.ωk+n/2 = -ωᵏ. T(n) = 2T(n/2) + Θ(n) = Θ(n lg n).ωₙ⁻¹ in place of ωₙ and a division by n.Chapter 31 covers number-theoretic algorithms: greatest common divisors, modular arithmetic, primality testing, and the RSA cryptosystem. It is the chapter where algorithms whose inputs are measured in bits rather than items come into their own, and where the hardness of factoring becomes a feature rather than a problem.