Polynomials and the FFT

Chapter 32 in CLR
Chapter 30 in CLRS
Goal

Fact 1: Straightforward method of adding two polynomials of degree n takes θ(n)

time.

Fact 2: Straightforward method of multiplying two polynomials of degree n takes

θ(n2) time.

Goal: Show how the Fast Fourier Transform (FFT) can reduce the time to
multiply polynomials to θ(nlogn).

2
Polynomials

A polynomial in the variable x over the field F is a function A(x) that can be
represented as
A(x) = ∑j=0,..,n-1 aj xj.

We call n a degree-bound of the polynomial, and the values a0,a1,*,an-1 are the
coefficients of the polynomial (they are drawn from the field F). A polynomial is
said to have degree k if its highest nonzero coefficient is ak.

3
Polynomial addition

If A(x) and B(x) are polynomials of degree bound n, we say that their sum is a
polynomial C(x), also of degree bound n, such that C(x) = A(x) + B(x) for all x in F.
That is, if
A(x) = ∑j=0,..,n-1 aj xj
and
B(x) = ∑j=0,..,n-1 bj xj
then
C(x) = ∑j=0,..,n-1 cj xj,
where cj = aj + bj for j = 0,1,*,n-1.

Clearly, we can calculate C(x) in θ(n) time.

4
Polynomial multiplication

If A(x) and B(x) are polynomials of degree bound n, we say that their product is a
polynomial C(x) (of degree bound 2n-1, but we may also say that 2n is a degree
bound), such that C(x) = A(x) B(x) for all x in F. That is, if
A(x) = ∑j=0,..,n-1 aj xj
and
B(x) = ∑j=0,..,n-1 bj xj
then
C(x) = ∑j=0,..,2n-2 cj xj,
where cj = ∑k=0,..,j ak bj-k for j = 0,1,*,2n-2. The resulting coefficient vector c is
also called the convolution of the input vectors a and b, and is denoted by
c=a ⊗ b.
Clearly, we can calculate C(x) in θ(n2) time.

5
32.1: Representation of polynomials

There are two common ways to represent polynomials:

• A coefficient representation of a polynomial A(x) = ∑j=0,..,n-1 aj xj of degree

bound n is the vector of coefficients (a0,a1,*,an-1).

Using this representation we easily multiply two polynomials of degree bound n in

O(n2) time.

• A point-value representation of a polynomial A(x) of degree bound n is a set of

n point-value pairs {(x0,y0),(x1,y1),*,(xn-1,yn-1)}, such that all the xk –s are distinct
and A(xk ) = yk for k = 0,1,*,n-1.

A polynomial has only one coefficient representation (given n), and many point-
value representations.

6
32.1: Representation of polynomials

Computing a point-value representation of a polynomial given its coefficient form

is straightforward, since all we have to do is to select n distinct points x0,x1,*,xn-1
and evaluate A(xk) for k = 0,1,*,n-1. This takes θ(n2) time “regularly”. Later we
shall see, how one can do this in θ(nlogn) time, if the xk –s are selected cleverly.

The inverse of evaluation – determining the coefficients form of a polynomial from

a point-value representation – is called interpolation. The following theorem, due
to Lagrange, shows that interpolation is well defined.

7
32.1: Representation of polynomials

Theorem 32.1 (Uniqueness of an interpolating polynomial):

For any set {(x0,y0),(x1,y1),*,(xn-1,yn-1)} of n point-value pairs, there is a unique

polynomial A(x) of degree bound n such that yk = A(xk) for k = 0,1,*,n-1.

8
32.1: Representation of polynomials

Proof:
The set of equations {yk = A(xk)} k=0,1,*,n-1 may be written in matrix form as

9
32.1: Representation of polynomials

and the existence of its solution is equivalent to the existence of the inverse of the
matrix of coefficients. This matrix on the left, denoted by V(x0,x1,*,xn-1) is the
Vandermonde matrix whose determinant is ∏j<k(xk-xj), and therefore it is invertible
if the xk –s are distinct. Therefore, the vector of coefficients, a, can be solved for
uniquely given the point value representation:
a = V(x0,x1,*,xn-1) -1 y. █

