• 沒有找到結果。

# Polynomial + Fast Fourier Transform

N/A
N/A
Protected

Share "Polynomial + Fast Fourier Transform"

Copied!
21
0
0

(1)

## Fast  Fourier  Transform

Michael  Tsai 2017/06/13

(2)

### 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

ajxj

0

1

2

2

n 1

n 1

0

1

n 1

k

(3)

• Using  it:

• Evaluation:

A(x) =

n 1X

j=0

ajxj

0

1

n 1

### )

Vector

A(x0) = a0 + x0(a1 + x0(a2 + · · · + x0(an 2 + x0(an 1)) . . . )) Horner’s  rule:

0

1

n 1

A(x) + B(x)

0

1

+ n 1

### )

(a0 + b0, a1 + b1, . . . , an 1 + bn 1) O(n)

O(n)

(4)

### Prob.:  Polynomial  Multiplication

• Example: A(x) = 6x3 + 7x2 10x + 9

B(x) = 2x3 + 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

ajxj and

B.x/ D

n!1

X

jD0

bjxj ; then

C.x/ D

n!1

X

jD0

cjxj ;

where cj D aj C bj for j D 0; 1; : : : ; n ! 1. For example, if we have the polynomials A.x/ D 6x3 C 7x2 ! 10x C 9 and B.x/ D !2x3 C 4x ! 5, then C.x/ D 4x3 C 7x2 ! 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 6x3 C 7x2 ! 10x C 9 and B.x/ D !2x3 C 4x ! 5 as follows:

6x3 C 7x2 ! 10x C 9

! 2x3 C 4x ! 5

! 30x3 ! 35x2 C 50x ! 45 24x4 C 28x3 ! 40x2 C 36x

! 12x6 ! 14x5 C 20x4 ! 18x3

! 12x6 ! 14x5 C 44x4 ! 20x3 ! 75x2 C 86x ! 45 Another way to express the product C.x/ is

C.x/ D

2nX!2 jD0

cjxj ; (30.1)

where cj D

Xj kD0

akbj!k : (30.2)

2

O(n) O(n) O(n)

(5)

### Prob.:  Polynomial  Multiplication

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

• Degree  bound:

• !

• Problem:  can  we  reduce  this  time  complexity?

C(x) =

2n 2X

j=0

cjxj cj =

Xj

k=0

akbj k

na + nb

2

(6)

### 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

ajxj

Degree-­‐bound:   n

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

yk = A(xk) k = 0, 1, . . . , n 1

A(xk) O(n)

A(xk) O(n2)

k

O(n log n)

(7)

### 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)

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

### x

k

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

yk = A(xk)

(8)

### Point-­‐Value  Representation

• Using  it:

since                                                                          ,

A(x) =

n 1X

j=0

ajxj

Degree-­‐bound:   n

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

yk = A(xk) k = 0, 1, . . . , n 1 Evaluation

Interpolation

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

{(x0, y0), (x1, y1), . . . , (xn 1, yn 1)} {(x0, y00), (x1, y10), . . . , (xn 1, yn 10 )} A:

B:

{(x0, y0 + y00), (x1, y1 + y10 ), . . . , (xn 1, yn 1 + yn 10 )}

C: O(n)

(9)

### Multiplication

• Multiplication:

since                                                                            ,

• Problem!

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

defined!

C(x) = A(x)B(x) C(xk) = A(xk)B(xk)

A:

B:

C:

{(x0, y0), (x1, y1), . . . , (xn 1, yn 1)} {(x0, y00), (x1, y10), . . . , (xn 1, yn 10 )}

{(x0, y0y00), (x1, y1y10), . . . , (xn 1, yn 1yn 10 )} n pairs

n pairs

n pairs

(10)

### 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(xk) = A(xk)B(xk)

A:

B:

C:

{(x0, y0y00 ), (x1, y1y10 ), . . . , (x2n 1, y2n 1y2n 10 )}

{(x0, y0), (x1, y1), . . . , (x2n 1, y2n 1)} {(x0, y00), (x1, y10), . . . , (x2n 1, y2n 10 )}

2n  pairs 2n  pairs

2n  pairs

### O(n)

(11)

904 Chapter 30 Polynomials and the FFT

a0; a1; : : : ; an!1

b0; b1; : : : ; bn!1

c0; c1; : : : ; c2n!2

Ordinary multiplication Time ‚.n2/ Evaluation

Time ‚.n lg n/

Time ‚.n lg n/

Interpolation

Pointwise multiplication Time ‚.n/

A.!2n0 /; B.!2n0 / A.!2n1 /; B.!2n1 /

A.!2n2n!1/; B.!2n2n!1/

::: :::

C.!2n0 / C.!2n1 /

C.!2n2n!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 !2nterms 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 !2n 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.

⇥(n2) ⇥(n2)

Can  we  improve  evaluation  and  interpolation  time  to

?

(12)

### 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:

n

iu

### e

2⇡ik/n k = 0, 1, . . . , n 1

(13)

### Example

30.2 The DFT and FFT 907

!1 1

i

!i

!80 D !88

!81

!82

!83

!84

!85

!86

!87

Figure 30.2 The values of !80; !81; : : : ; !87 in the complex plane, where !8 D e2! i=8 is the prin- cipal 8th root of unity.

!n D e2! i=n (30.6)

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

The n complex nth roots of unity,

