Efficient Polynomial Multiplication Techniques and Algorithms

CSE 417: Algorithms and Computational Complexity Winter 2001 Lecture 14 Instructor: Paul Beame

Multiplying Faster • On the first HW you analyzed our usual algorithm for multiplying numbers • Q(n2) time • We can do better! • We’ll describe the basic ideas by multiplying polynomials rather than integers • Advantage is we don’t get confused by worrying about carries at first

Note on Polynomials • These are just formal sequences of coefficients so when we show something multiplied by xk it just means shifted k places to the left

Polynomial Multiplication • Given: • Degree m-1 polynomials P and Q • P = a0 +a1 x+a2 x2 + … +am-2xm-2 +am-1xm-1 • Q = b0+b1 x+b2 x2+ … +bm-2xm-2 +bm-1xm-1 • Compute: • Degree 2m-2 Polynomial PQ • PQ =a0b0 + (a0b1+a1b0) x + (a0b2+a1b1 +a2b0) x2 +...+ (am-2bm-1+am-1bm-2) x2m-3+ am-1bm-1x2m-2 • Obvious Algorithm: • Compute all aibj and collect terms • Q (n2)time

Naive Divide and Conquer • Assume m=2k • P = (a0+ a1 x + a2 x2 + ... + ak-1 xk-1) + (ak + ak+1 x +… + am-2xk-2 + am-1xk-1) xk = P0 + P1xk • Q = Q0 + Q1xk • PQ = (P0+P1xk)(Q0+Q1xk) = P0Q0 + (P1Q0+P0Q1)xk + P1Q1x2k • 4 sub-problems of size k=m/2 plus linear combining • T(m)=4T(m/2)+cm • Solution T(m) = O(m2)

Karatsuba’s Algorithm • A better way to compute the terms • Compute • P0Q0 • P1Q1 • (P0+P1)(Q0+Q1) which is P0Q0+P1Q0+P0Q1+P1Q1 • Then • P0Q1+P1Q0 = (P0+P1)(Q0+Q1) - P0Q0- P1Q1 • 3 sub-problems of size m/2 plus O(m) work • T(m) = 3 T(m/2) + cm • T(m) = O(ma) where a = log23 = 1.59...

Karatsuba’s Algorithm Alternative • Compute • A0=P0Q0 • A1=(P0+P1)(Q0+Q1), i.e. P0Q0+P1Q0+P0Q1+P1Q1 • A-1=(P0 - P1)(Q0 - Q1), i.e. P0Q0 - P1Q0 -P0Q1+P1Q1 • Then • P1Q1= (A1+A-1)/2 - A0 • P0Q1+P1Q0 = (A1 - A-1)/2 • 3 sub-problems of size m/2 plus O(m) work • T(m) = 3 T(m/2) + cm • T(m) = O(ma) where a = log23 = 1.59...

What Karatsuba’s Algorithm did • For y=xk we wanted to compute P(y)Q(y)=(P0+P1 y)(Q0+Q1y) • We evaluated • P(0)= P0Q(0)= Q0 • P(1) = P0+P1 and Q(1)= Q0+Q1 • P(-1) = P0 - P1 and Q(-1)= Q0 -Q1 • We multiplied P(0)Q(0),P(1)Q(1), P(-1)Q(-1) • We then used these 3 values to figure out what the degree 2 polynomial P(y)Q(y) was • Interpolation

Interpolation Given set of values at 5 points

Interpolation Find unique degree 4 polynomial going through these points

Multiplying Polynomials by Evaluation & Interpolation • Any degree n-1 polynomial R(y) is determined by R(y0), ... R(yn-1) for any n distinct y0,...,yn-1 • To compute PQ(assume degree at most n-1) • Evaluate P(y0),..., P(yn-1) • Evaluate Q(y0),...,Q(yn-1) • Multiply values P(yi)Q(yi) for i=0,..., n-1 • Interpolate to recover PQ

Multiplying Polynomials by Evaluation & Interpolation ordinary polynomial multiplication Q(n2) P: a0,a1,...,am-1 Q: b0,b1,...,bm-1 R:c0,c1,...,cn-1 evaluation at y0,...,yn-1 O(?) interpolation from y0,...,yn-1 O(?) R(y0)=P(y0)Q(y0) R(y1)=P(y1)Q(y1) ... R(yn-1)=P(yn-1)Q(yn-1) P(y0),Q(y0) P(y1),Q(y1) ... P(yn-1),Q(yn-1) point-wise multiplication of numbersO(n)

