6.006- Introduction to Algorithms Lecture 22 Prof. Patrick JailletOutline • “Numerics” - algorithms for operations on large numbers – high precision – cryptography, simulations, etc • Today we will see: – irrationals – large number operations: • multiplication • division 4.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006… 2.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571470109559971605970274534596862014728517418640889198609552329230484308714321450839762603627995251407989687253396546331808829640620615258352395054745750287759961729835575220337531857011354374603 ...Computing to lots of digits ... why? € h11√21. 414 213 562 373 095 048 801 688 724 209 698 078 569 671 875 376 948 073 176 679 ...Computing to lots of digits ... why? • high precision may be needed in some applications • consider Dijkstra for paths between lattice points on plane: – lengths have form – where • is > ? € hii∑€ hai bi hi € hi2= ai2+ bi2€ 1 + 40 + 60€ 12 + 17 + 56€ 1 + 40 + 60 = 15.07052201275€ 12 + 17 + 56 = 15.07052201430Computing to lots of digits ... why? • geometry problem • BD =1 • what is AD? € hABCD1000,000,000,0001€ AD = AC − CD = 500,000,000,000 − 500,000,000,0002−1Computing ; Heron’s method • Iterative approach; also called the Babylonian method • can be seen to be equivalent to the Newton method (more later) • y0 =h; x0 =1 • y1 =(x0 +y0 )/2; x1 =h/(y1) • in general – yi+1 =(xi +yi )/2 – xi+1 =h/(yi+1) € h(c. 10–70 AD)Computing ; Heron’s method • y0 =h; x0 =1 • y1 =(x0 +y0 )/2; x1 =h/(y1) • in general – yi+1 =(xi +yi )/2 – xi+1 =h/(yi+1) • converges quickly € h(c. 10–70 AD) x y 0 1.000000000000000 2.000000000000000 1 1.333333333333330 1.500000000000000 2 1.411764705882350 1.416666666666670 3 1.414211438474870 1.414215686274510 4 1.414213562371500 1.414213562374690 5 1.414213562373090 1.414213562373090Computing ; Newton’s method • Newton’s method: Iterative approach to solving f(x)=0 • Iterative step: – find a line y= f(xi)+f’(xi)(x- xi) tangent to f(x) at xi – set xi+1 to the solution of f(xi)+f’(xi)(x- xi) =0 => xi+1=xi-f(xi)/f’(xi) -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 xi xi+1 f € hfor computing f(x) = x2 – h xi+1= xi - (xi2 -h) /2xi= (xi +h/xi )/2 € hComputing ; Newton’s method • xi+1= xi - (xi2 -h) /2xi= (xi +h/xi )/2 • convergence? € h0 1.000000000000000 1 1.500000000000000 2 1.416666666666670 3 1.414215686274510 4 1.414213562374690 5 1.414213562373090Large number addition • Given: two positive n-digit numbers x,y with radix r (e.g., r=2, r=10) • Goal: compute x+y, using only operations on single digits • Algorithm: keep adding digits from right to left, keeping track of carryovers • Time = O(n) xn-1 xn-2 …x1 x0 yn-1 yn-2 …y1 y0 +Large number multiplication • Given: two positive n-digit numbers x,y with radix r • Goal: compute x*y, using only operations on single digits • Ideas ? xn-1 xn-2 …x1 x0 yn-1 yn-2 …y1 y0 *Divide and conquer: attempt I • Split each number into “high” and “low” parts, each n/2 digits long – x=xHrn/2+xL – y=yHrn/2+yL • We have x*y = (xHrn/2+xL)*(yHrn/2+yL) = xH*yH rn +(xH*yL+xL*yH)rn/2 + xL*yL • Algorithm for mult(x,y): a= mult(xH,yH ), b= mult(xH,yL ), c= mult(xL,yH ), d= mult(xL,yL ) return add( lshift(a,n), lshift(b,n/2), lshift(c,n/2), d) xn-1 xn-2 …x1 x0 yn-1 yn-2 …y1 y0Attempt I: analysis • Algorithm: a= mult(xH,yH ), b= mult(xH,yL ), c= mult(xL,yH ), d= mult(xL,yL ) return add( lshift(a,n), lshift(b,n/2), lshift(c,n/2), d) • Running time T(n) ? – recurrence • T(n)=4 T(n/2) +O(n) – this solves to • T(n)=O(n2) • How to improve on this ?Attempt II • Need to compute x*y = xH*yH rn +(xH*yL+xL*yH)rn/2 + xL*yL using fewer (than 4) recursive multiplications ? • Here is how to do it (Karatsuba’62) a = xH*yH d = xL*yL e = (xH+xL) * (yH+yL) - a - d = xH*yH +(xH*yL+xL*yH) + xL*yL - a - d = xH*yL+xL*yH x*y = a rn + e rn/2 + d • This leads to – T(n)=3 T(n/2) +O(n) which solves to T(n)=O(nlog 3) = O(n1.58496… )Division 1/b • Inversion: – given positive n-digit number b (radix r) – compute 1/b i.e., for some R=rk, compute ⎣R/b⎦ • Iterative algorithm: – start from some x0 – keep computing xi+1 from xi – hope prove this converges to the answerDivision ... again ... Newton’s method • Newton’s method: Iterative approach to solving f(x)=0 • Iterative step: – find a line y= f(xi)+f’(xi)(x- xi) tangent to f(x) at xi – set xi+1 to the solution of f(xi)+f’(xi)(x- xi) =0 => xi+1=xi-f(xi)/f’(xi) -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 xi xi+1 f here we have (next slide) f(x)=1/x-b/R xi+1=2xi - xi2b/RDivision • want to solve f(x)=1/x-b/R=0 • we have f’(x)=-1/x2 • iterative step xi+1=xi-f(xi)/f’(xi) yields xi+1=xi+ xi2 (1/xi -b/R) i.e., xi+1=2xi - xi2b/R • convergence ?Convergence of xi+1=2xi - xi2b/R • Assume xi = R/b (1+ei), ei =error • Assumptions: – |ei| is small – Ignore the round-off errors caused by the “/R” operation • How does each iteration affect ei ? xi+1=2xi - xi2b/R = 2R/b (1+ei) - [R/b (1+ei)]2 b/R = R/b[2+2ei -1 -2ei - ei2] = R/b[1- ei2] = R/b[1+ ei+1], where ei+1=- ei2 • We have |ei+1|=|ei|2 - “quadratic”
View Full Document
Unlocking...