Transcript The Fast Fourier Transform and Applications to Multiplication Analysis of Algorithms Prepared by John Reif, Ph.D.

The Fast Fourier Transform and
Applications to Multiplication
Analysis of Algorithms
Prepared by
John Reif, Ph.D.
Topics and Readings:
- The Fast Fourier Transform
• Reading Selection:
• CLR, Chapter 30
Advanced Material :
- Using FFT to solve other Multipoint Evaluation
Problems
- Applications to Multiplication
Nth Roots of Unity
• Assume Commutative Ring (R,+,·, 0,1)
•  is principal nth root of unity if
– i  1 for i = 1, …, n-1
– n = 1, and
n-1
jp

  0 for 1  p  n
j=0
• Example:
 e
2 i/n
for complex numbers
Example of nth Root of Unity for
Complex Numbers
e
2 i/8
is the 8th root of unity
Fourier Matrix
1
1 1



n 1

1 

M n ( ) = 1  2
 2( n 1) 





n 1
( n 1)( n 1) 

1




so M( )ij =  ij for 0  i, j<n
 a0 


given a = 

a 
 n-1 
Discrete Fourier Transform
Input a column n-vector a = (a0, …, an-1)T
Output an n-vector which is the product of
the Fourier matrix times the input vector
DFTn (a) = M( ) x a
 f0 
 
=   where
f 
 n-1 
n-1
fi =
ik
a

 k
k=0
Inverse Fourier Transform
DFTn-1 (a) = M( )-1 x a
1 -ij
Theorem M( ) = 
n
proof We must show M( )  M( )-1 = I
-1
ij
1
n
n-1
ik -kj

 
k=0
1
=
n
n-1
k(i-j)


k=0
0 if i-j  0
= 
1 if i-j = 0
n-1
using identity   kp  0, for 1  p < n
k=0
Fourier Transform is Polynomial
Evaluation at the Roots of Unity
Input a column n-vector a = (a0, …, an-1)T
Output an n-vector (f0, …, fn-1)T which are the
values polynomial f(x)at the n roots of unity
 f0 
 
DFTn (a) =   where
f 
 n-1 
i
f i = f( ) and
n-1
f(x) =  a j  x j
j=0
Fast Fourier Transform
• Viewed as Evaluation Problem:
naïve algorithm takes n2 ops
• Divide and Conquer gives FFT
with O(n log n) ops for n a power of 2
• Key Idea:
• If  is nth root of unity
then 2 is n/2th root of unity
• So can reduce the problem to two
subproblems of size n/2
Algorithm FFTn
•
Input a = (a0, …, an-1)T, n a power of 2
[1] If n=1 then ouput
T
 '

'
[2]  f 0 ,..., f n   FFTn ((a0 , a2 ,..., an 2 )T )
1
2 
2

T
 "
" 
T
 f 0 ,..., f n 1   FFTn ((a1 , a3 ,..., an 1 ) )
2 
2

n
[3] For i=0, ...,  1 do fi  fi'   i fi"
2
f n  fi'   i fi"
i+
[4] Output (f 0 , f1 , ..., f n-1 )
2
FFT Circuit
(also known as Butterfly Network)
• Total Recursion depth = log n
• Communication Distance 2d at depth d
f i = a 0 + a1 i + a 2 2i +...+a n-1 (n-1)i
f i = f i' +  i f i" where
i(n-2)
2
2
f i' = a 0 + a 2 ( 2 )i + a 4 ( 2 ) 2i +...+ a n-2 ( )
i(n-2)
2
2
f i" = a1 + a 3 ( 2 )i +...+ a n-1 ( )
 ' 
 a0 
f 0 


a

  M ( 2 )  2   DFT ((a , a ,..., a )T )
n
n
0
2
n2




2
2
f n' 


 2 1 
 an  2 
 " 
 a1 
f 0 
a 

  M ( 2 )  3   DFT (( a , a ,..., a )T )
n
n
1
3
n 1




2
2
f n" 


 2 1 
 an 1 
