Polynomial + Fast Fourier Transform

Polynomial  +  

Fast  Fourier  Transform  

Michael  Tsai 2017/06/13

Polynomials

• Coefficients:

• Degree:  highest  order  term  with  nonzero  coefficient   (k  if  highest  nonzero  term  is )

• Degree-­‐bound:  any  integer  strictly  greater  than  the   degree

A(x) =

n 1X

j=0

a_jx^j

A(x) = a

+ a

x + a

x

+ · · · + a

x

a

, a

, . . . , a

a

Coefficient  Representation

• Using  it:

• Evaluation:

• Addition:

A(x) =

n 1X

j=0

a_jx^j

(a

, a

, . . . , a

)

Vector

A(x₀) = a₀ + x₀(a₁ + x₀(a₂ + · · · + x0(a^{n 2} + x₀(a_{n 1})) . . . )) Horner’s  rule:

(a

, a

, . . . , a

)

A(x) + B(x)

(b

, b

, . . . , b

)

(a₀ + b₀, a₁ + b₁, . . . , a_{n 1} + b_{n 1}) O(n)

O(n)

Prob.:  Polynomial  Multiplication

• Example: A(x) = 6x³ + 7x² 10x + 9

B(x) = 2x³ + 4x 5 C(x) = A(x)B(x)

Chapter 30 Polynomials and the FFT 899

mial C.x/, also of degree-bound n, such that C.x/ D A.x/ C B.x/ for all x in the underlying field. That is, if

A.x/ D

n!1

X

jD0

a_jx^j and

B.x/ D

n!1

X

jD0

bjx^j ; then

C.x/ D

n!1

X

jD0

cjx^j ;

where c^j D a^j C b^j for j D 0; 1; : : : ; n ! 1. For example, if we have the polynomials A.x/ D 6x³ C 7x² ! 10x C 9 and B.x/ D !2x³ C 4x ! 5, then C.x/ D 4x³ C 7x² ! 6x C 4.

For polynomial multiplication, if A.x/ and B.x/ are polynomials of degree- bound n, their product C.x/ is a polynomial of degree-bound 2n ! 1 such that C.x/ D A.x/B.x/ for all x in the underlying field. You probably have multi- plied polynomials before, by multiplying each term in A.x/ by each term in B.x/

and then combining terms with equal powers. For example, we can multiply A.x/ D 6x³ C 7x² ! 10x C 9 and B.x/ D !2x³ C 4x ! 5 as follows:

6x³ C 7x² ! 10x C 9

! 2x³ C 4x ! 5

! 30x³ ! 35x² C 50x ! 45 24x⁴ C 28x³ ! 40x² C 36x

! 12x⁶ ! 14x⁵ C 20x⁴ ! 18x³

! 12x⁶ ! 14x⁵ C 44x⁴ ! 20x³ ! 75x² C 86x ! 45 Another way to express the product C.x/ is

C.x/ D

2nX!2 jD0

c_jx^j ; (30.1)

where c_j D

Xj kD0

a_kb_j!k : (30.2)

O(n

)

O(n) O(n) O(n)^…

Prob.:  Polynomial  Multiplication

• Degree(C)=degree(A)+degree(B)

• Degree  bound:

• !

• Problem:  can  we  reduce  this  time  complexity?

C(x) =

2n 2X

j=0

c_jx^j c_j =

Xj

k=0

a_kb_{j k}

n_a + n_b

O(n

)

Point-­‐Value  Representation

• Computing  a  point-­‐value  representation:

• Calculate  each                                  using  Horner’s  rule  takes                

• Thus  calculate  all  n                              values  take

• If  we  choose  the  points                      wisely,  we  can  reduce   this  to                                            !

A(x) =

n 1X

j=0

a_jx^j

Degree-­‐bound:   n

{(x⁰, y₀), (x₁, y₁), . . . , (x_{n 1}, y_{n 1})}

