The Fast Fourier Transform

Sixteen multiplications for one product

Multiply two polynomials:

$$A(x)=5+2x+7x^2+3x^3,$$

$$B(x)=4+6x^2+x^3.$$

Store each as its list of coefficients, lowest power first: $a$ = [5, 2, 7, 3] and $b$ = [4, 0, 6, 1]. The product has degree 6, and its coefficient of $x^k$ collects every pair of terms whose powers add up to $k$:

$$c_k=\sum _{i+j=k}{a}_i{b}_j.$$

This list $c$ is called the convolution of $a$ and $b$, written $c=a*b$. The schoolbook method computes it by multiplying every $a_i$ by every $b_j$.

import cmath, math, random

def schoolbook(a, b, count):
    """The convolution c = a * b; count[0] gains one per coefficient product."""
    c = [0] * (len(a) + len(b) - 1)
    for i, x in enumerate(a):
        for j, y in enumerate(b):
            count[0] += 1
            c[i + j] += x * y
    return c

A, B = [5, 2, 7, 3], [4, 0, 6, 1]
count = [0]
C = schoolbook(A, B, count)
assert C == [20, 8, 58, 29, 44, 25, 3] and count[0] == 16
assert C[3] == 5 * 1 + 2 * 6 + 7 * 0 + 3 * 4 == 29

Two lists of $n$ coefficients take $n^2$ multiplications. That is 16 here and ${10}^{10}$ at $n=100,000$. The same sum appears whenever numbers are multiplied digit by digit: the integer 3725 is $A(10)$ with its digits as coefficients, so multiplying two 100,000-digit integers is a convolution of their digit lists followed by carrying. This module computes a convolution of length $n$ in $O(n\log n)$ arithmetic operations.

Values multiply; coefficients don't

A polynomial can also be described by its values. Evaluate $A$ and $B$ at $x=2$: $A(2)=61$ and $B(2)=36$, so $AB(2)=61·36=2196$, one multiplication. Values multiply pointwise.

And values determine a polynomial. If two polynomials of degree less than $N$ agree at $N$ distinct points, their difference has degree less than $N$ but $N$ roots, so it is the zero polynomial. (That uses the fact that a nonzero polynomial of degree $d$ has at most $d$ roots, which holds over the complex numbers and over the integers modulo a prime; Polynomials in Algorithm Design proves it and builds much more on it.) $AB$ has degree 6, so 7 of its values fix it. That suggests a plan:

1. evaluate $A$ and $B$ at $N\ge 7$ points,
2. multiply the values pointwise, $N$ multiplications,
3. interpolate: recover the 7 coefficients of $AB$ from its $N$ values.

For two inputs of length $n$, the product has $2n-1$ coefficients, so $N\ge 2n-1$.

def horner(a, x):
    v = 0
    for coef in reversed(a):
        v = v * x + coef
    return v

assert horner(A, 2) == 61 and horner(B, 2) == 36
assert horner(C, 2) == 61 * 36 == 2196

The plan moves the waste instead of removing it. Horner's rule costs $n$ multiplications per point, so evaluating at $N$ arbitrary points costs about $Nn$, and interpolating back from arbitrary points costs $\Theta (N^2)$ by Lagrange's formula. Each point is evaluated from scratch: nothing computed for $x=2$ helps with $x=3$. The rest of this module chooses the $N$ points so that their evaluations share almost all of their work.

The model

We count arithmetic operations on exact values of a number system in which the points we need exist: the complex numbers, treated as exact, or the integers modulo a prime $p$. The operation we count is multiplication. Every multiplication in the algorithm comes with one addition and one subtraction, so additions change the total by a constant factor. The algorithm flips no coins, and every bound below is worst case over all inputs, for a length $N$ that is a power of two (shorter inputs are padded with zeros, which at most doubles $N$). What happens when complex numbers are rounded to floating point is a separate question, answered in its own section.

Roots of unity and the even/odd split

