Fast Fourier Transform

Karatsuba cut integer multiplication below quadratic by trading one multiplication for a few additions, reaching $\Theta (n^{{\log }_{2}3})\approx \Theta (n^{1.585})$. This lesson goes further, to $\Theta (n\log n)$, with the Fast Fourier Transform. The route is indirect. The obstacle to fast multiplication is the coefficient representation of a polynomial, in which the product is a convolution costing $\Theta (n^2)$. If instead we represent a polynomial by its values at enough points, the product is computed pointwise in linear time. The FFT converts between the two representations in $\Theta (n\log n)$, and the points it chooses — the complex roots of unity — are what make that conversion fast.

The problem: polynomial multiplication is convolution

A polynomial of degree $n-1$ is given by its coefficient vector $a=(a_0,\dots ,a_{n-1})$, meaning $A(x)={\sum }_{j=0}^{n-1}{a}_j{x}^j$. The product $C(x)=A(x)B(x)$ has coefficients

This sum is the convolution $a\ast b$. Computed directly, each of the $\sim 2n$ output coefficients is a sum over up to $n$ products, so the cost is $\Theta (n^2)$ — the schoolbook method, and exactly the cost of grade-school long multiplication when $a$ and $b$ are the digit vectors of two integers. Convolution is everywhere (signal filtering, string matching, probability of sums), so a faster algorithm pays off far beyond polynomials.

A polynomial is its values

A degree-$(n-1)$ polynomial is uniquely determined by its values at any $n$ distinct points — the interpolation fact that $n$ points fix a degree-$(n-1)$ curve (two points fix a line, three fix a parabola). So besides the coefficient vector there is a second, equally faithful description: the point-value representation, a list of pairs $\{(x_0,A(x_0)),\dots ,(x_{n-1},A(x_{n-1}))\}$.

In this representation, multiplication is cheap. If $A$ and $B$ are both sampled at the same points $x_i$, then the product $C=AB$ is sampled there by a single multiply per point:

That is $\Theta (n)$ work, not $\Theta (n^2)$. The catch is that the product has degree up to $2n-2$, so $2n-1$ points are needed to determine it; we must therefore sample $A$ and $B$ at (at least) $2n$ points up front, padding their coefficient vectors with zeros to length $2n$. The strategy is now a three-step detour:

