J2ee Web Hosting, Perl Scripting And Developing Blog

350 ARITHMETIC 4.5.3 With a slightly different choice of p, we can obtain tighter bounds (see ex- ercise 21). In fact, Wirsing went much further in his paper, proving that Fn(z) = lg(1 + Z) + (-x) *(Z) + O(Z(1 -Z)(X -o.o31)n), (38) where X = 0.30366 30028 98732 65860.. . = /3,3,2,2,3,13,1,174,1,1,1,2,2,2,1,1,1,2,2,1,. . . / (39) is a fundamental constant (apparently unrelated to more familiar constants), and where 9 is an interesting function that is analytic in the entire complex plane except for the negative real axis from -1 to –00. Wirsing s function satisfies 3(O) = Q(1) = 0, @ (O) < 0, and SQ = --Xk; thus by (35) it satisfies the identity Q(z) -qz+ 1) = ;Q -!-. ( 1+z > Furthermore, Wirsing demonstrated that a(-E+$)=ch- logN+O(l) asN-+co, where c is a constant and 72 = T(u, V) is the number of iterations when Euclid s algorithm is applied to the integers u > v > 0. A complete solution to Gauss s problem was found a few years later by K. I. Babenko [Doklady Akad. Nauk SSSR 238 (1978), 1021-10241, who used powerful techniques of functional analysis to prove that for all 0 5 z 5 1, 72 2 1. Here 1×21 > [X3) 2 1X4( > .. ., and each @j(z) is an analytic function in the complex plane except for a cut at [-00, -11. The function 92 is Wirsing s k, and Xz = –X, while X3 = 0.1009, X4 = -0.0408, Xs = -0.0355, Xs = 0.0128. Babenko also established further properties of the eigenvalues Xj, proving in particular that they are exponentially small as j + 00, and that the sum for j 2 k in (42) is bounded by (n2/6)]_Xk]n-1 min(s, 1 - z). [Further information appears in a paper by Babenko and Iur ev, Doklady Akad. Nauk SSSR 240 (1978), 1273-1276.1 F rom continuous to discrete. We have now derived results about the prob- ability distributions for continued fractions when X is a real number uniformly distributed in the interval [ 0,l). But a real number is rational with probability zero (almost all numbers are irrational), so these results do not apply directly to Euclid s algorithm. Before we can apply Theorem W to our problem, some technicalities must be overcome. Consider the following observation based on elementary measure theory:

45.3 ANALYSIS OF EUCLID S ALGORITHM 349 Thus everything hinges on our being able to prove that Un produces small function values, where U is the linear operator defined in (29). Note that U is a positive operator, in the sense that Ucp(x) 2 0 for all x if (o(x) 2 0 for all x. It follows that U is order-preserving: If (01(x) 5 ps(x) for all x then we have V(oi(x) 5 Upz(x) for all x. One way to exploit this property is to find a function yo for which we can calculate Up exactly and to use constant multiples of this function to bound the ones that we are really interested in. First let us look for a function g such that Tg is easy to compute. If we consider functions defined for all x > 0, instead of only on [0, 11, it is easy to remove the summation from (25) by observing that SG(x+l)-SG(x) = G(&)-;A&&) =$$+W (35) when G is continuous. Since T((l + x)G ) = (1 + x)(SG) , it follows (see exercise 20) that –Tdx) Tdx+ 1) = -&)g(&>* (36) x+1 (& x+2 If we set Tg(x) = l/(x + l), we find that the corresponding value of g(x) is 1 -j-x -l/(1 + x). Let p(x) = g (x) = 1 + l/(1 + x)~, so that Up(x) = -(Tg) (x) = l/(1 + x)~; th is is the function (o we have been looking for. For this choice of cp we have 2 2 p(x)/Up(x) = (1 + x)~ + 1 5 5 for 0 5 x 5 1, hence i(P I UP I 3P. By the positivity of U and p we can apply U to this inequality again, obtaining &p < 4Up 2 U2p < 3Up 5 $(p; and after n - 1 applications we have 5- (0 5 unp 2 2-Q (37) for this particular (o. Let x(x) = &(x) = 1 be the constant function; then for 0 5 x 2 1 we have ix < (o 5 2x, hence It follows by (33) that s #(ln2)25-n 5 (-l)nRK(x) 2 +j(ln2)22-n, for 0 < 2 5 1; hence by (30) and (34) we have proved the following result: Theorem W. The distribution F,(x) equals lg(1 + x) + 0(2- ) as n –+ 00. In fact, F,(x)-lg(1 -j-x) lies between &(-l)nf15-n(ln(l+x))(ln2/(1+x)) and ~(-l)n+12-n(ln(l + x))(ln2/(1+ x)),, for 0 2 x 5 1. 1

348 ARITHMETIC 4.5.3 Continuing, we see that if g has a bounded first derivative, we can differen- tiate term by term to show that Tg does also: (TgY(x) x)2 &) = - kFl ((k + 1 + - (L;x;Jg( - -+k;x,2g (i&)) =- (g(&)-g(k+:+x)) + (k + x):(:; 1 + 2) There is consequently a third linear operator, U, such that (Tg) = -U(g ), namely = dt UP(X)c k v(t) (k + 1 + x) k>l + (k+x)i(:;l+x) . (29) What is the relevance of all this to our problem? Well, if we set Fn(x) = w + x) + q@ + xl), (30) fn(x) = cl+ x) F:,(x) = &Cl + qdl+ x,)), (31) we have f;(x) = RE(Ml + x))/((ln 2) P + x)); (32) the effect of the lg(1 + x) term disappears, after these transformations. F urther- more since Fn = PFo we have fn = Tnfo and f; = (-l)nU fb. Both F, and fn have bounded derivatives, by induction on 72. Thus (32) becomes (-l)nRK(lg(l + x)) = (1 + x)(ln 2)2 U fb(x). (33) Now Fe(x) = x, fo(x) = 1 + x, and fb(x) is the constant function 1. We are going to show that the operator U takes the constant function into a function with very small values, hence (Rx(x)1 must be very small for 0 5 x 5 1. Finally we can clinch the argument by showing that R,(x) itself is small: Since we have R,(O) = Rn(l) = 1, it follows from a well-known interpolation formula (cf. exercise 4.6.4-15 with x0 = 0, x1 = x, x2 = 1) that R,(x) = - ( ; ) R;@(x)) (34) for some function c(x), where 0 5 E(x) 5 1 when 0 2 x 2 1.

4.5.3 ANALYSIS OF EUCLID S ALGORITHM 347 Gauss s problem, namely to find the asymptotic behavior of Fn(x) -lg(1 + x), was not really resolved until 1974, when Eduard Wirsing published a beautiful analysis of the situation [Acta Arithmetica 24 (1974), 507-5281. We shall study the simplest aspects of Wirsing s approach here, since his method is an instructive use of linear operators. If G is any function of x defined for 0 5 x 5 1, let SG be the function defined by -(25) SG(x) = kg (G(i) - G( &))- Thus, S is an operator that changes one function into another. In particular, by (23) we have F,+l(x) = SF,(x), hence F, = SnFo. (26) (In this discussion F, stands for a distribution function, not for a Fibonacci number.) Note that S is a linear operator ; i.e., S(cG) = c(SG) for all con- stants c, and S(G1 + Ga) = SG1 + SG2. Now if G has a bounded first derivative, we can differentiate (25) term by term to show that (SG) (x)= k-l (k ; x)2 G( &): (27) - hence SG also has a bounded first derivative. (Term-by-term differentiation of a convergent series is justified when the series of derivatives is uniformly convergent; cf. K. Knopp, Theory and Application of Infinite series (Glasgow: Blackie, 1951), $47.) Let H = SG, and let g(x) = (1 + x)G (x), h(x) = (1 + x)H (x). It follows that h(x)=c -1+x (1+&o) k>l (k+xY - In other words, h = Tg, where T is the linear operator defined by (28)

346 ARITHMETIC 4.5.3 We might conjecture, for example, that F(i) = lg($) z 0.58496; let us see how close F,( 4) comes to this value for small n. We have F (i) = 4, and F&g= f - &+;-++..- 2 f =2 ( Li+L!+… 1= 2(1 -ln2) z 0.6137; E&H = & ; ($p -& + & -L5m+2 + * * -> - 2 1 = -L+L… m2 2 3 4 C-C TTL>l 4 1 1 +… -c m>l ii ( 2m(2m + 2) -3m(3m + 2) > - = i7rz(1 -ln2)-C %, ?7L>l where S, = 1/(4m + 4) -1/(9m + 6) + 1/(16m + -. . . . Using the values of H, for fractional z found in Table 3 of Appendix B, we find that s1= A, S2 = 4 -In 2, s3 = # -T/(2&, etc.; a numerical evaluation yields Fz(i) ==: 0.5748. Although F1(z) = Hz, it is clear that F,(z) is difficult to calculate exactly when n is large. The distributions F,(z) were first studied by K. F. Gauss, who thought of the problem in 1800. His notebook for that year lists various recurrence relations and gives a brief table of values, including the four-place value for F2( 4) that has just been mentioned. After performing these calculations, Gauss wrote, Tam complicate evadunt, ut nulla spes superesse videatur, i.e., They come out so complicated that no hope appears to be left. Twelve years later, Gauss wrote a letter to Laplace in which he posed the problem as one he could not resolve to his satisfaction. He said, I found by very simple reasoning that, for n infinite, F,(s) = log(1 + x)/log 2. But the efforts that I made since then in my inquiries to assign F,(s) -log(1 + x)/log2 for very large but not infinite values of n were fruitless. He never published his very simple reasoning, and it is not completely clear that he had found a rigorous proof. More than 100 years went by before a proof was finally published, by R. 0. Kuz min [Atti de1 Congress0 internazionale dei matematici 6 (Bologna, 1928), 83-891, who showed that F,(z) = lg(1 + X) + O(eCAfi) for some positive constant A. The error term was improved to O(ewAn) by Paul L&y shortly afterward [Bull. Sot. Math. de fiance 57 (1929), 178-194]*; but *An exposition of L&y s interesting proof appeared in the first edition of this book.

4.5.3 ANALYSIS OF EUCLID S ALGORITHM 345 value of n mod k is less than the average value of i k, in the range 1 5 k 5 n: En— 1 1441 + 1 n c( 2 > 1141n = 1-g n+O(logn) (21) ( > (cf. exercise 4.5.2-10(c)). This is only about .1775n, not .25n; so the value of nmod k tends to be smaller than the above model predicts, and Euclid s algorithm works faster than we might expect. A continuous model. The behavior of Euclid s algorithm with v = N is essen- tially determined by the behavior of the regular continued fraction process when x = O/N, l/N, . . . , (N -1)/N. Assuming that N is very large, we are led naturally to a study of regular continued fractions when X is a random real number uniformly distributed in [ 0,l). Therefore let us define the distribution function F,(z) = probability that X, 5 5, for 0 5 5 5 1, (22) given a uniform distribution of X = X0. By the definition of regular continued fractions, we have Fe(z) = z, and F,+i(z) = c probability that (k 2 l/Xn 5 k + x) k>l = c probability that (l/(k + CT) 5 X, 5 l/k) k>l = c (%(1/k) -%(ll(k + 4)). (23) k>l If the distributions Fe(s), Fi(z), . . . defined by these formulas approach a limiting distribution F,(z) = F(z), we will have F(z) = c (W/k) -J (ll(k + 4)). k>l One function that satisfies this relation is F(z) = log,(l + z), for any base b > 1; see exercise 19. The further condition F(1) = 1 implies that we should take b = 2. Thus it is reasonable to make a guess that F(z) = lg(1 + z), and that F,(z) approaches this behavior.

344 ARITHMETIC 4.5.3 An approximate model. Now that we know the maximum number of division steps that can occur, let us attempt to find the average number. Let T(m, n) be the number of division steps that occur when u = m and v = n are input to Euclid s algorithm. Thus T(m, 0) = 0; T(m, n) = 1 + T(n, m mod n) ifn 2 1. (18) Let T, be the average number of division steps when v = n and when u is chosen at random; since only the value of umodv affects the algorithm after the first division step, we may write T, = k c T(k,n). (19) O 1. (20) (This approximation is analogous to the lattice-point model used to investigate Algorithm B in Section 4.5.2.) The recurrence (20) is readily solved by the use of generating functions. A more direct way to solve it, analogous to our solution of the lattice-point model, is by noting that Sn+1 = 1+ -& (So + Sl + * * * + sn-1 + Sn) = 1+ & (n(% -1) + Sn) = Sn + &; hence S, is 1 + JJ + . . . + 6 = H,, a harmonic number. The approximation T12 M S, now suggests that T, z Inn + O(1). Comparison of this approximation with tables of the true value of T, show, however, that Inn is too large; T, does not grow this fast. One way to account for the fact that this approximation is too pessimistic is to observe that the average

4.5.3 ANALYSIS OF EUCLID S ALGORITHM 343 by (5)1 _ = . . ,A,) 2, &n–1(&, (16) U Qn(Al,Aa,. . . ,A,) This formula holds also if Euclid s algorithm is applied for u < U, when Al = 0. Furthermore, because of (8), Qn-1(A2, . . . , A,) and Qn(A1, Aa, . . , A,) are rela- tively prime, and the fraction on the right-hand side of (16) is in lowest terms; therefore u=Qn(Al,Aa ,..., A&, 2, = &n-l&, . . . , A&, (17) where d = gcd(u, v). The worst case. We can now apply these observations to determine the behavior of Euclid s algorithm in the worst case, or in other words to give an upper bound on the number of division steps. The worst case occurs when the inputs are consecutive Fibonacci numbers: Theorem F (G. Lamk, 1845). For n 2 1, let u and u be integers with u > w > 0 such that Euclid s algorithm applied to u and v requires exactly n division steps, and such that u is as small as possible satisfying these conditions. Then u = F n+2 and w = F,+l. Proof. By (17), we must have u = Qn(A1,A2,. . . ,A,)d, where Al, A2, . . . , &, and d are positive integers and A, 2 2. Since Qn is a polynomial with nonnegative coefficients, involving all of the variables, the minimum value is achieved only when Al = 1, . . . , A,-1 = 1, A, = 2, d = 1. Putting these values in (17) yields the desired result. 1 (This theorem has the historical claim of being the first practical application of the Fibonacci sequence; since then many other applications of Fibonacci numbers to algorithms and to the study of algorithms have been discovered.) As a consequence of Theorem F we have an important corollary: Corollary L. If 0 2 u, w < N, the number of division steps required when Algorithm 4.5.2A is applied to u and u is at most [log+(& N)] -2. Proof. By Theorem F, the maximum number of steps, n, occurs when u = F, and v = F%+l, where n is as large as possible with F,+l < N. (The first division step in this case merely interchanges u and u when n > 1.) Since F,+, < N, we have @-l/v% < N (see Eq. 1.2.8-15), so n+ 1 < log+(& N). This completes the proof. 1 Note that log+(& N) is approximately 2.078 In N + 1.672 M 4.785 log,, N + 1.672. See exercises 31 and 36 for extensions of Theorem F.

342 ARITHMETIC 4.5.3 some familiar real numbers, we find, for example, that z&3 = /3,1,1,1,2/; &=I 1, 1,9,2,2,3,2,2,9,1,2,1,9,2,2,3,2,2,9,1,2,1,9,2,2,3,2,2,9,1, . . . /; i% = 1 + /3,1,5,1,1,4,1,1,8,1,14,1,10,2,1,4,12,2,3,2,1,3,4,1,1,2,14,3,. . . /; T = 3 + /7,15,1,292,1,1,1,2,1,3,1,14,2,1,1,2,2,2,2,1,84,2,1,1,15,3,13,. . . /; e = 2 + /l, 2,1,1,4,1,1,6,1,1,8,1,1,10,1,1,12,1,1,14,1,1,16,1,1,18,1,. . . /; 7 = /I, 1,2,1,2,1,4,3,13,5,1,1,8,1,2,4,1,1,40,1,11,3,7,1,7,1,1,5,. . . /; f#J= 1 + /I, 1, 1, 1, 1, 1, 1, 1, 1, . . . /. (13) The numbers AI, AZ, . . . are called the partial quotients of X. Note the regular pattern that appears in the partial quotients for $@9, I$, and e; the reasons for this behavior are discussed in exercises 12 and 16. There is no apparent pattern in the partial quotients for fi, 7r, or 7. It is interesting to note that the ancient Greeks first definition of real numbers, once they had discovered the existence of irrationals, was essentially stated in terms of infinite continued fractions. (Later they adopted the suggestion of Eudoxus that 5 = y should be defined instead as Z < T if and only if y < r, for all rational r. ) See 0. Becker, Quellen und Studien zur Geschichte Math., A&on., Physik (B) 2 (1933), 311-333. When X is a rational number, the regular continued fraction corresponds in a natural way to Euclid s algorithm. Let us assume that X = v/u, where u > 21 2 0. The regular continued fraction process starts with X0 = X; let us define U, = u, Vo = u. Assuming that X, = l&/U, # 0, (10) becomes A n+1 = LG./KLl, X n+l = Un/Vn -&+I = (U, modV,)/V,. (14) Therefore, if we define u n+1 –V,, Vn+l = U,modV,, (15) the condition X, = V,/U,, holds throughout the process. Furthermore, (15) is precisely the transformation made on the variables u and ZJ in Euclid s algorithm (see Algorithm 4.5.2A, step A2). For example, since & = /3,1,1,1,2/, we know that Euclid s algorithm applied to u = 29 and v = 8 will require exactly five division steps, and the quotients lu/vj in step A2 will be successively 3, 1, 1, 1, and 2. Note that the last partial quotient A, must be 2 or more when X, = 0 and n 2 1, since X,-r is less than unity. From this correspondence with Euclid s algorithm we can see that the regular continued fraction for X terminates at some step with X, = 0 if and only if X is rational; for it is obvious that X, cannot be zero if X is irrational, and, conversely, we know that Euclid s algorithm always terminates. If the partial quotients obtained during Euclid s algorithm are Al, AZ, . . . , A,, then we have,

4.5.3 ANALYSIS OF EUCLID S ALGORITHM 341 Thus the &-polynomials are intimately related to continued fractions. Every real number X in the range 0 5 X < 1 has a regular continued fraction defined as follows: Let X0 = X, and for all n 2 0 such that X, # 0 let A n+l = 11/-&l, X,+1 = l/-G –A,+,. (10) If X, = 0, the quantities A,+1 and X,+1 are not defined, and the regular continued fraction for X is /Al, . . . , A,/. If X, # 0, this definition guarantees that 0 2 Xntl < 1, so each of the A s is a positive integer. The definition (10) clearly implies that x=x,= l 1 AI +X1 =Ai+l/(As+Xs) = hence X=/Al,… ,&-I,& f-G/ (11) for all n 2 1, whenever X, is defined. In particular, we have X = /Al, . . . , An/ when X, = 0. If X, # 0, the number X always lies between /Al,. . . , An/ and I&,.. . , A, + l/, since by (7) the quantity qn = Qn(A1, . . . , A, +Xn) increases monotonically from Qn(A1, . . . , &) up to Qn(A1, . . . , A, + 1) as X, increases from 0 to 1, and by (9) the continued fraction increases or decreases when qn increases, according as n is even or odd. In fact, IX-/Al,…,&/1 = l/Al ,…, &+X,/-/A, ,…, An/l = l/Al,… , Ant l/-G/ -/AI,. . . , An/l Q&b,. . . ,A,, l/-G) &n-1(4 . . . , An) = Qn+l(& . . . , An, l/X) -Q&h,. . . , A,) = l/(&&b,. . . ,An)Qn+l(Al,. . . ,A,, l/X)) 5 ll(Qn(Al, . . . ,4JQn+l(Al,. . . , An, &+I)) (12) by (5), (7), (8), and (10). Therefore /AI,. . . ,An/ is an extremely close approx- imation to X, unless n is small. If X is irrational, it is impossible to have X, = 0 for any n, so the regular continued fraction expansion in this case is an infinite continued fraction /Al, AZ, As, . . . /. The value of an infinite continued fraction is defined to be lim /Al,&…,AJ, n+cc and from the inequality (12) it is clear that this limit equals X. The regular continued fraction expansion of real numbers has several prop- erties analogous to the representation of numbers in the decimal system. If we use the formulas above to compute the regular continued fraction expansions of