Understanding Fast Fourier Transform from scratch — to solve Polynomial Multiplication.

Aiswarya Prakasan
Feb 23, 2017 · 13 min read
E.g., multiply (x^2 + 2x + 3)(2x^2 + 5) = 

2x^2 + 0x + 5
1x^2 + 2x + 3
-------------
6 0 15
4 0 10
2 0 5
----------------------------
2x^4 + 4x^3 + 11x^2 +10x+15
ci = a0*bi + a1*bi-1 + ... + ai*b0.
ζ(f⊗g) = ζ(f).ζ(g)
A set of n pairs {(x0, y0),(x1, y1),...,(xn, yn)} such that 
• for all i != j, xi != xj. ie, the points are unique.
• for every k, yk = A(xk);
If 𝐴 and 𝐵 are of degree-bound 'n', then 𝐶 is of degree-bound '2n'.
Need to start with “extended” point-value forms for 𝐴 and 𝐵 consisting of
2𝑛 point-value pairs each. –

𝐴: (𝑥0, 𝑦0) , (𝑥1, 𝑦1) , … ,( 𝑥2𝑛−1, 𝑦2𝑛−1)
𝐵: (𝑥0, 𝑦0') , (𝑥1, 𝑦1 ′) , … , (𝑥2𝑛−1, 𝑦2𝑛−1 ′)
𝐶: (𝑥0, 𝑦0𝑦0 ′) ,( 𝑥1, 𝑦1𝑦1 ′) , … ,( 𝑥2𝑛−1, 𝑦2𝑛−1𝑦2𝑛−1 ′)
Let m = 2n-1. [so degree of C is less than m]
1. Pick m points x_0, x_1, ..., x_{m-1} chosen cleverly.
2. Evaluate A at each of the points: A(x_0),..., A(x_{m-1}).
3. Same for B.
4. Now compute C(x_0),..., C(x_{m-1}), where C is A(x)*B(x)
5. Interpolate to get the coefficients of C.
  polar coordinates : re^iθ
cartesian coordinates :(rcosθ, rsinθ)
As we will see, the fast Fourier transform algorithm cleverly makes 
use of the following properties about ωn:
𝜔𝑛^n =1
𝜔𝑛^n+𝑘 = 𝜔𝑛^𝑘
𝜔𝑛^n/2 = -1
𝜔𝑛^(n/2+𝑘) = -𝜔𝑛^
𝑘
The 2nd half of points are the negative of the 1st half   (***)
Fast Fourier Transform (FFT) The problem of evaluating 𝐴(𝑥) at 𝜔𝑛^0 , 𝜔𝑛^1 , … , 𝜔𝑛^𝑛−1 reduces to 
1. evaluating the degree-bound 𝑛/2 polynomials 𝐴even(𝑥) and 𝐴odd(𝑥) at the points (𝜔𝑛^0)^2 ,(𝜔𝑛^1)^2 , … , (𝜔𝑛^𝑛−1)^2.
2. combining the results by 𝐴(𝑥) = 𝐴even(𝑥2) + 𝑥𝐴odd(𝑥2).

Why bother?
– The list (𝜔𝑛^0)^2 ,(𝜔𝑛^1)^2 , … , (𝜔𝑛^𝑛−1)^2 does not contain
𝑛 distinct values, but 𝑛/2 complex 𝑛/2-th roots of unity.
– Polynomials 𝐴even and 𝐴odd are recursively evaluated at the 𝑛/2 complex 𝑛/2-th roots of unity.
– Subproblems have exactly the same form as the original problem, but are half the size
Algorithm
Here is the general algorithm in pseudo-C:

Let A be array of length m, w be primitive mth root of unity.
Goal: produce DFT F(A): evaluation of A at 1, w, w^2,...,w^{m-1}.
FFT(A, m, w)
{
if (m==1) return vector (a_0)
else {
A_even = (a_0, a_2, ..., a_{m-2})
A_odd = (a_1, a_3, ..., a_{m-1})
F_even = FFT(A_even, m/2, w^2)//w^2 is a m/2-th root of unity
F_odd = FFT(A_odd, m/2, w^2)
F = new vector of length m
x = 1
for (j=0; j < m/2; ++j) {
F[j] = F_even[j] + x*F_odd[j]
F[j+m/2] = F_even[j] - x*F_odd[j]
x = x * w
}
return F
}
if ω = a+ib   => ω^-1 = a-ib
if ω = e^2Πi/k => ω^-1 = e^-2Πi/k
So, the final algorithm is:

Let F_A = FFT(A, m, ω) // time O(n log n)
Let F_B = FFT(B, m, ω) // time O(n log n)
For i=1 to m, let F_C[i] = F_A[i]*F_B[i] // time O(n)
Output C = 1/m * FFT(F_C, m, ω-1). // time O(n log n)

Aiswarya Prakasan

Written by

Software engineer, algorithm enthusiast.

Welcome to a place where words matter. On Medium, smart voices and original ideas take center stage - with no ads in sight. Watch
Follow all the topics you care about, and we’ll deliver the best stories for you to your homepage and inbox. Explore
Get unlimited access to the best stories on Medium — and support writers while you’re at it. Just $5/month. Upgrade