Algorithms COMP3121/9101
THE UNIVERSITY OF
NEW SOUTH WALES
Algorithms
COMP3121/9101
Aleks Ignjatović
School of Computer Science and Engineering
University of New South Wales
5. THE FAST FOURIER TRANSFORM
(not examinable material)
COMP3121/9101 1 / 33
Our strategy to multiply polynomials fast:
Given two polynomials of degree at most n,
PA(x) = Anx
n
+ . . .+A0; PB(x) = Bnx
n
+ . . .+B0
1 convert them into value representation at 2n + 1 distinct points
x0, x1, . . . , x2n:
PA(x)↔ {(x0, PA(x0)), (x1, PA(x1)), . . . , (x2n, PA(x2n))}
PB(x)↔ {(x0, PB(x0)), (x1, PB(x1)), . . . , (x2n, PB(x2n))}
2 multiply them point by point using 2n + 1 multiplications:{
(x0, PA(x0)PB(x0)︸ ︷︷ ︸
PC(x0)
), (x1, PA(x1)PB(x1)︸ ︷︷ ︸
PC(x1)
), . . . , (x2n, PA(x2n)PB(x2n)︸ ︷︷ ︸
PC(x2n)
)
}
3 Convert such value representation of PC(x) to its coefficient form
PC(x) = C2nx
2n + C2n−1x
2n−1 + . . . + C1x + C0;
COMP3121/9101 2 / 33
Our strategy to multiply polynomials fast:
Given two polynomials of degree at most n,
PA(x) = Anx
n
+ . . .+A0; PB(x) = Bnx
n
+ . . .+B0
1 convert them into value representation at 2n + 1 distinct points
x0, x1, . . . , x2n:
PA(x)↔ {(x0, PA(x0)), (x1, PA(x1)), . . . , (x2n, PA(x2n))}
PB(x)↔ {(x0, PB(x0)), (x1, PB(x1)), . . . , (x2n, PB(x2n))}
2 multiply them point by point using 2n + 1 multiplications:{
(x0, PA(x0)PB(x0)︸ ︷︷ ︸
PC(x0)
), (x1, PA(x1)PB(x1)︸ ︷︷ ︸
PC(x1)
), . . . , (x2n, PA(x2n)PB(x2n)︸ ︷︷ ︸
PC(x2n)
)
}
3 Convert such value representation of PC(x) to its coefficient form
PC(x) = C2nx
2n + C2n−1x
2n−1 + . . . + C1x + C0;
COMP3121/9101 2 / 33
Our strategy to multiply polynomials fast:
Given two polynomials of degree at most n,
PA(x) = Anx
n
+ . . .+A0; PB(x) = Bnx
n
+ . . .+B0
1 convert them into value representation at 2n + 1 distinct points
x0, x1, . . . , x2n:
PA(x)↔ {(x0, PA(x0)), (x1, PA(x1)), . . . , (x2n, PA(x2n))}
PB(x)↔ {(x0, PB(x0)), (x1, PB(x1)), . . . , (x2n, PB(x2n))}
2 multiply them point by point using 2n + 1 multiplications:{
(x0, PA(x0)PB(x0)︸ ︷︷ ︸
PC(x0)
), (x1, PA(x1)PB(x1)︸ ︷︷ ︸
PC(x1)
), . . . , (x2n, PA(x2n)PB(x2n)︸ ︷︷ ︸
PC(x2n)
)
}
3 Convert such value representation of PC(x) to its coefficient form
PC(x) = C2nx
2n + C2n−1x
2n−1 + . . . + C1x + C0;
COMP3121/9101 2 / 33
Our strategy to multiply polynomials fast:
Given two polynomials of degree at most n,
PA(x) = Anx
n
+ . . .+A0; PB(x) = Bnx
n
+ . . .+B0
1 convert them into value representation at 2n + 1 distinct points
x0, x1, . . . , x2n:
PA(x)↔ {(x0, PA(x0)), (x1, PA(x1)), . . . , (x2n, PA(x2n))}
PB(x)↔ {(x0, PB(x0)), (x1, PB(x1)), . . . , (x2n, PB(x2n))}
2 multiply them point by point using 2n + 1 multiplications:{
(x0, PA(x0)PB(x0)︸ ︷︷ ︸
PC(x0)
), (x1, PA(x1)PB(x1)︸ ︷︷ ︸
PC(x1)
), . . . , (x2n, PA(x2n)PB(x2n)︸ ︷︷ ︸
PC(x2n)
)
}
3 Convert such value representation of PC(x) to its coefficient form
PC(x) = C2nx
2n + C2n−1x
2n−1 + . . . + C1x + C0;
COMP3121/9101 2 / 33
Our strategy to multiply polynomials fast:
Given two polynomials of degree at most n,
PA(x) = Anx
n
+ . . .+A0; PB(x) = Bnx
n
+ . . .+B0
1 convert them into value representation at 2n + 1 distinct points
x0, x1, . . . , x2n:
PA(x)↔ {(x0, PA(x0)), (x1, PA(x1)), . . . , (x2n, PA(x2n))}
PB(x)↔ {(x0, PB(x0)), (x1, PB(x1)), . . . , (x2n, PB(x2n))}
2 multiply them point by point using 2n + 1 multiplications:{
(x0, PA(x0)PB(x0)︸ ︷︷ ︸
PC(x0)
), (x1, PA(x1)PB(x1)︸ ︷︷ ︸
PC(x1)
), . . . , (x2n, PA(x2n)PB(x2n)︸ ︷︷ ︸
PC(x2n)
)
}
3 Convert such value representation of PC(x) to its coefficient form
PC(x) = C2nx
2n + C2n−1x
2n−1 + . . . + C1x + C0;
COMP3121/9101 2 / 33
Our strategy to multiply polynomials fast:
So, we need 2n+ 1 values of PA(xi) and PB(xi), 0 ≤ i ≤ 2n.
If we use 2n+ 1 integers which are the smallest by their absolute value, i.e.,
−n,−(n− 1), . . . ,−1, 0, 1, . . . , n− 1, n
among the values which we need to compute is
PA(n) = A0 +A1n+ . . .+An−1n
n−1 +Ann
n
We saw that the trouble is that, as the degree n of the polynomials PA(x) and
PB(x) increases, the value of n
n 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.
Key Question: What values should we take for x0, . . . , x2n to avoid
“explosion” of size when we evaluate xni while computing
PA(xi) = A0 +A1x+ . . .+Anx
n
i ?
COMP3121/9101 3 / 33
Our strategy to multiply polynomials fast:
So, we need 2n+ 1 values of PA(xi) and PB(xi), 0 ≤ i ≤ 2n.
If we use 2n+ 1 integers which are the smallest by their absolute value, i.e.,
−n,−(n− 1), . . . ,−1, 0, 1, . . . , n− 1, n
among the values which we need to compute is
PA(n) = A0 +A1n+ . . .+An−1n
n−1 +Ann
n
We saw that the trouble is that, as the degree n of the polynomials PA(x) and
PB(x) increases, the value of n
n 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.
Key Question: What values should we take for x0, . . . , x2n to avoid
“explosion” of size when we evaluate xni while computing
PA(xi) = A0 +A1x+ . . .+Anx
n
i ?
COMP3121/9101 3 / 33
Our strategy to multiply polynomials fast:
So, we need 2n+ 1 values of PA(xi) and PB(xi), 0 ≤ i ≤ 2n.
If we use 2n+ 1 integers which are the smallest by their absolute value, i.e.,
−n,−(n− 1), . . . ,−1, 0, 1, . . . , n− 1, n
among the values which we need to compute is
PA(n) = A0 +A1n+ . . .+An−1n
n−1 +Ann
n
We saw that the trouble is that, as the degree n of the polynomials PA(x) and
PB(x) increases, the value of n
n 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.
Key Question: What values should we take for x0, . . . , x2n to avoid
“explosion” of size when we evaluate xni while computing
PA(xi) = A0 +A1x+ . . .+Anx
n
i ?
COMP3121/9101 3 / 33
Our strategy to multiply polynomials fast:
So, we need 2n+ 1 values of PA(xi) and PB(xi), 0 ≤ i ≤ 2n.
If we use 2n+ 1 integers which are the smallest by their absolute value, i.e.,
−n,−(n− 1), . . . ,−1, 0, 1, . . . , n− 1, n
among the values which we need to compute is
PA(n) = A0 +A1n+ . . .+An−1n
n−1 +Ann
n
We saw that the trouble is that, as the degree n of the polynomials PA(x) and
PB(x) increases, the value of n
n 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.
Key Question: What values should we take for x0, . . . , x2n to avoid
“explosion” of size when we evaluate xni while computing
PA(xi) = A0 +A1x+ . . .+Anx
n
i ?
COMP3121/9101 3 / 33
Complex numbers revisited
Complex numbers z = a+ ib can be represented using their modulus |z| =
√
a2 + b2
and their argument, arg(z), which is an angle taking values in (−π, π] and satisfying:
z = |z|ei arg(z) = |z|(cos arg(z) + i sin arg(z)),
see figure below.
!
b!
a!
|z|!
arg!(z)!
z!=!a!+!i!b!
COMP3121/9101 4 / 33
Complex numbers revisited
Recall that
z
n
=
(
|z|ei arg(z)
)n
= |z|nei n arg(z) = |z|n(cos(n arg(z)) + i sin(n arg(z))),
see the figure.
!
z!
z2!
|z|2!
2!arg!(z)!
COMP3121/9101 5 / 33
Complex roots of unity
Roots of unity of order n are complex numbers which satisfy zn = 1.
If zn = |z|n(cos(n arg(z)) + i sin(n arg(z))) = 1 then
|z| = 1 and n arg(z) is an integer multiple of 2π;
Thus, n arg(z) = 2π k, i.e., arg(z) =
2π k
n
We denote ωn = e
i 2π/n; such ωn is called a primitive root of unity of order n.
!
2!/!!
!! = !!!!!!/!!
|!!| = 1!
1!
ω16
(ω16)2
(ω16)15
(ω16)3
(ω16)4 (ω16)5
1=(ω16)0 (ω16)8
Roots of unity of order 16
A root of unity ω of order n is primitive if all other roots of unity of the same
order can be obtained as its powers ωk.
COMP3121/9101 6 / 33
Complex roots of unity
Roots of unity of order n are complex numbers which satisfy zn = 1.
If zn = |z|n(cos(n arg(z)) + i sin(n arg(z))) = 1 then
|z| = 1 and n arg(z) is an integer multiple of 2π;
Thus, n arg(z) = 2π k, i.e., arg(z) =
2π k
n
We denote ωn = e
i 2π/n; such ωn is called a primitive root of unity of order n.
!
2!/!!
!! = !!!!!!/!!
|!!| = 1!
1!
ω16
(ω16)2
(ω16)15
(ω16)3
(ω16)4 (ω16)5
1=(ω16)0 (ω16)8
Roots of unity of order 16
A root of unity ω of order n is primitive if all other roots of unity of the same
order can be obtained as its powers ωk.
COMP3121/9101 6 / 33
Complex roots of unity
Roots of unity of order n are complex numbers which satisfy zn = 1.
If zn = |z|n(cos(n arg(z)) + i sin(n arg(z))) = 1 then
|z| = 1 and n arg(z) is an integer multiple of 2π;
Thus, n arg(z) = 2π k, i.e., arg(z) =
2π k
n
We denote ωn = e
i 2π/n; such ωn is called a primitive root of unity of order n.
!
2!/!!
!! = !!!!!!/!!
|!!| = 1!
1!
ω16
(ω16)2
(ω16)15
(ω16)3
(ω16)4 (ω16)5
1=(ω16)0 (ω16)8
Roots of unity of order 16
A root of unity ω of order n is primitive if all other roots of unity of the same
order can be obtained as its powers ωk.
COMP3121/9101 6 / 33
Complex roots of unity
Roots of unity of order n are complex numbers which satisfy zn = 1.
If zn = |z|n(cos(n arg(z)) + i sin(n arg(z))) = 1 then
|z| = 1 and n arg(z) is an integer multiple of 2π;
Thus, n arg(z) = 2π k, i.e., arg(z) =
2π k
n
We denote ωn = e
i 2π/n; such ωn is called a primitive root of unity of order n.
!
2!/!!
!! = !!!!!!/!!
|!!| = 1!
1!
ω16
(ω16)2
(ω16)15
(ω16)3
(ω16)4 (ω16)5
1=(ω16)0 (ω16)8
Roots of unity of order 16
A root of unity ω of order n is primitive if all other roots of unity of the same
order can be obtained as its powers ωk.
COMP3121/9101 6 / 33
Complex roots of unity
Roots of unity of order n are complex numbers which satisfy zn = 1.
If zn = |z|n(cos(n arg(z)) + i sin(n arg(z))) = 1 then
|z| = 1 and n arg(z) is an integer multiple of 2π;
Thus, n arg(z) = 2π k, i.e., arg(z) =
2π k
n
We denote ωn = e
i 2π/n; such ωn is called a primitive root of unity of order n.
!
2!/!!
!! = !!!!!!/!!
|!!| = 1!
1!
ω16
(ω16)2
(ω16)15
(ω16)3
(ω16)4 (ω16)5
1=(ω16)0 (ω16)8
Roots of unity of order 16
A root of unity ω of order n is primitive if all other roots of unity of the same
order can be obtained as its powers ωk.
COMP3121/9101 6 / 33
Complex roots of unity
For ωn = e
i 2π/n and for all k such that 0 ≤ k ≤ n− 1,
((ωn)
k
)
n
= (ωn)
nk
= ((ωn)
n
)
k
= 1
k
= 1.
Thus, ωkn = (ωn)
k is also a root of unity (and it can be shown that it is
primitive just in case k is relatively prime with n).
Since ωkn are roots of unity for k = 0, 1, . . . , n− 1 and there are at most n
distinct roots of unity of order n (i.e., solutions to the equation xn − 1 = 0) we
conclude that every root of unity of order n must be of the form ωkn.
For the product of any two roots of unity ωkn and ω
m
n of the same order we
have ωknω
m
n = ω
k+m
n .
If k +m ≥ n then k +m = n+ l for l = (k +m) mod n and we have
ωknω
m
n = ω
k+m
n = ω
n+l
n = ω
n
nω
l
n = 1 · ωln = ωln where 0 ≤ l < n.
Thus, the product of any two roots of unity of the same order is just another
root of unity of the same order.
COMP3121/9101 7 / 33
Complex roots of unity
For ωn = e
i 2π/n and for all k such that 0 ≤ k ≤ n− 1,
((ωn)
k
)
n
= (ωn)
nk
= ((ωn)
n
)
k
= 1
k
= 1.
Thus, ωkn = (ωn)
k is also a root of unity (and it can be shown that it is
primitive just in case k is relatively prime with n).
Since ωkn are roots of unity for k = 0, 1, . . . , n− 1 and there are at most n
distinct roots of unity of order n (i.e., solutions to the equation xn − 1 = 0) we
conclude that every root of unity of order n must be of the form ωkn.
For the product of any two roots of unity ωkn and ω
m
n of the same order we
have ωknω
m
n = ω
k+m
n .
If k +m ≥ n then k +m = n+ l for l = (k +m) mod n and we have
ωknω
m
n = ω
k+m
n = ω
n+l
n = ω
n
nω
l
n = 1 · ωln = ωln where 0 ≤ l < n.
Thus, the product of any two roots of unity of the same order is just another
root of unity of the same order.
COMP3121/9101 7 / 33
Complex roots of unity
For ωn = e
i 2π/n and for all k such that 0 ≤ k ≤ n− 1,
((ωn)
k
)
n
= (ωn)
nk
= ((ωn)
n
)
k
= 1
k
= 1.
Thus, ωkn = (ωn)
k is also a root of unity (and it can be shown that it is
primitive just in case k is relatively prime with n).
Since ωkn are roots of unity for k = 0, 1, . . . , n− 1 and there are at most n
distinct roots of unity of order n (i.e., solutions to the equation xn − 1 = 0) we
conclude that every root of unity of order n must be of the form ωkn.
For the product of any two roots of unity ωkn and ω
m
n of the same order we
have ωknω
m
n = ω
k+m
n .
If k +m ≥ n then k +m = n+ l for l = (k +m) mod n and we have
ωknω
m
n = ω
k+m
n = ω
n+l
n = ω
n
nω
l
n = 1 · ωln = ωln where 0 ≤ l < n.
Thus, the product of any two roots of unity of the same order is just another
root of unity of the same order.
COMP3121/9101 7 / 33
Complex roots of unity
For ωn = e
i 2π/n and for all k such that 0 ≤ k ≤ n− 1,
((ωn)
k
)
n
= (ωn)
nk
= ((ωn)
n
)
k
= 1
k
= 1.
Thus, ωkn = (ωn)
k is also a root of unity (and it can be shown that it is
primitive just in case k is relatively prime with n).
Since ωkn are roots of unity for k = 0, 1, . . . , n− 1 and there are at most n
distinct roots of unity of order n (i.e., solutions to the equation xn − 1 = 0) we
conclude that every root of unity of order n must be of the form ωkn.
For the product of any two roots of unity ωkn and ω
m
n of the same order we
have ωknω
m
n = ω
k+m
n .
If k +m ≥ n then k +m = n+ l for l = (k +m) mod n and we have
ωknω
m
n = ω
k+m
n = ω
n+l
n = ω
n
nω
l
n = 1 · ωln = ωln where 0 ≤ l < n.
Thus, the product of any two roots of unity of the same order is just another
root of unity of the same order.
COMP3121/9101 7 / 33
Complex roots of unity
For ωn = e
i 2π/n and for all k such that 0 ≤ k ≤ n− 1,
((ωn)
k
)
n
= (ωn)
nk
= ((ωn)
n
)
k
= 1
k
= 1.
Thus, ωkn = (ωn)
k is also a root of unity (and it can be shown that it is
primitive just in case k is relatively prime with n).
Since ωkn are roots of unity for k = 0, 1, . . . , n− 1 and there are at most n
distinct roots of unity of order n (i.e., solutions to the equation xn − 1 = 0) we
conclude that every root of unity of order n must be of the form ωkn.
For the product of any two roots of unity ωkn and ω
m
n of the same order we
have ωknω
m
n = ω
k+m
n .
If k +m ≥ n then k +m = n+ l for l = (k +m) mod n and we have
ωknω
m
n = ω
k+m
n = ω
n+l
n = ω
n
nω
l
n = 1 · ωln = ωln where 0 ≤ l < n.
Thus, the product of any two roots of unity of the same order is just another
root of unity of the same order.
COMP3121/9101 7 / 33
Complex roots of unity
For ωn = e
i 2π/n and for all k such that 0 ≤ k ≤ n− 1,
((ωn)
k
)
n
= (ωn)
nk
= ((ωn)
n
)
k
= 1
k
= 1.
Thus, ωkn = (ωn)
k is also a root of unity (and it can be shown that it is
primitive just in case k is relatively prime with n).
Since ωkn are roots of unity for k = 0, 1, . . . , n− 1 and there are at most n
distinct roots of unity of order n (i.e., solutions to the equation xn − 1 = 0) we
conclude that every root of unity of order n must be of the form ωkn.
For the product of any two roots of unity ωkn and ω
m
n of the same order we
have ωknω
m
n = ω
k+m
n .
If k +m ≥ n then k +m = n+ l for l = (k +m) mod n and we have
ωknω
m
n = ω
k+m
n = ω
n+l
n = ω
n
nω
l
n = 1 · ωln = ωln where 0 ≤ l < n.
Thus, the product of any two roots of unity of the same order is just another
root of unity of the same order.
COMP3121/9101 7 / 33
Complex roots of unity
So in the set of all roots of unity of order n, i.e., {1, ωn, ω2n, . . . , ωn−1n } we can
multiply any two elements or raise an element to any power without going out
of this set.
Note that this is not true for addition, i.e., the sum of two roots of unity is
NOT another root of unity!
A most important property of the roots of unity is:
The Cancelation Lemma: ωkmkn = ω
m
n for all integers k,m, n.
Proof:
ω
km
kn = (ωkn)
km
= (e
i 2π
kn )
km
= e
i 2πkm
kn = e
i 2πm
n = (e
i 2π
n )
m
= ω
m
n
Thus, in particular, (ωk2n)
2 = ω2k2n = (ω
2
2n)
k = ωkn.
So the squares of the roots of unity of order 2n are just the roots of unity of
order n.
COMP3121/9101 8 / 33
Complex roots of unity
So in the set of all roots of unity of order n, i.e., {1, ωn, ω2n, . . . , ωn−1n } we can
multiply any two elements or raise an element to any power without going out
of this set.
Note that this is not true for addition, i.e., the sum of two roots of unity is
NOT another root of unity!
A most important property of the roots of unity is:
The Cancelation Lemma: ωkmkn = ω
m
n for all integers k,m, n.
Proof:
ω
km
kn = (ωkn)
km
= (e
i 2π
kn )
km
= e
i 2πkm
kn = e
i 2πm
n = (e
i 2π
n )
m
= ω
m
n
Thus, in particular, (ωk2n)
2 = ω2k2n = (ω
2
2n)
k = ωkn.
So the squares of the roots of unity of order 2n are just the roots of unity of
order n.
COMP3121/9101 8 / 33
Complex roots of unity
So in the set of all roots of unity of order n, i.e., {1, ωn, ω2n, . . . , ωn−1n } we can
multiply any two elements or raise an element to any power without going out
of this set.
Note that this is not true for addition, i.e., the sum of two roots of unity is
NOT another root of unity!
A most important property of the roots of unity is:
The Cancelation Lemma: ωkmkn = ω
m
n for all integers k,m, n.
Proof:
ω
km
kn = (ωkn)
km
= (e
i 2π
kn )
km
= e
i 2πkm
kn = e
i 2πm
n = (e
i 2π
n )
m
= ω
m
n
Thus, in particular, (ωk2n)
2 = ω2k2n = (ω
2
2n)
k = ωkn.
So the squares of the roots of unity of order 2n are just the roots of unity of
order n.
COMP3121/9101 8 / 33
Complex roots of unity
So in the set of all roots of unity of order n, i.e., {1, ωn, ω2n, . . . , ωn−1n } we can
multiply any two elements or raise an element to any power without going out
of this set.
Note that this is not true for addition, i.e., the sum of two roots of unity is
NOT another root of unity!
A most important property of the roots of unity is:
The Cancelation Lemma: ωkmkn = ω
m
n for all integers k,m, n.
Proof:
ω
km
kn = (ωkn)
km
= (e
i 2π
kn )
km
= e
i 2πkm
kn = e
i 2πm
n = (e
i 2π
n )
m
= ω
m
n
Thus, in particular, (ωk2n)
2 = ω2k2n = (ω
2
2n)
k = ωkn.
So the squares of the roots of unity of order 2n are just the roots of unity of
order n.
COMP3121/9101 8 / 33
Complex roots of unity
So in the set of all roots of unity of order n, i.e., {1, ωn, ω2n, . . . , ωn−1n } we can
multiply any two elements or raise an element to any power without going out
of this set.
Note that this is not true for addition, i.e., the sum of two roots of unity is
NOT another root of unity!
A most important property of the roots of unity is:
The Cancelation Lemma: ωkmkn = ω
m
n for all integers k,m, n.
Proof:
ω
km
kn = (ωkn)
km
= (e
i 2π
kn )
km
= e
i 2πkm
kn = e
i 2πm
n = (e
i 2π
n )
m
= ω
m
n
Thus, in particular, (ωk2n)
2 = ω2k2n = (ω
2
2n)
k = ωkn.
So the squares of the roots of unity of order 2n are just the roots of unity of
order n.
COMP3121/9101 8 / 33
Complex roots of unity
So in the set of all roots of unity of order n, i.e., {1, ωn, ω2n, . . . , ωn−1n } we can
multiply any two elements or raise an element to any power without going out
of this set.
Note that this is not true for addition, i.e., the sum of two roots of unity is
NOT another root of unity!
A most important property of the roots of unity is:
The Cancelation Lemma: ωkmkn = ω
m
n for all integers k,m, n.
Proof:
ω
km
kn = (ωkn)
km
= (e
i 2π
kn )
km
= e
i 2πkm
kn = e
i 2πm
n = (e
i 2π
n )
m
= ω
m
n
Thus, in particular, (ωk2n)
2 = ω2k2n = (ω
2
2n)
k = ωkn.
So the squares of the roots of unity of order 2n are just the roots of unity of
order n.
COMP3121/9101 8 / 33
The Discrete Fourier Transform
Let A = 〈A0, A1, . . . , An−1〉 be a sequence of n real or complex numbers.
We can form the corresponding polynomial PA(x) =
∑n−1
j=0
Ajx
j .
We can evaluate it at all complex roots of unity of order n, i.e., we compute
PA(ω
k
n) for all 0 ≤ k ≤ n− 1.
The sequence of values 〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉, is called the
Discrete Fourier Transform (DFT) of the sequence
A = 〈A0, A1, . . . , An−1〉.
The value PA(ω
k
n) is usually denoted by Âk and the sequence of values
〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉 is usually denoted by
 = 〈Â0, Â1, . . . , Ân−1〉.
