J2ee Web Hosting, Perl Scripting And Developing Blog

340 ARITHMETIC 4.5.3 Thus, for example, 1 1 -x2 /Xl/ = -7 /Xl, x2/ = (3) x1 21+1/x2 -x1×2+1 If n = 0, the symbol /xl,…, x,/ is taken to mean 0. Let us also define the polynomials Qn(2rr x2, . . . ,x,) of n variables, for n 2 0, by the rule 1, if n = 0; &7x(x1,22,. . . = XI, ifn=l; , xn) ~&n–1(x2,. . . ,xn)+&n–2(x3,…,xn), ifn > 1. (4) Thus &2(x1, x2) = xix2 + 1, Qs(xr, x2, xs) = ~1~2×3 +si +x3, etc. In general, as noted by L. Euler in the eighteenth century, Qn(xr , x2,. . . , 2,) is the sum of all terms obtainable by starting with ~1×2 . . . xn and deleting zero or more non- overlapping pairs of consecutive variables xjsj+l; there are F,+i such terms. The polynomials defined in (4) are called continuants. The basic property of the Q-polynomials is that /x1, x2,-.. , xn/ = &n–1(x2,. . . , xn)/Qn(xl, 22,. . . , xn), n > 1. (5) This can be proved by induction, since it implies that xo+/~1,…txn/ = Qn+l(xo, . . . ,xJ; XI,. . . ,xn)/Qn(xl, hence 1~0, xl,. . . , xn/ is the reciprocal of the latter quantity. The Q-polynomials are symmetrical in the sense that Qn(xl, x2,. . . , xn) = Qn(xn, . . . , x2, xl). (6) This follows from Euler s observation above, and as a consequence we have Qn(xl, . . . , xn) = xnQn-~(x~, . . . , G-I) + &n–2(x1,. . . , xn-2) (7) for n > 1. The Q-polynomials also satisfy the important identity Qn(xl, . . . , xn)Qn(xz,. . . , x,+1) -Qn+l(xl,. . . ,x,+1)&n-1(52,. . . , xn) = (-l) , n 2 1. (8) (See exercise 4.) The latter equation in connection with (5) implies that /xl,…,xn/= - 1 -L+L&…+(, :,; , 4041 q1q2 n n where qk = Qk(Zr,. . . , xk). (9)

4.53 ANALYSIS OF EUCLID S ALGORITHM 339 31. [AK% ] Use Euclid s algorithm to find a simple formula for gcd(2 -1,2 -1) when m and n are nonnegative integers. 32. [A&$3] Can the upper bound O(n ) in exercise 30 be decreased, if another algo- rithm for calculating the greatest common divisor is used? 33. [M@] Analyze V. C. Harris s binary Euclidean algorithm. b 34. [M3Z] (R. W. Gosper.) Demonstrate how to modify Algorithm B for large num- bers, using ideas analogous to those in Algorithm L. b 35. [MZ?] (V. R. Pratt.) Extend Algorithm B to an Algorithm Y that is analogous to Algorithm X. 36. [Hi!4491 Find a rigorous proof that Brent s model describes the asymptotic be- havior of Algorithm B. *4.5.3. Analysis of Euclid s Algorithm The execution time of Euclid s algorithm depends on T, the number of times the division step A2 is performed. (See Algorithm 4.5.2A and Program 4.5.2A.) The quantity T is also an important factor in the running time of other algorithms, such as the evaluation of functions satisfying a reciprocity formula (see Section 3.3.3). We shall see in this section that the mathematical analysis of this quantity T is interesting and instructive. Relation to continued fractions. Euclid s algorithm is intimately connected with continued fractions, which are expressions of the form bl = h/(al+bz/(az+h/(. . .l(an-l+bnl4. . .I,). bz al + b3 a2 + (1) bn %x-1+- an Continued fractions have a beautiful theory that is the subject of several books. [See, for example, 0. Perron, Die Lehre von den Kettenbriichen, 3rd ed. (Stutt- gart: Teubner, 1954), 2 ~01s.; A. Khinchin, Continued Fractions, tr. by Peter Wynn (Groningen: P. Noordhoff, 1963); H. S. Wall, Analytic Theory of Con- tinued Fractions (New York: Van Nostrand, 1948); and see also J. Tropfke, Geschichte der Elementar-Mathematik 6 (Berlin: Gruyter, 1924), 74-84, for the early history of the subject.] It is necessary to limit ourselves to a comparatively brief treatment of the theory here, studying only those aspects that give us more insight into the behavior of Euclid s algorithm. The continued fractions of primary interest to us are those in which all of the b s in (1) are equal to unity. For convenience in notation, let us define /51,~2,***, xn/ = l/(x1 + l/(x2 + I/(. . . (X7x-l + l/h) * .I,>. (2)

338 ARITHMETIC 45.2 20. [MZ6] Let X = [lggcd( lL,2))]. Sh ow that the lattice-point model gives X = 1 with probability 4, X = 2 with probability &, X = 3 with probability ha, etc., plus correction terms that go rapidly to zero as u and w approach infinity; hence the average value of X is approximately 5. [Hint: Consider the relation between the probability of a path from (m,n) to (k + 1, k + 1) and a corresponding path from (m -k,n -k) to (1, 11.1 21. [HMZ6] Let C,, and D,, be the average number of subtraction and shift cycles, respectively, in Algorithm B, when u and v are odd, Llg uJ = m, Llg w] = n. Show that for fixed n, C,, = $m + O(1) and D,, = m + O(1) as m + co. 22. [.%I] Prove Eq. (44). 23. [MB] Show that if C,, = cym + pn + 7 for some constants o, /3, and y, then c (N -m)(N -n)2m+n-2Cmn = 2 (#(a + @IV + O(l)), l

4.5.2 THE GREATEST COMMON DIVISOR 337 b 10. (HM.Z4] Let qn be the number of ordered pairs of integers (u, w) lying in the range 1 5 2~, v 2 n such that gcd(u, V) = 1. The object of this exercise is to prove that we have lim,,, q,/n2 = 6/7?, thereby establishing Theorem D. a) Use the principle of inclusion and exclusion (Section 1.3.3) to show that qn = n2 -ClnImJ2 + C ln/plp2j2-. . . , PI Pl

336 ARITHMETIC 4.5.2 This completes our analysis of the average values of C and D. The other three quantities appearing in the running time of Algorithm B are rather easily analyzed; see exercises 6, 7, and 8. Thus we know approximately how Algorithm B behaves on the average. Let us now consider a worst case analysis: What values of u and v are in some sense the hardest to handle? If we assume as before that LlguJ = m and 1kvJ = n, let us try to find (u, v) that make the algorithm run most slowly. In view of the fact that the subtractions take somewhat longer than the shifts, when the auxiliary bookkeeping is considered, this question may be rephrased by asking for the inputs u and v that require most subtractions. The answer is somewhat surprising; the maximum value of C is exactly max(m, n) + 1, (46) although the lattice-point model would predict that substantially higher values of C are possible (see exercise 26). The derivation of the worst case (46) is quite interesting, so it has been left as an amusing problem for the reader to work out by himself (see exercises 27 and 28). EXERCISES 1. [A&21] How can (8), (9), (lo), (11) and (12) be derived easily from (6) and (7)? 2. [A4.86] Given that u divides zlizl~ . . . v,, prove that u divides gcd(u, VI) gcd(u, ~2). . . gcd(u, v,). 3. [Mz%] Show that the number of ordered pairs of positive integers (u, v) such that lcm(u, v) = n is the number of divisors of n2. 4. [A&Z] Given positive integers u and v, show that there are divisors u of u and v of v such that gcd(u , v ) = 1 and u v = lcm(u, v). b 5. [AL%] Invent an algorithm (analogous to Algorithm B) for calculating the greatest common divisor of two integers based on their balanced ternary representation. Dem- onstrate your algorithm by applying it to the calculation of gcd(40902,24140). 6. [MZ?] Given that u and v are random positive integers, find the mean and the standard deviation of the quantity A that enters into the timing of Program B. (This is the number of right shifts applied to both u and v during the preparatory phase.) 7. [A4,20] Analyze the quantity B that enters into the timing of Program B. b 8. [MEY] Show that in Program B, the average value of E is approximately equal to +C,,,, where C,,, is the average value of C. 9. [18] Using Algorithm B and hand calculation, find gcd(31408,2718). Also find integers m and n such that 31408m + 2718n = gcd(31408,2718), using Algorithm X.

4.5.2 THE GREATEST COMMON DIVISOR 335 These experiments showed a rather small standard deviation from the observed average values. The coefficients $ and 1 of m in (42) and (43) can be verified rigorously without using the lattice-point approximation (see exercise 21); so the error in the lattice-point model is apparently in the coefficient of n, which is too high. The above calculations have been made under the assumption that u and v are odd and in the ranges 2m 5 u < 2m+1 and 2n 2 v < 2n+1. If we assume instead that 1~ and v are to be any integers, independently and uniformly distributed over the ranges 12 IL< 2N, 1 2 v < 2N, then we can calculate the average values of C and D from the data already given; in fact, if C,, denotes the average value of C under our earlier assumptions, exercise 22 shows that we have (zN -1)2C = N2Coo + 2N c (N -n)2n-1C,0 l 1, or by 2d/(X-1 -1) if X < 1. Th e random variable X has a limiting distribution that makes it possible to analyze the average value of the ratio by which max(u, v) decreases at each iteration; see exercise 25. Numerical calculations show that this maximum decreases by p = 0.705971246102 bits per iteration; the agreement with experiment is so good that Brent s constant p must be the true value of the number 0.70 in (45), and we should replace 0.203 by 0.206 in (43). [See Algorithms and Complexity, ed. by J. F. Daub (New York: Academic Press, 1976), 321-355.1

334 ARITHMETIC 4.5.2 Table 1 NUMBER OF SUBTRACTIONS IN ALGORITHM B n n 0 12 3 4 5 0 12 3 4 5 0 1.00 2.00 2.50 3.00 3.50 4.00 1.00 2.00 2.50 3.00 3.50 4.00 0 1 2.00 1.00 2.50 3.00 3.50 4.00 2.00 1.00 3.00 2.75 3.63 3.94 1 2 2.50 2.50 2.25 3.38 3.88 4.38 2.50 3.00 2.00 3.50 3.88 4.25 2 3 3.00 3.00 3.38 3.25 4.22 4.72 3.00 2.75 3.50 2.88 4.13 4.34 3 4 3.50 3.50 3.88 4.22 4.25 5.10 3.50 3.63 3.88 4.13 3.94 4.80 4 5 4.00 4.00 4.38 4.72 5.10 5.19 4.00 3.94 4.25 4.34 4.80 4.60 5 m Predicted by model Actual average values m execution of Algorithm B, with u and v both odd, the corresponding path in the lattice model must satisfy the relation number of shifts + number of diagonal points f 2Llg gcd(u, v)] = m + n, since we are assuming that T in (33) is always either 0 or 1. The average value of [lggcd(u, w)J predicted by the lattice-point model is approximately 2 (see exercise 20). Therefore we have, for the total number of shifts, D=m+n-(in+&++&,)–* (41) =m+jn–g-Zj3Smn, when m 2 n, plus terms that decrease rapidly to zero for large n. To summarize the most important facts we have derived from the lattice- point model, we have shown that if u and w are odd and if [lg u] = m, [lg v] = n, then the quantities C and D that are the critical factors in the running time of Program B will have average values given by C = $m + i$ + O(l), D = m -I- i$ -I- O(l), m >_ 72. (42) But the model that we have used to derive (42) is only a pessimistic approxima- tion to the true behavior; Table 1 compares the true average values of C, com- puted by actually running Algorithm B with all possible inputs, to the values predicted by the lattice-point model, for small m and n. The lattice model is completely accurate when m or n is zero, but it tends to be less accurate when (m - n( is small and min(m, n) is large. When m = n = 9, the lattice-point model gives C = 8.78, compared to the true value C = 7.58. Empirical tests of Algorithm B with several million random inputs and with various values of m, n in the range 29 5 m,n 5 37 indicate that the actual average behavior of the algorithm is given by CR5 irn + 0.203n+ 1.9 - 0.4(0.6)m-n, m 2 n. (43) DZ m + 0.41n -0.5 - 0.7(0.6)m-n,

4.5.2 THE GREATEST COMMON DIVISOR 333 Thus the double recurrence (34) can be transformed into the single recurrence relation in (37). Use of the generating function G(z) = & + Alz + A2z2 + . . . now transforms (37) into the equation (1-$z-&z2)G(z) = ao+aiz+a2z2+~ 1–+$1–z/2 A?-+ (l-z/2)2 7 (38) where a~, al, and a2 are constants that can be determined by the values of Ao, Al, and Aa. Since 1-$z + az2 = (1+ az)(l -z), we can express G(z) by the method of partial fractions in the form Tedious manipulations produce the values of these constants bo, . . . , bs, and thus the coefficients of G(z) are determined. We finally obtain the solution + (-:,Y-& -+@ + &c + ++bo) + &c&o; A,, = +mc + n($a + +A + (&a + $b + & + @oo) + 2- (44 + (-$,+&a -@ + $c + @bo), m > n. (39) With these elaborate calculations done, we can readily determine the be- havior of the lattice-point model. Assume that the inputs u and v to the algo- rithm are odd, and let m = LlguJ, n = [lg vJ. The average number of subtrac- tion cycles, namely the quantity C in the analysis of Program B, is obtained by setting a = 1, b = 0, c = 1, A00 = 1 in the recurrence (34). By (39) we see that (for m 2: n) the lattice model predicts subtraction cycles, plus terms that rapidly go to zero as n approaches infinity. The average number of times that gcd(u, v) = 1 is obtained by setting a = b = c = 0, A00 = 1; this gives the probability that u and v are relatively prime, approximately 2. Actually, since u and v are assumed to be odd, they should be relatively prime with probability 8/.rr2 (see exercise 13), so this reflects the degree of inaccuracy of the lattice-point model. The average number of times that a path from (m, n) goes through one of the diagonal points (m , m ) for some m 2 1 is obtained by setting a = 1, b = c = &O = 0 in (34); so we find that this quantity is approximately m n. *n + & + &L,, when 2 Now we can determine the average number of shifts, i.e., the number of times that step B3 is performed. (This is the quantity D in Program B.) In any

332 ARITHMETIC 4.5.2 An analysis of the lattice-point model can be carried out by solving the following rather complicated set of double recurrence relations: 1 &, = a + 5&qm–2) + … -I LA,o &, ifm 2 1; 2m-l + A mn = c + iAcrn-~jn + … + &&AI, + &hn, if m > 72 2 0; A mn - &V2, if n > m 2 0. (34) The problem is to solve for A,, in terms of m, n, and the parameters a, b, c, and Aoo. This is an interesting set of recurrence equations, which have an interesting solution. First we observe that if 0 5 rz < m, A(,+,),= c+ c 2-kA(m+,-k), + 2-m&n l 0. (35) Now let A,,, = A,,. If m > 0, we have A(m+qm= c+ c 2-k-+n+wc)m+ 2-m-1&m l_ 2, (37)

4.5.2 THE GREATEST COMMON DIVISOR 331 not on the actual values of u and v. Let us call this approximation a lattice-point model, since we will say that we are at the point (m, n) when Llg ~1 = m and [lg v] = n. From point (m, n) the algorithm takes us to (m , n) if u > U, or to (m,n ) if u < U, or terminates if u = v. For example, the calculation starting with u = 20451 and u = 6035 begins at the point (14,12), then goes to (9,12), P!ll), (%9), (4,9), (4,5), (4,4), and terminates. In line with the comments of the preceding paragraph, we will make the following assumptions about the probability that we reach a given point just after point (m, n): Case 1, m > n. Case 2, m < 12. Next point Probability Next point Probability Cm - Ln) % (m,n-1) a (m -2, n) a (m,n-2) a (rn,l) w-1. . h 0) w -1 Case 3, m = n > 0. Next point Probability Cm - 2, n), Cm, n -2) a, a Cm - 3, n), (3 n -3) QJ R (0, ni im, 0) c:,i,~id,m terminate w-1 For example, from points (5,3) the lattice-point model would go to points (4,3), (3,3), (2,3), (1,3), (0,3) with the respective probabilities 4, a, i, &, &; from (4,4) it would go to (2,4), (1,4), (0,4), (4,2), (4, I), (4,0), or would terminate, with the respective probabilities a, A, &, $, 4, &, i. When m and n are both 0, the formulas above do not apply; the algorithm always terminates in such a case, since m = n = 0 implies that u = v = 1. The detailed calculations of exercise 18 show that this lattice-point model is somewhat pessimistic. In fact, when m > 3 the actual probability that Algorithm B goes from (m, m) to one of the two points (m -2, m) or (m, m -2) is equal to i, although we have assumed the value 3; the algorithm actually goes from (m, m) to (m-3, m) or (m, m-3) with probability $, not a; and it actually goes from (m + 1, m) to (m, m) with probability 4, not 4. The probabilities in the model are nearly correct when Irn -n( is large, but when Irn -nl is small the model predicts slower convergence than is actually obtained. In spite of the fact that our model is not a completely faithful representation of Algorithm B, it has one great virtue: It can be completely analyzed! Furthermore, empirical experiments with Algorithm B show that the behavior predicted by the lattice- point model is analogous to the true behavior.