10
32.1: Representation of polynomials

The proof describes an algorithm via a solution of a set of linear equations. This
may be done in time O(n3).

A faster algorithm for n-point interpolation is based on Lagrange’s formula:

A(x) = ∑ k=0,1,*,n-1 yk (∏ j≠k (x-xj) / ∏ j≠k (xk-xj)).

An easy exercise shows that using this formula, one may calculate the
coefficients in time θ(n2).

11
32.1: Representation of polynomials

The point-value representation is quite convenient for many operations on

polynomials. For addition, if C(x) = A(x) + B(x), then C(xk) = A(xk) + B(xk), for any
point xk.

More precisely, if we have a point-value representation for A, {(x0,y0),(x1,y1),*,

(xn-1,yn-1)}, and for B, {(x0,y’0),(x1,y’1),*,(xn-1,y’n-1)} (same xk –s!), then a point-value
representation for C is {(x0,y0+y’0),(x1,y1+y’1),*,(xn-1,yn-1+y’n-1)}, and it is
calculated in time θ(n).

12
32.1: Representation of polynomials

Similarly, for multiplication, if C(x) = A(x)B(x), then C(xk) = A(xk)B(xk), for any point
xk. However, this time, n points will not suffice, since the degree bound for C is
now 2n. So now we need to represent A via {(x0,y0),(x1,y1),*,(x2n-1,y2n-1)}, and B
via {(x0,y’0),(x1,y’1),*,(x2n-1,y’2n-1)} (again, same xk–s!), and the point-value
representation for C is {(x0,y0y’0),(x1,y1y’1),*,(x2n-1,y2n-1y’2n-1)} (calculated in time
θ(n)).

13
32.1: Representation of polynomials

Can we use the linear-time multiplication method for polynomials in point-value

form to expedite polynomial multiplication in coefficient form? The answer hinges
on our ability to convert a polynomial quickly from one form to another (and back).
A key observation is that we may choose the evaluation points in any useful way.
It turns out that if we use “complex roots of unity” as evaluation points, we can
produce a point-value by taking the Discrete Fourier Transform (DFT) of a
coefficient vector. The inverse operation turns out to be the “inverse DFT”,
yielding the coefficient vector. FFT performs both the DFT and the inverse DFT in
time θ(nlogn).

The following figure explains it all.

14
32.1: Representation of polynomials

15
32.1: Representation of polynomials

Given the FFT, we have the following θ(nlogn)-time procedure for multiplying two
polynomials A(x) and B(x) of degree-bound n, where the input and output
representations are in coefficient form. We assume that n is a power of 2; this
requirement can always be met by adding high-order zero coefficients.

Step 1: Double degree-bound: Create coefficient representation of A(x) and

B(x) as degree-bound 2n polynomials by adding n high-order zero coefficients to
each.
Step 2: Evaluate: Compute point-value representations of A(x) and B(x) of length
2n through two applications of FFT of order 2n. These representations contains
the values of the two polynomials at the 2n-th roots of unity.

16
32.1: Representation of polynomials

Step 3: Pointwise multiply: Compute a point-value representation for the

polynomial C(x) = A(x)B(x) by multiplying these values together pointwise. This
representation contains the value C(x) at each 2n-th root of unity.
Step 4: Interpolate: Create the coefficient representation of the polynomial C(x)
through a single application of an FFT on 2n point-value pairs to compute the
inverse DFT.

Steps (1) & (3) take time θ(n), and steps (2) & (4) take time θ(nlogn). Therefore,
once we show how to use FFT, we will have proven the following theorem:

17
32.1: Representation of polynomials

Theorem 32.2:

The product of two polynomials of degree bound n can be computed in time

θ(nlogn), with both input and output representations in coefficient forms.

18
32.2: The DFT and FFT – mathematical background

Complex roots of unity

