Complex roots of unity

A complex nth root of unity is a complex number w such that

wⁿ = 1.

There are exactly n complex nth roots of unity: e^πik/n for k = 0, 1,..., n - 1. To interpret this formula, we use the definition of the exponential of a complex number:

e^iu = cos(u) + i sin(u).

Figure 30.2 shows that 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

Figure 30.2: The values of in the complex plane, where w₈ = e^πi/8 is the principal 8th root of unity.

(30.6)

is called the principal nth root of unity; all of the other complex nth roots of unity are powers of w_n.^[2]

The n complex nth roots of unity,

form a group under multiplication (see Section 31.3). This group has the same structure as the additive group (Z_n, +) modulo n, since implies that . Similarly, . Essential properties of the complex nth roots of unity are given in the following lemmas.

Lemma 30.3: (Cancellation lemma)

For any integers n ≥ 0, k ≥ 0, and d > 0,

(30.7)

Proof The lemma follows directly from equation (30.6), since

Corollary 30.4

For any even integer n > 0,

Proof The proof is left as Exercise 30.2-1.

Lemma 30.5: (Halving lemma)

If n > 0 is even, then the squares of the n complex nth roots of unity are the n/2 complex (n/2)th roots of unity.

Proof By the cancellation lemma, we have , for any nonnegative integer k. Note that if we square all of the complex nth roots of unity, then each (n/2)th root of unity is obtained exactly twice, since

Thus, and have the same square. This property can also be proved using Corollary 30.4, since implies , and thus .

As we shall see, the halving lemma is essential to our divide-and-conquer approach for converting between coefficient and point-value representations of polynomials, since it guarantees that the recursive subproblems are only half as large.

Lemma 30.6: (Summation lemma)

For any integer n ≥ 1 and nonnegative integer k not divisible by n,

Proof Equation (A.5) applies to complex values as well as to reals, and so we have

Requiring that k not be divisible by n ensures that the denominator is not 0, since only when k is divisible by n.

The DFT

Recall that we wish to evaluate a polynomial

of degree-bound n at (that is, at the n complex nth roots of unity).^[3] Without loss of generality, 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.^[4] We assume that A is given in coefficient form: a = (a₀, a₁,..., a_n-1). Let us define the results y_k, for k = 0, 1,..., n - 1, by

(30.8)

The vector y = (y₀, y₁,..., y_n-1) is the Discrete Fourier Transform (DFT) of the coefficient vector a = (a₀, a₁,..., a_n-1). We also write y = DFT_n(a).

The 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 DFT_n(a)in time Θ(n lg n), as opposed to the Θ(n²) 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 polynomials A^[0](x) and A^[1](x) of degree-bound n/2:

A^[0](x) = a₀ + a₂x + a₄x² + ··· + a_n-2x^n/2-1,

A^[1](x) = a₁ + a₃x + a₅x² + ··· + a_n-1x^n/2-1.

Note that A^[0] contains all the even-index coefficients of A (the binary representation of the index ends in 0) and A^[1] contains all the odd-index coefficients (the binary representation of the index ends in 1). It follows that

(30.9)

so that the problem of evaluating A(x) at reduces to

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

   (30.10)

   and then

2. combining the results according to equation (30.9).

By the halving lemma, the list of values (30.10) 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] and A^[1] of degree-bound n/2 are recursively evaluated at the n/2 complex (n/2)th roots of unity. These subproblems have exactly the same form as the original problem, but are half the size. We have now successfully divided an n-element DFT_n computation into two n/2-element DFT_n/2 computations. This decomposition is the basis for the following recursive FFT algorithm, which computes the DFT of an n-element vector a = (a₀, a₁,..., a_n-1), where n is a power of 2.

RECURSIVE-FFT(a)
 1  n ← length[a]     ▹ n is a power of 2.
 2  if n = 1
 3     then return a
 4  wn ← eπi/n
 5  w ← 1
 6  a[0] ← (a0, a2,..., an-2)
 7  a[1] ← (a1, a3,..., an-1)
 8  y[0] ← RECURSIVE-FFT(a[0])
 9  y[1] RECURSIVE-FFT(a[1])
10  for k ← 0 to n/2 - 1
11       do 
12          
13          w ← w wn
14  return y       ▹ y is assumed to be a column vector.

The RECURSIVE-FFT procedure works as follows. Lines 2 -3 represent the basis of the recursion; the DFT of one element is the element itself, since in this case