A number $\omega$ is a primitive $N$-th root of unity if ${\omega }^N=1$ and ${\omega }^j\ne 1$ for $0<j<N$. Over the complex numbers, $\omega =e^{2\pi i/N}$, and its powers ${\omega }^0,{\omega }^1,\dots ,{\omega }^{N-1}$ are $N$ equally spaced points on the unit circle. The discrete Fourier transform (DFT) of a coefficient list $(a_0,\dots ,a_{N-1})$ is the list of values

$$(A({\omega }^0),A({\omega }^1),\dots ,A({\omega }^{N-1})).$$

Two facts about these points do all the work, for $N$ even:

• Halving lemma. ${\omega }^2$ is a primitive $(N/2)$-th root of unity. So squaring the $N$ points gives the $N/2$ points of the half-size problem, each twice.
• Opposite points. ${\omega }^{N/2}=-1$, so ${\omega }^{k+N/2}=-{\omega }^k$.

The second fact is where primitivity matters: $({\omega }^{N/2}{)}^2=1$ gives ${\omega }^{N/2}=\pm 1$, and $+1$ is ruled out because $N/2<N$.

Now split $A$ by the parity of its powers. Let $E$ hold the even-indexed coefficients and $O$ the odd-indexed ones, each a polynomial with $N/2$ coefficients, so that

$$A(x)=E(x^2)+xO(x^2).$$

At $x={\omega }^k$ and at $x={\omega }^{k+N/2}=-{\omega }^k$ the square is the same point ${\omega }^{2k}$. So for $0\le k<N/2$:

$$A({\omega }^k)=E({\omega }^{2k})+{\omega }^{k}O({\omega }^{2k}),$$

$$A({\omega }^{k+N/2})=E({\omega }^{2k})-{\omega }^{k}O({\omega }^{2k}).$$

The values $E({\omega }^{2k})$ for $k<N/2$ are exactly the DFT of $E$ at the $(N/2)$-th roots, and likewise for $O$. One product ${\omega }^k·O({\omega }^{2k})$, one addition and one subtraction give two outputs. That pair of lines is called a butterfly.

Recursing on $E$ and $O$ gives the fast Fourier transform (FFT), in the form published by Cooley and Tukey in 1965. Here it is as a recursion, with the root chosen by a helper so that the same code runs over the complex numbers or modulo a prime $p$.