The DFT Â of a sequence A can be computed VERY FAST using a
divide-and-conquer algorithm called the Fast Fourier Transform.
COMP3121/9101 9 / 33
The Discrete Fourier Transform
Let A = 〈A0, A1, . . . , An−1〉 be a sequence of n real or complex numbers.
We can form the corresponding polynomial PA(x) =
∑n−1
j=0
Ajx
j .
We can evaluate it at all complex roots of unity of order n, i.e., we compute
PA(ω
k
n) for all 0 ≤ k ≤ n− 1.
The sequence of values 〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉, is called the
Discrete Fourier Transform (DFT) of the sequence
A = 〈A0, A1, . . . , An−1〉.
The value PA(ω
k
n) is usually denoted by Âk and the sequence of values
〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉 is usually denoted by
 = 〈Â0, Â1, . . . , Ân−1〉.
The DFT Â of a sequence A can be computed VERY FAST using a
divide-and-conquer algorithm called the Fast Fourier Transform.
COMP3121/9101 9 / 33
The Discrete Fourier Transform
Let A = 〈A0, A1, . . . , An−1〉 be a sequence of n real or complex numbers.
We can form the corresponding polynomial PA(x) =
∑n−1
j=0
Ajx
j .
We can evaluate it at all complex roots of unity of order n, i.e., we compute
PA(ω
k
n) for all 0 ≤ k ≤ n− 1.
The sequence of values 〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉, is called the
Discrete Fourier Transform (DFT) of the sequence
A = 〈A0, A1, . . . , An−1〉.
The value PA(ω
k
n) is usually denoted by Âk and the sequence of values
〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉 is usually denoted by
 = 〈Â0, Â1, . . . , Ân−1〉.