y_k = A(x_k) k = 0, 1, . . . , n 1

A(x_k) O(n)

A(x_k) O(n²)

x

O(n log n)

Interpolation

• Interpolation is  the  inverse  operation  of  evaluation

• Interpolation is  well-­‐defined when  the  

interpolating  polynomial  have  a  degree-­‐bound ==  

given  number  of  point-­‐value  pairs

Theorem:    

For  any  set                                                                                                                    of  n  point-­‐value  

pairs  such  that  all  the                  values  are  distinct,  there  is  a  unique   polynomial  A(x)  of  degree-­‐bound  n  such  that                                                      for

.        (see  p.  902  in  Cormen for  the  proof)

{(x⁰, y₀), (x₁, y₁), . . . , (x_{n 1}, y_{n 1})}

x

k = 0, 1, . . . , n 1

y_k = A(x_k)

Point-­‐Value  Representation

• Using  it:

• Addition:

since                                                                          ,    

A(x) =

n 1X

j=0

a_jx^j

Degree-­‐bound:   n

{(x⁰, y0), (x1, y1), . . . , (xn 1, yn 1)}

y_k = A(x_k) k = 0, 1, . . . , n 1 Evaluation

Interpolation

C(x) = A(x) + B(x) C(x_k) = A(x_k) + B(x_k)

{(x⁰, y₀), (x₁, y₁), . . . , (x_{n 1}, y_{n 1})} {(x⁰, y₀⁰), (x₁, y₁⁰), . . . , (x_{n 1}, y_{n 1}⁰ )} A:

B:

{(x⁰, y0 + y₀⁰), (x1, y1 + y₁⁰ ), . . . , (xn 1, yn 1 + y_{n 1}⁰ )}

C: O(n)

Point-­‐Value  Representation:  

Multiplication

• Multiplication:

since                                                                            ,  

• Problem!  

We  need  2n  point-­‐value  pairs  so  that  C(x)  is  well-­‐

defined!

C(x) = A(x)B(x) C(x_k) = A(x_k)B(x_k)

A:

B:

C:

{(x⁰, y0), (x1, y1), . . . , (xn 1, yn 1)} {(x⁰, y₀⁰), (x₁, y₁⁰), . . . , (x_{n 1}, y_{n 1}⁰ )}

{(x⁰, y₀y₀⁰), (x₁, y₁y₁⁰), . . . , (x_{n 1}, y_{n 1}y_{n 1}⁰ )} n pairs

n pairs

n pairs

Point-­‐Value  Representation:  

Multiplication

• Multiplication:

since                                                                            ,  

• Solution:  Extend  A  and  B  to  2n  point-­‐value  pairs (add  n  zero  coefficient  high-­‐order  terms)

• C(x)  is  now  well-­‐defined  with  2n  point-­‐value  pairs C(x) = A(x)B(x) C(x_k) = A(x_k)B(x_k)

A:

B:

C:

{(x⁰, y₀y₀⁰ ), (x₁, y₁y₁⁰ ), . . . , (x_{2n 1}, y_{2n 1}y_{2n 1}⁰ )}

{(x⁰, y0), (x1, y1), . . . , (x2n 1, y2n 1)} {(x⁰, y₀⁰), (x₁, y₁⁰), . . . , (x_{2n 1}, y_{2n 1}⁰ )}

2n  pairs 2n  pairs

2n  pairs

O(n)

904 Chapter 30 Polynomials and the FFT

a0; a1; : : : ; an!1

b0; b1; : : : ; bn!1

c0; c1; : : : ; c2n!2

Ordinary multiplication Time ‚.n²/ Evaluation

Time ‚.n lg n/

Time ‚.n lg n/

Interpolation

Pointwise multiplication Time ‚.n/

A.!_2n⁰ /; B.!_2n⁰ / A.!_2n¹ /; B.!_2n¹ /

A.!_2n^2n!1/; B.!_2n^2n!1/