def root(N, k, invert, p):
    """omega^k (or omega^-k) for a primitive N-th root omega: complex, or mod p."""
    if p is None:
        return cmath.exp((-2j if invert else 2j) * math.pi * k / N)
    w = pow(3, (p - 1) // N, p)                  # why 3 works: see "Exact or approximate"
    return pow(w, (N - k) % N if invert else k, p)

def fft(a, invert=False, p=None, count=None):
    """The DFT of a (len(a) a power of two); count[0] gains one per butterfly product."""
    N = len(a)
    if N == 1:
        return list(a)
    E = fft(a[0::2], invert, p, count)
    O = fft(a[1::2], invert, p, count)
    out = [0] * N
    for k in range(N // 2):
        t = root(N, k, invert, p) * O[k]
        if count is not None:
            count[0] += 1
        out[k], out[k + N // 2] = E[k] + t, E[k] - t
        if p is not None:
            out[k], out[k + N // 2] = out[k] % p, out[k + N // 2] % p
    return out

Tracing $A$ with eight points

Pad $a$ to length 8 (the product needs $N\ge 7$) and split it all the way down. Each level separates the even and odd positions of every piece.

The leaf order is not an accident. Write the index in three bits. The first split sorts by the last bit, the second by the middle bit, the third by the first bit, so the leaf in position $j$ holds the coefficient whose index is $j$ with its bits reversed: position 1 is 001, which holds index 100 = 4. An iterative FFT therefore permutes its input into this bit-reversed order once, then combines neighbouring blocks of length 2, 4, 8, … in place.

def leaf_order(idx):
    """The leaves of the even/odd recursion on a list of indices, left to right."""
    if len(idx) == 1:
        return idx
    return leaf_order(idx[0::2]) + leaf_order(idx[1::2])

assert (A + [0] * 4)[0::2] == [5, 7, 0, 0] and (A + [0] * 4)[1::2] == [2, 3, 0, 0]
assert leaf_order(list(range(8))) == [0, 4, 2, 6, 1, 5, 3, 7]
assert all(int(format(j, "03b")[::-1], 2) == leaf for j, leaf in enumerate(leaf_order(list(range(8)))))

Now the top butterflies. $E(y)=5+7y$ and $O(y)=2+3y$ (the zeros drop out), evaluated at the fourth roots $1,i,-1,-i$, which are ${\omega }^{2k}$ for $\omega =e^{2\pi i/8}$ and $k=0,1,2,3$. Each row is one butterfly, giving outputs $k$ and $k+4$.

Row 2 is worth doing by hand: ${\omega }^2=i$, so $A(i)=E(-1)+iO(-1)=-2-i$, and indeed $5+2i-7-3i=-2-i$. Row 0 is the pair $17,7$ from the prediction above.

E4 = [5 + 7 * y for y in (1, 1j, -1, -1j)]
O4 = [2 + 3 * y for y in (1, 1j, -1, -1j)]
assert E4 == [12, 5 + 7j, -2, 5 - 7j] and O4 == [5, 2 + 3j, -1, 2 - 3j]
w8 = cmath.exp(2j * math.pi / 8)
vals = fft(A + [0] * 4)
for k in range(8):
    assert abs(vals[k] - horner(A, w8 ** k)) < 1e-9                  # the DFT is evaluation
    assert abs(vals[k] - (E4[k % 4] + w8 ** k * O4[k % 4])) < 1e-9   # ... joined by butterflies
assert abs(vals[2] - (-2 - 1j)) < 1e-12 and abs(vals[0] - 17) < 1e-12 and abs(vals[4] - 7) < 1e-12
assert [round(v.real, 2) for v in vals[1:4:2]] == [4.29, 5.71]
assert [round(vals[1].imag, 2), round(vals[5].imag, 2)] == [10.54, 3.46]

Why it is correct, and how to go back

Correctness of the transform is induction on $N$. For $N=1$ the DFT of $(a_0)$ is $A({\omega }^0)=a_0$. For larger $N$, the recursive calls return the DFTs of $E$ and $O$ at the $(N/2)$-th roots (by the halving lemma, ${\omega }^2$ is the right root for them), and the two butterfly lines are the identity $A(x)=E(x^2)+xO(x^2)$ at $x=\pm {\omega }^k$. The step that must not be skipped is ${\omega }^{N/2}=-1$: with a root that is not primitive, output $k+N/2$ is not the negated partner of output $k$, and the halves do not share their work.

The inverse is the same algorithm. Write the DFT as a matrix $V$ with entries $V_{jk}={\omega }^{jk}$. The product of $V$ with the matrix of ${\omega }^{-1}$ has entries

$$\sum _{k=0}^{N-1}{\omega }^{(j-l)k}.$$

If $j=l$ every term is 1 and the sum is $N$. Otherwise the terms form a geometric series with ratio $r={\omega }^{j-l}\ne 1$ (primitivity again), and its sum is $(r^N-1)/(r-1)=0$ because $r^N=1$. So the inverse of $V$ is $\frac{1}{N}$ times the matrix of ${\omega }^{-1}$: interpolate by running the FFT with ${\omega }^{-1}$ and dividing by $N$.

The product follows. Evaluation respects multiplication, so the pointwise product of the two DFTs is the DFT of the product polynomial, provided that polynomial has fewer than $N$ coefficients. Here is the whole algorithm, over the complex numbers (rounded at the end) or modulo $p$:

def multiply(a, b, p=None, count=None):
    m = len(a) + len(b) - 1
    N = 1
    while N < m:
        N *= 2
    fa = fft(a + [0] * (N - len(a)), False, p, count)
    fb = fft(b + [0] * (N - len(b)), False, p, count)
    if count is not None:
        count[0] += 2 * N                          # N pointwise products, N scalings
    prod = [x * y if p is None else x * y % p for x, y in zip(fa, fb)]
    back = fft(prod, True, p, count)
    if p is None:
        return [round((v / N).real) for v in back[:m]]
    n_inv = pow(N, p - 2, p)                       # 1/N mod p, by Fermat's little theorem
    return [v * n_inv % p for v in back[:m]]

P = 998244353
assert multiply(A, B) == C and multiply(A, B, P) == C

Two mistakes produce plausible wrong answers. With too few points, $N=4$ for our pair, the product's seven coefficients do not fit: $x^4=1$ at every fourth root of unity, so $x^4$ folds onto $x^0$, $x^5$ onto $x^1$, and so on. The result is the cyclic convolution $c_k+c_{k+4}$, [64, 33, 61, 29]. Forgetting to divide by $N$ returns every coefficient eight times too large.

fa, fb = fft(A, False, P), fft(B, False, P)                    # N = 4: too few points
wrapped = [v * pow(4, P - 2, P) % P for v in fft([x * y % P for x, y in zip(fa, fb)], True, P)]
assert wrapped == [64, 33, 61, 29] == [C[k] + (C[k + 4] if k + 4 < 7 else 0) for k in range(4)]

Complexity

Assume the required roots, the twiddle factors ${\omega }^k$, are precomputed; we count multiplications of transform values by those factors. Under this convention one length-$N$ radix-2 transform performs exactly $\frac{N}{2}{\log }_{2}N$ such multiplications, as we now show. (The code above computes each factor with root, which the counter does not see.) Building the twiddle tables also costs arithmetic, but does not change the $O(N\log N)$ bound.

Let $T(N)$ count these multiplications for one transform. The two recursive calls cost $T(N/2)$ each and the butterflies $N/2$, so $T(N)=2T(N/2)+N/2$ with $T(1)=0$. The recursion has ${\log }_{2}N$ levels, each doing $N/2$ butterflies, so exactly

$$T(N)=\frac{N}{2}{\log }_{2}N.$$

Compare Selection in Linear Time: median of medians satisfies $T(n)\le T(n/5)+T(7n/10)+cn$, whose subproblems add up to $9n/10$, so the level sums shrink geometrically and the total is linear. Here the subproblems add up to exactly $N$, every level costs the same, and the $\log N$ factor stays.

A product makes two forward transforms, $N$ pointwise products, one inverse transform and $N$ scalings. With $N\ge 2n-1$ a power of two, and in the same convention, that is

$$3·\frac{N}{2}{\log }_{2}N+2N$$

multiplications.

At $n=100,000$ the next power of two above $2n-1$ is $N=2^{18}$, and the count is 7,602,176, against ${10}^{10}$ for the schoolbook method. For small $n$ the schoolbook method wins: at $n=16$ it needs 256 multiplications and the transform 304. From $n=32$ on (1,024 against 704) the transform is cheaper. Space is $O(N)$.

def product_cost(n):
    N = 1
    while N < 2 * n - 1:
        N *= 2
    return 3 * (N // 2) * (N.bit_length() - 1) + 2 * N

for n in (16, 32, 64, 1024):
    count = [0]
    x = [random.Random(n).randrange(1000) for _ in range(n)]
    assert multiply(x, x, P, count) == schoolbook(x, x, [0])
    assert count[0] == product_cost(n)
assert (product_cost(16), product_cost(32)) == (304, 704)
assert product_cost(100_000) == 7_602_176

Is $N\log N$ the end? For the DFT, $O(N\log N)$ is the best known algorithm. No lower bound above linear is known in general arithmetic models. Morgenstern (1973) proved $\Omega (N\log N)$ for a restricted class: linear algorithms whose constants are bounded in absolute value. So the FFT is optimal among algorithms of that kind, and nobody knows whether something else could be faster.

Exact or approximate

Modulo a prime. The integers mod $p$ have a primitive $N$-th root of unity exactly when $N$ divides $p-1$. The prime $p=998244353$ is chosen for that: $p-1=119·2^{23}$, and 3 generates all of the nonzero residues, so $\omega =3^{(p-1)/N}$ is a primitive $N$-th root for every power of two $N\le 2^{23}$. The transform modulo $p$ (often called the number-theoretic transform, NTT) makes no rounding error at all. It returns each $c_{k}modp$, which is $c_k$ itself whenever $0\le c_k<p$. Past that it is silently wrong as an integer: ${40000}^2$ comes back as 601,755,647. Larger coefficients need two or three primes combined by the Chinese remainder theorem, or smaller digits.

assert P - 1 == 119 * 2 ** 23 and 119 == 7 * 17
assert all(pow(3, (P - 1) // q, P) != 1 for q in (2, 7, 17))      # 3 generates mod P
w = pow(3, (P - 1) // 8, P)
assert pow(w, 8, P) == 1 and pow(w, 4, P) == P - 1                 # primitive 8th root
assert multiply([40000], [40000], P) == [601_755_647] == [40000 ** 2 - P]

Big integers. Take $90417263\times 58820419$. Their digit lists, least significant first, convolve to 15 sums of which the largest is 190, far below $p$, so the transform mod $p$ gives them exactly. Carrying from the lowest digit up turns the sums into the product.

def big_product(x, y):
    digits = lambda z: [int(ch) for ch in str(z)[::-1]]
    sums = multiply(digits(x), digits(y), P)
    out, carry = [], 0
    for s in sums:
        carry += s
        out.append(carry % 10)
        carry //= 10
    while carry:
        out.append(carry % 10)
        carry //= 10
    return int("".join(map(str, reversed(out)))), sums

value, sums = big_product(90417263, 58820419)
assert sums == [27, 57, 36, 89, 30, 101, 84, 190, 113, 118, 75, 55, 92, 72, 45]
assert value == 90417263 * 58820419 == 5318381294493197

Over the complex numbers in floating point. A double has a 53-bit significand, every butterfly rounds, and rounding the final values to integers is right only while the accumulated error stays below $1/2$. That error grows with the size of the coefficients and with the number of levels. Having every exact coefficient below $2^{53}$ is not enough. Here are two random lists of 8,192 entries below ${10}^6$: every coefficient of the product is below $2.1·{10}^{15}$, exactly representable as a double, and the rounded complex product still gets thousands of them wrong.

def exact_product(a, b, bits=64):
    """The exact reference: pack each list into one Python integer and multiply."""
    ia = sum(v << (bits * i) for i, v in enumerate(a))
    ib = sum(v << (bits * i) for i, v in enumerate(b))
    prod, mask = ia * ib, (1 << bits) - 1
    return [(prod >> (bits * i)) & mask for i in range(len(a) + len(b) - 1)]

rng = random.Random(2025)
wrong = {}
for n in (512, 2048, 8192):
    x = [rng.randrange(10 ** 6) for _ in range(n)]
    y = [rng.randrange(10 ** 6) for _ in range(n)]
    exact = exact_product(x, y)
    assert max(exact) < 2.1e15 < 2 ** 53
    wrong[n] = sum(g != e for g, e in zip(multiply(x, y), exact))
assert wrong == {512: 0, 2048: 0, 8192: 2239}

The failure is not a bug in the code; the same lists multiplied mod $p$ would still be wrong, but for the other reason (these coefficients exceed $p$). A floating-point FFT can be made safe by splitting numbers into small enough pieces that a proven error bound stays below $1/2$.

Measure the claim

Multiplications for two lists of $n$ coefficients, counted by the code above (the transform column matches $3·\frac{N}{2}{\log }_{2}N+2N$ exactly):

The last column falls towards $3/2$, since the ratio is exactly $3/2+2/{\log }_{2}N$: the growth is $N\log N$ with constant $3/2$, while the schoolbook column grows sixteenfold for every fourfold step in $n$. Rounded complex products of random entries below ${10}^6$ were exact at $n=512$ and $n=2,048$ and wrong in 2,239 of 16,383 coefficients at $n=8,192$ (one seed; a measurement, not a theorem).

assert [round(product_cost(n) / (2 * n * math.log2(2 * n)), 2) for n in (16, 32, 64, 1024, 4096)] \
    == [1.90, 1.83, 1.79, 1.68, 1.65]
assert all(abs(product_cost(n) / (2 * n * math.log2(2 * n)) - (1.5 + 2 / math.log2(2 * n))) < 1e-12
           for n in (16, 4096))
assert product_cost(4096) == 176_128

In real software

GMP, the GNU multiple-precision library, multiplies large integers with a sequence of algorithms by size: schoolbook, Karatsuba, several Toom–Cook variants, and at the largest sizes an FFT in the style of Schönhage and Strassen. That FFT works exactly, modulo $2^{N^{\prime }}+1$, not in floating point. CPython's int does not use an FFT: once both factors have more than 70 of its internal 30-bit digits (140 when squaring), it switches from the schoolbook method to Karatsuba's $O(n^{1.59})$ algorithm. NumPy's numpy.fft computes floating-point complex DFTs with the pocketfft library, so a product computed through it rounds, as ours did.

A problem that looks different

Two lists each hold 100,000 integers between 0 and 1,000,000. For every target $t$ you want to know how many pairs, one number from each list, add up to $t$. A double loop makes ${10}^{10}$ additions. Is there a polynomial hiding in this question? (Not solved here; Convolutions and Their Applications takes up questions of this kind, and the lab's last problem is a different one.)

Practise

In the lab you evaluate a polynomial by splitting it once and joining the halves with butterflies, watching each pair of outputs appear; predict and then trace where every coefficient lands when a list is split all the way down; write the iterative transform modulo a prime, its inverse and an exact product; count the multiplications of a product and derive the constant yourself; and finish with a problem that doesn't say what it is.

Recap

• You can now: multiply polynomials and big integers in $O(n\log n)$ arithmetic operations by evaluating at roots of unity, multiplying pointwise and interpolating with ${\omega }^{-1}$ and $1/N$; explain the halving lemma, the butterfly and bit-reversed order; and say when the result is exact.
• Invariant: a length-$N$ DFT is the two length-$N/2$ DFTs of the even and odd coefficients, joined by butterflies $E({\omega }^{2k})\pm {\omega }^{k}O({\omega }^{2k})$.
• Complexity achieved: exactly $\frac{N}{2}{\log }_{2}N$ multiplications per transform and $3·\frac{N}{2}{\log }_{2}N+2N$ per product ($N\ge 2n-1$ a power of two), worst case, counting operations on exact complex numbers or residues mod $p$; against $n^2$ for the schoolbook method.
• Failure mode: too few points (the product wraps around into a cyclic convolution), a forgotten $1/N$, trusting the result mod $p$ when coefficients reach $p$, or trusting rounded floating point on large inputs.
• In real software: GMP's largest multiplications use a Schönhage–Strassen FFT modulo $2^{N^{\prime }}+1$; CPython's int stops at Karatsuba; numpy.fft is floating point (pocketfft).
• Next: Convolutions and Their Applications asks what else a convolution computes, starting from the question above; Polynomials in Algorithm Design returns to polynomials, and to why their values are such good fingerprints.

Check yourself

1. Why do the $N$-th roots of unity make the recursion work, when evaluating at $0,1,\dots ,N-1$ would not? Point to the line of the argument that needs ${\omega }^{N/2}=-1$.
2. Replace the transform mod $p$ by complex floating point and round at the end. When does it break, and what would you change to make it safe for 100,000-digit integers?
3. You need the product of a polynomial with 100,000 coefficients and one with 5. Schoolbook or transform? Count both.

Topic list modeled on the public schedule of CMU 15-451/651 Algorithm Design and Analysis (Fall 2025). All lecture text, proofs, examples, code, and exercises are original to SciMigo; no CMU course materials are reproduced. This course is independent and not affiliated with or endorsed by Carnegie Mellon University.