The DFT Â of a sequence A can be computed VERY FAST using a
divide-and-conquer algorithm called the Fast Fourier Transform.
COMP3121/9101 9 / 33
The Discrete Fourier Transform
Let A = 〈A0, A1, . . . , An−1〉 be a sequence of n real or complex numbers.
We can form the corresponding polynomial PA(x) =
∑n−1
j=0
Ajx
j .
We can evaluate it at all complex roots of unity of order n, i.e., we compute
PA(ω
k
n) for all 0 ≤ k ≤ n− 1.
The sequence of values 〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉, is called the
Discrete Fourier Transform (DFT) of the sequence
A = 〈A0, A1, . . . , An−1〉.
The value PA(ω
k
n) is usually denoted by Âk and the sequence of values
〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉 is usually denoted by
 = 〈Â0, Â1, . . . , Ân−1〉.
The DFT Â of a sequence A can be computed VERY FAST using a
divide-and-conquer algorithm called the Fast Fourier Transform.
COMP3121/9101 9 / 33
The Discrete Fourier Transform
Let A = 〈A0, A1, . . . , An−1〉 be a sequence of n real or complex numbers.
We can form the corresponding polynomial PA(x) =
∑n−1
j=0
Ajx
j .
We can evaluate it at all complex roots of unity of order n, i.e., we compute
PA(ω
k
n) for all 0 ≤ k ≤ n− 1.
The sequence of values 〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉, is called the
Discrete Fourier Transform (DFT) of the sequence
A = 〈A0, A1, . . . , An−1〉.
The value PA(ω
k
n) is usually denoted by Âk and the sequence of values
〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉 is usually denoted by
 = 〈Â0, Â1, . . . , Ân−1〉.
The DFT Â of a sequence A can be computed VERY FAST using a
divide-and-conquer algorithm called the Fast Fourier Transform.
COMP3121/9101 9 / 33
The Discrete Fourier Transform
Let A = 〈A0, A1, . . . , An−1〉 be a sequence of n real or complex numbers.
We can form the corresponding polynomial PA(x) =
∑n−1
j=0
Ajx
j .
We can evaluate it at all complex roots of unity of order n, i.e., we compute
PA(ω
k
n) for all 0 ≤ k ≤ n− 1.
The sequence of values 〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉, is called the
Discrete Fourier Transform (DFT) of the sequence
A = 〈A0, A1, . . . , An−1〉.
The value PA(ω
k
n) is usually denoted by Âk and the sequence of values
〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ωn−1n )〉 is usually denoted by
 = 〈Â0, Â1, . . . , Ân−1〉.