::: :::

C.!_2n⁰ / C.!_2n¹ /

C.!_2n^2n!1/

Coefficient

Point-value representations representations

Figure 30.1 A graphical outline of an efficient polynomial-multiplication process. Representations on the top are in coefficient form, while those on the bottom are in point-value form. The arrows from left to right correspond to the multiplication operation. The !²ⁿterms are complex .2n/th roots of unity.

on whether we can convert a polynomial quickly from coefficient form to point- value form (evaluate) and vice versa (interpolate).

We can use any points we want as evaluation points, but by choosing the eval- uation points carefully, we can convert between representations in only ‚.n lg n/

time. As we shall see in Section 30.2, if we choose “complex roots of unity” as the evaluation points, we can produce a point-value representation by taking the discrete Fourier transform (or DFT) of a coefficient vector. We can perform the inverse operation, interpolation, by taking the “inverse DFT” of point-value pairs, yielding a coefficient vector. Section 30.2 will show how the FFT accomplishes the DFT and inverse DFT operations in ‚.n lg n/ time.

Figure 30.1 shows this strategy graphically. One minor detail concerns degree- bounds. The product of two polynomials of degree-bound n is a polynomial of degree-bound 2n. Before evaluating the input polynomials A and B, therefore, we first double their degree-bounds to 2n by adding n high-order coefficients of 0.

Because the vectors have 2n elements, we use “complex .2n/th roots of unity,”

which are denoted by the !²ⁿ terms in Figure 30.1.

Given the FFT, we have the following ‚.n lg n/-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; we can always meet this requirement by adding high-order zero coefficients.

1. Double degree-bound: Create coefficient representations of A.x/ and B.x/ as degree-bound 2n polynomials by adding n high-order zero coefficients to each.

⇥(n²) ⇥(n²)

Can  we  improve  evaluation  and  interpolation  time  to

?  

⇥(n log n)

Complex  Roots  of  Unity

• A  complex  n-­‐th root  of  unity  is  a  complex  number   such  that

• Exponential  of  a  complex  number:

• There  are  exactly  n  complex  n-­‐th roots  of  unity:

!

= 1

!

e

= cos(u) + i sin(u)

e

Example

30.2 The DFT and FFT 907

!1 1

i

!i

!₈⁰ D !8⁸

!₈¹

!₈²

!₈³

!₈⁴

!₈⁵

!₈⁶

!₈⁷

Figure 30.2 The values of !8⁰; !₈¹; : : : ; !₈⁷ in the complex plane, where !⁸ D e^{2! i=8} is the prin- cipal 8th root of unity.

!n D e^{2! i=n} (30.6)

is the principal nth root of unity;² all other complex nth roots of unity are powers of !ⁿ.

The n complex nth roots of unity,

!_n⁰; !_n¹; : : : ; !_n^n!1 ;

form a group under multiplication (see Section 31.3). This group has the same structure as the additive group .Zⁿ; C/ modulo n, since !nⁿ D !n⁰ D 1 implies that

!_n^j!_n^k D !n^jCk D !n^{.jCk/ mod n}. Similarly, !n^!1 D !n^n!1. The following lemmas furnish some essential properties of the complex nth roots of unity.

Lemma 30.3 (Cancellation lemma)

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

!_{d n}^{d k} D !n^k : (30.7)

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

!_{d n}^{d k} D !

e^{2! i=d n}"^{d k}

D !

e^{2! i=n}"^k D !n^k :

2Many other authors define !ⁿ differently: !ⁿ D e^!2!i=n. This alternative definition tends to be used for signal-processing applications. The underlying mathematics is substantially the same with either definition of !ⁿ.

Principle  n-­‐th root  of  unity:

n  complex  n-­‐th roots  of  unity:

!

= e

!

, !

, . . . , !

Discrete  Fourier  Transform

• Evaluate  A(x)  of  degree-­‐bound  n  at

• The  vector  

