3.1 ANSWERS TO EXERCISES 519 where Q(m) is defined in Section 1.2.11.3, Eq. (2). By Eq. (25) in that section, the probability is approximately 6&7% z 1.25/fi. The chance of Algorithm K converging as it did is only about one in 80000; the author was decidedly unlucky. But see exercise 15 for further comments on the colossalness. 12. c XP(b,X)=;(1+3(1-;)+6(1-;)(+2-)+..j l
518 ANSWERS TO EXERCISES 3.1 Notes: Brent has also considered a more general method in which the successive values of Y = X,% satisfy ni = 0, ni+l = 1 + [pni] where p is any number greater than 1. He showed that the best choice of p, approximately 2.4771, saves about 3 percent of the iterations by comparison with p = 2. (See exercise 4.5.4-4.) The method in part (b) has a serious deficiency, however, since it might generate a lot of nonrandom numbers before shutting off. For example, we might have a particularly bad case such as X = 1, p = 2k. A method based on Floyd s idea in exercise 6(b), namely one that maintains Y = X sn and X = X, for n = 0, 1, 2, 9 will require a few more function evaluations than Brent s method, but it will stop before any number has been output twice. On the other hand if f is unknown (e.g., if we are receiving the values Xc, Xi, . . on-line from an outside source) or if f is difficult to apply, the following cycle detection algorithm due to R. W. Gosper will be preferable: Maintain an auxiliary table To, Tl, I T,, where m = Llg n] when receiving X,. Initially TO + X0. For n = 1, 2, .I compare X, with each of TO, . . . , Tllgnl; if no match is found, set T,(,) +-X,, where e(n) = max{ e 1 2e divides n + 1 }. But if a match X, = Tk is found, then X = n -max{ m ] m < n and e(m) = k}. This procedure does not stop before generating a number twice, but at most [3x] of the X s will have occurred more than once. [MIT AI Laboratory Memo 239 (Feb. 29, 1972), Hack 132.1 See also R. Sedgewick and T. G. Szymanski, Proc. ACM Symp. Th. Comp. 11 (1979) 74-80. 8. (a,b) OO,OO,. . . [62 starting values]; lO,lO,. . [19]; 60,60,. [15]; 50,50,. . [l]; 24,57,24,57,. . . [3]. (c) 42 or 69; these both lead to a set of fifteen distinct values, namely (42 or 69), 76, 77, 92, 46, 11, 12, 14, 19, 36, 29, 84, 05, 02, 00. 9. Since X < b , we have X2 < bzn, and the middle square is [X2/bnj 5 X2/b . If X > 0, then X2/bn < Xb lb = X. 10. If X = abn, the next number of the sequence has the same form; it is (a2 mod b )b . If a is a multiple of b or of all the prime factors of b, the sequence will soon degenerate to zero; if not, the sequence will degenerate into a cycle of numbers having the same general form as X. Further facts about the middle-square method have been found by B. Jansson, Random Number Generators (Stockholm: Almqvist & Wiksell, 1966), Section 3A. Numerologists will be interested to learn that the number 3792 is self-reproducing in the four-digit middle-square method, since 3792 = 14379264; furthermore (as Jansson has observed), it is self-reproducing in another sense, too, since its prime factorization is 3 79.24! 11. The probability that X = 1 and p = 0 is the probability that Xi = Xc, namely l/m. The probability that X = 1,~ = 1, or X = 2,~ = 0, is the probability that XI # Xc and that XZ has a certain value, so it is (1 -l/m)(l/m). Similarly, the probability that the sequence has any given ,U and X is a function only of p + X, namely P(P,1) = ; l
3.1 ANSWERS TO EXERCISES 517 (c) The markings on the faces will slightly bias the die, but for practical purposes this method is quite satisfactory (and it has been used by the author in the preparation of several examples in this set of books). See Math. Comp. 15 (1961), 94-95, for further discussion of these dice. (d) (This is a hard question thrown in purposely as a surprise.) The number is not quite uniformly random. If the average number of emissions per minute is m, the probability that the counter registers k is e- mk/k! (the Poisson distribution); so the digit 0 is selected with probability e- c,,, mlok /(lOk)!, etc. The units digit will be even with probability e- cash m = 4 + ieZ2n, and this is never equal to i (although the error is negligibly small when m is large). (e) Okay, provided that the time since the previous digit selected in this way is random. However, there is possible bias in borderline cases. (f,g) No, people usually think of certain digits (like 7) with higher probability. (h) Okay; your assignment of numbers to the horses had probability & of assigning a given digit to the winning horse. 2. The number of such sequences is the multinomial coefficient 1000OOO!/(lOOOOO!) ; the probability is this number divided by 101oooooo, the total number of sequences of a million digits. By Stirling s approximation we find that the probability is close to 1/(167r410z2fi) x 2.55 X 10-26, about one chance in 4 X 10z5. 3. 3040504030. 4. Step Kll can be entered only from step KlO or step K2, and in either case we find it impossible for X to be zero by a simple argument. If X could be zero at that point, the algorithm would not terminate. 5. Since only lOlo ten-digit numbers are possible, some value of X must be repeated during the first 10 + 1 steps; and as soon as a value is repeated, the sequence continues to repeat its past behavior. 6. (a) Arguing as in the previous exercise, the sequence must eventually repeat a value; let this repetition occur for the first time at step p + X, where X,+X = X,. (This condition defines p and X.) We have 0 2 p < m, 0 < X 5 m, p + X 5 m. The values p = 0, X = m are attained iff f is a cyclic permutation; and p = m -1, X = 1 occurs, e.g., if X0 = 0, f(z) = z $ 1 for 2 < m -1, and f(m -1) = m -1. (b) We have, for r 2 n, X, = X, iff r -n is a multiple of X and n 2 p. Hence X2, = X, iff n is a multiple of X and n 2 p. The desired results now follow immediately. [Note: This is essentially a proof of the familiar mathematical result that the powers of an element in a finite semigroup include a unique idempotent element: take Xc = a, f(x) = ax.1 (c) Once n has been found, generate X, and X,+, for i > 0 until first finding x = Xntz; then p = i. If none of the values of Xn+i for 0 27 5 p is equal to X,, it follows that X = n, otherwise X is the smallest such i. 7. (a) The least n > 0 such that n -(e(n) -1) is a multiple of X and e(n) -1 2 p is n = 2r gmax(~f1,X)1 -1 + X. [This may be compared with the least n > 0 such that Xz,, = X,, namely X6,0 + p + X - 1 - ((p + X - 1) mod X).] (b) Start with X = Y = Xc, k = m = 1. (At key places in this algorithm we will have X = Xzm-k-1, Y = X,-l, and m = e(2m -k).) To generate the next random number, do the following steps: Set X + f(X) and k + k -1. If X = Y, stop (the period length X is equal to m-k). Otherwise if k = 0, set Y + X, m + 2m, k + m. Output X.
ANSWERS TO EXERCISES This branch of mathematics is the only one, r believe, in which good writers frequently get results entirely erroneous. It may be doubted if there is a single extensive treatise on probabilities in existence which does not contain solutions absolutely indefensible. -C. S. PEIRCE, in Popular Science Monthly (1878) NOTES ON THE EXERCISES 1. An average problem for a mathematically inclined reader. 3. See W. J. LeVeque, Topics in Number Theory 2 (Reading, Mass.: Addison-Wesley, 1956), Chapter 3. [Note: One of the men who read a preliminary draft of the manuscript of this book reported that he had discovered a truly remarkable proof, which the margin of his copy was too small to contain.] SECTION 3.1 1. (a) This will usually fail, since round telephone numbers are often selected by the telephone user when possible. In some communities, telephone numbers are perhaps assigned randomly. But it would be a mistake in any case to try to get several successive random numbers from the same page, since the same telephone number is often listed several times in a row. (b) But do you use the left-hand page or the right-hand page? Say, use the left- hand page number, divide by 2, and take the units digit. The total number of pages should be a multiple of 20; but even so, this method will have some bias. 516
4.7 MANIPULATION OF POWER SERIES 515 12. [MZU] Find a connection between polynomial division and power series division: Given polynomials U(X) and w(z) of respective degrees m and n over a field, show how to find the polynomials q(x), r(z) such that U(X) = q(z)v(z) + r(z) and deg(r) < n, using only operations on power series. 13. [A&7] (Rational function approximation.) It is occasionally desirable to find polynomials whose quotient has the same initial terms as a given power series. For example, if W(z) = 1 + z + 32 + 7z3 f . , there are essentially four different ways to express W(z) as w~(z)/w~(z) + O(z4) w h ere wi(z) and ws(z) are polynomials with deg(wi) + deg(wz) < 4: (1 + z + 3z2 + 7z3) / 1 = 1 + z + 3z2 + 7z3 + o.z4 + . ) (3 -42 + 222) / (3 -72) = 1 + z + 3z2 + 72s + Yz + . , (1 -2) / (1 -22 -z ) = 1 + z + 3z2 + 7z3 + 17z4 + , l/(1-z-zz2-z+=1+Z+3z2+723+15z4+…. Rational functions of this kind are commonly called Pad6 approximations, since they were studied extensively by H. E. Pade [Annales Scient. de J&oJe Normale SupCrieure (3) 9 (1892), Sl-S93]. Show that all Pad6 approximations W(z) = w~(z)/w~(z) + O(zN) with deg(wi) f deg(ws) < N can be obtained by applying an extended Euclidean algorithm to the polynomials zN and WO+WI~+…+WN-IZ~- ; and design an all-integer algorithm for the case that each W, is an integer. [Hint: See exercise 4.6.1-26.1 b 14. [HMXI] Fill in the details of Brent and Traub s method for calculating U (z) when U(z) = z + Uk zk + . . , using (27) and (28). And it shall be, when thou hast made an end of reading this book, that thou shalt bind a stone to it, and cast it Into the midst of Euphrates. -Jeremiah 51:63
514 ARITHMETIC 4.7 Algebraic functions. The coefficients of each power series W(z) that satisfies a general equation of the form An(z)W(z)n + . . . + Al(z)Wz) + Adz) = 0, (29) where each A,(z) is a polynomial, can be computed efficiently by using methods due to H. T. Kung and J. F. Traub; see JACM 25 (1978), 245-260. EXERCISES 1. [M10] The text explains how to divide V(z) by V(z) when VO # 0; how should the division be done when VO = O? 2. [Zoo] If the coefficients of U(z) and V(z) are integers and Vo # 0, find a recurrence relation for the integers Vt+ W,, where W, is defined by (3). How would you use this for power series division? 3. [Ml51 Does formula (9) give the right results when cx = O? When (Y = l? b 4. [HA&?3] Show that simple modifications of (9) can be used to calculate evcz) and ln(1 + V(z)), when V(z) = VIZ + vZz2 + … . 5. [MOO] What happens when a power series is reverted twice; i.e., if the output of Algorithm L or T is reverted again? b 6. [k%] (H. T. Kung.) Apply Newton s method to the computation of W(z) = l/V(z), when V(0) # 0, by finding the power series root of the equation f(z) = 0, where f(z) = z-l -V(z). 7. [MB ] Use Lagrange s inversion formula (12) to find a simple expression for the coefficient W, in the reversion of z = t -P. b 8. (A&Z] Lagrange s inversion formula can be generalized as follows: If W(z) = wlz+wzz2+~~~ = Glt + Gzt + G3t3 + . . . = G(t), where z and t are related by Eq. (lo), then W, = U,-l/n where Uo +U~z+Uzz2 +… = (Gl +2G2z+3G3z2 +…)(l +vZz+~z2 +…)- . (Equation (12) is the special case G1 = 1, Gz = Gs = . . . = 0. This equation can be proved, for example, by using tree-enumeration formulas as in exercise 2.3.4.4-33.) Extend Algorithm L so that it obtains the coefficients WI, W2, . . . in this more general situation, without substantially increasing its running time. 9. [11] Find the values of T,, computed by Algorithm T as it determines the first five coefficients in the reversion of z = t -t2. 10. [A&%] Given that y = x + alxoL+l + u2xa+ + . .. , cy # 0, show how to compute the coefficients in the expansion x = y + bzy21a + b3y3/ + . b 11. [M,Z5] (Composition of power series.) Let U(z)=U0+Ulz+U2z2+~~~ and V(z)=Vlz+~z2+V3z3+~~~; design an algorithm that computes the first N coefficients of U(V(z)).
4.7 MANIPULATION OF POWER SERIES 513 The running time T(n) of this procedure satisfies T(2n) = 2T(n) + C(n), (26) where C(n) is the time to compute R(z), I@( z ), and g(z). The function C(n) is dominated by the time to compute V(U(z)) modulo z2n, and C(n) presumably grows faster than order nl+ ; therefore the solution r(n) to (26) will be of order C(n). For example, if C(n) = cn3 we have T(n) z $cn3; or if C(n) is O(nlogn)3/2 using fast composition, we have T(n) = O(nlogn)3/2. The procedure breaks down when W(0) = 1 and S(0) # 0, so we need to investigate when this can happen. It is easy to prove by induction on n that the solution of (22) by the Brent-Traub method entails consideration of exactly n subproblems, in which the coefficient of V(z) on the right-hand side takes the respective values W(z)(z/U(z))j for 0 5 j < n in some order. If W(0) = U1 and if U1 is not a root of unity, we therefore have W(0) = 1 only when j = 1; the procedure will fail in this case only if (22) has no solution for n = 2. The Schrijder function for U can therefore be found by solving (22) for n = 2, 4, 8, 16, . . , with W(z) = U 1 and S(z) = 0, whenever VI is nonzero and not a root of unity. If Ul = 1, there is no Schrader function unless U(z) = z. But Brent and Traub have found a fast way to compute P(z) even when Ul = 1, by making use of a function V(z) such that V(U(z)) = U (z)V(z). (27) If two functions U(z) and o(z) both satisfy (27), for the same V, it is easy to check that their composition U(o(z)) does too; therefore all iterates of U(z) are solutions of (27). Suppose that U(z) = z + Uk zk + Uk+lzk+l + ... where k 2 2 and uk # 0. Then it can be shown that there is a unique power series of the form V(z) = zk + Vk+lzk+ + Vk+2zk+2 +. . . satisfying (27). Conversely if such a function V(z) is given, and if k 2 2 and uk are given, then there is a unique power series of the form U(z) = z + uk zk + uk+lz + +. . . satisfying (27); the desired iterate P(Z) is the unique power series P(z) satisfying V(P(z)) = P (z)V(z) (28) such that P(z) = z + nukzk + ... . Both V(z) and P(z) can be found by appropriate algorithms (see exercise 14). If UI is a kth root of unity, but not equal to 1, the same method can be applied to the function Uk(z) = z + ... , and Uk(z) can be found from U(z) by doing l(k) composition operations (cf. Section 4.6.3). We can also handle the case VI = 0: If U(z) = Ukzk + Uk+lzk+ + ... where k > 2 and uk # 0, the idea is to find a solution to the equation V(U(z)) = ukV(qk; then Un(z) = ~-1(u)ck -W(k-1) V(Z)~ ). Finally, if U(z) = UO + UIZ + … where UO # 0, let cr be a fixed point such that U(a) = e, and let o(z) = U(a + z) -a = zU (a) + zV ((Y)/2! + . . . ; then P(Z) = Un(z -0) + a. Further det#ails can be found in Brent, and Traub s paper [SIAM J. Computing 9 (1980), 54-661.
512 ARITHMETIC 4.7 We obtained the simple state of affairs in (20) by starting with V and u, then defining U. But in practice we generally want to go the other way: Starting with some given function U, we want to find V and u such that (19) holds, i.e., such that V(U(z)) = uV(z). (21) Such a function V is called the Schriider function of U, because it was introduced by Ernst Schriider in Math. Annalen 3 (1871), 296-322. Let us now look at the problem of finding the Schrijder function V(z) = z + V2z2 +. . . of a given power series U = Ulz $ Uzz2 + a . . . Clearly u = VI if (21) is to hold. Expanding (21) with u = U1 and equating coefficients of z leads to a sequence of equations that begins u34 + u, = UIV2 up3 + 2u1uzvz + u, = UlV3 u; vi + 3u;u&s + 2ulusv2 + u;v2 + u, = ulb -, and so on. Clearly there is no solution when VI = 0 (unless trivially U2 = U, = . . . = 0); otherwise there is a unique solution unless VI is a root of unity. We might have expected that something funny would happen when Uy = 1, since Eq. (20) tells us that U (z) = z if the Schrijder function exists. For the moment let us assume that VI is nonzero and not a root of unity; then the Schrijder function does exist, and the next question is how to compute it without doing too much work. The Collowing procedure has been suggested by R. P. Brent and J. F. Traub. Equation (21) leads to subproblems of a similar but more complicated form, so we set ourselves a more general task whose subtasks have the same form: Let us try to find V(z) = V. + Viz + … + Vn-l~n-l such that V(W) = Wz)V(z) + S(z) + O(z L (22) W(z), SC z , and n, where n is a power of 2 and U(0) = 0. If n = 1 we simply let Vi = S(O)/(l -W(O)), with Vi = 1 if S(0) = 0 and W(0) = 1. Furthermore it is possible to go from n to 2n: First we find R(z) such that given u(z), ) V(U(z)) = W(z)V(z) + S(z) -z R(z) + O(z2 ). (23) Then we compute Wz) = W4(zlW)) , 3(z) = R(z)(z/U(z))n (24) and find P(z) = V, + Vntlz + . . . + V2n–1~n-1 such that P(U(z)) = vv(z)P(z) + 3(z) + O(P). (25) It follows that V*(U(z)) = W(z)V*(z) + S(Z) + O(z2 ) as desired, where V*(z) is equal to V(z) + z V(z).
4.7 MANIPULATION OF POWER SERIES 511 The running time for this algorithm to obtain the coefficients up to iV = 2K is T(N), where T(2N) = T(N) + (time to do step N4) + O(N). (17) Straightforward algorithms for composition and division in step N4 will take order N3 steps, so Algorithm N will run slower than Algorithm T. However, Brent and Kung have found a way to do the required composition of power series with O(N log iV)3/2 arithmetic operations, and exercise 6 gives an even faster algorithm for division; hence (17) shows that power series reversion can be achieved by doing only O(Nlog N)3/2 operations as N + co. (On the other hand the constant of proportionality is such that N must be really large before Algorithms L and T lose out to this high-speed method.) Historical note: J. N. Bramhall and M. A. Chapple published the first O(n ) method for power series reversion in CACM 4 (1961), 317-318, 503. It was an off-line algorithm whose running time was approximately the same as that of Algorithms L and T. Iteration of series. If we want to study the behavior of an iterative process X, t f(~~-~), we are interested in studying the n-fold composition of a given function f with itself, namely 2, = f(f( . . . f(zo) . . . )). Let us define f (x) = x and f (z) = f(fnY1(,)), so that f +w = frn(Y(4) (18) for all integers m, 72 2 0. It also makes sense to define jn(,) wh& n is a negative integer, namely to let f and f- be inverse functions such that 2 = fn(fPn(,)); then (18) holds for all integers m and n. Reversion of series is essentially the operation of finding the inverse function f- (z); for example, Eqs. (10) and (11) essentially state that z = V(W(.z)) and that t = W(V(t)), so w = v-l. Suppose we are given two power series V(z) = z + V2z2 + . . . and W(Z) = z + w2.2 + … such that W = V-l. Let u be any nonzero constant, and consider the function U(z) = W(uV(z)). (19) It is easy to see that U(U@)) = W(u2V(z)), and in general that U (z) = W(u V(2)) (20) for all integers 72. Therefore we have a simple expression for the nth iterate Un, which can be calculated with roughly the same amount of work for all n. Furthermore, we can even use (20) to define U for non-integer values of n; the half iterate U1i2, for example, is a function such that U1/2(U1/2(~)) = U(Z). (There are two such functions U112, obtained by using fi and -fi as the value of u l .)
510 ARITHMETIC 4.7 Equation (16) explains the mechanism of this algorithm, which is due to Henry C. Thacher, Jr. [CACM 9 (1966), 10-111. The running time is essentially the same as Algorithm L, but considerably more storage space is required. An example of this algorithm is worked out in exercise 9. Still another approach to power series reversion has been proposed by R. P. Brent and H. T. Kung [JACM 25 (1978), 581-5951, based on the fact that stand- ard iterative procedures used to find roots of equations over the real numbers can also be applied to equations over power series. In particular, we can con- sider Newton s method for computing approximations to a real number t such that f(t) = 0, given a function f that is well-behaved near t: If z is a good approximation to t, then 4(z) = z - f(~)/f (~) will be even better, for if we write z = t + E we have f(x) = f(t) + ef (t) + O(e2), f (s) = f (t) + O(t); consequently 4(5) = t + t -(0 + ef (t) + O(?))/(f (t) + O(E)) = t + O(e2). Applying this idea to power series, let f(z) = V(s) -V(Z), where U and V are the power series in Eq. (14). We wish to find the power series t in z such that f(t) = 0. Let z = W1 z +. . . + W,-l zn- = t + O(F) be an approximation to t of order n; then @(z) = z - f(~)/f (z) will be an approximation of order 2n, since the assumptions of Newton s method hold for this f and t. In other words, we can use the following procedure: Algorithm N (General power series reversion by Newton s method). This semi- on-line algorithm inputs the values of U, and V, in (14) for 2k 5 n < 2k+1 and then outputs the values of W, in (15) for 2k 5 n < 2kf1, thereby producing its answers in batches of 2k at a time, for k = 0, 1, 2, . . . , K. Nl. [Initialize.] Set N c 1. (We will have N = ak.) Input the first coefficients VI and VI (where VI = l), and set WI e VI. N2. [Output.] Output W, for N 5 n < 2N. N3. [Input.] Set N e 2N. If N > 2K, the algorithm terminates; otherwise input the values U, and V, for N 5 n < 2N. N4. [Newtonian step.] Use an algorithm for power series composition (see exer- cise 11) to evaluate the coefficients Qj and Rj (0 < j < N) in the power series u,z+ + U2jl,-$2N-1 -v(wlz + -+wN-lzN- ) = RozN + RlzN+l + . . . + RN-1z2N- + O(Z~~), v (w,z + + wj,T-lz N-1) = Qo + QIZ + . . . + QN-~z~- + O(zN), where V(x) = 2 + V2s2 + . .. and V (x) = 1 + 2Vzz + . . . Then set WN, . . . . W2N–1 to the coefficients in the power series Ro+Rlz+. * .+RN-~z N-l = WN + + W&X-12 N-1 + O(ZN) Qo+Qlz+. . .+QN-Iz -~ and return to step N2. m