Coefficient Representation vs. Point-Value Representation
A polynomial $A(x) = \sum_{j=0}^{n-1} a_j x^j$ of degree $< n$ can be represented in two dual ways:
- Coefficient Form: The vector $(a_0, a_1, \dots, a_{n-1})$. Adding two polynomials takes $O(n)$, but multiplying $C(x) = A(x)B(x)$ requires discrete convolution $c_k = \sum_{j=0}^{k} a_j b_{k-j}$ in $\Theta(n^2)$ time.
- Point-Value Form: A set of $n$ distinct sample points $\{(x_0, y_0), \dots, (x_{n-1}, y_{n-1})\}$ where $y_k = A(x_k)$. By the Fundamental Theorem of Algebra, $n$ point-value pairs uniquely determine a polynomial of degree $< n$. In point-value form, multiplying $C(x) = A(x)B(x)$ (padded to $N \ge 2n - 1$ points) is simply pointwise multiplication $C(x_k) = A(x_k) \cdot B(x_k)$ in $\Theta(N)$ time!
The Collapsing Property of Complex $N$-th Roots of Unity
How can we evaluate $A(x)$ at $N$ points in $O(N \log N)$ time instead of $O(N^2)$? Group the terms of $A(x)$ into even-indexed and odd-indexed coefficients:
$$A(x) = \underbrace{(a_0 + a_2 x^2 + a_4 x^4 + \dots)}_{A_{\text{even}}(x^2)} + x \underbrace{(a_1 + a_3 x^2 + a_5 x^4 + \dots)}_{A_{\text{odd}}(x^2)}$$Notice that both $A_{\text{even}}(u)$ and $A_{\text{odd}}(u)$ are polynomials of half the size ($N/2$) evaluated at $u = x^2$. If we choose evaluation points that come in opposite pairs $\pm x$, then $(+x)^2 = (-x)^2 = x^2$, so evaluating $A$ at $N$ points reduces to evaluating $A_{\text{even}}$ and $A_{\text{odd}}$ at only $N/2$ squared points! For this +/- symmetry to hold recursively across all $\log_2 N$ levels (where squares of squares must also be +/- pairs), we must choose the complex $N$-th roots of unity:
$$\omega_N^k = e^{2\pi i k / N} = \cos\!\left(\frac{2\pi k}{N}\right) + i \sin\!\left(\frac{2\pi k}{N}\right), \qquad k = 0, 1, \dots, N-1$$Inspect the Radix-2 Cooley-Tukey FFT butterfly y_k = E_k + ω_n^k O_k, Karatsuba 3-multiplication algebraic identity, and CDQ Divide-and-Conquer left-to-right contribution transfer.
Invariant & State
| Transform / Convolution Method | Time Complexity | Arithmetic Domain | Numerical Precision / Constraints |
|---|---|---|---|
| Direct Convolution / Karatsuba | $\Theta(N^2)$ / $\Theta(N^{1.585})$ | Exact Integers / Reals | Exact; slow for $N > 5000$ |
| Cooley-Tukey Radix-2 FFT | $\Theta(N \log N)$ | Complex double ($mathbb{C}$) |
~53-bit mantissa; rounding error for coefficients $> 10^{14}$ |
| Number-Theoretic Transform (NTT) | $\Theta(N \log N)$ | $\mathbb{Z}_p$ ($p = 998244353 = 119 \cdot 2^{23} + 1$) | 100% Exact Modular Arithmetic (Primitive Root $g = 3$) |
For any even $N$ and $k \in [0, N/2)$:
- Halving Lemma: $\left(\omega_N^k\right)^2 = \left(\omega_N^{k + N/2}\right)^2 = e^{2\pi i k / (N/2)} = \omega_{N/2}^k$. Squaring all $N$ roots of unity collapses them into the $N/2$-th roots of unity twice over.
- Reflection Lemma: $\omega_N^{k + N/2} = -\omega_N^k$.
Combining these gives the Radix-2 Butterfly Operation that computes both $y_k$ and $y_{k + N/2}$ from $u = y^{[\text{even}]}_k$ and $v = \omega_N^k \cdot y^{[\text{odd}]}_k$ in $O(1)$:
$$y_k = u + v, \qquad y_{k + N/2} = u - v$$And because the Vandermonde matrix $V_N[j, k] = \omega_N^{jk}$ satisfies $V_N^{-1}[j, k] = \frac{1}{N} \omega_N^{-jk}$, the Inverse FFT (IFFT) uses the exact same algorithm with $\omega_N^{-1} = e^{-2\pi i / N}$ followed by dividing each output by $N$!
A recursive FFT allocates temporary arrays at every level to separate even indices ($0, 2, 4, 6$) and odd indices ($1, 3, 5, 7$). Observe what happens to the binary representation of index $i$ as we repeatedly partition by parity (least significant bit first, then second bit, etc.): an index $i$ ends up at the position given by reversing its $\log_2 N$ binary bits! By permuting the array in-place via a single $O(N)$ bit-reversal pass upfront, we can execute all $\log_2 N$ butterfly stages purely in-place and bottom-up ($len = 2, 4, 8, \dots, N$) with zero heap allocations and sequential SIMD memory strides.
Production C++20 Implementation: Iterative In-Place FFT & Polynomial Convolution
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <numbers>
#include <vector>
using cd = std::complex<double>;
// In-place Iterative Radix-2 Cooley-Tukey FFT (invert = false) and IFFT (invert = true)
void fft_inplace(std::vector<cd>& a, bool invert) noexcept {
const std::size_t n = a.size();
// 1. O(N) In-place Bit-Reversal Permutation
for (std::size_t i = 1, j = 0; i < n; ++i) {
std::size_t bit = n >> 1;
for (; j & bit; bit >>= 1) {
j ^= bit;
}
j ^= bit;
if (i < j) std::swap(a[i], a[j]);
}
// 2. Bottom-Up Butterfly Stages: len = 2, 4, 8, ..., n
for (std::size_t len = 2; len <= n; len <<= 1) {
const double ang = 2.0 * std::numbers::pi / static_cast<double>(len) * (invert ? -1.0 : 1.0);
const cd wlen(std::cos(ang), std::sin(ang));
const std::size_t half = len >> 1;
for (std::size_t i = 0; i < n; i += len) {
cd w(1.0, 0.0);
for (std::size_t j = 0; j < half; ++j) {
const cd u = a[i + j];
const cd v = a[i + j + half] * w;
a[i + j] = u + v;
a[i + j + half] = u - v;
w *= wlen;
}
}
}
// 3. Normalize by 1/N for Inverse FFT
if (invert) {
const double inv_n = 1.0 / static_cast<double>(n);
for (cd& x : a) x *= inv_n;
}
}
// Convolves two integer polynomials A and B in O(N log N) time
std::vector<std::int64_t> multiply_polynomials_fft(const std::vector<std::int64_t>& A,
const std::vector<std::int64_t>& B) {
if (A.empty() || B.empty()) return {};
const std::size_t result_len = A.size() + B.size() - 1;
std::size_t n = 1;
while (n < result_len) n <<= 1;
std::vector<cd> fa(n, cd(0.0, 0.0)), fb(n, cd(0.0, 0.0));
for (std::size_t i = 0; i < A.size(); ++i) fa[i] = static_cast<double>(A[i]);
for (std::size_t i = 0; i < B.size(); ++i) fb[i] = static_cast<double>(B[i]);
fft_inplace(fa, false);
fft_inplace(fb, false);
for (std::size_t i = 0; i < n; ++i) {
fa[i] *= fb[i];
}
fft_inplace(fa, true);
std::vector<std::int64_t> result(result_len);
for (std::size_t i = 0; i < result_len; ++i) {
result[i] = static_cast<std::int64_t>(std::llround(fa[i].real()));
}
return result;
}Two bugs frequently corrupt FFT convolutions:
- Under-padding (Circular Convolution): Multiplying degree-$d_1$ and degree-$d_2$ polynomials produces degree $d_1 + d_2$ ($d_1 + d_2 + 1$ coefficients). If you only pad $N$ to $\max(d_1, d_2)$ instead of the next power of two $\ge d_1 + d_2 + 1$, terms with $k \ge N$ alias mod $N$ ($x^N \equiv 1$) and corrupt the low coefficients!
- Twiddle Factor Drift & Integer Overflow: Repeatedly multiplying
w *= wlenover $N = 2^{20}$ iterations accumulates IEEE-754 rounding error. For exact integer convolutions where coefficients exceed $2^{53} \approx 9 \times 10^{15}$, use the Number-Theoretic Transform (NTT) modulo a prime $p = c \cdot 2^k + 1$ (such as $998244353 = 119 \cdot 2^{23} + 1$, where $3^{(p-1)/N} \bmod p$ serves as an exact integer $N$-th root of unity).