Lines 6 -7 define the coefficient vectors for the polynomials A^[0] and A^[1]. Lines 4, 5, and 13 guarantee that w is updated properly so that whenever lines 11 -12 are executed, . (Keeping a running value of w from iteration to iteration saves time over computing from scratch each time through the for loop.) Lines 8-9 perform the recursive DFT_n/2 computations, setting, for k = 0, 1,..., n/2 - 1,

or, since by the cancellation lemma,

Lines 11-12 combine the results of the recursive DFT_n/2 calculations. For y₀, y₁,..., y_n/2-1, line 11 yields

For y_n/2, y_n/2+1,..., y_n-1, letting k = 0, 1,..., n/2 - 1, line 12 yields

Thus, the vector y returned by RECURSIVE-FFT is indeed the DFT of the input vector a.

Within the for loop of lines 10 -13, each value is multiplied by , for k = 0, 1,..., n/2 - 1. The product is both added to and subtracted from . Because each factor is used in both its positive and negative forms, the factors are known as twiddle factors.

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

Thus, we can evaluate a polynomial of degree-bound n at the complex nth roots of unity in time Θ(n lg n) using the Fast Fourier Transform.

Interpolation at the complex roots of unity

We now complete the polynomial multiplication scheme by showing how to interpolate the complex roots of unity by a polynomial, which enables us to convert from point-value form back to coefficient form. We interpolate by writing the DFT as a matrix equation and then looking at the form of the matrix inverse.

From equation (30.4), we can write the DFT as the matrix product y = V_na, where V_n is a Vandermonde matrix containing the appropriate powers of w_n:

The (k, j) entry of V_n is , for j, k = 0, 1,..., n - 1, and the exponents of the entries of V_n form a multiplication table.

For the inverse operation, which we write as , we proceed by multiplying y by the matrix , the inverse of V_n.

Theorem 30.7

For j, k = 0, 1,..., n - 1, the (j, k) entry of is .

Proof We show that , the n × n identity matrix. Consider the (j, j') entry of :

This summation equals 1 if j' = j, and it is 0 otherwise by the summation lemma (Lemma 30.6). Note that we rely on -(n - 1) < j' - j < n - 1, so that j' - j is not divisible by n, in order for the summation lemma to apply.

Given the inverse matrix , we have that is given by

(30.11)

for j = 0, 1,..., n - 1. By comparing equations (30.8) and (30.11), we see that by modifying the FFT algorithm to switch the roles of a and y, replace w_n by , and divide each element of the result by n, we compute the inverse DFT (see Exercise 30.2-4). Thus, can be computed in Θ(n lg n) time as well.

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 Θ(n lg n). In the context of polynomial multiplication, we have shown the following.

Theorem 30.8: (Convolution theorem)

For any two vectors a and b of length n, where n is a power of 2,

where the vectors a and b are padded with 0's to length 2n and · denotes the componentwise product of two 2n-element vectors.

Exercises 30.2-1

Prove Corollary 30.4.

Exercises 30.2-2

Compute the DFT of the vector (0, 1, 2, 3).

Exercises 30.2-3

Do Exercise 30.1-1 by using the Θ(n lg n)-time scheme.

Exercises 30.2-4

Write pseudocode to compute in Θ(n lg n) time.

Exercises 30.2-5

Describe the generalization of the FFT procedure to the case in which n is a power of 3. Give a recurrence for the running time, and solve the recurrence.

Exercises 30.2-6: ⋆

Suppose that instead of performing an n-element FFT over the field of complex numbers (where n is even), we use the ring Z_m of integers modulo m, where m = 2^tn/2 + 1 and t is an arbitrary positive integer. Use w = 2^t instead of w_n as a principal nth root of unity, modulo m. Prove that the DFT and the inverse DFT are well defined in this system.

Exercises 30.2-7

Given a list of values z₀, z₁,..., z_n-1 (possibly with repetitions), show how to find the coefficients of a polynomial P(x) of degree-bound n + 1 that has zeros only at z₀, z₁,..., z_n-1 (possibly with repetitions). Your procedure should run in time O(n lg² n). (Hint: The polynomial P(x) has a zero at z_j if and only if P(x) is a multiple of (x - z_j).)

Exercises 30.2-8: ⋆

The chirp transform of a vector a = (a₀, a₁,..., a_n-1) is the vector y = (y₀, y₁,..., y_n-1), where and z is any complex number. The DFT is therefore a special case of the chirp transform, obtained by taking z = w_n. Prove that the chirp transform can be evaluated in time O(n lg n) for any complex number z. (Hint: Use the equation

to view the chirp transform as a convolution.)