I want to explain how the Fast Fourier Transform (FFT) lets me multiply polynomials — and, more generally, transform between coefficient and value representations of a polynomial — far faster than the naive approach. Multiplying two polynomials of degree n the obvious way costs me $O(n^2)$ operations, but with FFT I can do it in $O(n \log n)$. I care about this algorithm because it sits at the intersection of algebra, signal processing, and algorithm design, and it shows up everywhere from audio processing to big-integer multiplication.
History and Background
I credit the modern, widely-used formulation of the FFT to James Cooley and John Tukey, who published their algorithm in 1965. However, I should note that similar ideas were discovered much earlier — Carl Friedrich Gauss described an equivalent method around 1805, long before the electronic computer existed, though his work went largely unnoticed for over a century. The Cooley-Tukey rediscovery in the context of digital computing transformed fields such as signal processing, telecommunications, and numerical analysis, and it is often cited as one of the most important algorithms of the 20th century.
Problem Statement
I frame the core problem as: given two polynomials A(x) and B(x) of degree less than n, each represented by their coefficients, compute their product C(x) = A(x) · B(x), a polynomial of degree less than 2n, as efficiently as possible. More generally, I want an efficient way to convert a polynomial between its coefficient representation and its value representation (its values at a specific set of points), since multiplication is easy in the value representation but hard in the coefficient representation.
Core Concepts
- Coefficient representation: a polynomial
A(x) = a_0 + a_1 x + ... + a_{n-1} x^{n-1}represented by the list of coefficients(a_0, ..., a_{n-1}). - Value representation: the same polynomial represented as a set of
npoint-value pairs(x_0, A(x_0)), ..., (x_{n-1}, A(x_{n-1})). - Discrete Fourier Transform (DFT): evaluating a polynomial at the
n-th roots of unity. - Roots of unity: complex numbers
ωsatisfyingω^n = 1; I use the principal root $\omega_n = e^{2\pi i / n}$. - Inverse DFT: the process of recovering coefficients from values at the roots of unity.
- Convolution: the coefficient sequence of the product polynomial is the convolution of the two input coefficient sequences.
How It Works
I break polynomial multiplication via FFT into three phases:
- Evaluation (DFT): I evaluate both input polynomials
AandBat then-th roots of unity, using the FFT algorithm, which does this in $O(n \log n)$ time instead of the naive $O(n^2)$. - Pointwise multiplication: I multiply the evaluated values together pointwise:
C(x_k) = A(x_k) · B(x_k)for each root of unityx_k. This step costs $O(n)$. - Interpolation (Inverse DFT): I convert the value representation of
Cback into coefficient form using the inverse FFT, again in $O(n \log n)$ time.
Before I start, I pad both polynomials with zero coefficients so their combined degree fits within a power of two n large enough to hold the product (n ≥ 2 · degree + 1), since the FFT’s divide-and-conquer structure works most naturally on power-of-two sizes.
Working Principle
I rely on the divide-and-conquer trick at the heart of FFT: I split a polynomial’s coefficients into even-indexed and odd-indexed subsets, forming two smaller polynomials of half the size. I recursively compute the DFT of each half, then combine the results using the special algebraic structure of roots of unity — in particular, the fact that $\omega_n^{k + n/2} = -\omega_n^k$, which lets me reuse each recursive result twice (once with a positive sign, once with a negative sign) instead of recomputing it. This halving-and-combining is what turns an $O(n^2)$ computation into an $O(n \log n)$ one.
Mathematical Foundation
I define the n-th principal root of unity as:
$$ \omega_n = e^{\frac{2\pi i}{n}} $$
The DFT of a coefficient vector $(a_0, \dots, a_{n-1})$ is the vector $(y_0, \dots, y_{n-1})$ where:
$$ y_k = \sum_{j=0}^{n-1} a_j , \omega_n^{jk}, \quad k = 0, 1, \dots, n-1 $$
The Cooley-Tukey recursive decomposition splits A(x) into even and odd terms:
$$ A(x) = A_{\text{even}}(x^2) + x \cdot A_{\text{odd}}(x^2) $$
which lets me write, for $k < n/2$:
$$ A(\omega_n^k) = A_{\text{even}}(\omega_{n/2}^k) + \omega_n^k \cdot A_{\text{odd}}(\omega_{n/2}^k) $$
$$ A(\omega_n^{k + n/2}) = A_{\text{even}}(\omega_{n/2}^k) – \omega_n^k \cdot A_{\text{odd}}(\omega_{n/2}^k) $$
The inverse DFT recovers the coefficients using the conjugate root of unity and a normalization factor:
$$ a_j = \frac{1}{n} \sum_{k=0}^{n-1} y_k , \omega_n^{-jk} $$
Diagrams
flowchart TD
A["Coefficient vectors A, B (size n)"] --> B["FFT: evaluate A at roots of unity"]
A --> C["FFT: evaluate B at roots of unity"]
B --> D["Pointwise multiply: C(x_k) = A(x_k) * B(x_k)"]
C --> D
D --> E["Inverse FFT: interpolate back to coefficients"]
E --> F["Coefficient vector of product polynomial C"]
Pseudocode
FFT(a) // a is a coefficient array of length n (power of 2)
n = length(a)
if n == 1
return a
w_n = e^(2*pi*i / n)
w = 1
a_even = FFT(a[0, 2, 4, ...])
a_odd = FFT(a[1, 3, 5, ...])
y = array of size n
for k = 0 to n/2 - 1
y[k] = a_even[k] + w * a_odd[k]
y[k + n/2] = a_even[k] - w * a_odd[k]
w = w * w_n
return y
MULTIPLY-POLYNOMIALS(A, B)
n = smallest power of 2 >= 2 * max(len(A), len(B))
pad A and B with zeros to length n
yA = FFT(A)
yB = FFT(B)
yC[k] = yA[k] * yB[k] for all k
C = INVERSE-FFT(yC) // divide every entry by n
return C
Step-by-Step Example
I multiply A(x) = 1 + 2x (coefficients [1, 2]) by B(x) = 3 + 4x (coefficients [3, 4]).
By hand, I know the answer should be:
$$ (1 + 2x)(3 + 4x) = 3 + 4x + 6x + 8x^2 = 3 + 10x + 8x^2 $$
so the coefficient vector should be [3, 10, 8].
Using FFT conceptually: I pad both to length n = 4 (the next power of 2 that fits a degree-2 result): A = [1, 2, 0, 0], B = [3, 4, 0, 0]. I evaluate both at the 4th roots of unity {1, i, -1, -i}, multiply the evaluations pointwise, and apply the inverse FFT. After interpolation, I recover the coefficient vector [3, 10, 8, 0], which matches my hand computation once I drop the trailing zero.
Time Complexity
- Best, average, and worst case: $O(n \log n)$ for the FFT and inverse FFT, since the divide-and-conquer recursion always splits evenly regardless of the input values — the running time depends only on
n, not on the specific coefficients. - Naive polynomial multiplication, by contrast, always costs $O(n^2)$, so FFT gives me an asymptotic improvement for large
n.
Space Complexity
I need $O(n)$ auxiliary space to hold the recursive arrays created during the FFT (for the even/odd splits) and $O(n)$ for the output arrays, so the overall space complexity is $O(n)$. An iterative, in-place FFT implementation (using bit-reversal permutation) can reduce constant factors but keeps the same asymptotic space bound.
Correctness Analysis
I justify correctness in two parts. First, the DFT/inverse-DFT pair is mathematically exact: I can prove using orthogonality of the roots of unity that applying the inverse DFT formula to the DFT values exactly recovers the original coefficients (subject to floating-point rounding in practice). Second, the recursive even/odd decomposition I use in the FFT is an exact algebraic identity — I am not approximating the DFT, I am just computing it in a smarter order — so the recursive FFT computes precisely the same values as the direct $O(n^2)$ DFT formula, just faster.
Advantages
- I get a dramatic asymptotic speedup, from $O(n^2)$ to $O(n \log n)$, for polynomial multiplication and convolution.
- It generalizes to many domains: signal processing, image processing, big-integer multiplication, and more.
- The same divide-and-conquer structure applies to related transforms (e.g., the Number Theoretic Transform for exact integer arithmetic).
Disadvantages
- Standard FFT implementations use floating-point complex numbers, introducing rounding error that I need to manage carefully (e.g., rounding results to the nearest integer for integer polynomial multiplication).
- The algorithm is conceptually more complex than naive multiplication, making it harder for me to implement correctly the first time.
- It requires padding to a power-of-two length, which can waste some computation for awkwardly-sized inputs (mitigated by more advanced mixed-radix FFTs).
Applications
- Digital signal processing: audio and image compression, filtering, and spectral analysis.
- Big-integer arithmetic: multiplying very large numbers (used in libraries like GMP) using FFT-based convolution.
- Polynomial and power series computations in computer algebra systems.
- Cryptography: some lattice-based cryptographic schemes use Number Theoretic Transforms, a variant of FFT over finite fields.
- String algorithms: certain pattern-matching problems (e.g., wildcard matching) can be solved efficiently using FFT-based convolution.
Implementation in C
#include <stdio.h>
#include <math.h>
#include <complex.h>
/* I implement a recursive, in-place-style FFT over complex numbers.
n must be a power of 2. sign = -1 for forward FFT, +1 for inverse FFT. */
void fft(double complex *a, int n, int sign) {
if (n == 1) {
return;
}
/* I split into even and odd indexed coefficients */
double complex even[n / 2], odd[n / 2];
for (int i = 0; i < n / 2; i++) {
even[i] = a[2 * i];
odd[i] = a[2 * i + 1];
}
/* I recursively transform each half */
fft(even, n / 2, sign);
fft(odd, n / 2, sign);
/* I combine the two halves using roots of unity */
for (int k = 0; k < n / 2; k++) {
double complex w = cexp(sign * 2.0 * I * M_PI * k / n);
double complex t = w * odd[k];
a[k] = even[k] + t;
a[k + n / 2] = even[k] - t;
}
}
/* I multiply two real-coefficient polynomials using FFT-based convolution */
void multiplyPolynomials(double *A, int lenA, double *B, int lenB, double *result) {
int resultSize = lenA + lenB - 1;
int n = 1;
while (n < resultSize) n <<= 1; /* I round up to the next power of 2 */
double complex fa[n], fb[n];
for (int i = 0; i < n; i++) {
fa[i] = (i < lenA) ? A[i] : 0.0;
fb[i] = (i < lenB) ? B[i] : 0.0;
}
fft(fa, n, -1);
fft(fb, n, -1);
for (int i = 0; i < n; i++) {
fa[i] = fa[i] * fb[i];
}
fft(fa, n, 1); /* inverse FFT */
for (int i = 0; i < resultSize; i++) {
result[i] = creal(fa[i]) / n; /* I normalize and take the real part */
}
}
int main(void) {
double A[] = {1, 2}; /* 1 + 2x */
double B[] = {3, 4}; /* 3 + 4x */
int resultSize = 2 + 2 - 1;
double result[resultSize];
multiplyPolynomials(A, 2, B, 2, result);
printf("Product coefficients: ");
for (int i = 0; i < resultSize; i++) {
printf("%.0f ", result[i]);
}
printf("\n");
return 0;
}
I wrote fft as a straightforward recursive implementation for clarity rather than maximum performance; a production-grade version would use an iterative, in-place implementation with bit-reversal permutation to avoid the overhead of recursive array splitting.
Sample Input and Output
Input: A(x) = 1 + 2x, B(x) = 3 + 4x
Output:
Product coefficients: 3 10 8
which matches 3 + 10x + 8x^2 exactly, as I computed by hand earlier.
Optimization Techniques
- Iterative FFT with bit-reversal: I can avoid recursion overhead by rearranging the input array according to bit-reversed indices and then processing it iteratively, level by level.
- Precomputed twiddle factors: I compute the roots of unity once and reuse them, rather than recomputing
cexpcalls repeatedly. - Number Theoretic Transform (NTT): for integer polynomial multiplication where I need exact results, I replace complex roots of unity with roots of unity in a finite field, avoiding floating-point rounding entirely.
- In-place computation: I can perform the FFT without allocating new arrays at each recursive level, reducing memory allocation overhead.
- SIMD/vectorization: for large-scale applications, I can vectorize the butterfly operations to exploit modern CPU instruction sets.
Common Mistakes
- Forgetting to pad polynomials to a power-of-two length before running FFT.
- Forgetting to divide by
nafter the inverse FFT, which leaves the result scaled incorrectly. - Using the wrong sign convention for forward vs. inverse FFT (mixing up $e^{+i\theta}$ and $e^{-i\theta}$).
- Not rounding the final real-valued output when I know the true answer should be an integer, letting floating-point noise show up in the result.
- Allocating variable-length arrays recursively without considering stack depth for very large
n, which can cause stack overflow in naive recursive implementations.
Further Reading
- Cormen, Leiserson, Rivest, Stein — Introduction to Algorithms, Chapter 30 (Polynomials and the FFT): https://mitpress.mit.edu/9780262046305/introduction-to-algorithms/
- Cooley, Tukey — “An Algorithm for the Machine Calculation of Complex Fourier Series,” Mathematics of Computation, 1965: https://www.ams.org/journals/mcom/1965-19-090/S0025-5718-1965-0178586-1/
- Press, Teukolsky, Vetterling, Flannery — Numerical Recipes: The Art of Scientific Computing: https://numerical.recipes/
- MIT OpenCourseWare — 6.046J Design and Analysis of Algorithms, FFT lecture notes: https://ocw.mit.edu/courses/6-046j-design-and-analysis-of-algorithms-spring-2015/
