Web And Email Hosting, Jsp Web Hosting Blog

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

3.3.4 THE SPECTRAL TEST 105 in dimensions 2 through 5, partly because it is easy to remember. Line 27 uses a random primitive root, modulo the prime 231 - 1, as multiplier. Line 28 is for the sequence X, = (271828183X,-1 -314159269X,-z) mod (231 -I), (39) which can be shown to have period length (231 -1)2 -1; it has been analyzed with the generalized spectral test of exercise 24. Theoretical upper bounds on pt, which can never be transcended for any m, are shown just below Table 1; it is known that every lattice with m points per unit volume has vt i 7t milt, (40) where yt takes the respective values (4/3) , 2l 7 2114, 23 1o, (64/3)1 12, 23/7, 2112 (41) for t = 2, . . . , 8. (See exercise 9 and J. W. S. Cassels, Inntroduction to the Geometry of Numbers (Berlin: Springer, 1959), p. 332.) These bounds hold for lattices generated by vectors with arbitrary real coordinates. For example, the optimum lattice for t = 2 is hexagonal, and it is generated by vectors of length 2/G that form two sides of an equilateral triangle. In three dimensions the optimum lattice is generated by vectors VI, V,, V3 that can be rotated into the form (w, ZI, -v), (w, -v,v), (-w, w, w), where w = l/a. *F. Relation to the serial test. In a series of important papers published during the 19 7Os, Harald Niederreiter has shown how to analyze the distribution of the t-dimensional vectors (1) by means of exponential sums. One of the main consequences of his theory is that the serial test in several dimensions will be passed by any generator that passes the spectral test, even when we consider only a sufficiently large part of the period instead of the whole period. We shall now turn briefly to a study of his interesting methods, in the case of linear congruential sequences (X0, a, c, m) of period length m. The first idea we need is the notion of discrepancy in t dimensions, a quantity that we shall define as the difference between the expected number and the actual number of t-dimensional vectors (z,, z,+~, . . . , s,ttYl) falling into a hyper- rectangular region, maximized over all such regions. To be precise, let (z,) be a sequence of integers in the range 0 5 x, < m. We define number of (x,, . . . , x,+t-l) in R for 0 5 n < N volume of R Djj) = max R N mt (42) where R ranges over all sets of points of the form R = { (Yl, * * * , Yt) I QIl I Yl < Pl, . . . , Qt I Yt < Pt >; (43) here cyj and & are integers in the range 0 5 CZ~ < pj 5 m, for 1 5 j 5 t. The volume of R is clearly (pl -ol). . . (& -at). To get the discrepancy DE , we imagine looking at all these sets R and finding the one with the greatest excess or deficiency of points (x,, . . . , Z,+t–l).

104 RANDOM NUMBERS 3.3.4 small that the numbers can hardly be called random; the vt values are terribly low. Line 6 is the generator discussed above; line 7 is a similar example, having an abnormally low value of p3. Line 8 shows a nonrandom multiplier for the same modulus m; all of its partial quotients are 1, 2, or 3. Such multipliers have been suggested by I. Borosh and H. Niederreiter because the Dedekind sums are likely to be especially small and because they produce best results in the two- dimensional serial test (cf. Section 3.3.3 and exercise 30). The particular example in line 8 has only one 3 as a partial quotient; there is no multiplier congruent to 1 modulo 20 whose partial quotients with respect to lOlo are only l s and 2 s. The generator in line 9 shows another multiplier chosen with malice aforethought, following a suggestion by A. G. Waterman that guarantees a reasonably high value of ~2 (see exercise 11). Lines 10 through 21 of Table 1 show further examples with m = 235, beginning with some random multipliers. The generators in lines 12 and 13 are reminders of the good old days-they were once used extensively since 0. Taussky first suggested them in the early 1950s. Lines 14 through 18 show various multipliers of maximum potency having only four l s in their binary representation. The point of having few l s is to replace multiplication by a few additions, but only line 16 comes near to being passable. Since these multipliers satisfy (a - 5)3 mod 2 35 = 0, all five of them achieve ~4 at the same point (zl,zz,zs,z4) = (-125, 75, -15,l). Another curiosity is the high value of ~3 following a very low ~2 in line 18 (see exercise . Lines 19 and 20 are respectively the Borosh-Niederreiter and Waterman multipliers for modulus 235; and line 21 was found by M. Lavaux and F. Janssens, in a computer search for spectrally good multipliers having a very high ~2. Lines 22 through 28 apply to System/370 and other machines with 32-bit words; in this case the comparatively small word size calls for comparatively greater care. Line 22 is, regrettably, the generator that has actually been used on such machines in most of the world s scientific computing centers for about a decade; its very name RANDU is enough to bring dismay into the eyes and stomachs of many computer scientists! The actual generator is defined by X0 odd, Xn+i = (65539X,) mod 231, (38) and exercise 20 indicates that 22g is the appropriate modulus for the spectral test. Since 9X, + 6Xn+z +Xn+z E 0 (modulo 231), the generator fails most three- dimensional criteria for randomness, and it should never have been used. Almost any multiplier -5 (modulo would be better. (A curious fact about RAND& noticed by R. W. Gosper, is that v4 = v5 = vs = ur = vs = vg = &i6, hence ps is a spectacular 11.98.) Lines 23 and 24 are the Borosh-Niederreiter and Waterman multipliers for modulus 232, lines 26 and 29 were found by Lavaux and Janssens, and line 30 (whose excellent multiplier 6364136223846793005 is too big to fit in the column) is due to C. E. Haynes. Line 25 was nominated by George Marsaglia as a candidate for the best of all multipliers, after a computer search

3.3.4 THE SPECTRAL TEST 103 k z k 3 kv4 ku5 lg& /& IL3 114 /15 /dLg Linl 4.5 4.5 4.5 4.5 4.4 2E5 564 0.01 0.34 4.62 1 7.0 7.0 7.0 7.0 4.0 2t6 3E4 0.04 4.66 2c3 2 17.5 1.3 1.0 1.0 1.0 3.14 2cg 2~~ 52 t* 3 15.7 10.0 7.4 5.0 5.0 0.27 0.13 0.11 0.01 0.21 4 4.0 2.5 1.9 1.3 1.0 3.36 2.69 3.78 1.81 1.29 5 16.0 10.0 8.0 5.4 4.5 1.44 0.44 1.92 0.07 0.08 6 16.0 9.0 8.3 5.9 5.6 1.35 0.06 4.69 0.35 6.98 7 16.7 10.7 7.8 6.3 4.8 3.39 1.75 1.20 1.39 0.28 8 16.5 11.1 7.0 6.4 5.2 2.89 4.15 0.14 2.04 1.25 9 16.8 11.2 8.2 6.8 5.8 1.24 1.70 1.12 2.79 3.81 10 17.2 10.7 8.6 6.7 5.8 2.10 0.55 3.15 1.85 3.72 11 17.5 10.7 8.4 6.8 5.6 3.03 0.61 1.85 2.99 1.73 12 17.2 11.6 8.7 6.3 5.7 2.02 4.02 4.03 0.40 2.62 13 13.7 11.5 7.2 7.0 5.3 0.02 2.79 0.07 5.61 0.65 14 13.7 11.2 7.2 5.5 5.1 0.02 1.50 0.07 0.03 0.22 15 16.8 11.2 7.2 6.9 5.5 1.12 1.67 0.07 3.13 1.26 16 16.5 10.4 7.2 6.5 5.8 0.75 0.30 0.07 0.82 4.39 17 11.5 11.5 7.2 7.0 5.3 8c4 2.95 0.07 5.53 0.50 18 17.5 9.7 7.2 6.6 5.7 2.95 0.07 0.07 1.37 2.83 19 17.5 11.4 8.3 5.8 5.7 3.33 2.44 1.30 0.08 3.52 20 17.6 11.6 8.8 7.1 5.6 3.60 3.56 5.91 7.38 2.03 21 14.5 3.4 3.4 3.4 3.4 3.14 3 8 63 0.02 22 16.0 10.2 6.9 6.1 4.9 3.16 1.73 0.26 2.02 0.89 23 16.1 10.5 7.7 6.1 4.7 3.41 2.92 2.32 1.81 0.35 24 16.0 10.5 7.8 6.4 4.0 3.10 2.91 3.20 5.01 0.02 25 16.1 10.6 8.0 6.0 5.0 3.61 3.45 4.66 1.31 1.35 26 15.2 9.9 7.6 5.9 5.1 2.10 1.66 3.14 1.69 3.60 27 31.0 20.2 15.6 11.8 9.8 3.14 1.49 0.44 0.69 0.66 28 24.1 16.0 12.0 9.1 7.9 3.60 3.92 5.27 0.97 3.82 29 31.5 21.3 16.0 12.7 10.4 1.50 3.68 4.52 4.02 1.76 30 upper bounds from (40): 3.63 5.92 9.87 14.89 23.87 Table 1 shows what sorts of values occur in typical sequences. Each line of the table considers a particular generator, and lists vt, pLt, and the number of bits of accuracy lg vt. Lines 1 through 4 show the generators that were the subject of Figs. 2 and 5 in Section 3.3.1. The generators in lines 1 and 2 suffer from too small a multiplier; a diagram like Fig. 8 will have a nearly vertical stripes when a is small. The terrible generator in line 3 has a good p2 but very poor ~3 and ~4; like nearly all generators of potency 2, it has v3 = & and ~4 = 2 (see exercise 3). Line 4 shows a random multiplier; this generator has satisfactorily passed numerous empirical tests for randomness, but it does not have especially high values of ~2, . . . , ps. In fact, the value of ~5 flunks our criterion. Line 5 shows the generator of Fig. 8. It passes the spectral test with very high-hying colors, when ~2 through PcLg are considered, but of course m is so