J2ee Web Hosting, Perl Scripting And Developing Blog

3.4.1 NUMERICAL DISTRIBUTIONS 115 and after these two instructions have been executed the desired integer will appear in register A. If a random integer between 1 and Ic is desired, we add one to this result. (The instruction INCA I would follow (l).) This method gives each integer with nearly equal probability. There is a slight error because the computer word size is finite (see exercise 2); but the error is quite negligible if Ic is small, for example if k/m < l/10000. In a more general situation we might want to give different weights to different integers. Suppose that the value X = x1 is to be obtained with probability pl, and X = xz with probability pa, . . . , X = xk with probability pk. We can generate a uniform number U and let Xl, if 0 5 U < pi; x2, ifn I U < p1 +p2; x= . (2) 1: xk, ifpl +p2 + +pk-1 5 u < 1. (Note that pl + p2 + . . . + pk = 1.) There is a best possible way to do the comparisons of U against various values of PI + ~2 + . . . + P,, as implied in (2); this situation is discussed in Section 2.3.4.5. Special cases can be handled by more efficient methods; for example, to obtain one of the eleven values 2, 3, . . . , 12 with the respective dice probabilities &, 6, . . . , &, . . . , &, 9, we could compute two independent random integers between 1 and 6 and add them together. However, none of the above techniques is really the fastest general method for selecting x1, . . . , xk with the correct probabilities. An ingenious way to do the trick has been discovered by A. J. Walker [Electronics Letters 10,8 (1974), 127-128; ACM Trans. Math. Software 3 (1976) 253-2561. Suppose we form kU and consider the integer part K = LkU] and fraction part V = (kU) mod 1 separately; for example, after the code (1) we will have K in register A and V in register X. Then we can always obtain the desired distribution by doing the operations if V < PK then X +- XK+~ otherwise X t Y,, (3) for some appropriate tables (PO, . . . , P&i) and (Yo, . . . , Yk-1). Exercise 7 shows how such tables can be computed in general. Walker s method is sometimes called the method of aliases. On a binary computer it is usually helpful to assume that k is a power of 2, so that multiplication can be replaced by shifting; this can be done without loss of generality by introducing additional x s that occur with probability zero. For example, let s consider dice again; suppose we want X = j to occur with the following 16 probabilities: j = 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Pj=oo~~&&?i&~~~~~o 0 0

114 RANDOMNUMBERS 3.4 3.4. OTHER TYPES OF RANDOM QUANTITIES WE HAVE NOW SEEN how to make a computer generate a sequence of numbers uo, Ul, u2, . . . that behaves as if each number were independently selected at random between zero and one with the uniform distribution. Applications of random numbers often call for other kinds of distributions, however; for example, if we want to make a random choice from among Ic alternatives, we want a random integer between 1 and k. If some simulation process calls for a random waiting time between occurrences of independent events, a random number with the exponential distribution is desired. Sometimes we don t even want random numbers-we want a random permutation (i.e., a random arrangement of n objects) or a random combination (i.e., a random choice of k objects from a collection of n). In principle, any of these other random quantities may be obtained from the uniform deviates UO, VI, U2, . . . . People have devised a number of important random tricks that may be used to perform these manipulations efficiently on a computer, and a study of these techniques also gives some insight into the proper use of random numbers in any Monte Carlo application. It is conceivable that someday somebody will invent a random number generator that produces one of these other random quantities directly, instead of getting it indirectly via the uniform distribution. But except for the random bit generator described in Section 3.2.2, no direct methods have so far proved to be practical. The discussion in the following section assumes the existence of a random sequence of uniformly distributed real numbers between zero and one. A new uniform deviate U is generated whenever we need it. These numbers are usually represented in a computer word with the decimal point assumed at the left. 3.4.1. Numerical Distributions This section summarizes the best techniques known for producing numbers from various important distributions. Many of the methods were originally suggested by John von Neumann in the early 195Os, and they have gradually been improved upon by other people, notably George Marsaglia, J. H. Ahrens, and U. Dieter. A. Random choices from a finite set. The simplest and most common type of distribution required in practice is a random integer. An integer between 0 and 7 can be extracted from three bits of U on a binary computer; in such a case, these bits should be extracted from the most significant (left-hand) part of the computer word, since the least significant bits produced by many random number generators are not sufficiently random. (See the discussion in Section 3.2.1.1.) In general, to get a random integer X between 0 and k -1, we can multiply by k, and let X = LkUj. On MIX, we would write LDA U MULK (1)

