CS 155, Programming Paradigms Fall 2014, SJSU important numeric algorithms

1. CS 155, Programming ParadigmsFall 2014, SJSUimportant numeric algorithms Jeff Smith

2. Applications of the greatest common divisor • An efficient algorithm for finding the greatest common divisor gcd(a,b) of two nonnegative integers a and b is useful for • representing fractions (exact noninteger values) • applications to public key encryption

3. About the greatest common divisor • It can be shown that gcd(a,b) is the smallest positive integer such that d = ax + by for integer values x and y. • cf. Theorem 31.2 of CLRS. • Sometimes gcd(a,b) is written just (a,b).

4. Divide-and-conquer and the GCD • If q is any integer, then any common divisor of a and b is a common divisor of b and a-qb. • If a >= b and q = a div b, then a-qb = a mod b, and gcd(a,b) = gcd(b, a mod b). • so we’ve reduced the computation of gcd(a,b) to a smaller problem of the same type • If b=0, we have an immediate solution • since then gcd(a,b) = a.

5. Euclid’s algorithm for the GCD • The resulting algorithm is known as Euclid’s algorithm • cf. p. 935 of CLRS • Note that when a<b, gcd(a,b) is reduced in one step to gcd(b,a) • so the asymptotic time complexity doesn’t depend on the relative sizes of a and b. • That the algorithm requires only Θ(log a) divisions is not hard to see.

6. Time complexity of Euclid’s algorithm • We may assume a≥b. • It’s enough to show that for every two steps, the length of the first argument is halved, • if a >= 2b, a single step reduces the first argument of the gcd function by at least half. • if instead a < 2b, then a mod b = a-b, so gcd(a,b) = gcd(b,a-b). Since a-b < a-a/2 = a/2, we have cut a at least in half in two steps

7. The worst case for Euclid’s algorithm • CLRS, p.936 gives a precise result • Only logα a divisions are required • where α is the golden ratio describing the growth rate of the Fibonacci sequence. • Intuition: all quotients are 1 in finding the gcd of two adjacent Fibonacci numbers • So the sequence of remainders should decrease as slowly as possible • but this case is still logarithmic

8. Reciprocals mod n • A residue mod n is an integer value strictly between 0 and n • i.e., a possible remainder when dividing by n • The multiplicative inverse (or reciprocal) mod n of a residue x is a residue y with xy = 1 mod n • the application to public key encryption depends on finding reciprocals mod n • division mod n can be defined in terms of reciprocals.

9. Existence of reciprocals mod n • Not all residues have reciprocals modulo all values of n. • e.g., there’s no x such that 2x mod 10 is 1 • A residue x has a reciprocal mod n iff gcd(x,n) = 1 • since gcd(x,n) = 1 holds iff • there are i and j with ix + jn = 1 iff • ix = 1 mod n

10. Extending Euclid’s algorithm • An algorithm for finding the reciprocal of x mod y can be obtained by extending Euclid’s algorithm. • Note that if gcd(a,b) = d and gcd(u,v) is a subproblem of gcd(a,b), then gcd(u,v) = d. • An extended Euclidean algorithm can maintain, bottom-up, values of x and y such that ux+vy= d = gcd(u,v).

11. The extended Euclidean algorithm • In the final subproblem, the arguments are d and 0, and clearly d = 1(d) + 0(0). • And if (u, v) and (v, u mod v) are successive arguments to Euclid’s algorithm, then there’s an integer q with • u mod v = u – qv, and since • d = v(x) + (u mod v)(y), for some x and y, then • d = v(x) + (u – qv)(y) = uy + (x-qy)v. • So the coefficients for u and v can be taken to be y and x-qy respectively

12. The final algorithm • The final version of the algorithm appears on page 937 of CLRS. • Note that the integer quotient of a and b is multiplied by the old y after the recursive call • this is the quotient we called q above • Given arguments a >= b, the extended algorithm halts after O(log a) divisions • just as in the original algorithm

13. Example: finding 11-1 mod 26 • Let n=a=26 and x=b=11. • We have (26, 11) = (11, 4) = (4, 3) = (3,1) = (1,0) = 1. • So 1 = 1(1) + 0(0), • = 3(0) + 1(1), since q=3 and 1 = 1-0q • = 4(1) + 3(-1), since q=1 and -1 = 0-1q • = 11(-1) + 4(3), since q=2 and 3 = 1-(-1)q • = 26(3) + 11(-7), since q=2 and -7=-1-3q • So the reciprocal is 19 (the residue for -7) • we can verify that 11(19) = 209 = 1 mod 26

14. Exponentiation by divide-and-conquer • When computing ab, many fewer than b-1 multiplications are required if the exponent is repeatedly halved. • We need a way of dealing with exponents that aren’t powers of 2. • So consider computing ab as • (ab/2)2 if b is even • a(a(b-1)/2)2 if b is odd

15. Bottom-up exponentiation • A bottom-up implementation is available. If b = 50 = 1100102, we can find a1 = a and • a3 = a*(a1)2 • a6= (a3)2 • a12 = (a6)2 • a25 = a*(a12)2 • a50 = (a25)2 • Here the 1 bits, reading from left to right, tell whether to multiply by a or not.

16. A drawback of divide-and-conquer exponentiation • It’s easy to show that this divide-and-conquer strategy (of CLRS, p. 957) uses Θ(log k) multiplications. • However in the new strategy, the factors are generally larger than in the naïve strategy • so it’s not clear that the new strategy is asymptotically better.