!n0; !n1; : : : ; !nn!1 ;

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

!nj!nk D !njCk D !n.jCk/ mod n. Similarly, !n!1 D !nn!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 nd k D !nk : (30.7)

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

!d nd k D !

e2! i=d n"d k

D !

e2! i=n"k D !nk :

2Many other authors define !n differently: !n 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 !n.

Principle  n-­‐th root  of  unity:

n  complex  n-­‐th roots  of  unity:

n

2⇡i/n

n0

n1

nn 1

(14)

### 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

ajxj

!n0, !n1, . . . , !nn 1

yk = A(!nk) =

n 1X

j=0

aj!nkj k = 0, 1, . . . , n 1

0

1

n 1

0

1

n 1

2

still

(15)

### Physical  Meaning  of  DFT

yk = A(!nk) =

n 1X

j=0

aj!nkj

aj = 1 n

n 1X

k=0

yk!n kj Inverse  DFT

Signal  at  different  frequencies   Weight  of  that

frequency

(16)

### 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) = a0 + a1x + a2x2 + a3x3 + · · · + an 2xn 2 + an 1xn 1

A[0](x) = a0 + a2x + a4x2 + · · · + an 2xn/2 1

A[1](x) = a1 + a3x + a5x2 + · · · + an 1xn/2 1

A(x) = A[0](x2) + xA[1](x2)

(17)

### Fast  Fourier  Transform

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

2. Combine  the  result  using  (ㄅ)

A[0](x) = a0 + a2x + a4x2 + · · · + an 2xn/2 1

A[1](x) = a1 + a3x + a5x2 + · · · + an 1xn/2 1

A(x) = A[0](x2) + xA[1](x2)

!n0, !n1, . . . , !nn 1

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

(!n0)2, (!n1)2, . . . , (!n 10 )2

(ㄅ)

(18)

### 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

(19)

### Fast  Fourier  Transform

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

2. Combine  the  result  using  (ㄅ)

A[0](x) = a0 + a2x + a4x2 + · · · + an 2xn/2 1

A[1](x) = a1 + a3x + a5x2 + · · · + an 1xn/2 1

A(x) = A[0](x2) + xA[1](x2)

!n0, !n1, . . . , !nn 1

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

(!n0)2, (!n1)2, . . . , (!n 10 )2

(ㄅ)

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

Combine  the  sub-­‐problem   solutions

(20)

### Pseudo-­‐code

30.2 The DFT and FFT 911

RECURSIVE-FFT.a/

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

3 return a

4 !n D e2! i=n

5 ! D 1

6 aŒ0" D .a0; a2; : : : ; an!2/ 7 aŒ1" D .a1; a3; : : : ; an!1/

8 yŒ0" D RECURSIVE-FFT.aŒ0"/

9 yŒ1" D RECURSIVE-FFT.aŒ1"/

10 for k D 0 to n=2 ! 1

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

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

13 ! D ! !n

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 y0 D a0 !10

D a0 " 1 D a0 :

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 !nk. (Keeping a running value of ! from iteration to iteration saves time over computing !nk from scratch each time through the for loop.) Lines 8–9 perform the recursive DFTn=2 computations, setting, for k D 0; 1; : : : ; n=2 ! 1,

ykŒ0" D AŒ0".!n=2k / ; ykŒ1" D AŒ1".!n=2k / ;

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

ykŒ1" D AŒ1".!n2k/ :

Divide:  2x    n/2 Conquer

Combine

Evaluate  at  !nk

O(n)

### O(n log n)

(21)

904 Chapter 30 Polynomials and the FFT

a0; a1; : : : ; an!1

b0; b1; : : : ; bn!1

c0; c1; : : : ; c2n!2

Ordinary multiplication Time ‚.n2/ Evaluation

Time ‚.n lg n/

Time ‚.n lg n/

Interpolation

Pointwise multiplication Time ‚.n/

A.!2n0 /; B.!2n0 / A.!2n1 /; B.!2n1 /

A.!2n2n!1/; B.!2n2n!1/

::: :::

C.!2n0 / C.!2n1 /

C.!2n2n!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 !2nterms 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 !2n 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.

⇥(n2)

Can  we  improve  evaluation  and  interpolation  time  to

?

### ⇥(n log n)

⇥(n log n)

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

1 As an aside, I don’t know if this is the best way of motivating the definition of the Fourier transform, but I don’t know a better way and most sources you’re likely to check

The notation T n (x) is used to represent the nth partial sum of this series and we can call it as it the nth-degree Taylor polynomial of f at a... In general, it can be shown

The sample point can be chosen to be any point in the subrectangle R ij , but if we choose it to be the upper right-hand corner of R ij [namely (x i , y j ), see Figure 3],

We are going to define the length of a general curve by first approximating it by a polygon and then taking a limit as the number of segments of the polygon is increased...

Since the assets in a pool are not affected by only one common factor, and each asset has different degrees of influence over that common factor, we generalize the one-factor

In Chapter 3, we transform the weighted bipartite matching problem to a traveling salesman problem (TSP) and apply the concepts of ant colony optimization (ACO) algorithm as a basis

A factorization method for reconstructing an impenetrable obstacle in a homogeneous medium (Helmholtz equation) using the spectral data of the far-field operator was developed

A factorization method for reconstructing an impenetrable obstacle in a homogeneous medium (Helmholtz equation) using the spectral data of the far-eld operator was developed