3.3.4 THE SPECTRAL TEST 113 b 24. [M.28] Generalize the spectral test to second-order sequences of the form X, = (ax,-1 + bX,-z)modp, having period length p2 -1. (Cf. Eq. 3.2.2-8.) How should Algorithm S be modified? 25. [HA&24] Let d be a divisor of m and let 0 5 9 < d. Prove that c r(k), summed over all 0 5 k < m such that Icmodd = 9, is at most (2/dr)ln(m/d) + O(1). (Here r(k) is defined in Eq. (46) when t = 1.) 26. [A&?$ Explain why the derivation of (53) leads to a similar bound on N-l for 0 < 9 < m. Where does the derivation of (53) break down when m = l? 27. [HA4391 (E. Hlawka, H. Niederreiter.) Let r(ui, . , urn) be the function defined in (46). Prove that c T(u~, . . . , it), summed over all 0 5 ~1,. , it < m such that (~1,. . . ,~t) # (0,. . ,O) and (15) holds, is 0((27rlgm)t ~~~~~~ where rmax is the maximum term ~(2~1,. . . , it) in the sum. b 28. [M.28] (H. Niederreiter.) Find an analog of Theorem N for the case m = prime, c = 0, a = primitive root modulo m, X0 $ 0 (modulo m). [Hint: Your exponential sums should involve 5 = e2Ti (m-1) as well as w.] Prove that in this case the average primitive root has discrepancy 0:)-r = O(t(log m) /#(m -l)), hence good primitive roots exist for all m. 29. [A& 11 Prove that quantity rmax of exercise 27 is never larger than &?/ut, unless UT 5 2(t -1). 30. [A4331 (S. K. Zaremba.) Prove that in two dimensions, rmax 5 m/max(al,. . . , a,), whereal,…, a, are the partial quotients obtained when Euclid s algorithm is applied to m and a. [Hint: We have a/m = /al,. , a,/ in the notation of Section 4.5.3; apply ex- ercise 4.5.3-42.1 31. [HM47] (I. Borosh.) Prove that for all sufficiently large m there exists a number a relatively prime to m such that all partial quotients of u/m are 2 3. Furthermore the set of all m satisfying this condition but with all partial quotients 5 2 has positive density.

112 RANDOM NUMBERS 3.3.4 13. [HM,%!] Lemma A uses the fact that U is nonsingular to prove that a positive definite quadratic form attains a definite, nonzero minimum value at nonzero integer points. Show that this hypothesis is necessary, by exhibiting a quadratic form (19) whose matrix of coefficients is singular, and for which the values of f(z~, . . . , Q) get arbitrarily near zero (but never reach it) at nonzero integer points (~1,. . . , Q). 14. [24] Perform Algorithm S by hand, for m = 100, a = 41, T = 3. b 15. [MZO] Let U be an integer vector satisfying (15). How many of the (t -l)- dimensional hyperplanes defined by U intersect the unit hypercube { (~1,. . . ,z,) ] 0 2 xj < 1 for 1 2 j 5 t }? (This is approximately the number of hyperplanes in the family that will suffice to cover LO.) 16. [MAO] (U. Dieter.) Show how to modify Algorithm S in order to calculate the minimum number A$ of parallel hyperplanes intersecting the unit hypercube as in exercise 15, over all U satisfying (15). [Hint: What are appropriate analogs to positive definite quadratic forms and to Lemma A?] 17. [ZO] Modify Algorithm S so that, in addition to computing the quantities vt, it outputs all integer vectors (~1,. . . , it) satisfying (15) such that u: + + . e + ZL~ = Y:, for 2 5 t 5 T. 18. [M30] (a) Let m = 2e, where e is even. By considering combinatorial matrices, i.e., matrices whose elements have the form y + x&j (cf. exercise 1.2.3-39), find 3 X 3 matrices of integers U and V satisfying (29) such that the transformation of step S5 does nothing for any j, but the corresponding values of zk in (31) are so huge that exhaustive search is out of the question. (The matrix U need not satisfy (28), we are interested here in arbitrary positive definite quadratic forms of determinant m.) (b) Although transformation (23) is of no use for the matrices constructed in (a), find another transformation that does produce a substantial reduction. F 19. [HMZ5] Suppose step S5 were changed slightly, so that a transformation with q = 1 would be performed when 2V,. q = I$ e V,. (Thus, q = [(vi. I,$ / I$. V,) + $1 in all cases.) Would it be possible for Algorithm S to get into an infinite loop? 20. [M21] Discuss how to carry out an appropriate spectral test for linear congruential sequences having c = 0, X0 odd, m = 2e, a mod 8 = 5. 21. [M80] (R. W. Gosper.) A certain application uses random numbers in batches of four, but throws away the second of each set. How can we study the grid structure of {&(X4% X4n+2, Xdn+s)}, given a linear congruential generator of period m = 2e? 22. [M46] What is the best upper bound on ps, given that ~2 is very near its maxi- mum value &?T? What is the best upper bound on ~2, given that ~3 is very near its maximum value +rfi? 23. [M46] Let U,, Vj be vectors of real numbers with Vi . Vj = &j for 1 5 i,j 5 t, and such that Ui . Ui = 1, 2]Uz . Ujl 5 1, 21Vz +l.$ < Vj . Vj for i # j. HOW large can !JI . VI be? (This question relates to the bounds in step S8, if both (23) and the transformation of exercise 18(b) fail to make any reductions. The maximum value known to be achievable is (n + 2)/3, which occurs when VI = II, Uj = $11 + &&Ij, Vl = I1 -(12 + * * . + In)/& Vj = 2Ij/fi, for 2 2 j 5 n, where (II,. . . , In) is the identity matrix; this construction is due to B. V. Alekseev [Matematicheskie Zametki, to appear] .)

3.3.4 THE SPECTRAL TEST 111 what happens when t = l?) 2. [Hi%?0] Let VI, , Vt be linearly independent vectors in t-space, let LO be the lattice of points defined by (lo), and let Vi, . . . , Ut be defined by (19). Prove that the maximum distance between (t -1)-dimensional hyperplanes, over all families of parallel hyperplanes that cover LO, is l/min{ f(zi, . . , Q) 1 (xi,. . , Q) # (0,. . . ,O)}, where f is defined in (17). 3. [M24] Determine ~3 and ~4 for all linear congruential generators of potency 2 and period length m. b 4. [ME?] Let ~11, ~12, usi, ~2s be elements of a 2 X 2 integer matrix such that ~11 $ a~12 = ~21 + ~2~22 E 0 (modulo m) and ~112122 -~21~12 = m. (a) Prove that all integer solutions (yi, 92) to the congruence y1 + ay2 E 0 (modulo m) have the form (~1,~s) = (51~11+222~21,51 (~12+222~22) for integer xl, 22. (b) If, in addition, 2~~11~214-~121122) 5 u~,+u& 5 u;,+&, prove that (yi, ys) = (~11, ~12) minimizes y: + yi over all nonzero solutions to the congruence. 5. [A4301 Prove that steps Sl through S3 of Algorithm S correctly perform the spectral test in two dimensions. [Hint: See exercise 4, and prove that (h + h)2 + (p + p) 2 h2 + p2 at the beginning of step S2.1 6. [A&N] Let ac, al, . . , at-1 be the partial quotients of a/m as defined in Section 3.3.3, and let A = rnaxc13 27r/(A + 1 + l/A). 7. [HiME??] Prove that question (a) and question (b) of the text have the same solution for real values of 41, . . . , q3-r, q3fl, . . , qt (cf. (24), (26)). 8. [AR61 Line 18 of Table 1 has a very low value of ~2, yet ~3 is quite satisfactory. What is the highest possible value of ~3 when ~2 = lo@ and m = lOlo? 9. [Hi!&%] (C. Hermite, 1846.) Let f(zi,. , xt) be a positive definite quadratic form, defined by the matrix U as in (17), and let 0 be the minimum value of f at nonzero integer points. Prove that ~9 5 ( $)(t-1)/2 1 det U12/t. [Hints: If W is any integer matrix of determinant 1, the matrix WU defines a form equivalent to f; and if S is any orthogonal matrix (i.e., S- = ST ), t he matrix US defines a form identically equal to f. Show that there is an equivalent form g whose minimum 9 occurs at (l,O, . ,O). Then prove the general result by induction on t, writing g(xi, . , xt) = B(XlfP2~2 +. . . + Ptxt) + h(x2,. . . , xt) where h is a positive definite quadratic form in t -1 variables.] 10. [M.28] Let (yi, ~2) be relatively prime integers such that yl +ay2 E 0 (modulo m) and Y: -I-Y; < m m. Show that there exist integers (Us, ~2) such that u1 +au2 F 0 (module m), 211~2-212~1 = m, 21~1~1 +UZYZ~ 5 min(zl:-t&, y:-l-&, and (UT + U$X (Y: + ~2) 2 m2. (Hence by exercise 4, V; = min(uf +u;, y:+yi).) b 11. [HM30] (Alan G. Waterman, 1974.) Invent a reasonably efficient procedure that computes multipliers a E 1 (modulo 4) for which there exists a relatively prime solution to the congruence yl + ay2 E 0 (modulo m) with y: + yz = mm-e, where E > 0 is as small as possible, given m = 2 . (By exercise 10, this choice of a will guarantee that ~2 2 m /(y: + yi) > m m, and there is a chance that ui will be near its optimum value mm. In practice we will compute several such multipliers having small c, choosing the one with best spectral values ~2, ~3, . . . .) 12. [HM.ZS] Prove, without geometrical handwaving, that any solution to the text s question (b) must also satisfy the set of equations (26).

110 RANDOM NUMBERS 3.3.4 in this case is extremely small in spite of the fact that there are parallelogram- shaped regions of area x l/fi containing no points (U,, Un+l). The fact that discrepancy can change so drastically when the points are rotated warns us that the serial test may not be as meaningful a measure of randomness as the rotation-invariant spectral test. G. Historical remarks. In 1959, while deriving upper bounds for the error in the evaluation of t-dimensional integrals by the Monte Carlo method, N. M. Korobov devised a way to rate the multiplier of a linear congruential sequence. His formula (which is rather complicated) is related to the spectral test since it is strongly influenced by small solutions to (15); but it is not quite the same. Korobov s test has been the subject of an extensive literature, surveyed by Kuipers and Niederreiter in Uniform Distribution of Sequences (New York: Wiley, 1974), 52.5. The spectral test was originally formulated by R. R. Coveyou and R. D. MacPherson [JACM 14 (196i ), lOO-1191, who introduced it in an interesting indirect way. Instead of working with the grid structure of successive points, they considered random number generators as sources of t-dimensional waves. The numbers ~xY +. . . + X; such that 21 + . . . + at- zt E 0 (modulo m) in their original treatment were the wave frequencies, or points in the spectrum defined by the random number generator, with low-frequency waves being the most damaging to randomness; hence the name spectral test. Coveyou and MacPherson introduced a procedure analogous to Algorithm S for performing their test, based on the principle of Lemma A. However, the original procedure (which used matrices UUT and VVT instead of U and V) dealt with extremely large numbers; the idea of working directly with U and V was independently suggested by F. Janssens and by U. Dieter. Several other authors pointed out that the spectral test could be understood in far more concrete terms; by introducing the study of the grid and lattice struc- tures corresponding to linear congruential sequences, the fundamental limitations on randomness became graphically clear. See G. Marsaglia, Proc. Nat. Acad. Sci. 61 (1968), 25-28; W. W. Wood, J. Chem. Phys. 48 (1968), 427; R. R. Coveyou, Studies in Applied Math. 3 (1969), 70-112; W. A. Beyer, R. B. Roof, and D. Williamson, Math. Comp. 25 (1971), 345-360; G. Marsaglia and W. A. Beyer, Applications of Number Theory to Numerical Analysis, ed. by S. K. Zaremba (New York: Academic Press, 1972), 249-285, 361-370. Harald Niederreiter s papers concerning the use of exponential sums to study the distribution of linear congruential sequences have appeared in Math. Comp. 26 (1972), 793-795; 28 (1974), 1117-1132; 30 (1976), 571-597; Advances in Math. 26 (1977), 99-181 [this is the most important paper of the series]; and Bull. Amer. Math. Sot. 84 (1978), 273-274: 957-1041 [this one summarizes the others and contains an extensive bibliography]. EXERCISES 1. [M10] To what does the spectral test reduce in one dimension? (In other words,

3.3.4 THE SPECTRAL TEST 109 In fact, the upper bound gets even smaller when q has a factor in common with m, since s and m/G generally become smaller. (See exercise 26.) We have now proved that the g(ul, . . . , ut) part of our upper bound (44) on the discrepancy is small, if N is large enough and if (ul,. . . , Q) does not satisfy the spectral test congruence (15). Exercise 27 proves that the f(ul, . . . , ut) part of our upper bound is small, when summed over all the nonzero vectors CR,…, ut) satisfying (15), provided that all such vectors are far away from (0,. . . ,O). Putting these results together leads to the following theorem of Niederreiter: Theorem N. Let (Xn) be a linear congruential sequence (X0, a, c, m) of period length m, and let s be the least positive integer such that us = 1 (modulo m). Let ut be the t-dimensional accuracy of (XJ, as determined by the spectral test. Then the t-dimensional discrepancy DN(t) determined by the first N values of (Xn), as defined in (42)) satisfies D(t) = o d%k4hwV N N ) + O( mFJ)t) + O((logm)t rmax); (54) ( Dg = O((logm)t rmrtx), (55) Here rmax is the maximum value of the quantity T(u~, . . . , Ut) defined in (46), taken over all nonaero integer vectors (ul, . . . , Q) satisfying (15). Proof. The first two 0 terms in (54) come from vectors (ul, . . . , Ut) in (44) that do not satisfy (15), since exercise 25 proves that f(zll, . . . , Ut) summed over all (Ul,… , ut) is O(((~/YT) lnm)t) and exercise 26 bounds each g(ul, . . . , Ut). (These terms are missing from (55) since g(ul, . . . , ut) = 0 in that case.) The remaining 0 term in (54) and (55) comes from nonzero vectors (~1,. . . , ut) that do satisfy (15), using the bound derived in exercise 27. (By examining this proof carefully, we could replace each 0 in these formulas by an explicit function of t.) I Eq. (55) relates to the serial test in t dimensions over the entire period, while Eq. (54) gives us useful information about the distribution of the first N generated values when N is less than m, provided that N is not too small. Note that (54) will guarantee low discrepancy only when s is sufficiently large, otherwise the m/G term will dominate. If m = p: . . . pF7 and gcd(a - 1, m) = fl P, ***P, fr , then s equals ~Tl-~l. . . pFr–fr (cf. Lemma 3.2.1.2P); thus, the largest values of s correspond to high potency. In the common case m = 2e and a G 5 (modulo , we have s = am, so0:) is O(fi(logm)t+l/N)+O((logm)tr,,,). It is not difficult to prove that T,,, 2 a/z+ unless vt is very small (see exercise 29); therefore Eq. (54) says in particular that the discrepancy will be low in t dimensions if the spectral test is passed and if N is somewhat larger than Jm (log m)t+l. In a sense Theorem N is almost too strong, for the result in exercise 30 shows that linear congruential sequences like those in lines 8, 19, and 23 of Table 1 have a discrepancy of order (logm)2/m in two dimensions. The discrepancy

108 RANDOM NUMBERS 3.3.4 where Now Ski = W-lkSko, so ISkl( = ISkO( for all I, and we can calculate this common value by further exponential-summery: tskOt2 = ; c lskl12 O O

3.3.4 THE SPECTRAL TEST 107 Both f and g can be further simplified in order to get a good upper bound on Ofi . We have when u # 0, and the sum is 5 1 when u = 0; hence f(u1,. . .7 ut) I T(%,. . . , Ut), (45) where 1 T(Ul,… 7%) = n (46) l

106 RANDOM NUMBERS 3.3.4 An upper bound for the discrepancy can be found by using exponential sums. Let w = earijrn be a primitive mth root of unity. If (Q, . . . , zt) and (yl, . . . , yt) are two vectors with all components in the range 0 1. xj, yj < m, we have W(51-Y1)211+…+(Zt-Yt)Ut = mt, if (XI,. . . , xt) =(~1,. . . , yt); c 0, if(xl,…,xt) # (YI,… ,yt). O