Note: f
'
n
i
2
'
i
=f , f
"
n
i
2
n
 f , i=0, ...,  1
2
"
i
n
1
2 2
But  n =1, so ( )
=  n  ( 2 )i   2i
n
for i=0, ...,  1
2
n
'
i
"
Thus, f i = f i +  f i for i=0, ...,  1
2
and f
i+
'
=
f
n
i +
n
2
 f i"
2
=f - f
'
i
i+
i
n
2 2
"
i
n
for i=0, ...,  1
2
n
2
since ( )   n  1, so  = -1
Operation Counts for FFT Algorithm
• Assume n = 2k
• # additions
Add(n) = 2· Add(n/2) + n
= n log n
• # multiplications
Mult(n) = 2· Mult(n/2) + n/2
= ½ n log n
• Total Time
O(n log n)
• Note in complex FFT,
# real ops is 5 n log n
Multipoint Polynomial Evaluation
• Input polynomial f ( x) 
n 1
i
a
x
 i
i 0
• Problem evaluate f(x) at x0, x1, …, xn-1
• Easy Cases:
FFT Case xi = i
= principal root of unity
Multipoint Polynomial Evaluation
(cont’d)
Summary of FFT:
method f ( x)  f '( y )  x f "( y )
where y  x
2
f '( x), f "( x) both degree halved
 needed to only evaluate at half as many points
Other Polynomial Evaluation
Problems Solved by FFT
Each costs O(n log n) time
•
Evaluate at points Xi = bai + d for i=0,…, n-1(Chirp
Transfom)
– Reduced to FFT
•
Single point evaluation of all derivatives of a polynomial
– Solve by reduction to above Chirp Transform of case 2)
•
Evaluate at points Xi = b(ai)2+ cai + d for i=0,…, n-1
– Solve by divide and conquer similar to FFT
Single Point Evaluation of all
Derivatives of Polynomial
• Input
n 1
f ( x)   ai xi
and point x0
i 0
• output
k
d
f ( x)
k
f ( x) 
x  x0 for k  0,..., n  1
dx
Single Point Evaluation of all
Derivatives of Polynomial (cont’d)
• Taylor Series Representation of
n 1
f ( x)   ci ( x  x0 )i
i 0
Then
f ( k ) ( x0 )  k !ck
 reduces to case of evaluation at points
xi  abi for i  0,..., n 1
• Solve this Chirp Transform problem by
reduction FFT
Advanced Material:
Further Applications of FFT
1) Convolution: Products and Powers of
Polynomials
• Used for for Integer Multiplication
Algorithms
• Also used for Filtering on infinite
input streams
2) Division and Inverse of Polynomials
3) Multipoint Evaluation and
Interpolation
Advanced Material: Products and
Powers of Polynomials
• Input vectors
a = (a0, a1, …, an-1)T
b = (b0, b1, …, bn-1)T
• Definition of Convolution c = a  b
n-1
Where
ci = a jbi-j
for i=0, …, 2n-1
j=0
define ak = bk = 0
if k< 0 or kn
Products and Powers of Polynomials
(cont’d)
• Convolution Theorem
a  b = FFT
-1
2n
 FFT2n (a)  FFT2n (b) 
• Application to Polynomial Products:
n-1
p(x) =  a i x i
i=0
n-1
q(x) =  b j x j
j=0
2n-2
n-1
i=0
j=0
p(x)  q(x) =  ci x i where c i =  a jb i-j
Products of m Polynomials
n-1
for k=1, ..., m let Pk (x) =  a k,i x i
i=0
m(n-1)
m
m
 P (x) =  c x , where c =  a
i
k
k=1
i
i=0
i
k,jk
k=1
 jk =1