The DFT Â of a sequence A can be computed VERY FAST using a
divide-and-conquer algorithm called the Fast Fourier Transform.
COMP3121/9101 9 / 33
New way for fast multiplication of polynomials
If we multiply a polynomial
PA(x) = A0 + . . . + An−1x
n−1
of degree n− 1 with a polynomial
PB(x) = B0 + . . . + Bm−1x
m−1
of degree m− 1 we get a polynomial
C(x) = PA(x)PB(x) = C0 + . . . + Cm+n−2x
m+n−2
of degree n− 1 + m− 1 = m + n− 2 with m + n− 1 coefficients.
To uniquely determine such a polynomial C(x) of degree m + n− 2 we
need m + n− 1 many values.
Thus, we will evaluate both PA(x) and PB(x) at all the roots of unity of
order n + m− 1 (instead of at −(n− 1), . . . ,−1, 0, 1, . . . ,m− 1 as we
would in Karatsuba’s method!)
COMP3121/9101 10 / 33
New way for fast multiplication of polynomials
If we multiply a polynomial
PA(x) = A0 + . . . + An−1x
n−1
of degree n− 1 with a polynomial
PB(x) = B0 + . . . + Bm−1x
m−1
of degree m− 1 we get a polynomial
C(x) = PA(x)PB(x) = C0 + . . . + Cm+n−2x
m+n−2
of degree n− 1 + m− 1 = m + n− 2 with m + n− 1 coefficients.
To uniquely determine such a polynomial C(x) of degree m + n− 2 we
need m + n− 1 many values.
Thus, we will evaluate both PA(x) and PB(x) at all the roots of unity of
order n + m− 1 (instead of at −(n− 1), . . . ,−1, 0, 1, . . . ,m− 1 as we
would in Karatsuba’s method!)
COMP3121/9101 10 / 33
New way for fast multiplication of polynomials
If we multiply a polynomial
PA(x) = A0 + . . . + An−1x
n−1
of degree n− 1 with a polynomial
PB(x) = B0 + . . . + Bm−1x
m−1
of degree m− 1 we get a polynomial
C(x) = PA(x)PB(x) = C0 + . . . + Cm+n−2x
m+n−2
of degree n− 1 + m− 1 = m + n− 2 with m + n− 1 coefficients.
To uniquely determine such a polynomial C(x) of degree m + n− 2 we
need m + n− 1 many values.
Thus, we will evaluate both PA(x) and PB(x) at all the roots of unity of
order n + m− 1 (instead of at −(n− 1), . . . ,−1, 0, 1, . . . ,m− 1 as we
would in Karatsuba’s method!)
COMP3121/9101 10 / 33
New way for fast multiplication of polynomials
Note that we defined the DFT of a sequence of length n as the values of the
corresponding polynomial of degree n− 1 at the n roots of unity of order n,
i.e., ωkn (0 ≤ k ≤ n− 1).
So the DFT of a sequence A is another sequence  of exactly the same length;
we do not have an operation which would evaluate a polynomial of degree
n− 1 at m roots of unity of order m 6= n.
For that reason, since we need n+m− 1 values of PC(x) = PA(x)PB(x), we
pad A with m− 1 zeros at the end, (A0, A1, . . . , An−1, 0, . . . , 0︸ ︷︷ ︸
m−1
) to make it of
length m+ n− 1, and similarly we pad B with n− 1 zeros at the end,
(B0, B1, . . . , Bm−1, 0, . . . , 0︸ ︷︷ ︸
n−1
) to also obtain a sequence of length n+m− 1.
Note that this does not change the associated polynomials because the added
higher powers have the corresponding coefficients equal to zero.
COMP3121/9101 11 / 33
New way for fast multiplication of polynomials
Note that we defined the DFT of a sequence of length n as the values of the
corresponding polynomial of degree n− 1 at the n roots of unity of order n,
i.e., ωkn (0 ≤ k ≤ n− 1).
So the DFT of a sequence A is another sequence  of exactly the same length;
we do not have an operation which would evaluate a polynomial of degree
n− 1 at m roots of unity of order m 6= n.
For that reason, since we need n+m− 1 values of PC(x) = PA(x)PB(x), we
pad A with m− 1 zeros at the end, (A0, A1, . . . , An−1, 0, . . . , 0︸ ︷︷ ︸
m−1
) to make it of
length m+ n− 1, and similarly we pad B with n− 1 zeros at the end,
(B0, B1, . . . , Bm−1, 0, . . . , 0︸ ︷︷ ︸
n−1
) to also obtain a sequence of length n+m− 1.
Note that this does not change the associated polynomials because the added
higher powers have the corresponding coefficients equal to zero.
COMP3121/9101 11 / 33
New way for fast multiplication of polynomials
Note that we defined the DFT of a sequence of length n as the values of the
corresponding polynomial of degree n− 1 at the n roots of unity of order n,
i.e., ωkn (0 ≤ k ≤ n− 1).
So the DFT of a sequence A is another sequence  of exactly the same length;
we do not have an operation which would evaluate a polynomial of degree
n− 1 at m roots of unity of order m 6= n.
For that reason, since we need n+m− 1 values of PC(x) = PA(x)PB(x), we
pad A with m− 1 zeros at the end, (A0, A1, . . . , An−1, 0, . . . , 0︸ ︷︷ ︸
m−1
) to make it of
length m+ n− 1, and similarly we pad B with n− 1 zeros at the end,
(B0, B1, . . . , Bm−1, 0, . . . , 0︸ ︷︷ ︸
n−1
) to also obtain a sequence of length n+m− 1.
Note that this does not change the associated polynomials because the added
higher powers have the corresponding coefficients equal to zero.
COMP3121/9101 11 / 33
New way for fast multiplication of polynomials
Note that we defined the DFT of a sequence of length n as the values of the
corresponding polynomial of degree n− 1 at the n roots of unity of order n,
i.e., ωkn (0 ≤ k ≤ n− 1).
So the DFT of a sequence A is another sequence  of exactly the same length;
we do not have an operation which would evaluate a polynomial of degree
n− 1 at m roots of unity of order m 6= n.
For that reason, since we need n+m− 1 values of PC(x) = PA(x)PB(x), we
pad A with m− 1 zeros at the end, (A0, A1, . . . , An−1, 0, . . . , 0︸ ︷︷ ︸
m−1
) to make it of
length m+ n− 1, and similarly we pad B with n− 1 zeros at the end,
(B0, B1, . . . , Bm−1, 0, . . . , 0︸ ︷︷ ︸
n−1
) to also obtain a sequence of length n+m− 1.
Note that this does not change the associated polynomials because the added
higher powers have the corresponding coefficients equal to zero.
COMP3121/9101 11 / 33
New way for fast multiplication of polynomials
We can now compute the DFTs of the two (0 padded) sequences:
DFT (〈A0, A1, . . . , An−1, 0, . . . , 0︸ ︷︷ ︸
m−1
〉) = 〈Â0, Â1, . . . , Ân+m−2)
and
DFT (〈B0, B1, . . . , Bm−1, 0, . . . , 0︸ ︷︷ ︸
n−1
〉) = 〈B̂0, B̂1, . . . , B̂n+m−2)
For each k we multiply the corresponding values Âk = PA(ω
k
n+m−1) and
B̂k = PB(ω
k
n+m−1), thus obtaining
Ĉk = ÂkB̂k = PA(ω
k
n+m−1)PB(ω
k
n+m−1) = PC(ω
k
n+m−1)
We then use the inverse transformation for DFT, called IDFT, to recover the
coefficients 〈C0, C1, . . . , Cn+m−1〉 of the product polynomial PC(x) from the
sequence 〈Ĉ0, Ĉ1, . . . , Ĉn+m−1〉 of its values Ck = PC(ωkn+m−1) at the roots of
unity of order n+m− 1.
COMP3121/9101 12 / 33
New way for fast multiplication of polynomials
We can now compute the DFTs of the two (0 padded) sequences:
DFT (〈A0, A1, . . . , An−1, 0, . . . , 0︸ ︷︷ ︸
m−1
〉) = 〈Â0, Â1, . . . , Ân+m−2)
and
DFT (〈B0, B1, . . . , Bm−1, 0, . . . , 0︸ ︷︷ ︸
n−1
〉) = 〈B̂0, B̂1, . . . , B̂n+m−2)
For each k we multiply the corresponding values Âk = PA(ω
k
n+m−1) and
B̂k = PB(ω
k
n+m−1), thus obtaining
Ĉk = ÂkB̂k = PA(ω
k
n+m−1)PB(ω
k
n+m−1) = PC(ω
k
n+m−1)
We then use the inverse transformation for DFT, called IDFT, to recover the
coefficients 〈C0, C1, . . . , Cn+m−1〉 of the product polynomial PC(x) from the
sequence 〈Ĉ0, Ĉ1, . . . , Ĉn+m−1〉 of its values Ck = PC(ωkn+m−1) at the roots of
unity of order n+m− 1.
COMP3121/9101 12 / 33
New way for fast multiplication of polynomials
We can now compute the DFTs of the two (0 padded) sequences:
DFT (〈A0, A1, . . . , An−1, 0, . . . , 0︸ ︷︷ ︸
m−1
〉) = 〈Â0, Â1, . . . , Ân+m−2)
and
DFT (〈B0, B1, . . . , Bm−1, 0, . . . , 0︸ ︷︷ ︸
n−1
〉) = 〈B̂0, B̂1, . . . , B̂n+m−2)
For each k we multiply the corresponding values Âk = PA(ω
k
n+m−1) and
B̂k = PB(ω
k
n+m−1), thus obtaining
Ĉk = ÂkB̂k = PA(ω
k
n+m−1)PB(ω
k
n+m−1) = PC(ω
k
n+m−1)
We then use the inverse transformation for DFT, called IDFT, to recover the
coefficients 〈C0, C1, . . . , Cn+m−1〉 of the product polynomial PC(x) from the
sequence 〈Ĉ0, Ĉ1, . . . , Ĉn+m−1〉 of its values Ck = PC(ωkn+m−1) at the roots of
unity of order n+m− 1.
COMP3121/9101 12 / 33
New way for fast multiplication of polynomials
PA(x) = A0 +. . .+ An−1x
n−1
+ 0 · xn +. . .+ 0 · xn+m−2;
PB(x) = B0 +. . .+ Bm−1x
m−1
+ 0 · xm +. . .+ 0 · xn+m−2
⇓ DFT ⇓ DFT
{PA(1), PA(ωn+m−1), . . . , PA(ω
n+m−2
n+m−1)}; {PB(1), PB(ωn+m−1), . . . , PB(ω
n+m−2
n+m−1)}
⇓ multiplication
{
PA(1)PB(1)︸ ︷︷ ︸
PC (1)
, PA(ωn+m−1)PB(ωn+m−1)︸ ︷︷ ︸
PC (ωn+m−1)
, . . . , PA(ω
n+m−2
n+m−1)PB(ω
n+m−2
n+m−1)︸ ︷︷ ︸
PC (ω
n+m−2
n+m−1)
}
⇓ IDFT
PC(x) = PA(x) · PB(x) =
n+m−2∑
j=0
Cjx
j
=
n+m−2∑
j=0
( j∑
i=0
AiBj−i
︸ ︷︷ ︸
Cj
)
x
j
COMP3121/9101 13 / 33
The Fast Fourier Transform (FFT)
Crucial fact: the values PA(ω
k
n) for all k such that 0 ≤ k < n can
be computed in O(n logn) time!
Note that a direct evaluation of a polynomial of degree n− 1 at n
roots of unity of order n would take n2 many multiplications, even
if we precompute all powers ωkmn , because we have to perform
multiplications Am · (ωkn)m for all 0 ≤ m ≤ n− 1 and all
0 ≤ k ≤ n− 1.
We can assume that n 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 n which is not a power of two the
smallest power of 2 larger or equal to n is smaller than 2n.
Hint: consider n in binary. How many bits does the nearest power
of two have?
COMP3121/9101 14 / 33
The Fast Fourier Transform (FFT)
Crucial fact: the values PA(ω
k
n) for all k such that 0 ≤ k < n can
be computed in O(n logn) time!
Note that a direct evaluation of a polynomial of degree n− 1 at n
roots of unity of order n would take n2 many multiplications, even
if we precompute all powers ωkmn , because we have to perform
multiplications Am · (ωkn)m for all 0 ≤ m ≤ n− 1 and all
0 ≤ k ≤ n− 1.
We can assume that n 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 n which is not a power of two the
smallest power of 2 larger or equal to n is smaller than 2n.
Hint: consider n in binary. How many bits does the nearest power
of two have?
COMP3121/9101 14 / 33
The Fast Fourier Transform (FFT)
Crucial fact: the values PA(ω
k
n) for all k such that 0 ≤ k < n can
be computed in O(n logn) time!
Note that a direct evaluation of a polynomial of degree n− 1 at n
roots of unity of order n would take n2 many multiplications, even
if we precompute all powers ωkmn , because we have to perform
multiplications Am · (ωkn)m for all 0 ≤ m ≤ n− 1 and all
0 ≤ k ≤ n− 1.
We can assume that n 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 n which is not a power of two the
smallest power of 2 larger or equal to n is smaller than 2n.
Hint: consider n in binary. How many bits does the nearest power
of two have?
COMP3121/9101 14 / 33
The Fast Fourier Transform (FFT)
Crucial fact: the values PA(ω
k
n) for all k such that 0 ≤ k < n can
be computed in O(n logn) time!
Note that a direct evaluation of a polynomial of degree n− 1 at n
roots of unity of order n would take n2 many multiplications, even
if we precompute all powers ωkmn , because we have to perform
multiplications Am · (ωkn)m for all 0 ≤ m ≤ n− 1 and all
0 ≤ k ≤ n− 1.
We can assume that n 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 n which is not a power of two the
smallest power of 2 larger or equal to n is smaller than 2n.
Hint: consider n in binary. How many bits does the nearest power
of two have?
COMP3121/9101 14 / 33
The Fast Fourier Transform (FFT)
Crucial fact: the values PA(ω
k
n) for all k such that 0 ≤ k < n can
be computed in O(n logn) time!
Note that a direct evaluation of a polynomial of degree n− 1 at n
roots of unity of order n would take n2 many multiplications, even
if we precompute all powers ωkmn , because we have to perform
multiplications Am · (ωkn)m for all 0 ≤ m ≤ n− 1 and all
0 ≤ k ≤ n− 1.
We can assume that n 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 n which is not a power of two the
smallest power of 2 larger or equal to n is smaller than 2n.
Hint: consider n in binary. How many bits does the nearest power
of two have?
COMP3121/9101 14 / 33
The Fast Fourier Transform (FFT)
Problem: Given a sequence A = 〈A0, A1, . . . , An〉 compute its DFT.
This amounts to finding values of PA(x) for all x = ω
k
n, 0 ≤ k ≤ n− 1.
The main idea of the FFT algorithm: divide-and-conquer by splitting the
polynomial PA(x) into the even powers and the odd powers:
PA(x) = (A0 + A2x
2 + A4x
4 + . . . + An−2x
n−2) + (A1x + A3x
3 + . . . + An−1x
n−1)
= A0 + A2x
2 + A4(x
2)2 + . . . + An−2(x
2)n/2−1
+ x
(
A1 + A3x
2 + A5(x
2)2 + . . . + An−1(x
2)n/2−1
)
Let us define A[0] = 〈A0, A2, A4, . . . An−2〉 and A[1] = 〈A1, A3, A5, . . . An−1〉; then
P
A[0]
(y) = A0 + A2y + A4y
2 + . . . + An−2y
n/2−1
P
A[1]
(y) = A1 + A3y + A5y
2 + . . . + An−1y
n/2−1
PA(x) = PA[0] (x
2) + xP
A[1]
(x2)
Note that the number of coefficients of the polynomials P
A[0]
(y) and P
A[1]
(y) is n/2
each, while the number of coefficients of the polynomial PA(x) is n.
COMP3121/9101 15 / 33
The Fast Fourier Transform (FFT)
Problem: Given a sequence A = 〈A0, A1, . . . , An〉 compute its DFT.
This amounts to finding values of PA(x) for all x = ω
k
n, 0 ≤ k ≤ n− 1.
The main idea of the FFT algorithm: divide-and-conquer by splitting the
polynomial PA(x) into the even powers and the odd powers:
PA(x) = (A0 + A2x
2 + A4x
4 + . . . + An−2x
n−2) + (A1x + A3x
3 + . . . + An−1x
n−1)
= A0 + A2x
2 + A4(x
2)2 + . . . + An−2(x
2)n/2−1
+ x
(
A1 + A3x
2 + A5(x
2)2 + . . . + An−1(x
2)n/2−1
)
Let us define A[0] = 〈A0, A2, A4, . . . An−2〉 and A[1] = 〈A1, A3, A5, . . . An−1〉; then
P
A[0]
(y) = A0 + A2y + A4y
2 + . . . + An−2y
n/2−1
P
A[1]
(y) = A1 + A3y + A5y
2 + . . . + An−1y
n/2−1
PA(x) = PA[0] (x
2) + xP
A[1]
(x2)
Note that the number of coefficients of the polynomials P
A[0]
(y) and P
A[1]
(y) is n/2
each, while the number of coefficients of the polynomial PA(x) is n.
COMP3121/9101 15 / 33
The Fast Fourier Transform (FFT)
Problem: Given a sequence A = 〈A0, A1, . . . , An〉 compute its DFT.
This amounts to finding values of PA(x) for all x = ω
k
n, 0 ≤ k ≤ n− 1.
The main idea of the FFT algorithm: divide-and-conquer by splitting the
polynomial PA(x) into the even powers and the odd powers:
PA(x) = (A0 + A2x
2 + A4x
4 + . . . + An−2x
n−2) + (A1x + A3x
3 + . . . + An−1x
n−1)
= A0 + A2x
2 + A4(x
2)2 + . . . + An−2(x
2)n/2−1
+ x
(
A1 + A3x
2 + A5(x
2)2 + . . . + An−1(x
2)n/2−1
)
Let us define A[0] = 〈A0, A2, A4, . . . An−2〉 and A[1] = 〈A1, A3, A5, . . . An−1〉; then
P
A[0]
(y) = A0 + A2y + A4y
2 + . . . + An−2y
n/2−1
P
A[1]
(y) = A1 + A3y + A5y
2 + . . . + An−1y
n/2−1
PA(x) = PA[0] (x
2) + xP
A[1]
(x2)
Note that the number of coefficients of the polynomials P
A[0]
(y) and P
A[1]
(y) is n/2
each, while the number of coefficients of the polynomial PA(x) is n.
COMP3121/9101 15 / 33
The Fast Fourier Transform (FFT)
Problem: Given a sequence A = 〈A0, A1, . . . , An〉 compute its DFT.
This amounts to finding values of PA(x) for all x = ω
k
n, 0 ≤ k ≤ n− 1.
The main idea of the FFT algorithm: divide-and-conquer by splitting the
polynomial PA(x) into the even powers and the odd powers:
PA(x) = (A0 + A2x
2 + A4x
4 + . . . + An−2x
n−2) + (A1x + A3x
3 + . . . + An−1x
n−1)
= A0 + A2x
2 + A4(x
2)2 + . . . + An−2(x
2)n/2−1
+ x
(
A1 + A3x
2 + A5(x
2)2 + . . . + An−1(x
2)n/2−1
)
Let us define A[0] = 〈A0, A2, A4, . . . An−2〉 and A[1] = 〈A1, A3, A5, . . . An−1〉; then
P
A[0]
(y) = A0 + A2y + A4y
2 + . . . + An−2y
n/2−1
P
A[1]
(y) = A1 + A3y + A5y
2 + . . . + An−1y
n/2−1
PA(x) = PA[0] (x
2) + xP
A[1]
(x2)
Note that the number of coefficients of the polynomials P
A[0]
(y) and P
A[1]
(y) is n/2
each, while the number of coefficients of the polynomial PA(x) is n.
COMP3121/9101 15 / 33
The Fast Fourier Transform (FFT)
Problem: Given a sequence A = 〈A0, A1, . . . , An〉 compute its DFT.
This amounts to finding values of PA(x) for all x = ω
k
n, 0 ≤ k ≤ n− 1.
The main idea of the FFT algorithm: divide-and-conquer by splitting the
polynomial PA(x) into the even powers and the odd powers:
PA(x) = (A0 + A2x
2 + A4x
4 + . . . + An−2x
n−2) + (A1x + A3x
3 + . . . + An−1x
n−1)
= A0 + A2x
2 + A4(x
2)2 + . . . + An−2(x
2)n/2−1
+ x
(
A1 + A3x
2 + A5(x
2)2 + . . . + An−1(x
2)n/2−1
)
Let us define A[0] = 〈A0, A2, A4, . . . An−2〉 and A[1] = 〈A1, A3, A5, . . . An−1〉; then
P
A[0]
(y) = A0 + A2y + A4y
2 + . . . + An−2y
n/2−1
P
A[1]
(y) = A1 + A3y + A5y
2 + . . . + An−1y
n/2−1
PA(x) = PA[0] (x
2) + xP
A[1]
(x2)
Note that the number of coefficients of the polynomials P
A[0]
(y) and P
A[1]
(y) is n/2
each, while the number of coefficients of the polynomial PA(x) is n.
COMP3121/9101 15 / 33
The Fast Fourier Transform (FFT)
Problem of size n:
Evaluate a polynomial with n coefficients at n many roots of unity.
Problem of size n/2:
Evaluate a polynomial with n/2 coefficients at n/2 many roots of unity.
We reduced evaluation of our polynomial PA(x) with n coefficients at inputs
x = ω0n, x = ω
1
n, x = ω
2
n, . . . , x = ω
n−1
n to evaluation of two polynomials PA[0] (y) and
P
A[1]
(y) each with n/2 coefficients, at points y = x2 for the same values of inputs x.
However, as x ranges through values {ω0n, ω1n, ω2n, . . . , ω
n−1
n }, the value of y = x2
ranges through {ω0
n/2
, ω1
n/2
, ω2
n/2
, . . . , ωn−1
n/2
}, and there are only n/2 distinct
such values.
Once we get these n/2 values of A[0](x2) and A[1](x2) we need n additional
multiplications with numbers ωkn to obtain the values of
PA(ω
k
n) = PA[0] ((ω
k
n)
2) + ωkn · PA[1] ((ω
k
n)
2)
= P
A[0]
(ωkn/2) + ω
k
n · PA[1] (ω
k
n/2).
Thus, we have reduced a problem of size n to two such problems of size n/2, plus a
linear overhead.
COMP3121/9101 16 / 33
The Fast Fourier Transform (FFT)
Problem of size n:
Evaluate a polynomial with n coefficients at n many roots of unity.
Problem of size n/2:
Evaluate a polynomial with n/2 coefficients at n/2 many roots of unity.
We reduced evaluation of our polynomial PA(x) with n coefficients at inputs
x = ω0n, x = ω
1
n, x = ω
2
n, . . . , x = ω
n−1
n to evaluation of two polynomials PA[0] (y) and
P
A[1]
(y) each with n/2 coefficients, at points y = x2 for the same values of inputs x.
However, as x ranges through values {ω0n, ω1n, ω2n, . . . , ω
n−1
n }, the value of y = x2
ranges through {ω0
n/2
, ω1
n/2
, ω2
n/2
, . . . , ωn−1
n/2
}, and there are only n/2 distinct
such values.
Once we get these n/2 values of A[0](x2) and A[1](x2) we need n additional
multiplications with numbers ωkn to obtain the values of
PA(ω
k
n) = PA[0] ((ω
k
n)
2) + ωkn · PA[1] ((ω
k
n)
2)
= P
A[0]
(ωkn/2) + ω
k
n · PA[1] (ω
k
n/2).
Thus, we have reduced a problem of size n to two such problems of size n/2, plus a
linear overhead.
COMP3121/9101 16 / 33
The Fast Fourier Transform (FFT)
Problem of size n:
Evaluate a polynomial with n coefficients at n many roots of unity.
Problem of size n/2:
Evaluate a polynomial with n/2 coefficients at n/2 many roots of unity.
We reduced evaluation of our polynomial PA(x) with n coefficients at inputs
x = ω0n, x = ω
1
n, x = ω
2
n, . . . , x = ω
n−1
n to evaluation of two polynomials PA[0] (y) and
P
A[1]
(y) each with n/2 coefficients, at points y = x2 for the same values of inputs x.
However, as x ranges through values {ω0n, ω1n, ω2n, . . . , ω
n−1
n }, the value of y = x2
ranges through {ω0
n/2
, ω1
n/2
, ω2
n/2
, . . . , ωn−1
n/2
}, and there are only n/2 distinct
such values.
Once we get these n/2 values of A[0](x2) and A[1](x2) we need n additional
multiplications with numbers ωkn to obtain the values of
PA(ω
k
n) = PA[0] ((ω
k
n)
2) + ωkn · PA[1] ((ω
k
n)
2)
= P
A[0]
(ωkn/2) + ω
k
n · PA[1] (ω
k
n/2).
Thus, we have reduced a problem of size n to two such problems of size n/2, plus a
linear overhead.
COMP3121/9101 16 / 33
The Fast Fourier Transform (FFT)
Problem of size n:
Evaluate a polynomial with n coefficients at n many roots of unity.
Problem of size n/2:
Evaluate a polynomial with n/2 coefficients at n/2 many roots of unity.
We reduced evaluation of our polynomial PA(x) with n coefficients at inputs
x = ω0n, x = ω
1
n, x = ω
2
n, . . . , x = ω
n−1
n to evaluation of two polynomials PA[0] (y) and
P
A[1]
(y) each with n/2 coefficients, at points y = x2 for the same values of inputs x.
However, as x ranges through values {ω0n, ω1n, ω2n, . . . , ω
n−1
n }, the value of y = x2
ranges through {ω0
n/2
, ω1
n/2
, ω2
n/2
, . . . , ωn−1
n/2
}, and there are only n/2 distinct
such values.
Once we get these n/2 values of A[0](x2) and A[1](x2) we need n additional
multiplications with numbers ωkn to obtain the values of
PA(ω
k
n) = PA[0] ((ω
k
n)
2) + ωkn · PA[1] ((ω
k
n)
2)
= P
A[0]
(ωkn/2) + ω
k
n · PA[1] (ω
k
n/2).
Thus, we have reduced a problem of size n to two such problems of size n/2, plus a
linear overhead.
COMP3121/9101 16 / 33
The Fast Fourier Transform (FFT)
Problem of size n:
Evaluate a polynomial with n coefficients at n many roots of unity.
Problem of size n/2:
Evaluate a polynomial with n/2 coefficients at n/2 many roots of unity.
We reduced evaluation of our polynomial PA(x) with n coefficients at inputs
x = ω0n, x = ω
1
n, x = ω
2
n, . . . , x = ω
n−1
n to evaluation of two polynomials PA[0] (y) and
P
A[1]
(y) each with n/2 coefficients, at points y = x2 for the same values of inputs x.
However, as x ranges through values {ω0n, ω1n, ω2n, . . . , ω
n−1
n }, the value of y = x2
ranges through {ω0
n/2
, ω1
n/2
, ω2
n/2
, . . . , ωn−1
n/2
}, and there are only n/2 distinct
such values.
Once we get these n/2 values of A[0](x2) and A[1](x2) we need n additional
multiplications with numbers ωkn to obtain the values of
PA(ω
k
n) = PA[0] ((ω
k
n)
2) + ωkn · PA[1] ((ω
k
n)
2)
= P
A[0]
(ωkn/2) + ω
k
n · PA[1] (ω
k
n/2).
Thus, we have reduced a problem of size n to two such problems of size n/2, plus a
linear overhead.
COMP3121/9101 16 / 33
The Fast Fourier Transform (FFT)
Problem of size n:
Evaluate a polynomial with n coefficients at n many roots of unity.
Problem of size n/2:
Evaluate a polynomial with n/2 coefficients at n/2 many roots of unity.
We reduced evaluation of our polynomial PA(x) with n coefficients at inputs
x = ω0n, x = ω
1
n, x = ω
2
n, . . . , x = ω
n−1
n to evaluation of two polynomials PA[0] (y) and
P
A[1]
(y) each with n/2 coefficients, at points y = x2 for the same values of inputs x.
However, as x ranges through values {ω0n, ω1n, ω2n, . . . , ω
n−1
n }, the value of y = x2
ranges through {ω0
n/2
, ω1
n/2
, ω2
n/2
, . . . , ωn−1
n/2
}, and there are only n/2 distinct
such values.
Once we get these n/2 values of A[0](x2) and A[1](x2) we need n additional
multiplications with numbers ωkn to obtain the values of
PA(ω
k
n) = PA[0] ((ω
k
n)
2) + ωkn · PA[1] ((ω
k
n)
2)
= P
A[0]
(ωkn/2) + ω
k
n · PA[1] (ω
k
n/2).
Thus, we have reduced a problem of size n to two such problems of size n/2, plus a
linear overhead.
COMP3121/9101 16 / 33
The Fast Fourier Transform (FFT) - a simplification
Note that by the Cancelation Lemma ω
n/2
n = ω
n/2
2n/2
= ω2 = −1.
Thus,
ω
k+n/2
n = ω
n/2
n ω
k
n = ω2ω
k
n = −ω
k
n;
We can now simplify evaluation of
PA(ω
k
n) = PA[0] (ω
k
n/2) + ω
k
nPA[1] (ω
k
n/2)
for n/2 ≤ k < n as follows:
let k = n/2 + m where 0 ≤ m < n/2; then
PA(ω
n/2+m
n ) = PA[0] (ω
n/2+m
n/2
) + ω
n/2+m
n PA[1] (ω
n/2+m
n/2
)
= P
A[0]
(ω
n/2
n/2
ωmn/2) + ω
n/2
n ω
m
n PA[1] (ω
n/2
n/2
ωmn/2)
= P
A[0]
(ωmn/2)− ω
m
n PA[1] (ω
m
n/2)
Compare this with PA(ω
m
n ) = PA[0] (ω
m
n/2
) + ωmn PA[1] (ω
m
n/2
) for 0 ≤ m < n/2.
COMP3121/9101 17 / 33
The Fast Fourier Transform (FFT) - a simplification
Note that by the Cancelation Lemma ω
n/2
n = ω
n/2
2n/2
= ω2 = −1.
Thus,
ω
k+n/2
n = ω
n/2
n ω
k
n = ω2ω
k
n = −ω
k
n;
We can now simplify evaluation of
PA(ω
k
n) = PA[0] (ω
k
n/2) + ω
k
nPA[1] (ω
k
n/2)
for n/2 ≤ k < n as follows:
let k = n/2 + m where 0 ≤ m < n/2; then
PA(ω
n/2+m
n ) = PA[0] (ω
n/2+m
n/2
) + ω
n/2+m
n PA[1] (ω
n/2+m
n/2
)
= P
A[0]
(ω
n/2
n/2
ωmn/2) + ω
n/2
n ω
m
n PA[1] (ω
n/2
n/2
ωmn/2)
= P
A[0]
(ωmn/2)− ω
m
n PA[1] (ω
m
n/2)
Compare this with PA(ω
m
n ) = PA[0] (ω
m
n/2
) + ωmn PA[1] (ω
m
n/2
) for 0 ≤ m < n/2.
COMP3121/9101 17 / 33
The Fast Fourier Transform (FFT) - a simplification
Note that by the Cancelation Lemma ω
n/2
n = ω
n/2
2n/2
= ω2 = −1.
Thus,
ω
k+n/2
n = ω
n/2
n ω
k
n = ω2ω
k
n = −ω
k
n;
We can now simplify evaluation of
PA(ω
k
n) = PA[0] (ω
k
n/2) + ω
k
nPA[1] (ω
k
n/2)
for n/2 ≤ k < n as follows:
let k = n/2 + m where 0 ≤ m < n/2; then
PA(ω
n/2+m
n ) = PA[0] (ω
n/2+m
n/2
) + ω
n/2+m
n PA[1] (ω
n/2+m
n/2
)
= P
A[0]
(ω
n/2
n/2
ωmn/2) + ω
n/2
n ω
m
n PA[1] (ω
n/2
n/2
ωmn/2)
= P
A[0]
(ωmn/2)− ω
m
n PA[1] (ω
m
n/2)
Compare this with PA(ω
m
n ) = PA[0] (ω
m
n/2
) + ωmn PA[1] (ω
m
n/2
) for 0 ≤ m < n/2.
COMP3121/9101 17 / 33
The Fast Fourier Transform (FFT) - a simplification
Note that by the Cancelation Lemma ω
n/2
n = ω
n/2
2n/2
= ω2 = −1.
Thus,
ω
k+n/2
n = ω
n/2
n ω
k
n = ω2ω
k
n = −ω
k
n;
We can now simplify evaluation of
PA(ω
k
n) = PA[0] (ω
k
n/2) + ω
k
nPA[1] (ω
k
n/2)
for n/2 ≤ k < n as follows:
let k = n/2 + m where 0 ≤ m < n/2; then
PA(ω
n/2+m
n ) = PA[0] (ω
n/2+m
n/2
) + ω
n/2+m
n PA[1] (ω
n/2+m
n/2
)
= P
A[0]
(ω
n/2
n/2
ωmn/2) + ω
n/2
n ω
m
n PA[1] (ω
n/2
n/2
ωmn/2)
= P
A[0]
(ωmn/2)− ω
m
n PA[1] (ω
m
n/2)
Compare this with PA(ω
m
n ) = PA[0] (ω
m
n/2
) + ωmn PA[1] (ω
m
n/2
) for 0 ≤ m < n/2.
COMP3121/9101 17 / 33
The Fast Fourier Transform (FFT) - a simplification
So we can replace evaluations of
PA(ω
k
n) = PA[0] (ω
k
n/2) + ω
k
nPA[1] (ω
k
n/2)
for k = 0 to k = n− 1
with such evaluations only for k = 0 to k = n/2− 1
and just let for k = 0 to k = n/2− 1
PA(ω
k
n) = PA[0] (ω
k
n/2) + ω
k
nPA[1] (ω
k
n/2)
PA(ω
n/2+k
n ) = PA[0] (ω
k
n/2)− ω
k
nPA[1] (ω
k
n/2)
We can now write a pseudo-code for our FFT algorithm:
COMP3121/9101 18 / 33
The Fast Fourier Transform (FFT) - a simplification
So we can replace evaluations of
PA(ω
k
n) = PA[0] (ω
k
n/2) + ω
k
nPA[1] (ω
k
n/2)
for k = 0 to k = n− 1
with such evaluations only for k = 0 to k = n/2− 1
and just let for k = 0 to k = n/2− 1
PA(ω
k
n) = PA[0] (ω
k
n/2) + ω
k
nPA[1] (ω
k
n/2)
PA(ω
n/2+k
n ) = PA[0] (ω
k
n/2)− ω
k
nPA[1] (ω
k
n/2)
We can now write a pseudo-code for our FFT algorithm:
COMP3121/9101 18 / 33
FFT algorithm
1: function FFT(A)
2: n← length[A]
3: if n = 1 then return A
4: else
5: A[0] ← (A0, A2, . . . An−2);
6: A[1] ← (A1, A3, . . . An−1);
7: y[0] ← FFT (A[0]);
8: y[1] ← FFT (A[1]);
9: ωn ← ei
2π
n ;
10: ω ← 1; % a variable to hold powers of ωn
11: for k = 0 to k = n/2− 1 do: % PA(ωkn)︸ ︷︷ ︸
yk
= P
A[0]
(ωkn/2)︸ ︷︷ ︸
y
[0]
k
+ωkn PA[1] (ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
12: yk ← y
[0]
k
+ ω · y[1]
k
;
13: yn/2+k ← y
[0]
k
− ω · y[1]
k
% PA(ω
n/2+k
n )︸ ︷︷ ︸
yn/2+k
= P
A[0]
(ωkn/2)︸ ︷︷ ︸
y
[0]
k
−ωkn PA[1] (ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
14: ω ← ω · ωn;
15: end for
16: return y
17: end if
18: end function
COMP3121/9101 19 / 33
How fast is the Fast Fourier Transform?
We have recursively reduced evaluation of a polynomial PA(x) with n
coefficients at n roots of unity of order n to evaluations of two polynomials
PA[0](y) and PA[1](y), each with n/2 coefficients, at n/2 many roots of unity of
order n/2.
Once we get these n/2 values of PA[0](y) and PA[1](y) we need n/2 additional
multiplications to obtain the values of
PA(ω
k
n)︸ ︷︷ ︸
yk
= PA[0](ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
+ω
k
n PA[1](ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(1)
and
PA(ω
n/2+k
n )︸ ︷︷ ︸
yn/2+k
= A
[0]
(ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
−ωknA
[1]
(ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(2)
for all 0 ≤ k < n/2.
Thus, we have reduced a problem of size n to two such problems of size n/2,
plus a linear overhead.
Consequently, our algorithm’s run time satisfies the recurrence
T (n) = 2T (n/2) + c n
The Master Theorem gives T (n) = Θ(n logn).
COMP3121/9101 20 / 33
How fast is the Fast Fourier Transform?
We have recursively reduced evaluation of a polynomial PA(x) with n
coefficients at n roots of unity of order n to evaluations of two polynomials
PA[0](y) and PA[1](y), each with n/2 coefficients, at n/2 many roots of unity of
order n/2.
Once we get these n/2 values of PA[0](y) and PA[1](y) we need n/2 additional
multiplications to obtain the values of
PA(ω
k
n)︸ ︷︷ ︸
yk
= PA[0](ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
+ω
k
n PA[1](ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(1)
and
PA(ω
n/2+k
n )︸ ︷︷ ︸
yn/2+k
= A
[0]
(ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
−ωknA
[1]
(ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(2)
for all 0 ≤ k < n/2.
Thus, we have reduced a problem of size n to two such problems of size n/2,
plus a linear overhead.
Consequently, our algorithm’s run time satisfies the recurrence
T (n) = 2T (n/2) + c n
The Master Theorem gives T (n) = Θ(n logn).
COMP3121/9101 20 / 33
How fast is the Fast Fourier Transform?
We have recursively reduced evaluation of a polynomial PA(x) with n
coefficients at n roots of unity of order n to evaluations of two polynomials
PA[0](y) and PA[1](y), each with n/2 coefficients, at n/2 many roots of unity of
order n/2.
Once we get these n/2 values of PA[0](y) and PA[1](y) we need n/2 additional
multiplications to obtain the values of
PA(ω
k
n)︸ ︷︷ ︸
yk
= PA[0](ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
+ω
k
n PA[1](ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(1)
and
PA(ω
n/2+k
n )︸ ︷︷ ︸
yn/2+k
= A
[0]
(ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
−ωknA
[1]
(ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(2)
for all 0 ≤ k < n/2.
Thus, we have reduced a problem of size n to two such problems of size n/2,
plus a linear overhead.
Consequently, our algorithm’s run time satisfies the recurrence
T (n) = 2T (n/2) + c n
The Master Theorem gives T (n) = Θ(n logn).
COMP3121/9101 20 / 33
How fast is the Fast Fourier Transform?
We have recursively reduced evaluation of a polynomial PA(x) with n
coefficients at n roots of unity of order n to evaluations of two polynomials
PA[0](y) and PA[1](y), each with n/2 coefficients, at n/2 many roots of unity of
order n/2.
Once we get these n/2 values of PA[0](y) and PA[1](y) we need n/2 additional
multiplications to obtain the values of
PA(ω
k
n)︸ ︷︷ ︸
yk
= PA[0](ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
+ω
k
n PA[1](ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(1)
and
PA(ω
n/2+k
n )︸ ︷︷ ︸
yn/2+k
= A
[0]
(ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
−ωknA
[1]
(ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(2)
for all 0 ≤ k < n/2.
Thus, we have reduced a problem of size n to two such problems of size n/2,
plus a linear overhead.
Consequently, our algorithm’s run time satisfies the recurrence
T (n) = 2T (n/2) + c n
The Master Theorem gives T (n) = Θ(n logn).
COMP3121/9101 20 / 33
How fast is the Fast Fourier Transform?
We have recursively reduced evaluation of a polynomial PA(x) with n
coefficients at n roots of unity of order n to evaluations of two polynomials
PA[0](y) and PA[1](y), each with n/2 coefficients, at n/2 many roots of unity of
order n/2.
Once we get these n/2 values of PA[0](y) and PA[1](y) we need n/2 additional
multiplications to obtain the values of
PA(ω
k
n)︸ ︷︷ ︸
yk
= PA[0](ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
+ω
k
n PA[1](ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(1)
and
PA(ω
n/2+k
n )︸ ︷︷ ︸
yn/2+k
= A
[0]
(ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
−ωknA
[1]
(ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(2)
for all 0 ≤ k < n/2.
Thus, we have reduced a problem of size n to two such problems of size n/2,
plus a linear overhead.
Consequently, our algorithm’s run time satisfies the recurrence
T (n) = 2T (n/2) + c n
The Master Theorem gives T (n) = Θ(n logn).
COMP3121/9101 20 / 33
How fast is the Fast Fourier Transform?
We have recursively reduced evaluation of a polynomial PA(x) with n
coefficients at n roots of unity of order n to evaluations of two polynomials
PA[0](y) and PA[1](y), each with n/2 coefficients, at n/2 many roots of unity of
order n/2.
Once we get these n/2 values of PA[0](y) and PA[1](y) we need n/2 additional
multiplications to obtain the values of
PA(ω
k
n)︸ ︷︷ ︸
yk
= PA[0](ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
+ω
k
n PA[1](ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(1)
and
PA(ω
n/2+k
n )︸ ︷︷ ︸
yn/2+k
= A
[0]
(ω
k
n/2)︸ ︷︷ ︸
y
[0]
k
−ωknA
[1]
(ω
k
n/2)︸ ︷︷ ︸
y
[1]
k
(2)
for all 0 ≤ k < n/2.
Thus, we have reduced a problem of size n to two such problems of size n/2,
plus a linear overhead.
Consequently, our algorithm’s run time satisfies the recurrence
T (n) = 2T (n/2) + c n
The Master Theorem gives T (n) = Θ(n logn).
COMP3121/9101 20 / 33
Matrix representation of polynomial evaluation
Evaluation of a polynomial PA(x) = A0 + A1x + . . . + An−1x
n−1 at roots of unity ωkn of
order n can be represented in the matrix form as follows:
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
A0
A1
A2
.
.
.
An−1
=
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
=
Â0
Â1
Â2
.
.
.
Ân−1
The FFT is just a method of replacing this matrix-vector multiplication taking n2 many
multiplications with an n log n procedure.
From PA(1) = PA(ω
0
n), PA(ωn), PA(ω
2
n), . . . , PA(ω
n−1
n ), we get the coefficients from
A0
A1
A2
.
.
.
An−1
=
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
−1
Â0
Â1
Â2
.
.
.
Ân−1
(3)
COMP3121/9101 21 / 33
Matrix representation of polynomial evaluation
Evaluation of a polynomial PA(x) = A0 + A1x + . . . + An−1x
n−1 at roots of unity ωkn of
order n can be represented in the matrix form as follows:
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
A0
A1
A2
.
.
.
An−1
=
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
=
Â0
Â1
Â2
.
.
.
Ân−1
The FFT is just a method of replacing this matrix-vector multiplication taking n2 many
multiplications with an n log n procedure.
From PA(1) = PA(ω
0
n), PA(ωn), PA(ω
2
n), . . . , PA(ω
n−1
n ), we get the coefficients from
A0
A1
A2
.
.
.
An−1
=
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
−1
Â0
Â1
Â2
.
.
.
Ân−1
(3)
COMP3121/9101 21 / 33
Matrix representation of polynomial evaluation
Evaluation of a polynomial PA(x) = A0 + A1x + . . . + An−1x
n−1 at roots of unity ωkn of
order n can be represented in the matrix form as follows:
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
A0
A1
A2
.
.
.
An−1
=
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
=
Â0
Â1
Â2
.
.
.
Ân−1
The FFT is just a method of replacing this matrix-vector multiplication taking n2 many
multiplications with an n log n procedure.
From PA(1) = PA(ω
0
n), PA(ωn), PA(ω
2
n), . . . , PA(ω
n−1
n ), we get the coefficients from
A0
A1
A2
.
.
.
An−1
=
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
−1
Â0
Â1
Â2
.
.
.
Ân−1
(3)
COMP3121/9101 21 / 33
Matrix representation of polynomial evaluation
Evaluation of a polynomial PA(x) = A0 + A1x + . . . + An−1x
n−1 at roots of unity ωkn of
order n can be represented in the matrix form as follows:
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
A0
A1
A2
.
.
.
An−1
=
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
=
Â0
Â1
Â2
.
.
.
Ân−1
The FFT is just a method of replacing this matrix-vector multiplication taking n2 many
multiplications with an n log n procedure.
From PA(1) = PA(ω
0
n), PA(ωn), PA(ω
2
n), . . . , PA(ω
n−1
n ), we get the coefficients from
A0
A1
A2
.
.
.
An−1
=
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
−1
Â0
Â1
Â2
.
.
.
Ân−1
(3)
COMP3121/9101 21 / 33
Another remarkable feature of the roots of unity:
To obtain the inverse of the above matrix, all we have to do is just change the signs of
the exponents and divide everything by n:
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
−1
=
1
n
1 1 1 . . . 1
1 ω−1n ω
−2
n . . . ω
−(n−1)
n
1 ω−2n ω
−2·2
n . . . ω
−2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ω
−(n−1)
n ω
−2(n−1)
n . . . ω
−(n−1)(n−1)
n
COMP3121/9101 22 / 33
To see this, note that if we compute the product
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
1 1 1 . . . 1
1 ω−1n ω
−2
n . . . ω
−(n−1)
n
1 ω−2n ω
−2·2
n . . . ω
−2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ω
−(n−1)
n ω
−2(n−1)
n . . . ω
−(n−1)(n−1)
n
the (i, j) entry in the product matrix is equal to a product of ith row and jth
column:
(
1 ωin ω
2·i
n . . . ω
i·(n−1)
n
)
1
ω−jn
ω−2jn
...
ω
−(n−1)j
n
=
n−1∑
k=0
ωikn ω
−jk
n =
n−1∑
k=0
ω(i−j)kn
COMP3121/9101 23 / 33
We now have two possibilities:
1 i = j: then
n−1∑
k=0
ω(i−j)kn =
n−1∑
k=0
ω0n =
n−1∑
k=0
1 = n;
2 i 6= j: then
∑n−1
k=0 ω
(i−j)k
n represents a sum of a geometric progression
with the ratio ωi−jn and thus
n−1∑
k=0
ω(i−j)kn =
1− ω(i−j)nn
1− ωi−jn
=
1− (ωnn)i−j
1− ωi−jn
=
1− 1
1− ωi−jn
= 0
So,
(
1 ωin ω
2·i
n . . . ω
i·(n−1)
n
)
1
ω−jn
ω−2jn
...
ω
−(n−1)j
n
=
n−1∑
k=0
ω(i−j)kn =
{
n if i = j
0 if i 6= j
(4)
COMP3121/9101 24 / 33
We now have two possibilities:
1 i = j: then
n−1∑
k=0
ω(i−j)kn =
n−1∑
k=0
ω0n =
n−1∑
k=0
1 = n;
2 i 6= j: then
∑n−1
k=0 ω
(i−j)k
n represents a sum of a geometric progression
with the ratio ωi−jn and thus
n−1∑
k=0
ω(i−j)kn =
1− ω(i−j)nn
1− ωi−jn
=
1− (ωnn)i−j
1− ωi−jn
=
1− 1
1− ωi−jn
= 0
So,
(
1 ωin ω
2·i
n . . . ω
i·(n−1)
n
)
1
ω−jn
ω−2jn
...
ω
−(n−1)j
n
=
n−1∑
k=0
ω(i−j)kn =
{
n if i = j
0 if i 6= j
(4)
COMP3121/9101 24 / 33
We now have two possibilities:
1 i = j: then
n−1∑
k=0
ω(i−j)kn =
n−1∑
k=0
ω0n =
n−1∑
k=0
1 = n;
2 i 6= j: then
∑n−1
k=0 ω
(i−j)k
n represents a sum of a geometric progression
with the ratio ωi−jn and thus
n−1∑
k=0
ω(i−j)kn =
1− ω(i−j)nn
1− ωi−jn
=
1− (ωnn)i−j
1− ωi−jn
=
1− 1
1− ωi−jn
= 0
So,
(
1 ωin ω
2·i
n . . . ω
i·(n−1)
n
)
1
ω−jn
ω−2jn
...
ω
−(n−1)j
n
=
n−1∑
k=0
ω(i−j)kn =
{
n if i = j
0 if i 6= j
(4)
COMP3121/9101 24 / 33
So we get:
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
1 1 1 . . . 1
1 ω−1n ω
−2
n . . . ω
−(n−1)
n
1 ω−2n ω
−2·2
n . . . ω
−2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ω
−(n−1)
n ω
−2(n−1)
n . . . ω
−(n−1)(n−1)
n
=
n 0 0 . . . 0
0 n 0 . . . 0
0 0 n . . . 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 0 . . . n
= n
1 0 0 . . . 0
0 1 0 . . . 0
0 0 1 . . . 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 0 . . . 1
i.e.,
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
−1
=
1
n
1 1 1 . . . 1
1 ω−1n ω
−2
n . . . ω
−(n−1)
n
1 ω−2n ω
−2·2
n . . . ω
−2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ω
−(n−1)
n ω
−2(n−1)
n . . . ω
−(n−1)(n−1)
n
COMP3121/9101 25 / 33
We now have
A0
A1
A2
.
.
.
An−1
=
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
−1
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
=
=
1
n
1 1 1 . . . 1
1 ω−1n ω
−2
n . . . ω
−(n−1)
n
1 ω−2n ω
−2·2
n . . . ω
−2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ω−(n−1)n ω
−2(n−1)
n . . . ω
−(n−1)(n−1)
n
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
This means that to covert from the values
〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ω
n−1
n )〉
which we denoted by 〈Â0, Â1, Â2, . . . , Ân−1〉 back to the coefficient form
PA(x) = A0 + A1x + A2x
2 + An−1x
n−1
we can use the same FFT algorithm with the only change that:
1 the root of unity ωn is replaced by ωn = e
−i 2π
n ,
2 the resulting output values are divided by n.
COMP3121/9101 26 / 33
We now have
A0
A1
A2
.
.
.
An−1
=
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
−1
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
=
=
1
n
1 1 1 . . . 1
1 ω−1n ω
−2
n . . . ω
−(n−1)
n
1 ω−2n ω
−2·2
n . . . ω
−2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ω−(n−1)n ω
−2(n−1)
n . . . ω
−(n−1)(n−1)
n
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
This means that to covert from the values
〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ω
n−1
n )〉
which we denoted by 〈Â0, Â1, Â2, . . . , Ân−1〉 back to the coefficient form
PA(x) = A0 + A1x + A2x
2 + An−1x
n−1
we can use the same FFT algorithm with the only change that:
1 the root of unity ωn is replaced by ωn = e
−i 2π
n ,
2 the resulting output values are divided by n.
COMP3121/9101 26 / 33
We now have
A0
A1
A2
.
.
.
An−1
=
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
−1
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
=
=
1
n
1 1 1 . . . 1
1 ω−1n ω
−2
n . . . ω
−(n−1)
n
1 ω−2n ω
−2·2
n . . . ω
−2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ω−(n−1)n ω
−2(n−1)
n . . . ω
−(n−1)(n−1)
n
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
This means that to covert from the values
〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ω
n−1
n )〉
which we denoted by 〈Â0, Â1, Â2, . . . , Ân−1〉 back to the coefficient form
PA(x) = A0 + A1x + A2x
2 + An−1x
n−1
we can use the same FFT algorithm with the only change that:
1 the root of unity ωn is replaced by ωn = e
−i 2π
n ,
2 the resulting output values are divided by n.
COMP3121/9101 26 / 33
We now have
A0
A1
A2
.
.
.
An−1
=
1 1 1 . . . 1
1 ωn ω
2
n . . . ω
n−1
n
1 ω2n ω
2·2
n . . . ω
2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ωn−1n ω
2(n−1)
n . . . ω
(n−1)(n−1)
n
−1
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
=
=
1
n
1 1 1 . . . 1
1 ω−1n ω
−2
n . . . ω
−(n−1)
n
1 ω−2n ω
−2·2
n . . . ω
−2·(n−1)
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 ω−(n−1)n ω
−2(n−1)
n . . . ω
−(n−1)(n−1)
n
PA(1)
PA(ωn)
PA(ω
2
n)
.
.
.
PA(ω
n−1
n )
This means that to covert from the values
〈PA(1), PA(ωn), PA(ω2n), . . . , PA(ω
n−1
n )〉
which we denoted by 〈Â0, Â1, Â2, . . . , Ân−1〉 back to the coefficient form
PA(x) = A0 + A1x + A2x
2 + An−1x
n−1
we can use the same FFT algorithm with the only change that:
1 the root of unity ωn is replaced by ωn = e
−i 2π
n ,
2 the resulting output values are divided by n.
COMP3121/9101 26 / 33
Inverse Fast Fourier Transform (IFFT):
1: function IFFT∗(Â)
2: n← length(Â)
3: if n = 1 then return Â
4: else
5: Â[0] ← (Â0, Â2, . . . Ân−2);
6: Â[1] ← (Â1, Â3, . . . Ân−1);
7: y[0] ← IFFT ∗(Â[0]);
8: y[1] ← IFFT ∗(Â[1]);
9: ωn ← e−i
2π
n ; ⇐ different from FFT
10: ω ← 1;
11: for k = 0 to k = n/2− 1 do;
12: yk ← y
[0]
k
+ ω · y[1]
k
;
13: yn/2+k ← y
[0]
k
− ω · y[1]
k
14: ω ← ω · ωn;
15: end for
16: return y;
17: end if
18: end function
1: function IFFT(Â) ⇐ different from FFT
2: return IFFT ∗(Â)/length(Â)
3: end function
COMP3121/9101 27 / 33
Important note:
Computer science books take the forward DFT operation to be the evaluation of the
corresponding polynomial at all roots of unity ωkn = cos
2πk
n
+ ı̇ sin 2πk
n
and the InverseDFT
to be the evaluation of the polynomial at the complex conjugates of the roots of unity, i.e.,
at ω−kn = cos
2πk
n
− ı̇ sin 2πk
n
.
However, Electrical engineering books do it just opposite, the direct DFT evaluates the
polynomial at ω−kn and the InverseDFT at ω
k
n!
While for the purpose of multiplying polynomials both choices are equally good, the choice
made by the electrical engineers is much better for all other purposes. We will explain this
in the Advanced Algorithms 4121 when we do the JPEG.
We did here only multiplication of polynomials, and did not apply it to multiplication of
large integers. This is possible to do but one has to be careful because roots of unity are
represented by floating point numbers so you have to show that if you do FFT with
sufficient precision you can round off the results and obtain correct integer values, but all of
this is tricky.
Earlier results along this line produced algorithms for multiplication of large integers which
operate in time n logn log(logn) but very recently David Harvey of the School of
Mathematics at UNSW came up with an algorithm to multiply large integers which runs in
time n logn.
COMP3121/9101 28 / 33
Important note:
Computer science books take the forward DFT operation to be the evaluation of the
corresponding polynomial at all roots of unity ωkn = cos
2πk
n
+ ı̇ sin 2πk
n
and the InverseDFT
to be the evaluation of the polynomial at the complex conjugates of the roots of unity, i.e.,
at ω−kn = cos
2πk
n
− ı̇ sin 2πk
n
.
However, Electrical engineering books do it just opposite, the direct DFT evaluates the
polynomial at ω−kn and the InverseDFT at ω
k
n!
While for the purpose of multiplying polynomials both choices are equally good, the choice
made by the electrical engineers is much better for all other purposes. We will explain this
in the Advanced Algorithms 4121 when we do the JPEG.
We did here only multiplication of polynomials, and did not apply it to multiplication of
large integers. This is possible to do but one has to be careful because roots of unity are
represented by floating point numbers so you have to show that if you do FFT with
sufficient precision you can round off the results and obtain correct integer values, but all of
this is tricky.
Earlier results along this line produced algorithms for multiplication of large integers which
operate in time n logn log(logn) but very recently David Harvey of the School of
Mathematics at UNSW came up with an algorithm to multiply large integers which runs in
time n logn.
COMP3121/9101 28 / 33
Important note:
Computer science books take the forward DFT operation to be the evaluation of the
corresponding polynomial at all roots of unity ωkn = cos
2πk
n
+ ı̇ sin 2πk
n
and the InverseDFT
to be the evaluation of the polynomial at the complex conjugates of the roots of unity, i.e.,
at ω−kn = cos
2πk
n
− ı̇ sin 2πk
n
.
However, Electrical engineering books do it just opposite, the direct DFT evaluates the
polynomial at ω−kn and the InverseDFT at ω
k
n!
While for the purpose of multiplying polynomials both choices are equally good, the choice
made by the electrical engineers is much better for all other purposes. We will explain this
in the Advanced Algorithms 4121 when we do the JPEG.
We did here only multiplication of polynomials, and did not apply it to multiplication of
large integers. This is possible to do but one has to be careful because roots of unity are
represented by floating point numbers so you have to show that if you do FFT with
sufficient precision you can round off the results and obtain correct integer values, but all of
this is tricky.
Earlier results along this line produced algorithms for multiplication of large integers which
operate in time n logn log(logn) but very recently David Harvey of the School of
Mathematics at UNSW came up with an algorithm to multiply large integers which runs in
time n logn.
COMP3121/9101 28 / 33
Important note:
Computer science books take the forward DFT operation to be the evaluation of the
corresponding polynomial at all roots of unity ωkn = cos
2πk
n
+ ı̇ sin 2πk
n
and the InverseDFT
to be the evaluation of the polynomial at the complex conjugates of the roots of unity, i.e.,
at ω−kn = cos
2πk
n
− ı̇ sin 2πk
n
.
However, Electrical engineering books do it just opposite, the direct DFT evaluates the
polynomial at ω−kn and the InverseDFT at ω
k
n!
While for the purpose of multiplying polynomials both choices are equally good, the choice
made by the electrical engineers is much better for all other purposes. We will explain this
in the Advanced Algorithms 4121 when we do the JPEG.
We did here only multiplication of polynomials, and did not apply it to multiplication of
large integers. This is possible to do but one has to be careful because roots of unity are
represented by floating point numbers so you have to show that if you do FFT with
sufficient precision you can round off the results and obtain correct integer values, but all of
this is tricky.
Earlier results along this line produced algorithms for multiplication of large integers which
operate in time n logn log(logn) but very recently David Harvey of the School of
Mathematics at UNSW came up with an algorithm to multiply large integers which runs in
time n logn.
COMP3121/9101 28 / 33
Important note:
Computer science books take the forward DFT operation to be the evaluation of the
corresponding polynomial at all roots of unity ωkn = cos
2πk
n
+ ı̇ sin 2πk
n
and the InverseDFT
to be the evaluation of the polynomial at the complex conjugates of the roots of unity, i.e.,
at ω−kn = cos
2πk
n
− ı̇ sin 2πk
n
.
However, Electrical engineering books do it just opposite, the direct DFT evaluates the
polynomial at ω−kn and the InverseDFT at ω
k
n!
While for the purpose of multiplying polynomials both choices are equally good, the choice
made by the electrical engineers is much better for all other purposes. We will explain this
in the Advanced Algorithms 4121 when we do the JPEG.
We did here only multiplication of polynomials, and did not apply it to multiplication of
large integers. This is possible to do but one has to be careful because roots of unity are
represented by floating point numbers so you have to show that if you do FFT with
sufficient precision you can round off the results and obtain correct integer values, but all of
this is tricky.
Earlier results along this line produced algorithms for multiplication of large integers which
operate in time n logn log(logn) but very recently David Harvey of the School of
Mathematics at UNSW came up with an algorithm to multiply large integers which runs in
time n logn.
COMP3121/9101 28 / 33
Back to fast multiplication of polynomials
PA(x) = A0 + A1x + . . . + An−1x
n−1
PB(x) = B0 + B1x + . . . + Bn−1x
n−1
⇓ DFT O(n log n) ⇓ DFT O(n log n)
{PA(1), PA(ω2n−1), PA(ω
2
2n−1), . . . , PA(ω
2n−2
2n−1)}; {PB(1), PB(ω2n−1), PB(ω
2
2n−1), . . . , PB(ω
2n−2
2n−1)}
⇓ multiplication O(n)
{PA(1)PB(1), PA(ω2n−1)PB(ω2n−1), . . . , PA(ω
2n−2
2n−1)PB(ω
2n−2
2n−1)}
⇓ IDFT O(n log n)
PC(x) =
2n−2∑
j=0
( j∑
i=0
AiBj−i
︸ ︷︷ ︸
Cj
)
x
j
=
2n−2∑
j=0
Cjx
j
= PA(x) · PB(x)
Thus, the product PC(x) = PA(x)PB(x) of two polynomials PA(x) and PB(x) can
be computed in time O(n logn).
COMP3121/9101 29 / 33
Back to fast multiplication of polynomials
PA(x) = A0 + A1x + . . . + An−1x
n−1
PB(x) = B0 + B1x + . . . + Bn−1x
n−1
⇓ DFT O(n log n) ⇓ DFT O(n log n)
{PA(1), PA(ω2n−1), PA(ω
2
2n−1), . . . , PA(ω
2n−2
2n−1)}; {PB(1), PB(ω2n−1), PB(ω
2
2n−1), . . . , PB(ω
2n−2
2n−1)}
⇓ multiplication O(n)
{PA(1)PB(1), PA(ω2n−1)PB(ω2n−1), . . . , PA(ω
2n−2
2n−1)PB(ω
2n−2
2n−1)}
⇓ IDFT O(n log n)
PC(x) =
2n−2∑
j=0
( j∑
i=0
AiBj−i
︸ ︷︷ ︸
Cj
)
x
j
=
2n−2∑
j=0
Cjx
j
= PA(x) · PB(x)
Thus, the product PC(x) = PA(x)PB(x) of two polynomials PA(x) and PB(x) can
be computed in time O(n logn).
COMP3121/9101 29 / 33
Computing the convolution C = A ∗B
A = 〈A0, A1, . . . , An−1〉 B = 〈B0, B1, . . . , Bn−1〉
⇓ O(n) ⇓ O(n)
PA(x) = A0 + A1x + . . . + An−1x
n−1
PB(x) = B0 + B1x + . . . + Bn−1x
n−1
⇓ DFT O(n log n) ⇓ DFT O(n log n)
{PA(1), PA(ω2n−1), PA(ω
2
2n−1), . . . , PA(ω
2n−2
2n−1)}; {PB(1), PB(ω2n−1), PB(ω
2
2n−1), . . . , PB(ω
2n−2
2n−1)}
⇓ multiplication O(n)
{PA(1)PB(1), PA(ω2n−1)PB(ω2n−1), . . . , PA(ω
2n−2
2n−1)PB(ω
2n−2
2n−1)}
⇓ IDFT O(n log n)
PC(x) =
2n−2∑
j=0
( j∑
i=0
AiBj−i
︸ ︷︷ ︸
Cj
)
x
j
⇓
C =
〈 j∑
i=0
AiBj−i
〉j=2n−2
j=0
Convolution C = A ∗B of sequences A and B is computed in time O(n logn).
COMP3121/9101 30 / 33
Computing the convolution C = A ∗B
A = 〈A0, A1, . . . , An−1〉 B = 〈B0, B1, . . . , Bn−1〉
⇓ O(n) ⇓ O(n)
PA(x) = A0 + A1x + . . . + An−1x
n−1
PB(x) = B0 + B1x + . . . + Bn−1x
n−1
⇓ DFT O(n log n) ⇓ DFT O(n log n)
{PA(1), PA(ω2n−1), PA(ω
2
2n−1), . . . , PA(ω
2n−2
2n−1)}; {PB(1), PB(ω2n−1), PB(ω
2
2n−1), . . . , PB(ω
2n−2
2n−1)}
⇓ multiplication O(n)
{PA(1)PB(1), PA(ω2n−1)PB(ω2n−1), . . . , PA(ω
2n−2
2n−1)PB(ω
2n−2
2n−1)}
⇓ IDFT O(n log n)
PC(x) =
2n−2∑
j=0
( j∑
i=0
AiBj−i
︸ ︷︷ ︸
Cj
)
x
j
⇓
C =
〈 j∑
i=0
AiBj−i
〉j=2n−2
j=0
Convolution C = A ∗B of sequences A and B is computed in time O(n logn).
COMP3121/9101 30 / 33
Visualizing Convolution C = A ∗B
C0 = A0B0
C1 = A0B1 + A1B0
Ck+m = AkBm + Ak+1Bm-1 + … +Ak+mB0
Cm+n = AnBm
0 0 0 0 A0 A1 A2 A3 A4 A5 A6 … Ak Ak+1 Ak+2 … Ak+m … … An-1 An 0 0 0 0
Bm Bm-1 … B1 B0
0 0 0 0 A0 A1 A2 A3 A4 A5 A6 … Ak Ak+1 Ak+2 … Ak+m An-1 An 0 0 0 0
Bm Bm-1 … B1 B0
0 0 0 0 A0 A1 A2 A3 A4 A5 A6 … Ak Ak+1 … Ak+m -1 Ak+m … … An-1 An 0 0 0 0
Bm Bm-1 … B1 B0
0 0 0 0 A0 A1 A2 A3 A4 A5 A6 … Ak Ak+1 Ak+2 … Ak+m -1 Ak+m … An-1 An 0 0 0 0
Bm Bm-1 … B1 B0
COMP3121/9101 31 / 33
An Exercise
Assume you are given a map of a straight sea shore of length 100n
meters as a sequence on 100n numbers such that Ai is the number of fish
between ith meter of the shore and (i + 1)th meter, 0 ≤ i ≤ 100n− 1.
You also have a net of length n meters but unfortunately it has holes in
it. Such a net is described as a sequence N of n ones and zeros, where
0’s denote where the holes are. If you throw such a net starting at meter
k and ending at meter k + n, then you will catch only the fish in one
meter stretches of the shore where the corresponding bit of the net is 1;
see the figure.
0 1 2 3 4 5 k k+1 . . . k+m n-1 n
C = Ak +Ak+1+0+ Ak+2+ 0 + Ak+4+…+0+ Ak+m-2+Ak+m-1
0 0 0 0 A0 A1 A2 A3 A4 . . . Ak Ak+1 Ak+2 Ak+3 Ak+4 … Ak+m—2 Ak+m -1 … An-2 An-1 0 0 0 0
1 1 0 1 0 1 1 0 1 1
Find the spot where you should place the left end of your net in order to
catch the largest possible number of fish using an algorithm which runs
in time O(n log n).
Hint: Let N ′ be the net sequence N in the reverse order; Compute A ∗B′
and look for the peak of that sequence.
COMP3121/9101 32 / 33
PUZZLE!!
On a circular highway there are n petrol stations, unevenly spaced, each
containing a different quantity of petrol. It is known that the total
quantity of petrol on all stations is enough to go around the highway
once, and that the tank of your car can hold enough fuel to make a trip
around the highway. Prove that there always exists a station among all
of the stations on the highway, such that if you take it as a starting point
and take the fuel from that station, you can continue to make a complete
round trip around the highway, never emptying your tank before
reaching the next station to refuel.
COMP3121/9101 33 / 33