Fast Fourier Transform, in depth: roots of unity, the butterfly, iterative in-place code, fast convolution, precision limits and where FFTs run in ML

Two ways to describe a polynomial

A polynomial A(x) = a0 + a1 x + ... + a(n-1) x^(n-1) can be stored as its n coefficients. It can equally be stored as its values at n distinct points: a polynomial of degree below n is uniquely determined by n samples. The two forms have opposite strengths. Multiplying is O(n squared) on coefficients but O(n) in point-value form: if you know A(x_k) and B(x_k) at the same points, then C(x_k) = A(x_k) * B(x_k).

So fast multiplication is a three-step plan: evaluate both polynomials at enough points, multiply pointwise, and interpolate back to coefficients. Naive evaluation at n points is O(n squared). The FFT makes evaluation and interpolation O(n log n) by choosing very special points: the complex n-th roots of unity.

The DFT is evaluation at roots of unity

Let w = exp(-2*pi*i/n), a primitive n-th root of unity. Its powers w^0, w^1, ..., w^(n-1) are n equally spaced points on the unit circle in the complex plane. The Discrete Fourier Transform of the coefficient vector a is the vector y with y[k] = sum over j of a[j] * w^(j*k). That is exactly A evaluated at w^k. This article uses NumPy's sign convention.

Two properties of these points make the FFT possible. The halving lemma: squaring the n-th roots of unity gives the (n/2)-th roots of unity, each twice. The cancellation property: w^(k + n/2) = -w^k. Halving keeps subproblems on the same kind of point set. Cancellation lets you reuse one multiplication for two outputs.

Divide and conquer: the even/odd split

Split A into its even-indexed and odd-indexed coefficients: E(x) = a0 + a2 x + a4 x^2 + ... and O(x) = a1 + a3 x + a5 x^2 + .... Then A(x) = E(x^2) + x * O(x^2). To evaluate A at all n roots, you need E and O at the squares of those roots, which by the halving lemma are just the n/2 roots of unity. That is two FFTs of half the size. Combining uses cancellation:

A(w^k)         = E(w^2k) + w^k * O(w^2k)
A(w^(k + n/2)) = E(w^2k) - w^k * O(w^2k)        for k = 0 .. n/2 - 1

Each pair of outputs costs one complex multiplication (by the twiddle factor w^k), one addition and one subtraction. This pattern is called a butterfly. The recurrence is T(n) = 2T(n/2) + O(n), which solves to O(n log n). The recursive version translates the math directly:

import cmath