A complex n-th root of unity is a complex number ω such that ωn=1.
There are exactly n complex n-th roots of unity - these are e2πik/n for
k=0,1,*,n-1. We may interpret this formula using a different representation
of complex numbers:
eiu = cos(u) + i sin(u).
The figure in the next slide shows the geometric picture: the n complex roots
of unity are equally spaced around the circle of unit radius centered at the
origin of the complex plane.
The value ωn = e2πi/n is called the principal n-th root of unity; all of the
other n-th roots of unity are powers of ωn.

19
32.2: The DFT and FFT – mathematical background

The n complex n-th roots of unity ωn0,ωn1, ,ωnn-1, form a multiplicative

group. This group is isomorphic to the additive group (Zn,+) modulo n, since
ωnn = ωn0 =1 implies ωnj ωnk = ωn(j+k) mod n. Similarly ωn-1 = ωnn-1.

Essential properties of the complex n-th roots of unity are given in the
following lemmas.

20
32.2: The DFT and FFT – mathematical background

Lemma 32.3 (Cancellation lemma): For any integers n ≥ 0, k ≥ 0, and

d > 0, ωdndk = ωnk.
Corollary 32.4: For any even integer n > 0, ωnn/2 = ω2 = -1.

Lemma 32.5 (Halving lemma): If n > 0 is even, then the squares of the n
complex n-th roots of unity are the n/2 complex (n/2)-th roots of unity. [This is
essential for the Divide and Conquer approach of FFT.]

Lemma 32.6 (Summation lemma): For any integer n > 1 and nonnegative
integer k not divisible by n,
∑j=0,..,n-1 (ωnk)j = 0.

21
32.2: The DFT and FFT

We wish to evaluate the polynomial A(x) = ∑j=0,..,n-1 aj xj of degree bound n* at

ωn0, ωn1, , ωnn-1, the n complex n-th roots of unity. Wlog we assume that n
is a power of 2, since a given degree-bound can always be raised – we can
always add new high-order zero coefficients as necessary. We assume that
A is given in coefficient form: a = (a0,a1,*,an-1).
Let us define the results yk, for k=0,..,n-1, by
yk = A(ωnk) = ∑j=0,..,n-1 aj ωnkj.

The vector y =(y0,y1,*,yn-1) is called the Discrete Fourier Transform (DFT)

of the coefficient vector a = (a0,a1,*,an-1). We also write y = DFTn(a).

*Actually it is 2n, since we double the degree bound. 22

32.2: The DFT and FFT

By using a method known as the Fast Fourier Transform (FFT), which

takes advantage of the special properties of the complex roots of unity. We
can compute DFTn(a) in time θ(nlogn), as opposed to the θ(n2) time of the
straightforward method.
The FFT method employs a divide-and-conquer strategy, using the even-
index and odd-index coefficients of A(x) separately to define the two new
degree-bound n/2 polynomials A[0](x) and A[1](x):
A[0](x) = a0 + a2x + a4x2 + * + an-2xn/2-1,
A[1](x) = a1 + a3x + a5x2 + * + an-1xn/2-1.
It follows that
A(x) = A[0](x2) + x A[1](x2),

23
32.2: The DFT and FFT

So that the problem of evaluating A(x) at ωn0, ωn1, , ωnn-1 reduces to

1. evaluating the degree-bound n/2 polynomials A[0](x) and A[1](x) at the

points (ωn0)2, (ωn1)2, *, (ωnn-1)2, and then

2. combining the results according to A(x) = A[0](x2) + x A[1](x2).

24
32.2: The DFT and FFT

By the halving lemma, the list of values (ωn0)2, (ωn1)2, *, (ωnn-1)2 consists
not of n distinct values but only of the n/2 complex (n/2)-th roots of unity,
with each root occurring exactly twice. Therefore, the polynomials A[0](x) and
A[1](x) of degree-bound n/2 are recursively evaluated at the n/2 complex
(n/2)-th roots of unity. These two sub-problems have exactly the same form
of the original problem, but are half the size. We have now successfully
divided an n-element DFTn computation into two n/2-element DFTn/2
computations.