is  the  discrete  Fourier  Transform  (DFT)  of   coefficient  vector    

A(x) =

n 1X

j=0

a_jx^j

!_n⁰, !_n¹, . . . , !_n^{n 1}

y_k = A(!_n^k) =

n 1X

j=0

a_j!_n^kj k = 0, 1, . . . , n 1

y = (y

, y

, . . . , y

)

a = (a

, a

, . . . , a

)

O(n

)

still

Physical  Meaning  of  DFT

y_k = A(!_n^k) =

n 1X

j=0

a_j!_n^kj

a_j = 1 n

n 1X

k=0

y_k!_n ^kj Inverse  DFT

Signal  at  different  frequencies   Weight  of  that

frequency

Fast  Fourier  Transform

• Taking  advantage  of  the  special  properties  of  the   complex  roots  of  unity,  we  can  compute  DFT  in

!

• Assumption:  n  is  a  power  of  2.

• Split  A(x)  into  two  parts:

⇥(n log n)

A(x) = a₀ + a₁x + a₂x² + a₃x³ + · · · + a_{n 2}x^{n 2} + a_{n 1}x^{n 1}

A^[0](x) = a₀ + a₂x + a₄x² + · · · + a_{n 2}x^{n/2 1}

A^[1](x) = a₁ + a₃x + a₅x² + · · · + a_{n 1}x^{n/2 1}

A(x) = A^[0](x²) + xA^[1](x²)

Fast  Fourier  Transform

• How  do  we  evaluate  A(x)  at                                                                              ? 1. Evaluate                                and                                  at

2. Combine  the  result  using  (ㄅ)

A^[0](x) = a₀ + a₂x + a₄x² + · · · + a_{n 2}x^{n/2 1}

A^[1](x) = a₁ + a₃x + a₅x² + · · · + a_{n 1}x^{n/2 1}

A(x) = A^[0](x²) + xA^[1](x²)

!_n⁰, !_n¹, . . . , !_n^{n 1}

A^[0](x) A^[1](x)

(!_n⁰)², (!_n¹)², . . . , (!_{n 1}⁰ )²

(ㄅ)

What  is  Divide-­‐and-­‐Conquer?

• When  dealing  with  a  problem:

1. Divide  the  problem  into  

smaller,  but  same  type  of,  problems

2. If  the  problem  is  small  enough  to  solve  (Conquer),

• then  solve  it

• Else  recursively  call  itself  to  solve  smaller  sub-­‐problems

3. Combine  the  solutions  of  smaller  sub-­‐problems  into   the  solution  of  the  original,  larger,  problem

18

Base  case Recursive  case

Fast  Fourier  Transform

• How  do  we  evaluate  A(x)  at                                                                              ? 1. Evaluate                                and                                  at

2. Combine  the  result  using  (ㄅ)

A^[0](x) = a₀ + a₂x + a₄x² + · · · + a_{n 2}x^{n/2 1}

A^[1](x) = a₁ + a₃x + a₅x² + · · · + a_{n 1}x^{n/2 1}

A(x) = A^[0](x²) + xA^[1](x²)

!_n⁰, !_n¹, . . . , !_n^{n 1}

A^[0](x) A^[1](x)

(!_n⁰)², (!_n¹)², . . . , (!_{n 1}⁰ )²

(ㄅ)

Divide  and  Conquer  2  n/2-­‐sized  problems

Combine  the  sub-­‐problem   solutions

Pseudo-­‐code

30.2 The DFT and FFT 911

R^ECURSIVE-FFT.a/

1 n D a:length // n is a power of 2 2 if n == 1

3 return a

4 !ⁿ D e^{2! i=n}

5 ! D 1

6 a^Œ0" D .a⁰; a₂; : : : ; a_n!2/ 7 a^Œ1" D .a¹; a₃; : : : ; a_n!1/

8 y^Œ0" D R^ECURSIVE-FFT.a^Œ0"/