Complex Numbers i2 = -1 i c+di a+bi e+fi To multiply complex numbers: 1. add angles 2. multiply lengths q+j j -1 q 1 e+fi = (a+bi)(c+di) a+bi =cosq+i sin q= eiq c+di =cosj+i sin j= eij e+fi =cos(q+j)+i sin (q+j)= ei(q+j) -i e2pi = 1 epi = -1

w2=i w3 w w4=-1 w0=1=w8 w7 w5 w6= -i Primitive n-th root of 1w=wn Let w= wn= ei 2p/n = cos(2p/n)+isin (2p/n) Properties: wn = 1. Any otherzs.t. zn = 1 hasz= wkfor somek<n. Ifnis evenw2=wn/2is a primitive n/2-th root of 1. wk+n/2= -wk i2= -1 e2pi = 1

Multiplying Polynomials by Fast Fourier Transform ordinary polynomial multiplication Q(n2) P: a0,a1,...,am-1 Q: b0,b1,...,bm-1 R:c0,c1,...,cn-1 evaluation at 1,w,...,wn-1 O(n log n) interpolation from 1,w,...,wn-1 O(n log n) P(1),Q(1) P(w),Q(w) P(w2),Q(w2) ... P(wn-1),Q(wn-1) R(1)=P(1)Q(1) R(w)=P(w)Q(w) R(w2)=P(w2)Q(w2) ... R(wn-1)=P(wn-1)Q(wn-1) point-wise multiplication of numbersO(n)

The key idea (since n is even) • P(w)=a0+a1w+a2w2+a3w3+a4w4+...+an-1wn-1 = a0 +a2w2+a4w4+...+ an-2wn-2+ a1w+a3w3 +a5w5 +...+an-1wn-1 = Peven(w2) + w Podd(w2) • P(-w)=a0 -a1w+a2w2 -a3w3+a4w4-... -an-1wn-1=a0 +a2w2+a4w4+...+ an-2wn-2-(a1w+a3w3 +a5w5 +...+an-1wn-1) = Peven(w2) - w Podd(w2)

The recursive idea for n a power of 2 • Also • Peven and Podd have degree n/2 • P(wk)=Peven(w2k)+wkPodd(w2k) • P(-wk)=Peven(w2k)-wkPodd(w2k) • Recursive Algorithm • Evaluate Peven at 1,w2,w4,...,wn-2 • Evaluate Podd at 1,w2,w4,...,wn-2 • Combine to compute P at 1,w,w2,...,wn/2-1 • Combine to compute P at -1,-w,-w2,...,-wn/2-1(i.e. atwn/2, wn/2+1 , wn/2+2,..., wn-1) w2 is an n/2-th root of 1 so problem is of same type but smaller size

Analysis and more • Run-time • T(n)=2T(n/2)+cn so T(n)=O(n log n) • So much for evaluation ... what about interpolation? • Given • r0=R(1), r1=R(w), r2=R(w2),..., rn-1=R(wn-1) • Compute • c0, c1,...,cn-1 s.t. R(x)=c0+c1x+...cn-1xn-1

Interpolation  Evaluation: strange but true • Weird fact: • If we define a new polynomial T(x) = r0 + r1x + r2x2 +...+ rn-1xn-1 where r0, r1, ... , rn-1 are the evaluations of R at 1, w, ... , wn-1 • Then ck=T(w-k)/n for k=0,...,n-1 • So... • evaluate T at 1,w-1,w-2,...,w-(n-1) then divide each by n to get the c0,...,cn-1 • w-1 behaves just like w did so the same O(n log n) evaluation algorithm applies !

Why this is called the discrete Fourier transform • Real Fourier series • Given a real valued function f defined on [0,2p] the Fourier series for f is given by f(x)=a0+a1cos(x)+a2cos(2x)+...+amcos(mx)+... where am= is the component of f of frequency m • In signal processing and data compression one ignores all but the components with large am and there aren’t many since

where fk=f(2kp/n) Why this is called the discrete Fourier transform • Complex Fourier series • Given a function f defined on [0,2p] the complex Fourier series for f is given by f(z)=b0+b1eiz+b2e2iz+...+bmemiz+... where bm= is the component of f of frequency m • If we discretize this integral using values at n equally spaced points between 0 and 2pwe get 2p/n apart just like interpolation!

Efficient Polynomial Multiplication Techniques and Algorithms