This decomposition is the basis for the following recursive FFT algorithm,
which computes the DFT of an n-element vector a = (a0,a1,*,an-1), where n
is a power of 2.
25
32.2: The DFT and FFT

RECURSIVE-FFT(a)
n ← length[a] ► n is a power of 2
if n = 1
then return a
ωn ← e2πi/n
ω←1
a[0] ← (a0,a2,*,an-2)
a[1] ← (a1,a3, *,an-1)
y[0] ← RECURSIVE-FFT(a[0])
y[1] ← RECURSIVE-FFT(a[1])
for k ← 0 to n/2 – 1
do yk ← yk[0] + ωyk[1]
yk+(n/2) ← yk[0] – ωyk[1]
ω ← ω ωn
return y ► y is assumed to be column vector

26
32.2: The DFT and FFT

The correctness of the algorithm is proved via easy calculations based on

the lemmas that were proved before.

To determine the running time of procedure RECURSIVE-FFT, we note that

exclusive of the recursive calls, each invocation takes time θ(n), where n is
the length of the input vector. The recurrence for the running time is
therefore

T(n) = 2T(n/2) + θ(n) = θ(nlogn).

Thus, we can evaluate a polynomial of degree-bound n at the complex n-th

roots of unity in time θ(nlogn) using the Fast Fourier Transform.

27
32.2: Interpolation at the complex roots of unity

To complete the polynomial multiplication scheme, we now show how to

interpolate the complex roots of unity by a polynomial. This enables us to
convert from point-value form back to coefficient form. We can easily
understand how to do it, if we write DFT as a matrix equation, and calculate
its inverse:

28
32.2: Interpolation at the complex roots of unity

So we have y = Vn a where Vn is the Vandermonde matrix containing the

appropriate powers ωn. The (k,j) entry of Vn is ωnkj for j, k = 0,1,*,n-1, and
the exponents of the entries of Vn form a multiplication table.
For the inverse operation, we write a = DFT-1n(y), and we calculate it by
using the inverse of the matrix Vn: a = Vn-1 y.

Theorem 32.7: For j, k = 0,1,*,n-1, The (j,k) entry of Vn-1 is ωn-kj /n.
Proof: Multiply.

Given the inverse matrix Vn-1, we have that DFT-1n(y) is given by

aj = 1/n ∑k=0,..,n-1 yk ωn-kj.

29
32.2: Interpolation at the complex roots of unity

We now compare:

aj = 1/n ∑k=0,..,n-1 yk ωn-kj yk = ∑j=0,..,n-1 aj ωnkj ,

and discover that by modifying the FFT algorithm to switch the roles of a
and y, and replace ωn by ωn-1 and divide each element of the result by n, we
compute the inverse DFT. Thus, DFT-1n can be computed in θ(nlogn) time as
well.

30
32.2: Interpolation at the complex roots of unity

Thus, by using the FFT and the inverse FFT, we can transform a polynomial
of degree bound n back and forth between its coefficient representation and
a point-value representation in time θ(nlogn). In the context of polynomial
multiplication, we have shown the following:

Theorem 32.7: For any two vectors a and b of length n, where n is a power
of 2,
a ⊗ b = DFT-12n(DFT2n(a) — DFT2n(b)),
where the vectors a and b are padded with 0’s to length 2n and “—” denotes
the componentwise product of two 2n-element vectors.

31
32.3: Implementation

32
What’s next?

Working over the complex numbers is problematic: the complex roots of

unity are (almost always) irrational in both their real and imaginary parts.
How does one handle them in the computer? Can we bypass this problem
in some applications?

Example: Suppose we wish to multiply 2 polynomials whose coefficients

are positive integers. We know that the resulting polynomial has the same
property. However, if we work through FFT, the rounding process will
probably not produce the exact results.

Idea: Work over Zm, for a suitable m.

33
Lemma: Let m = 2tn/2 + 1 where t is a positive integer, and n is an even
positive number. Then ω = 2t is an n-th root of unity in Zm, and we may use
it instead of ωn to perform FFT over Zm.

34