9 y^Œ1" D R^ECURSIVE-FFT.a^Œ1"/

10 for k D 0 to n=2 ! 1

11 y_k D yk^Œ0" C ! yk^Œ1"

12 y_kC.n=2/ D yk^Œ0" ! ! yk^Œ1"

13 ! D ! !ⁿ

14 return y // y is assumed to be a column vector

The R^ECURSIVE-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 y₀ D a⁰ !₁⁰

D a⁰ " 1 D a⁰ :

Lines 6–7 define the coefficient vectors for the polynomials A^Œ0" and A^Œ1". Lines 4, 5, and 13 guarantee that ! is updated properly so that whenever lines 11–12 are executed, we have ! D !n^k. (Keeping a running value of ! from iteration to iteration saves time over computing !n^k from scratch each time through the for loop.) Lines 8–9 perform the recursive DFTⁿ⁼² computations, setting, for k D 0; 1; : : : ; n=2 ! 1,

y_k^Œ0" D A^Œ0".!_n=2^k / ; y_k^Œ1" D A^Œ1".!_n=2^k / ;

or, since !n=2^k D !n^2k by the cancellation lemma, y_k^Œ0" D A^Œ0".!_n^2k/ ;

y_k^Œ1" D A^Œ1".!_n^2k/ :

Divide:  2x    n/2 Conquer

Combine

Evaluate  at  !_n^k

O(n)

O(n log n)

904 Chapter 30 Polynomials and the FFT

a0; a1; : : : ; an!1

b0; b1; : : : ; bn!1

c0; c1; : : : ; c2n!2

Ordinary multiplication Time ‚.n²/ Evaluation

Time ‚.n lg n/

Time ‚.n lg n/

Interpolation

Pointwise multiplication Time ‚.n/

A.!_2n⁰ /; B.!_2n⁰ / A.!_2n¹ /; B.!_2n¹ /

A.!_2n^2n!1/; B.!_2n^2n!1/

::: :::

C.!_2n⁰ / C.!_2n¹ /

C.!_2n^2n!1/

Coefficient

Point-value representations representations

Figure 30.1 A graphical outline of an efficient polynomial-multiplication process. Representations on the top are in coefficient form, while those on the bottom are in point-value form. The arrows from left to right correspond to the multiplication operation. The !²ⁿterms are complex .2n/th roots of unity.

on whether we can convert a polynomial quickly from coefficient form to point- value form (evaluate) and vice versa (interpolate).

We can use any points we want as evaluation points, but by choosing the eval- uation points carefully, we can convert between representations in only ‚.n lg n/

time. As we shall see in Section 30.2, if we choose “complex roots of unity” as the evaluation points, we can produce a point-value representation by taking the discrete Fourier transform (or DFT) of a coefficient vector. We can perform the inverse operation, interpolation, by taking the “inverse DFT” of point-value pairs, yielding a coefficient vector. Section 30.2 will show how the FFT accomplishes the DFT and inverse DFT operations in ‚.n lg n/ time.

Figure 30.1 shows this strategy graphically. One minor detail concerns degree- bounds. The product of two polynomials of degree-bound n is a polynomial of degree-bound 2n. Before evaluating the input polynomials A and B, therefore, we first double their degree-bounds to 2n by adding n high-order coefficients of 0.

Because the vectors have 2n elements, we use “complex .2n/th roots of unity,”

which are denoted by the !²ⁿ terms in Figure 30.1.

Given the FFT, we have the following ‚.n lg n/-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; we can always meet this requirement by adding high-order zero coefficients.

1. Double degree-bound: Create coefficient representations of A.x/ and B.x/ as degree-bound 2n polynomials by adding n high-order zero coefficients to each.

⇥(n²)

Can  we  improve  evaluation  and  interpolation  time  to

?  

⇥(n log n)

⇥(n log n)

How  about  interpolation??      (p.  912  on  Cormen)