Fast Fourier transform

The fast Fourier transform (FFT) computes the discrete Fourier transform of $n$ values in $O (n lo g n)$ arithmetic operations, where the definition takes $O (n^{2})$ . The discrete Fourier transform evaluates a polynomial at the $n$ complex $n$ -th roots of unity, and the Cooley–Tukey algorithm does it by splitting the coefficients into those at even and odd positions, transforming each half recursively and combining the two with $n /2$ butterfly operations^[1]. The transform turns convolution into pointwise multiplication, so the FFT multiplies polynomials and large integers in close to linear time^[3]. Over a finite field with the right roots of unity the same algorithm is exact, the number-theoretic transform^[6]. The idea goes back to Gauss; Cooley and Tukey's 1965 paper made it the standard algorithm of signal processing^[2].

§ 01

The discrete Fourier transform

Let $ω_{n} = e^{- 2 π i / n}$ , a primitive $n$ -th root of unity. The discrete Fourier transform of $a_{0}, \dots, a_{n - 1}$ is

A_{k} = j = 0 \sum n - 1 a_{j} ω_{n}^{j k}, k = 0, \dots, n - 1

which is the value of the polynomial $a (x) = \sum_{j} a_{j} x^{j}$ at $x = ω_{n}^{k}$ . The inverse transform uses $ω_{n}^{- 1}$ and divides by $n$ :

a_{j} = \frac{1}{n} k = 0 \sum n - 1 A_{k} ω_{n}^{- j k}

