Cooley-Tukey FFT

Purpose: The recursive divide-and-conquer scheme that computes a Discrete Fourier Transform in $O(n\log n)$. When people say “the FFT”, they almost always mean radix-2 Cooley-Tukey.

The Decimation

Split the input by index parity:

$$X_k=\sum _{j=0}^{n-1}{x}_j{\omega }^{jk}=\underset{\text{even}}{\underbrace{\sum _j{x}_{2j}{\omega }^{2jk}}}+{\omega }^k\underset{\text{odd}}{\underbrace{\sum _j{x}_{2j+1}{\omega }^{2jk}}}$$

Since ${\omega }^2$ is a primitive $(n/2)$-th root of unity, both sums are DFTs of half the size. Writing $E_k$ and $O_k$ for them:

$X_k=E_k+{\omega }^k{O}_k,X_{k+n/2}=E_k-{\omega }^k{O}_k.$

That pair is the butterfly, and it is why one DFT of size $n$ costs two of size $n/2$ plus $O(n)$:

$T(n)=2T(n/2)+O(n)=O(n\log n).$

Iterative implementation

The recursion permutes the input into bit-reversed order. Doing that permutation up front turns the algorithm into three tidy loops with no recursion:

using cd = complex<double>;
 
void fft(vector<cd>& a, bool invert) {
    int n = a.size();
    for (int i = 1, j = 0; i < n; i++) {          // bit-reversal permutation
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit;
        if (i < j) swap(a[i], a[j]);
    }
    for (int len = 2; len <= n; len <<= 1) {      // butterfly levels
        double ang = 2 * M_PI / len * (invert ? -1 : 1);
        cd wlen(cos(ang), sin(ang));
        for (int i = 0; i < n; i += len) {
            cd w(1);
            for (int j = 0; j < len / 2; j++) {
                cd u = a[i + j], v = a[i + j + len / 2] * w;
                a[i + j] = u + v;
                a[i + j + len / 2] = u - v;
                w *= wlen;
            }
        }
    }
    if (invert) for (cd& x : a) x /= n;
}

Complexity

• Time: $O(n\log n)$, with $\frac{n}{2}\log n$ butterflies
• Space: $O(n)$, fully in place after the bit-reversal

Radix and variants

Variant	Idea
Radix-2 DIT	split by index parity (the code above)
Radix-2 DIF	split by output parity; no bit-reversal on input, but on output
Radix-4 / split-radix	fewer multiplications; ~20% faster
Good-Thomas	$n=n_1{n}_2$ with $\gcd (n_1,n_2)=1$; no twiddle factors at all
Rader	prime $n$, via a convolution over the multiplicative group
Bluestein	arbitrary $n$, via a chirp-z convolution

Precision

Doubles give roughly ${\log }_{2}n+{\log }_2(\max |a_i|\cdot \max |b_i|)$ bits of error growth. For products of integers up to $10^9$ with $n=10^6$, plain complex FFT will lose the low bits. Use NTT over a prime modulus, or split each coefficient into two 15-bit halves and do four transforms.

Variants / Use Cases

• FFT and NTT — the topic page with the full contest toolkit
• Fast Fourier Transform — the overview page
• NTT — exact integer version, no precision loss
• Convolution — the reason you want an FFT at all
• Karatsuba / Toom-Cook — the sub-quadratic methods used below the FFT crossover