17. An important special case • But for exponentiation mod M, we may assume the factors are bounded in size • this special case is relevant for public-key encryption • In this special case, the length of each factor is linear in the length of M. • So the time required for exponentiation mod M is quadratic in the length of M • and linear in the length of the exponent

18. Representations for polynomial multiplication • Sometimes algorithms can be made faster by changing the representation of their input. • One example is polynomial multiplication. • Here we take the standard representation of the polynomial p(x) = Σaixi to be the sequence (ai) of coefficients.

19. A new representation for polynomial multiplication • But p(x) = Σaixi could be represented instead in terms of its values on d+1 distinct points • where d is the degree of the polynomial • that is, we could represent p(x) in terms of the values {(xj, p(xj)} • Note that for d=1, we are representing a linear equation by by two points • that is, two points determine a line!

20. Multiplication in the new representation • It’s easy to find the product of two polynomials in the new representation • we just multiply pointwise! • That is, the product of the polynomials p and q represented by {(xj, yj)} and {(xj, zj)} would have representation {(xj, yjzj)}. • Warning: we need to have enough points {xj}to uniquely represent the product • its degree is the sum of the degrees of p and q

21. Changing representation • To profit from the new representation, we need to be able to covert to it and from it. • Converting to the new representation is polynomial evaluation • given {xj} and p(x) = Σaixi , we want {p(xj)} • Converting from the new representation is polynomial interpolation • given {(xj, p(xj))}, we want {ai}

22. Improving evaluation and interpolation • To profit from the new representation, we need efficient algorithms for evaluation and interpolation • The standard general algorithms are too slow in each case to give an overall improvement for multiplication

23. General algorithms for evaluation and interpolation • Horner’s rule gives a Θ(n) algorithm for polynomial evaluation • but n evaluations would require Θ(n2) time • Given {(xi, yi)}, we can interpolate by first finding polynomials pj with value 1 on xj and 0 on the other xi. • each pj is of the form cjΠ(x-xi) over all i≠j • the desired polynomial would be Σyjpj • but again this would take time nΩ(n), or Ω(n2)

24. Divide-and-conquer to the rescue • It turns out that we can profit in each case from a divide-and-conquer approach • if we are careful about the choice of the xj • Consider evaluation of x3 + 2x2 + 3x + 4 at points 1 and -1. • This gives a value of 1 + 2 + 3 + 4 in the first case and -1 + 2 + -3 + 4 in the second. • Since -1 is a root of 1, much of the work in these two computations can be shared.

25. Nth roots of 1 • To generalize this to more than 2 points, we need to find and use nth roots of 1 for n>2. • there are n such nth roots that we can use • but they are complex numbers • all are of the form ωk for some specific ω • We must have n larger than the degree of any polynomial we need to represent. • The relevant properties of complex numbers and ω, are covered in the “slide” complex.html on the class web site.

26. The (discrete) Fast Fourier transform • For any polynomial p = Σaixi of degree n-1, let peven = Σa2ixi and podd = Σa2i+1xi. • Then both polynomials have degree n/2-1, • p(x) = peven(x2) + xpodd(x2), • and p(-x) = peven(x2) – xpodd(x2). • This reduces an evaluation of size n at 2 points to 2 evaluations of size n/2 at 1 point. • note: if x = ωk, then x2 =ω2k is an n/2th root of 1

27. Implementation of the FFT algorithm • This gives a divide-and-conquer algorithm: the (discrete) Fast Fourier transform (FFT). • There are several ways of implementing it. • one complication we’ll see soon: the FFT can be used for interpolation as well as evaluation • this can make it tricky for arithmetic operations to locate their operands quickly • Examples of two implementations are given in fft-examples.html.

28. Analysis of the discrete Fast Fourier transform • The FFT algorithm can be described as tree processing in at least two different ways • The first way gives the code of CLRS, p. 911 • The second way corresponds to my text • cf. Figures 4.45 and 4.46 of Sec. 4.17 of my text • In each case, processing time is Θ(n log n) • since there are Θ(log n) levels of the tree • and each level takes time Θ(n) to process • So the FFT has time complexity Θ(n log n) .

29. The FFT for interpolation • The FFT is very nearly its own inverse! • if we treat the FFT as an operator mapping a sequence (ai) to another sequence (yi) • So interpolation using the FFT is very much like evaluation.

30. Interpolation and repeated FFT application • Consider applying the FFT to its own output • (yj) = (p(ωj)) = (Σaiωji). • The resulting sequence has kth term • Σyjωjk = ΣΣai ωjkωji = ΣaiΣωj(k+i) = nan-k • since the last sum is n if i+k=n and 0 otherwise • So the inverse FFT needs only to • divide all values by n, and • reverse the order of this output sequence • perhaps by replacing k by –k

31. Another interpretation of the FFT • Suppose that we think of the FFT as an algorithm that transforms a column vector a into a column vector y • Here, the FFT computes y as F∙a, where F[i,j] = ωij • And the inverse FFT computes a as F-1∙y • where F-1[i,j] = ω-ij/n

32. Summary • The point of the FFT algorithm for us is mostly as an example of the advantages of a change of representation • although Fourier transforms have many applications • Pointers to applications of Fourier transforms are given at the beginning of Chapter 30, CLRS.