because $\sum_{k} ω_{n}^{k (j - j^{'})}$ is $n$ when $j = j^{'}$ and 0 otherwise. Computed from the definition, the transform takes $n^{2}$ multiplications.

§ 02

The Cooley–Tukey algorithm

For $n$ a power of two, split the polynomial by the parity of the exponents, $a (x) = e (x^{2}) + x o (x^{2})$ , where $e$ and $o$ have the even- and odd-indexed coefficients. The squares of the $n$ -th roots of unity are the $n /2$ -th roots of unity, each occurring twice, since $ω_{n}^{2} = ω_{n /2}$ and $ω_{n}^{n /2} = - 1$ . So with $E$ and $O$ the transforms of length $n /2$ of $e$ and $o$ :

A_{k} = E_{k} + ω_{n}^{k} O_{k} A_{k + n /2} = E_{k} - ω_{n}^{k} O_{k} 0 \leq k < n /2

Each pair of outputs costs one multiplication by $ω_{n}^{k}$ , the twiddle factor, and an addition and a subtraction: a butterfly. The cost satisfies $T (n) = 2 T (n /2) + O (n)$ , so $T (n) = O (n lo g n)$ , with exactly $\frac{n}{2} lo g_{2} n$ twiddle multiplications^[3].

The butterflies of an 8-point FFT, computed in place. The inputs are in bit-reversed order; each stage combines transforms of size 2, 4 and 8, and each crossing pair of lines is one butterfly. Red lines carry the twiddle factor with a minus sign.

The transform from the definition, the recursive FFT, the inverse, and polynomial multiplication.

(* The discrete Fourier transform, directly in O(n^2) and by the
   recursive radix-2 Cooley-Tukey FFT in O(n log n). *)
open Complex
 
let omega n k = polar 1. (-2. *. Float.pi *. float k /. float n)  (* e^(-2 pi i k / n) *)
 
let mults = ref 0
let ( *! ) a b = incr mults; mul a b
 
let dft a =
  let n = Array.length a in
  Array.init n (fun k -> let s = ref zero in Array.iteri (fun j x -> s := add !s (x *! omega n (j * k mod n))) a; !s)
 
(* Split into even and odd positions, transform each half, and combine
   with the butterflies X_k = E_k + w^k O_k, X_{k+n/2} = E_k - w^k O_k. *)
let rec fft a =
  let n = Array.length a in
  if n = 1 then [| a.(0) |]
  else begin
    let e = fft (Array.init (n / 2) (fun i -> a.(2 * i))) and o = fft (Array.init (n / 2) (fun i -> a.((2 * i) + 1))) in
    let x = Array.make n zero in
    for k = 0 to (n / 2) - 1 do
      let t = omega n k *! o.(k) in
      x.(k) <- add e.(k) t;
      x.(k + (n / 2)) <- sub e.(k) t
    done;
    x
  end
 
(* The inverse transform is the forward one on conjugates, scaled by 1/n. *)
let ifft x = let n = float (Array.length x) in Array.map (fun z -> div (conj z) { re = n; im = 0. }) (fft (Array.map conj x))
 
(* Multiplying polynomials: transform, multiply pointwise, transform back. *)
let poly_mul a b =
  let len = Array.length a + Array.length b - 1 in
  let n = let rec up n = if n >= len then n else up (2 * n) in up 1 in
  let pad p = Array.init n (fun i -> if i < Array.length p then { re = float p.(i); im = 0. } else zero) in
  let c = ifft (Array.map2 mul (fft (pad a)) (fft (pad b))) in
  Array.init len (fun i -> Float.to_int (Float.round c.(i).re))
 
let schoolbook a b =
  let r = Array.make (Array.length a + Array.length b - 1) 0 in
  Array.iteri (fun i x -> Array.iteri (fun j y -> r.(i + j) <- r.(i + j) + (x * y)) b) a; r

Running it.

     n    mults (DFT)    mults (FFT)   max |diff|
    16            256             32     2.33e-15
    64           4096            192     6.28e-15
   256          65536           1024     4.01e-14
  1024        1048576           5120     1.15e-13
ifft (fft a) = a to within 5.00e-16 for n = 4096
product of two degree-1999 polynomials with coefficients below 1000 equals schoolbook: true

The two transforms agree to within rounding error, and the number of multiplications drops from $n^{2}$ to $\frac{n}{2} lo g_{2} n$ : about two hundred times fewer at $n = 1024$ . Unrolling the recursion gives the in-place form in the figure, which first permutes the input into bit-reversed order and then runs $lo g_{2} n$ stages of butterflies over the array^[3].

Other factorizations

The splitting works for any factorization $n = n_{1} n_{2}$ : the transform of length $n$ becomes $n_{2}$ transforms of length $n_{1}$ , twiddle multiplications, and $n_{1}$ transforms of length $n_{2}$ . Radix-4 and split-radix variants reduce the operation count further, and when $n_{1}$ and $n_{2}$ are coprime, Good's prime-factor algorithm needs no twiddle factors at all^[10]. Libraries such as FFTW choose among these at run time by measuring which is fastest on the machine^[7].

§ 03

Convolution and multiplication

The product of two polynomials has as coefficients the convolution of their coefficient sequences, $c_{k} = \sum_{i + j = k} a_{i} b_{j}$ . Evaluating at the roots of unity turns the product into a pointwise one, the convolution theorem^[3]:

DFT (a * b) = DFT (a) \cdot DFT (b) ⟹ a * b = DFT^{- 1} (DFT (a) \cdot DFT (b))

with both sequences padded with zeros to a length $n$ at least the length of the product, so the cyclic convolution computed by the transform does not wrap around. Two polynomials of degree below $n /2$ are multiplied in $O (n lo g n)$ operations, against $O (n^{2})$ by the schoolbook method and $O (n^{1.585})$ by Karatsuba. With floating-point arithmetic the result has to be rounded to integers, which is exact only while the rounding errors stay below one half; for large inputs or large coefficients this limits the method^[11].

Exact transforms

The transform only needs a primitive $n$ -th root of unity and the inverse of $n$ , so it works in any ring that has them. Modulo a prime $p$ with $n ∣ p - 1$ , the multiplicative group is cyclic of order $p - 1$ and contains a primitive $n$ -th root of unity^[6]. The prime $998244353 = 119 \cdot 2^{23} + 1$ has roots of unity of every order $2^{k}$ up to $2^{23}$ , and the transform over it, the number-theoretic transform, computes convolutions exactly as long as the true coefficients are below $p$ :

An iterative, in-place number-theoretic transform, and multiplication of polynomials and decimal integers with it.

(* The number-theoretic transform: the same algorithm over Z/pZ with
   p = 998244353 = 119 * 2^23 + 1, where 3 generates the multiplicative
   group, so there are n-th roots of unity for every n = 2^k up to 2^23.
   Arithmetic is exact. This version is iterative and in place. *)
let p = 998_244_353
let rec pw b e = if e = 0 then 1 else let h = pw (b * b mod p) (e / 2) in if e land 1 = 1 then h * b mod p else h
 
let ntt a invert =
  let n = Array.length a in
  (* Reorder by bit-reversed index, so the butterflies can work in place. *)
  let j = ref 0 in
  for i = 1 to n - 1 do
    let bit = ref (n lsr 1) in
    while !j land !bit <> 0 do j := !j lxor !bit; bit := !bit lsr 1 done;
    j := !j lor !bit;
    if i < !j then (let t = a.(i) in a.(i) <- a.(!j); a.(!j) <- t)
  done;
  let len = ref 2 in
  while !len <= n do
    let w = pw 3 ((p - 1) / !len) in
    let w = if invert then pw w (p - 2) else w in
    let i = ref 0 in
    while !i < n do
      let wk = ref 1 in
      for k = 0 to (!len / 2) - 1 do
        let u = a.(!i + k) and v = a.(!i + k + (!len / 2)) * !wk mod p in
        a.(!i + k) <- (u + v) mod p;
        a.(!i + k + (!len / 2)) <- (u - v + p) mod p;
        wk := !wk * w mod p
      done;
      i := !i + !len
    done;
    len := 2 * !len
  done;
  if invert then (let ni = pw n (p - 2) in Array.iteri (fun i x -> a.(i) <- x * ni mod p) a)
 
let mul a b =
  let len = Array.length a + Array.length b - 1 in
  let n = let rec up n = if n >= len then n else up (2 * n) in up 1 in
  let fa = Array.init n (fun i -> if i < Array.length a then a.(i) else 0)
  and fb = Array.init n (fun i -> if i < Array.length b then b.(i) else 0) in
  ntt fa false; ntt fb false;
  let c = Array.map2 (fun x y -> x * y mod p) fa fb in
  ntt c true;
  Array.sub c 0 len
 
(* Big integers as arrays of decimal digits, least significant first:
   multiply the digit polynomials, then propagate carries. *)
let big_mul x y =
  let digits s = Array.init (String.length s) (fun i -> Char.code s.[String.length s - 1 - i] - 48) in
  let c = mul (digits x) (digits y) in
  let buf = Buffer.create (Array.length c + 1) and carry = ref 0 in
  let out = Array.map (fun d -> let t = d + !carry in carry := t / 10; t mod 10) c in
  let rest = ref [] in
  while !carry > 0 do rest := (!carry mod 10) :: !rest; carry := !carry / 10 done;
  List.iter (fun d -> Buffer.add_char buf (Char.chr (48 + d))) !rest;
  for i = Array.length out - 1 downto 0 do Buffer.add_char buf (Char.chr (48 + out.(i))) done;
  let s = Buffer.contents buf in
  let k = ref 0 in while !k < String.length s - 1 && s.[!k] = '0' do incr k done;
  String.sub s !k (String.length s - !k)

Running it; the integer product agrees with Python's.

degree-49999 product: coefficients 0, 1, 12345, 49999, 75000, 99998 checked against the convolution sum, all equal: true
123456789123456789123456789123456789123456789123456789123456
  * 987654321987654321987654321987654321987654321
  = 121932631356500531591068431825636332060204232172839501172838721791646821557078921322511021087943120853376

When the coefficients may be larger, the convolution is computed modulo several such primes and the results combined by the Chinese remainder theorem, as in the modular algorithms of computer algebra^[8].

§ 04

Integer multiplication

Writing integers as polynomials in a base $B$ , with their digits as coefficients, and evaluating the product at $B$ by propagating carries, reduces integer multiplication to polynomial multiplication. Schönhage and Strassen performed the transform over the ring $Z / (2^{k} + 1)$ , in which 2 is a root of unity and the twiddle multiplications are shifts, and multiplied $n$ -bit integers in $O (n lo g n lo g lo g n)$ bit operations^[4]. Harvey and van der Hoeven reached $O (n lo g n)$ in 2019, published in 2021, the bound Schönhage and Strassen had conjectured to be optimal^[5]. Arbitrary-precision libraries such as GMP switch from schoolbook to Karatsuba, Toom–Cook and finally FFT-based multiplication as the operands grow.

§ 05

Applications

In signal processing the transform computes the frequency spectrum of a sampled signal, and fast convolution implements digital filters and correlation. Image and audio compression use its real-valued relative, the discrete cosine transform. Spectral methods for differential equations differentiate in the frequency domain. In computer algebra, fast polynomial multiplication gives fast division through Newton iteration, fast multipoint evaluation and interpolation, and fast greatest common divisors^[8].

§ 06

History

Gauss used the method in 1805 to interpolate the orbits of asteroids, in work published after his death; Heideman, Johnson and Burrus traced this and several later rediscoveries^[2]. Good described the prime-factor algorithm in 1958^[10]. Cooley and Tukey published the radix-2 algorithm in 1965, when digital computers made it immediately useful^[1], and Gentleman and Sande analysed its variants in 1966^[9]. Pollard introduced transforms over finite fields in 1971^[6], the year Schönhage and Strassen applied the FFT to integer multiplication^[4].