• Generalized Convolution Theorem
a1  a 2  ...  a m =
FFTn-1m  FFTn m (a1 )  FFTn m (a 2 )  ...  FFTn m (a m ) 
Wrapped Convolutions
• a = (a0, a1, …, an-1)T , b = (b0, b1, …,
bn-1)T
• Positive wrapped convolution is
c = (c0, c1, …, cn-1)T
i
n-1
j=0
j=i+1
ci = a jbi-j +  a jbn+i-j
• Negative wrapped convolution is
d = (d0, d1, …, dn-1)T
i
n-1
j=0
j=i+1
di = a jbi-j -  a jbn+i-j
Application of Wrapped Convolution
to Modular Polynomial Products
n-1
p(x) =  a i x
i
i=0
n-1
q(x) =  b j x
j
j=0
p(x)  q(x) mod(x +1)
n
n-1
=  d i x since x = -1 mod(x +1)
i
i=0
n
n
Computing Positive Wrapped
Convolution
• Let  = principal nth root of unity
• Assume n has multiplicative inverse,
Theorem
-1
n
c = FFT
 FFTn (a)  FFTn (b) 
is the positive wrapped convolution of
n-vectors a and b.
Computing Negative Wrapped
Convolution
• Also

ˆd = FFT-1 FFT (a)
ˆ
ˆ

FFT
(b)
n
n
n

is the negatively wrapped convolution of
n-vectors a and b
where

bˆ =  b , b , ...,
aˆ = a 0 , a1 , ..., a n-1
n-1
0
1
n-1


b n-1
T
T
and 2 =  = principal nth root of unity
Integer Multiplication by Polynomial
Product (solved via FFT)
• Input n bit integers a,b
define polynomials degree k = n/L
k-1
a(x) = a i x ,
i
0  ai  2
L
i=0
k-1
b(x) = bi x i ,
0  bi  2L
i=0
L
L
so a = a(2 ), b = b(2 )
Integer Multiplication by Polynomial
Product (cont’d)
• Idea
1) Compute c(x) = a(x)· b(x)
by convolution
2) Evaluate c(2L) = a· b
Integer Multiplication Algorithms using
Reduction to Polynomial Product
• Pollard Mult Algorithm

2
O(n(logn) )(loglogn) use L = logn
• Karp Mult Algorithm
O(n(logn)2 ) use L = n
• Schönage-Strassen Mult Algorithm
O(n(logn)(loglogn)) use L = n
and wrapped convolution
Pollard Multiplication Algorithm
• n = kL, L = 1 + log k
1) Choose primes P1, P2, P3 where
P1  P2  P3  4  k
3
and Pi = αi  2 + 1, αi = O(1)
L
2) Compute C(x) by convolution over
finite field Zpi for i =1,2,3
(requires k mults on 2L bit integers)
Pollard Multiplication Algorithm
(cont’d)
3) Evaluate C(2L)
• Time Bounds
recursive mults
FFT
T(n) = 3kT(2L) + O(k logk)  O(L)
2
= 3kT (2(1 + logk)) + O(k(log k) )
 O(n(log n) 2 (log log n) ) for any  > 0
Korp Multiplication Algorithm
n = 2s = kL
 2s
if s even
2
k =  (s-1)
2 2 else

1) Compute C(x) modulo k by
convolution
2) Compute C(x) modulo (22L+1) by
convolution
3) Compute C(x) coefficients from C(x)
mod k, C(x) mod (22L+1) by Chinese
remaindering
Korp Multiplication Algorithm
(cont’d)
4) Compute C(2L)
• Time
recursive mults
FFT
T(n) = 2kT(2L) + O(k logk)O(L)
= 2 nT (2 n ) + O(n log n)
 O(n(log n)2 )
Schönage-Strassen Multiplication
Algorithm
(2’) Compute C(x) mod (xk+1) modulo
(22L+1) by wrapped convolution
 requires only k recursive mults on
2L bit numbers
• Time
recursive mults
FFT
T(n) = kT(2L) + O(k logk)O(L)
= nT (2 n ) + O(n log n)
 O(n log n)(log log n)
Still Open Problem: How Fast Can
You Multiply Integers?
• Can you mult n bit integers in
O(n log n) time?
The Fast Fourier Transform and
Applications to Multiplication
Analysis of Algorithms
Prepared by
John Reif, Ph.D.