THE UNIVERSITY OF NEW SOUTH WALES
5. THE FAST FOURIER TRANSFORM
Raveen de Silva,
office: K17 202
Course Admin: ,
School of Computer Science and Engineering UNSW Sydney
Term 3, 2021
Table of Contents
1. Revision
2. The Fast Fourier Transform
3. Puzzle
Our strategy to multiply polynomials fast
Given two polynomials of degree at most n, PA(x)=Anxn+…+A0 andPB(x)=Bnxn+…+B0,
we want to find PC(x) = PA(x)PB(x).
1. Evaluate PA(x) and PB(x) at 2n + 1 distinct points
x0,x1,…,x2n
2. Multiply them point by point using 2n + 1 multiplications of large numbers, to get 2n + 1 values of PC (x )
3. Reconstruct the coefficients of PC(x), i.e.
PC (x) = C2nx2n + C2n−1x2n−1 + . . . + C1x + C0.
Our strategy to multiply polynomials fast
PA(x) evaluate PA(x0) · · · PA(x2n) ××
PB(x) evaluate PB(x0) ··· PB(x2n)
PC (x) reconstruct PC (x0) · · · PC (x2n)
Our strategy to multiply polynomials fast
Previously, we chose x0, . . . , x2n to be the 2n + 1 integers which are the smallest by their absolute value, i.e.,
−n, −(n − 1), . . . , −1, 0, 1, . . . , n − 1, n. However, this choice requires us to compute values such as
PA(n)=A0 +A1n+…+An−1nn−1 +Annn
As the degree n of the polynomials PA(x) and PB(x) increases, the value of nn increases very fast and causes rapid increase of the computational complexity of the algorithm for polynomial multiplication which we used in the generalised Karatsuba algorithm.
Our strategy to multiply polynomials fast
Question
What values should we take for x0, . . . , x2n to avoid “explosion” of size when we evaluate xin while computing
P A ( x i ) = A 0 + A 1 x + . . . + A n x in ?
Answer
We would like |xin| = 1, but this can’t be achieved with only real numbers . . .
Complex numbers revisited
Let z = a+ib be a complex number, where a,b ∈ R. Then we can
define the
√ modulus |z|, a positive real number r such that r = a2 + b2
and
argument argz, an angle θ ∈ (−π,π] such that a = r cosθ
and b = r sin θ. Then we have
z = a + ib
= r (cos θ + i sin θ) = r eiθ.
Complex numbers revisited
b
r
z = a + ib = r eiθ
θ
a
Complex numbers revisited
To multiply two complex numbers, we multiply the moduli and add the arguments.
Theorem
Ifz=reiθ andw=seiφ,thenzw=(rs)ei(θ+φ).
Corollary
In particular, if z = r eiθ then zn = rn ei (nθ).
Complex numbers revisited
z2 r2
z
θ θ
Complex roots of unity
We are interested in the solutions of zm = 1, known as roots of unity of order m.
Solving rm ei (mθ) = 1, we see that
rm = 1, i.e. r = 1
and
mθ = 2π k i.e. θ = 2π k . m
Letωm =ei2π/m. Thenz=ωmk,i.e.,allrootsofunityof order m can be written as powers of ωm. We say that ωm is a primitive root of unity of order m.
Note that there are only m distinct values here, as ωm = ωm0 .
Complex roots of unity
1
ωm = ei 2π/m
2π/m
Complex roots of unity
ω92 ω93
ω91
ω94
ω95
ω96
ω90 ω98
ω97
Complex roots of unity
For ωm = ei 2π/m and for all k such that 0 ≤ k ≤ m − 1, ((ωm)k)m =(ωm)mk =((ωm)m)k =1k =1.
Thus, ωmk = (ωm)k is also a root of unity.
Theorem
ωmk is a primitive root of unity of order n if and only if gcd(m,k) = 1.
Complex roots of unity
The product of two roots of unity of order m is given by ωiωj =ωi+j,
mmm
which is itself a root of unity of the same order. So the set of all roots of unity of order m, i.e.,
{1,ω ,ω2,…,ωm−1}isclosedundermultiplication(andby mmm
extension, under taking powers).
Note that this is not true for addition, i.e., the sum of two roots of unity is NOT another root of unity!
Note
This is an example of a group.
Complex roots of unity
Cancellation Lemma
For all positive integers k, m and integers l, ωkl = ωl . km m
Proof
kl i 2π kl i 2πkl i 2πl i 2π l l ωkm=ekm =ekm =em =em =ωm.
In particular, (ωk )2 = ω2k = (ω2 )k = ωk , i.e. the squares 2m 2m 2m m
of the roots of unity of order 2m are just the roots of unity of order m.
Table of Contents
1. Revision
2. The Fast Fourier Transform
3. Puzzle
A new way to multiply polynomials fast
Given two polynomials of degree at most n, PA(x)=Anxn+…+A0 andPB(x)=Bnxn+…+B0,
we want to find PC(x) = PA(x)PB(x).
1. Evaluate PA(x) and PB(x) at m = 2n + 1 distinct points
ω0 ,ω1 ,…,ωm−1 mmm
2. Multiply them point by point using 2n + 1 multiplications of large numbers, to get 2n + 1 values of PC (x )
3. Reconstruct the coefficients of PC(x), i.e.
PC (x) = C2nx2n + C2n−1x2n−1 + . . . + C1x + C0.
Our strategy to multiply polynomials fast
P (x) evaluate P ω0 ··· P ωm−1 AAmAm
××
P (x) evaluate P ω0 ··· P ωm−1
BBmBm
P (x) P ω0 ··· P ωm−1 C reconstruct C m C m
The Discrete Fourier Transform
Definition
Let A = ⟨A0,A1,…,Am−1⟩ be a sequence of m real or complex
numbers, and let the corresponding polynomial be
PA(x) = m−1 Ajxj. j=0
LetAk =PA(ωmk)forall0≤k≤m−1,thevaluesofPA atthe roots of unity of order m.
The sequence of values A = ⟨A0, A1, . . . , Am−1⟩ is called the Discrete Fourier Transform (DFT) of A.
The Discrete Fourier Transform
Our polynomials PA(x) and PB(x) have degree n, so the corresponding sequences of coefficients have only n + 1 terms.
However, we defined the DFT of a sequence as having the same length as the original sequence, and we must obtain the values at all m roots of unity of order m.
This is easily rectified by padding the sequences A and B with zeros at the end, corresponding to the terms
0xn+1 + 0xn+2 + . . . + 0x2n, since m = 2n + 1.
The DFT A of a sequence A can be computed VERY FAST using a divide-and-conquer algorithm called the Fast Fourier Transform.
Our strategy to multiply polynomials fast
We can now compute the DFTs of the two (0 padded) sequences:
DFT(⟨A0,A1,…,An,0,…,0⟩) = ⟨A0,A1,…,Am−1⟩
n
and
DFT(⟨B0,B1,…,Bn,0,…,0⟩) = ⟨B0,B1,…,Bm−1⟩
n
For each k we multiply the corresponding values Ak = PA(ωmk ) and Bk = PB(ωmk ), thus obtaining
Ck = AkBk = PA(ωmk )PB(ωmk ) = PC(ωmk )
We then use the inverse transformation for DFT, called IDFT,
to recover the coefficients ⟨C0, C1, . . . , Cm−1⟩ of the product
polynomial PC (x) from the sequence ⟨C0, C1, . . . , Cm−1⟩ of its
Our strategy to multiply polynomials fast
⟨A0,A1,…,An,0,…,0⟩ DFT
n
⟨B0,B1,…,Bn,0,…,0⟩ DFT
A0 ··· Am−1 ××
B0 ···Bm−1
C ··· C
0 m−1
n
⟨C0,C1,…,Cm−1⟩
IDFT
The Fast Fourier Transform (FFT)
Question
How long does it take to compute the DFT of a sequence A, i.e. find the values
PA(ωmk)=A0+A1ωmk +A2(ωmk)2+…+Am−1(ωmk)m−1 for all k such that 0 ≤ k < m?
Answer
Even if we were to precompute all powers of ωm (recall, there are only m distinct values), there are m multiplications required for each of m values PA(ωmk ).
However, we can use divide-and-conquer to find all the required values in O(n log n) time!
The Fast Fourier Transform (FFT)
We can assume that m is a power of 2 - otherwise we can pad PA(x) with zero coefficients until its number of coefficients becomes equal to the nearest power of 2.
Exercise
Show that for every m which is not a power of two the smallest power of 2 larger or equal to m is smaller than 2m.
Hint
Consider m in binary. How many bits does the nearest power of two have?
The Fast Fourier Transform (FFT)
Problem
Given a sequence A = ⟨A0, A1, . . . , Am−1⟩, compute its DFT, i.e. findthevaluesPA(x)atx=ωmk forallksuchthat0≤k