AlgoVerse

Fast Fourier transform (FFT)

Evaluate a polynomial at the n complex roots of unity by divide and conquer in O(n log n) (and back, with the inverse FFT); the basis of fast polynomial multiplication.

We want the polynomial A(x) = a_0 + a_1·x + … + a_{n−1}·x^{n−1} at n points. If the points are the n-th roots of unity ω^0, ω^1, …, ω^{n−1} (with ω = e^{2πi/n}), a trick works:

Read the full explanation ↓
Input:
coeffs 1×4 · inverse=0PDF export

A(x) = 1 + 2x + 3x² + 4x³ at the four 4th roots of unity (1, i, −1, −i): 10, −2−2i, −2, −2+2i.

  1. import cmath
  4. def transform(coeffs, inverse):
  5. a = coeffs[0]
  6. n = len(a)
  7. if n & (n - 1) != 0:
  8. return None
  9. y = fft(a, inverse)
  10. if inverse:
  11. y = [v / n for v in y]
  12. rows = [(k, show(v)) for k, v in enumerate(y)]
  13. return [show(v) for v in y]
  16. def fft(a, inverse):
  17. n = len(a)
  18. if n == 1:
  19. return [complex(a[0])]
  20. sign = -1 if inverse else 1
  21. w_n = cmath.exp(sign * 2j * cmath.pi / n) # the principal n-th root of unity
  22. even = fft(a[0::2], inverse)
  23. odd = fft(a[1::2], inverse)
  24. y = [0] * n
  25. w = 1
  26. for k in range(n // 2):
  27. t = w * odd[k]
  28. y[k] = even[k] + t
  29. y[k + n // 2] = even[k] - t
  30. w = w * w_n
  31. return y
  34. def show(z):
  35. # Rounded for reading, e.g. "2-1.414i".
  36. re, im = round(z.real, 3) + 0.0, round(z.imag, 3) + 0.0
  37. return f"{re:g}" if im == 0 else f"{re:g}{im:+g}i"

Invariant: A(ω^k) = A_even(ω^(2k)) + ω^k · A_odd(ω^(2k)), and the squares of the n-th roots of unity are the (n/2)-th roots; so each level halves the work.

Start

n—
k—

The table is still empty.

Recursion tree

Filled orange = running call · Orange outline = still on the stack · Grey = finished

Local variables

inverse0

Orange = changed in this step

coeffs list

[1, 2, 3, 4]

Call stack

  1. transform()line 4
Step 1 of 101
Black = butterfly · Orange = key stepbase: 0/4 · butterfly: 0/4 · done: 0/1
    Explanation

    The idea

    We want the polynomial A(x) = a_0 + a_1·x + … + a_{n−1}·x^{n−1} at n points. If the points are the n-th roots of unity ω^0, ω^1, …, ω^{n−1} (with ω = e^{2πi/n}), a trick works:

    A(x) = A_even(x²) + x · A_odd(x²)

    and the squares of the n-th roots are the (n/2)-th roots. So we solve two half-size problems and combine them with butterflies:

    y[k] = y⁰[k] + ω^k·y¹[k]      y[k + n/2] = y⁰[k] − ω^k·y¹[k]

    T(n) = 2T(n/2) + O(n) = O(n log n), against O(n²) to compute the n values one by one.

    The inverse FFT

    Going back from values to coefficients is the same FFT with ω⁻¹ instead of ω, then dividing by n (the parameter inverse = 1).

    What to watch in the trace

    • The call tree is binary (even and odd); the last table holds the rounded values.