def fft(a):
    """Recursive radix-2 DFT: y[k] = sum_j a[j] * exp(-2*pi*i*j*k/n). len(a) must be a power of 2."""
    n = len(a)
    if n == 1:
        return list(a)
    even = fft(a[0::2])                  # A_even evaluated at the (n/2)-th roots of unity
    odd = fft(a[1::2])                   # A_odd  evaluated at the same points
    y = [0j] * n
    for k in range(n // 2):
        w = cmath.exp(-2j * cmath.pi * k / n)   # twiddle factor w_n^k
        t = w * odd[k]
        y[k] = even[k] + t               # A(w^k)         = E(w^2k) + w^k O(w^2k)
        y[k + n // 2] = even[k] - t      # A(w^(k + n/2)) = E(w^2k) - w^k O(w^2k)
    return y

Iterative, in place, with bit reversal

The recursive version allocates new lists at every level and makes about n function calls. Production code unrolls the recursion. If you follow the even/odd splits down to single elements, the leaf order is the bit reversal of each index: with n = 8, index 1 (001) lands at position 4 (100), and index 3 (011) at position 6 (110). So first permute the array into bit-reversed order, then run the combine steps bottom-up: block length 2, then 4, then up to n. Each stage overwrites values in place.

import math

def fft_inplace(a, invert=False):
    """Iterative in-place radix-2 FFT on a list of complex numbers."""
    n = len(a)
    assert n and n & (n - 1) == 0, "length must be a power of two"
    # 1. bit-reversal permutation: element at index i moves to reverse_bits(i)
    j = 0
    for i in range(1, n):
        bit = n >> 1
        while j & bit:
            j ^= bit
            bit >>= 1
        j |= bit
        if i < j:
            a[i], a[j] = a[j], a[i]
    # 2. butterflies, block length 2, 4, 8, ..., n
    length = 2
    while length <= n:
        ang = (2 * math.pi / length) * (1 if invert else -1)
        w_len = complex(math.cos(ang), math.sin(ang))
        half = length // 2
        for start in range(0, n, length):
            w = 1 + 0j
            for k in range(half):
                u = a[start + k]
                v = a[start + k + half] * w
                a[start + k] = u + v
                a[start + k + half] = u - v
                w *= w_len
        length <<= 1
    if invert:
        for i in range(n):
            a[i] /= n
    return a

Production code precomputes twiddle factors into a table instead of accumulating w *= w_len, which compounds rounding error.

The inverse transform

Interpolation, going from values back to coefficients, is almost the same computation. The inverse DFT is a[j] = (1/n) * sum over k of y[k] * w^(-j*k). Use the conjugate root and divide by n. That is why fft_inplace takes an invert flag that flips the sign of the angle and scales at the end. Either way, interpolation also costs O(n log n), so the three-step plan is now O(n log n) end to end.

Fast convolution and the padding rule

Multiplying polynomials is the same as convolving their coefficient sequences: c[k] = sum over i of p[i] * q[k - i]. The convolution theorem says the DFT turns convolution into pointwise multiplication. There is a catch, and it causes most FFT bugs. The DFT of length n computes a circular convolution: indices wrap around modulo n. A linear product of lengths len(p) and len(q) has len(p) + len(q) - 1 coefficients. If n is smaller than that, the high coefficients wrap around and are silently added to the low ones. So pad both inputs with zeros to a length n of at least len(p) + len(q) - 1 (rounded up to a power of two for radix-2), and the circular result equals the linear one.

def multiply(p, q):
    """Coefficients of p(x) * q(x) for integer coefficient lists."""
    need = len(p) + len(q) - 1           # length of the LINEAR convolution
    n = 1
    while n < need:
        n <<= 1                          # pad to a power of two >= need
    fa = [complex(x) for x in p] + [0j] * (n - len(p))
    fb = [complex(x) for x in q] + [0j] * (n - len(q))
    fft_inplace(fa)
    fft_inplace(fb)
    fc = [x * y for x, y in zip(fa, fb)] # pointwise product = convolution
    fft_inplace(fc, invert=True)
    out = [round(z.real) for z in fc[:need]]
    worst = max(abs(z.real - r) for z, r in zip(fc[:need], out))
    if worst > 0.25:                     # rounding is no longer trustworthy
        raise ArithmeticError(f"FFT rounding error {worst:.3f}; split coefficients or use NTT")
    return out

assert multiply([1, 2, 3], [4, 5]) == [4, 13, 22, 15]

A worked example by hand

Multiply 1 + 2x + 3x^2 by 4 + 5x. The product has 3 + 2 - 1 = 4 coefficients, so n = 4 needs no extra padding: a = [1, 2, 3, 0] and b = [4, 5, 0, 0]. With n = 4, w = exp(-2*pi*i/4) = -i, so the evaluation points are 1, -i, -1 and i.

Now interpolate with the conjugate points 1, i, -1, -i and divide by 4. For c0: (54 + (-18 + 2i) + (-2) + (-18 - 2i)) / 4 = 16 / 4 = 4. For c1, multiply each product by i^k: 54 + (-18 + 2i)i + (-2)(-1) + (-18 - 2i)(-i) = 54 - 2 - 18i + 2 - 2 + 18i = 52, giving 13. The same steps give c2 = 88 / 4 = 22 and c3 = 60 / 4 = 15. The product is 4 + 13x + 22x^2 + 15x^3, which matches direct expansion: 1*4 = 4, 1*5 + 2*4 = 13, 2*5 + 3*4 = 22, 3*5 = 15. Notice that the imaginary parts cancel exactly. With real inputs the true result is real, and in floating point you take the real part and round.

Precision: when rounding stops being safe

Floating-point FFTs are not exact. Using double precision with 53-bit mantissas, the error in each output grows roughly with the size of the coefficients and slowly with log n. Rounding to the nearest integer recovers the exact product only while that error stays well below 0.5. Keep the largest exact output coefficient many orders of magnitude below 2^53, and do not trust a fixed threshold. Measure the distance from each output to the nearest integer, as multiply does, and fail loudly if it gets close to 0.5.

For big integers, split each number into small limbs, for example base 10^4 or 2^15, so the convolution sums stay small. Run the FFT, round, then propagate carries from the lowest limb upward. For guaranteed exactness, the number-theoretic transform runs the same butterflies in modular arithmetic with roots of unity modulo a prime, and has no rounding at all.

Performance in practice

• Real inputs. If the input is real, the output is conjugate-symmetric, so half of it is redundant. numpy.fft.rfft computes only the n/2 + 1 unique outputs, roughly halving the work.
• Sizes. Radix-2 needs powers of two, but library FFTs use mixed radices (2, 3, 5, 7) and special algorithms for prime lengths. A size with only small prime factors is fast; a large prime length is much slower. Pad to a fast size: scipy.fft.next_fast_len finds one.
• Memory, not arithmetic. FFTs are usually limited by memory bandwidth, because every stage sweeps the whole array with strided access. Libraries such as FFTW benchmark several plans and keep the fastest; reuse plans for repeated sizes.
• When not to bother. For short filters, direct convolution wins. Measure the crossover on your hardware; overlap-add handles streams in fixed-size blocks.
• Complexity. The Big-O article explains why n log n beats n squared at scale.

Where FFTs show up in machine learning

FNet replaced the self-attention sublayer of a Transformer encoder with a two-dimensional Fourier transform over the sequence and hidden dimensions, keeping the real part, as a cheap token mixer. Long-convolution sequence models, such as Hyena, apply filters as long as the input sequence and compute them with FFT convolution, which is what makes those filters affordable. GPU libraries such as cuFFT provide the transforms, and the same rules apply: pad correctly, and watch precision, since half precision makes rounding error far worse.

Testing an FFT implementation

FFT bugs return plausible numbers. Test against an O(n squared) reference on random inputs of awkward sizes, compare with a trusted library on complex data, and check the round trip (ifft(fft(x)) == x within tolerance).

import random
import numpy as np

def naive(p, q):
    out = [0] * (len(p) + len(q) - 1)
    for i, x in enumerate(p):
        for j, y in enumerate(q):
            out[i + j] += x * y
    return out

for trial in range(500):
    p = [random.randint(-1000, 1000) for _ in range(random.randint(1, 300))]
    q = [random.randint(-1000, 1000) for _ in range(random.randint(1, 300))]
    assert multiply(p, q) == naive(p, q)

x = np.random.randn(64) + 1j * np.random.randn(64)
assert np.allclose(fft_inplace(list(x)), np.fft.fft(x))   # same sign convention as NumPy

Add length 1, all zeros and a single impulse (whose transform is constant), and undersize the padding on purpose to confirm the test catches circular aliasing.

What to do next

1. Type in the recursive fft and check it against numpy.fft.fft on random vectors of length 8, 64 and 1024.
2. Implement fft_inplace, print the bit-reversal permutation for n = 16, and confirm the outputs match the recursive version.
3. Write multiply, reproduce the 4 + 13x + 22x^2 + 15x^3 example, and run the 500-trial randomized test against the naive product.
4. Deliberately pad to len(p) + len(q) - 2 and watch the test fail, so you recognise circular aliasing when you see it.
5. Multiply two 10,000-digit integers using base-10^4 limbs, then increase the limb size until the rounding guard fires.