Evaluate, multiply pointwise, interpolate. The middle step is linear, so the total cost depends on the two end steps: evaluating a polynomial at $2n$ points, and recovering coefficients from $2n$ values. Naively each is $\Theta (n^2)$ (evaluating at one point by Horner's rule is $\Theta (n)$, times $2n$ points), so a poor choice of points gains nothing. For one special set of points, both steps drop to $\Theta (n\log n)$.

The roots of unity

The special points are the $n$-th roots of unity: the $n$ complex numbers $\omega$ with ${\omega }^n=1$. They are equally spaced around the unit circle in the complex plane,

where ${\omega }_n$ is the principal $n$-th root and its powers ${\omega }_n^0,{\omega }_n^1,\dots ,{\omega }_n^{n-1}$ enumerate all $n$ of them. Geometrically, multiplying by ${\omega }_n$ rotates a point by $2\pi /n$ radians; $n$ such rotations come full circle back to the start, which is the statement ${\omega }_n^n=1$.

Two algebraic properties of these points make the FFT work.

The $n$ evaluation points, when squared, become only $n/2$ distinct points, and they pair up as $\pm$ of each other. These two facts turn one size-$n$ evaluation into two size-$(n/2)$ evaluations.

The DFT and the Cooley–Tukey split

Evaluating $A$ at all $n$ roots of unity is, by definition, the Discrete Fourier Transform of its coefficient vector.

Computed as written, that is $n$ evaluations of an $n$-term sum: $\Theta (n^2)$. The FFT computes the same vector in $\Theta (n\log n)$ by a divide-and-conquer split on the parity of the coefficient index. Separate $A$ into its even- and odd-indexed coefficients and factor an $x$ out of the odd part:

where $A_{\text{even}}$ and $A_{\text{odd}}$ are each polynomials of half the degree, in the variable $x^2$. Repeating the split recursively builds a tree: at each level the input is sorted by parity into two halves, dark for the even indices and light for the odd, until single coefficients sit at the leaves.

Now evaluate at a root ${\omega }_n^k$. By the halving lemma, $({\omega }_n^k{)}^2={\omega }_{n/2}^k$, so

The right-hand side asks for $A_{\text{even}}$ and $A_{\text{odd}}$ evaluated at the $\frac{n}{2}$-th roots of unity — two DFTs of size $n/2$. And here the $\pm$ pairing pays off: the point ${\omega }_n^{k+n/2}=-{\omega }_n^k$ squares to the same ${\omega }_{n/2}^k$, so the two half-size results serve both $k$ and $k+n/2$ at once, with only a sign flip on the odd part:

These paired equations are the butterfly: two outputs assembled from two inputs by one multiply (${\omega }_n^k$ times the odd value, the twiddle factor) and one add/subtract. A full level of $n/2$ butterflies combines two half-DFTs into one full DFT in $\Theta (n)$ work.

Recursing this split bottoms out at size-$1$ DFTs (a single coefficient evaluates to itself), and at every level above it does $\Theta (n)$ combining work, giving the recurrence

This is the mergesort recurrence. By the Master Theorem ($a=2$, $b=2$, $f(n)=\Theta (n)$, the balanced case ${\log }_{b}a=1$), it solves to

import cmath
from collections.abc import Sequence

def next_power_of_two(value: int) -> int:
  """
    Smallest power of two that is >= `value` (and >= 1).\n
  """
  size: int = 1
  while size < value:
    size <<= 1
  return size

def naive_dft(coefficients: Sequence[complex]) -> list[complex]:
  """
    The Discrete Fourier Transform computed straight from the definition,\n
    hat_a[k] = sum_j a[j] * omega^(k j), in O(n^2). Serves as the slow but\n
    obviously-correct reference the fast transform is checked against.\n
  """
  length: int = len(coefficients)
  if length == 0:
    return []

  # the principal n-th root of unity all frequencies are sampled against.
  principal_root: complex = cmath.exp(2j * cmath.pi / length)

  # for each frequency k, sum a[j] * omega^(k j) over every coefficient.
  transform: list[complex] = []
  for frequency in range(length):
    total: complex = sum(
      coefficient * principal_root ** (frequency * index)
      for index, coefficient in enumerate(coefficients)
    )
    transform.append(total)

  return transform

def _transform(coefficients: list[complex], inverse: bool) -> list[complex]:
  """
    The shared recursive butterfly engine. `inverse` flips the sign of the\n
    root's angle (conjugate roots); the 1/n scaling is applied once by the\n
    public `inverse_fft` wrapper, not here, so the recursion stays uniform.\n
  """
  # a size-1 DFT is the identity: one coefficient evaluates to itself.
  length: int = len(coefficients)
  if length == 1:
    return list(coefficients)

  # recurse on the even- and odd-indexed halves.
  even_half: list[complex] = _transform(coefficients[0::2], inverse)
  odd_half: list[complex] = _transform(coefficients[1::2], inverse)

  # principal root for this level; angle flips sign when inverting.
  angle: float = (-2.0 if inverse else 2.0) * cmath.pi / length
  principal_root: complex = cmath.exp(1j * angle)
  root: complex = 1 + 0j

  # combine halves with butterflies, each filling the k and k + n/2 outputs.
  transform: list[complex] = [0j for _ in range(length)]
  half: int = length // 2
  for frequency in range(half):
    twiddle: complex = root * odd_half[frequency]
    transform[frequency] = even_half[frequency] + twiddle
    transform[frequency + half] = even_half[frequency] - twiddle
    root *= principal_root

  return transform

def _pad_to_power_of_two(coefficients: Sequence[complex]) -> list[complex]:
  """
    Copy `coefficients` into a fresh list zero-padded up to the next power\n
    of two, the length the radix-2 split requires.\n
  """
  # copy into a zero-filled buffer sized to the next power of two.
  target: int = next_power_of_two(max(1, len(coefficients)))
  padded: list[complex] = [0j for _ in range(target)]
  for index, coefficient in enumerate(coefficients):
    padded[index] = coefficient

  return padded

def fft(coefficients: Sequence[complex]) -> list[complex]:
  """
    Forward FFT: evaluate the polynomial with these coefficients at the N-th\n
    roots of unity, where N is the next power of two >= len(coefficients).\n
    The input is zero-padded as needed, so the output has length N.\n
  """
  return _transform(_pad_to_power_of_two(coefficients), inverse=False)

def inverse_fft(values: Sequence[complex]) -> list[complex]:
  """
    Inverse FFT: recover coefficients from values sampled at the roots of\n
    unity. Runs the same butterflies on the conjugate roots, then divides\n
    every output by N. `values` must already have power-of-two length.\n
  """
  length: int = len(values)
  if length == 0:
    return []

  # invert on the conjugate roots, then divide every output by N.
  transform: list[complex] = _transform(list(values), inverse=True)
  return [entry / length for entry in transform]

The algorithm requires $n$ to be a power of two so the halving is exact; pad the coefficient vector with zeros up to the next power of two (and to at least $2n$, for the product's degree) before transforming. A non-recursive, in-place variant first permutes the inputs into bit-reversed order and then runs the butterflies bottom-up, which is what production libraries implement, but the recursion above is the cleanest statement of why it is $\Theta (n\log n)$.

Inverting the transform

Evaluation is half the round trip; we still must turn the $2n$ product values back into coefficients — the inverse DFT. The DFT is a linear map, multiplication of $a$ by the matrix $V$ with $V_{kj}={\omega }_n^{kj}$ (a Vandermonde matrix at the roots of unity). Its inverse has a simple form.

The consequence is that the inverse transform is the same algorithm: replace ${\omega }_n$ by its conjugate ${\omega }_n^{-1}=e^{-2\pi i/n}$, run $\text{FFT}$ unchanged, and divide every output by $n$. We never invert a matrix or solve a linear system — interpolation at the roots of unity is just another FFT.

Putting it together: fast multiplication

The three-step plan now has an $\Theta (n\log n)$ cost on every leg. Pad $a$ and $b$ to length $N$, the next power of two $\ge 2n$; transform both; multiply the two value vectors pointwise; inverse-transform. Three FFTs plus one linear pointwise multiply:

from collections.abc import Sequence

from fft import fft, inverse_fft, next_power_of_two

def naive_convolution(left: Sequence[float], right: Sequence[float]) -> list[float]:
  """
    The schoolbook convolution c[k] = sum_j left[j] * right[k - j], in\n
    O(n m). Exact reference for the FFT-based product.\n
  """
  if not left or not right:
    return []

  # accumulate every term left[j] * right[i] into output slot j + i.
  result: list[float] = [0.0 for _ in range(len(left) + len(right) - 1)]
  for left_index, left_value in enumerate(left):
    for right_index, right_value in enumerate(right):
      result[left_index + right_index] += left_value * right_value

  return result

def poly_multiply(left: Sequence[float], right: Sequence[float]) -> list[float]:
  """
    Multiply two real coefficient vectors (lowest degree first) via three\n
    FFTs and a pointwise product. The result has length len(left) +\n
    len(right) - 1; the imaginary parts of the inverse transform are rounding\n
    noise and are discarded.\n
  """
  if not left or not right:
    return []

  result_length: int = len(left) + len(right) - 1

  # pad to a power of two large enough to hold the product's degree.
  size: int = next_power_of_two(result_length)
  left_padded: list[complex] = [complex(value) for value in left] + \
    [0j for _ in range(size - len(left))]
  right_padded: list[complex] = [complex(value) for value in right] + \
    [0j for _ in range(size - len(right))]

  left_values: list[complex] = fft(left_padded)
  right_values: list[complex] = fft(right_padded)

  product_values: list[complex] = [
    left_value * right_value
    for left_value, right_value in zip(left_values, right_values)
  ]

  coefficients: list[complex] = inverse_fft(product_values)
  return [coefficient.real for coefficient in coefficients[:result_length]]

def poly_multiply_round(
  left: Sequence[int],
  right: Sequence[int],
) -> list[int]:
  """
    Integer-coefficient product: multiply via FFT, then round each output to\n
    the nearest integer. Correct as long as the coefficients stay well below\n
    the float precision limit.\n
  """
  return [round(value) for value in poly_multiply(left, right)]

Worked example (a size-$4$ transform). Take $A(x)=1+2x$, coefficient vector $a=(1,2,0,0)$ padded to length $4$. The $4$-th roots of unity are ${\omega }_4^0=1,{\omega }_4^1=i,{\omega }_4^2=-1,{\omega }_4^3=-i$. Splitting by parity, $a_{\text{even}}=(1,0)$ gives $E=(1,1)$ (a constant evaluates to itself at both $2$nd roots), and $a_{\text{odd}}=(2,0)$ gives $O=(2,2)$. The butterflies with twiddles ${\omega }_4^0=1$ and ${\omega }_4^1=i$ assemble

Check directly: $A(1)=3$, $A(i)=1+2i$, $A(-1)=-1$, $A(-i)=1-2i$ — the transform is exactly $A$ evaluated at the four roots. To square $A$ (compute $(1+2x{)}^2=1+4x+4x^2$), multiply the transform pointwise — $\hat{c}=(9,(1+2i{)}^2,1,(1-2i{)}^2)=(9,-3+4i,1,-3-4i)$ — then run the inverse FFT (conjugate roots, divide by $4$), recovering $c=(1,4,4,0)$, the coefficients of $1+4x+4x^2$.

This is the asymptotic payoff that beats both schoolbook $\Theta (n^2)$ and Karatsuba's $\Theta (n^{1.585})$. The application that started us off, big-integer multiplication, is the same algorithm: an $n$-digit integer in base $b$ is the polynomial $\sum d_j{b}^j$ evaluated at $x=b$, so multiplying two big integers is multiplying their digit polynomials, then carrying the over-large coefficients — the problem behind Multiply Strings, and at large enough sizes it is what bignum libraries do; the same FFT pipeline (with carries) is the route from Karatsuba's $\Theta (n^{{\log }_{2}3})$ down toward the near-linear $\Theta (n\log n)$ that the divide-and-conquer lesson flagged as the next rung. Add Two Numbers shares the carry mechanics on a per-digit scale.

from collections.abc import Sequence

# a classic NTT-friendly prime: 998244353 = 119 * 2^23 + 1, so it admits
# primitive 2^k-th roots of unity for every k up to 23. 3 is a primitive root.
MODULUS: int = 998244353
PRIMITIVE_ROOT: int = 3

def _next_power_of_two(value: int) -> int:
  """
    Smallest power of two that is >= `value` (and >= 1).\n
  """
  size: int = 1
  while size < value:
    size <<= 1
  return size

def _transform(coefficients: list[int], inverse: bool) -> list[int]:
  """
    Iterative radix-2 NTT in place over a working copy. `inverse` uses the\n
    modular inverse of the root; the 1/n scaling is left to the caller. The\n
    length of `coefficients` must already be a power of two.\n
  """
  length: int = len(coefficients)
  values: list[int] = list(coefficients)

  # bit-reversal permutation so the butterflies can run bottom-up in place.
  position: int = 0
  for index in range(1, length):
    bit: int = length >> 1
    while position & bit:
      position ^= bit
      bit >>= 1
    position ^= bit
    if index < position:
      values[index], values[position] = values[position], values[index]

  # merge blocks bottom-up, doubling the block size each pass.
  block: int = 2
  while block <= length:

    # the primitive block-th root of unity in Z_p (inverted when inverting).
    exponent: int = (MODULUS - 1) // block
    root: int = pow(PRIMITIVE_ROOT, exponent, MODULUS)
    if inverse:
      root = pow(root, MODULUS - 2, MODULUS)

    # butterfly every block, walking the twiddle through its powers of root.
    for start in range(0, length, block):
      twiddle: int = 1
      half: int = block // 2
      for offset in range(half):
        even: int = values[start + offset]
        odd: int = values[start + offset + half] * twiddle % MODULUS
        values[start + offset] = (even + odd) % MODULUS
        values[start + offset + half] = (even - odd) % MODULUS
        twiddle = twiddle * root % MODULUS

    block <<= 1

  return values

def ntt(coefficients: Sequence[int]) -> list[int]:
  """
    Forward NTT modulo `MODULUS`: evaluate at the n-th roots of unity in Z_p,\n
    where n is the next power of two >= len(coefficients). Input is reduced\n
    mod p and zero-padded to length n.\n
  """
  size: int = _next_power_of_two(max(1, len(coefficients)))
  padded: list[int] = [value % MODULUS for value in coefficients]
  padded += [0 for _ in range(size - len(padded))]
  return _transform(padded, inverse=False)

def inverse_ntt(values: Sequence[int]) -> list[int]:
  """
    Inverse NTT: recover coefficients mod `MODULUS`. Runs the transform with\n
    the inverse root, then scales by the modular inverse of n. `values` must\n
    have power-of-two length.\n
  """
  length: int = len(values)
  if length == 0:
    return []
  transformed: list[int] = _transform(list(values), inverse=True)
  scale: int = pow(length, MODULUS - 2, MODULUS)
  return [entry * scale % MODULUS for entry in transformed]

def ntt_convolution(
  left: Sequence[int],
  right: Sequence[int],
) -> list[int]:
  """
    Exact convolution of two integer vectors modulo `MODULUS`, in\n
    O(n log n). The result has length len(left) + len(right) - 1, each entry\n
    reduced mod p.\n
  """
  if not left or not right:
    return []

  result_length: int = len(left) + len(right) - 1
  size: int = _next_power_of_two(result_length)

  left_padded: list[int] = [value % MODULUS for value in left] + \
    [0 for _ in range(size - len(left))]
  right_padded: list[int] = [value % MODULUS for value in right] + \
    [0 for _ in range(size - len(right))]

  left_values: list[int] = _transform(left_padded, inverse=False)
  right_values: list[int] = _transform(right_padded, inverse=False)

  product_values: list[int] = [
    left_value * right_value % MODULUS
    for left_value, right_value in zip(left_values, right_values)
  ]

  coefficients: list[int] = inverse_ntt(product_values)
  return coefficients[:result_length]

from poly_multiply import poly_multiply_round

# base for the digit decomposition; powers of ten keep the conversion to and
# from decimal strings trivial while staying small enough for float precision.
_BASE: int = 10

def _digits_of(magnitude: int) -> list[int]:
  """
    Least-significant-first base-`_BASE` digit list of a non-negative int.\n
  """
  if magnitude == 0:
    return [0]

  # peel off base-_BASE digits from the bottom up.
  digits: list[int] = []
  while magnitude > 0:
    digits.append(magnitude % _BASE)
    magnitude //= _BASE

  return digits

def _carry(raw_coefficients: list[int]) -> int:
  """
    Fold a list of (possibly over-large) base-`_BASE` coefficients, lowest\n
    first, into a single integer by propagating carries.\n
  """
  value: int = 0
  place: int = 1
  carry: int = 0

  # absorb each coefficient plus the running carry into one digit place.
  for coefficient in raw_coefficients:
    total: int = coefficient + carry
    value += (total % _BASE) * place
    carry = total // _BASE
    place *= _BASE

  # flush any carry left above the most significant coefficient.
  while carry > 0:
    value += (carry % _BASE) * place
    carry //= _BASE
    place *= _BASE

  return value

def big_integer_multiply(first: int, second: int) -> int:
  """
    The product `first * second`, computed by convolving the two integers'\n
    digit vectors with the FFT and then carrying. Handles negatives and\n
    zero; the result matches ordinary integer multiplication exactly.\n
  """
  if first == 0 or second == 0:
    return 0

  sign: int = 1 if (first < 0) == (second < 0) else -1

  first_digits: list[int] = _digits_of(abs(first))
  second_digits: list[int] = _digits_of(abs(second))

  raw_product: list[int] = poly_multiply_round(first_digits, second_digits)
  magnitude: int = _carry(raw_product)
  return sign * magnitude

The FFT beyond powers of two

Cooley and Tukey's 1965 paper is one of the most-cited in all of computing, and the ideas around it stretch from signal processing to the fastest multiplication algorithm known.

Not just powers of two. The recursion above needs $n$ to be a power of two. Real FFT libraries handle any length by mixed-radix splitting (factor $n=n_1{n}_2$ and do $n_1$ transforms of size $n_2$ then $n_2$ of size $n_1$), and for a prime length they use Bluestein's algorithm, which re-expresses the size-$n$ DFT as a convolution and pads that to a convenient power of two.² So the $\Theta (n\log n)$ bound holds for every $n$, not only the padded powers of two — a detail that matters when you cannot afford to round the transform size up.

The FFT was not new. Cooley and Tukey rediscovered and popularized the algorithm, but Gauss had written down the same even/odd recursion around 1805 to interpolate asteroid orbits — before Fourier's own work was published.³ The 1965 paper's real impact was timing: it arrived just as digital computers made $\Theta (n\log n)$ spectral analysis suddenly practical, launching modern digital signal processing.

Toward optimal multiplication. FFT-based multiplication leads to the asymptotically fastest integer multiplication known. Schönhage–Strassen (1971) runs the transform over a ring of the form $\mathbb{Z}/(2^m+1)\mathbb{Z}$ — no floating-point, no NTT-prime constraints — to multiply $n$-bit integers in $O(n\log n\log \log n)$.⁴ The $\log \log n$ factor stood for nearly fifty years until Harvey and van der Hoeven (2021) reached the conjectured optimum $O(n\log n)$.⁵ These are galactic algorithms — the crossover point is enormous — but Schönhage–Strassen is genuinely used by GMP for numbers beyond a few tens of thousands of digits.

Beyond multiplication. The same transform accelerates string matching (match a pattern against a text by convolving indicator vectors, catching wildcards), the computation of the autocorrelation and power spectrum of a signal, large-scale polynomial interpolation, and the multiplication step inside the Fiduccia/Kitamasa recurrence solver from the matrix-exponentiation lesson. All of these reduce to convolution, which the FFT makes cheap.

Takeaways

• Multiplying polynomials (or big integers) is convolution, $\Theta (n^2)$ in the coefficient representation; but in the point-value representation the product is a $\Theta (n)$ pointwise multiply at shared sample points.
• The bottleneck is converting between the two forms. Sampling at the $n$-th roots of unity makes both conversions cheap because squaring halves the root set and the second half is the negation of the first.
• The DFT evaluates $A$ at all roots of unity; the Cooley–Tukey FFT splits on even/odd coefficient parity into two half-size DFTs combined by butterflies, giving $T(n)=2T(n/2)+\Theta (n)=\Theta (n\log n)$ by the Master Theorem.
• The inverse FFT is the same algorithm with conjugate roots and a $1/n$ scaling, because the DFT matrix's inverse is $\frac{1}{n}\overline{V}$ (a geometric-series argument).
• Three FFTs plus a pointwise multiply give $\Theta (n\log n)$ polynomial and big-integer multiplication; the NTT does the same exactly modulo a prime.