Web And Email Hosting, Jsp Web Hosting Blog

3.3.4 THE SPECTRAL TEST 95 vl,…, Vt such that For example, in the special form (16) that arises in the spectral test, we have Ul = ( m, 070,. . . ,017 Vl =+$,a$2 ,…) at- ), u2 = ( –a, 1, 0, . . . , O), v2 = (0, 1, 0,. . . , O), us = ( -u2, 0, 1, . . . ) O), v, = (O,O, 1,. . . , O), (20) . . . . . . . . . . . . ut = (-d-l, 0, 0, . . . , l), vt = (O,O, 0,. . . , 1). These V, are precisely the vectors (8), (9) that we used to define our original lattice Lo. As the reader may well suspect, this is not a coincidence-indeed, if we had begun with an arbitrary lattice LO, defined by any set of linearly independent vectors VI, . . ,V,, the argument we have used above can be generalized to show that the maximum separation between hyperplanes in a covering family is equivalent to minimizing (17), where the coefficients uij are defined by (19). (See exercise 2.) Our first step in minimizing (18) is to reduce it to a finite problem, i.e., to show that we won t need to test infinitely many vectors (~1, . . . , Q) to find the minimum. This is where the vectors VI, . . . , Vt come in handy; we have xk = (xl& + +xtut) vk, and Cauchy s inequality tells us that ((xlul + + xtut) vk)2 2 j-(x1,. . . > xt)(Vk vk). Hence we have derived a useful upper bound on each coordinate xk: Lemma A. Let (xl,…, q) be a nonaero vector that minimizes (18) and let (Yl,.. . , yt) be any nonzero integer vector. Then 2: 2 (vk vk/k)f(Yl , . . , Yt ), for 1 < k 5 t. (21) In particular, letting yz = Sij for all i, 2; 5 (vk vk)(uy u,), for 15 j,k 5 t. 1 (22) Lemma A reduces the problem to a finite search, but the right-hand side of (21) is usually much too large to make an exhaustive search feasible; we need at least one more idea. On such occasions, an old maxim provides sound advice: If you can t solve a problem as it is stated, change it into a simpler problem that

94 RANDOM NUMBERS 3.3.4 hence the minimum in (12) occurs when each xj = uj/(uf + .. . + uz); the distance between neighboring hyperplanes is 1 dm = l/length(U). (14I In other words, the quantity ut we seek is precisely the length of the shortest ve&or U that defines a family of hyperplanes { X. U = q 1 integer q } containing all the elements of Lo. Such a vector U = (~1,. . . , ut) must be nonzero, and it must satisfy V.U = integer for all V in Lo. In particular, since the points (1, 0, . . . , 0), (0, 1, . . . , 0), . . . . (0,O )…) 1) are all in Lo, all of the uj must be integers. Furthermore since VI is in Lo, we must have &(ul + au2 + . . . + atP1ut) = integer, i.e., Ul + au2 + . . . + &lut s 0 (modulo m). (15) Conversely, any nonzero integer vector U = (ul, . . . , ut) satisfying (15) defines a family of hyperplanes with the required properties, since all of LO will be covered: (y1Vl + . . + + yt&) . U will be an integer for all integers yl, . . . , yt. We have proved that u:= min {uf+-..+uF Iu1+uu2+..-+cP1ut =O(modulom)} (~1,..-,~t)#(O,…,O) =(%I,…,$; (O,…,0)(( mx~-ux~-u2×3-*. . -ut-1xt)2+x;+x;+~~~ +xt ). (16) C. Deriving a computational method. We have now reduced the spectral test to the problem of finding the minimum value (16); but how on earth can we determine that minimum value in a reasonable amount of time? A brute-force search is out of the question, since m is very large in cases of practical interest. It will be interesting and probably more useful if we develop a computational method for solving an even more general problem: Find the minimum value of the quantity f(Xl,… , xt) = (%1X1 + . . . + Utlxt)2 + . . . + (%X1 + * *. + %xt)2 (17) over all nonzero integer vectors (xl,. . . , xt), given any nonsingular matrix of oefficients U = (uij). The expression (17) is called a positive definite quadratic . form m t variables. Since U is nonsingular, (17) cannot be zero unless the xj are all zero. Let us write VI, .,., Ut for the rows of U. Then (17) may be written f(Xl,.-. , xt) = (Xl& + . . . + XtUt) . (Xl& + … $-x&s), (18) the square of the length of the vector x1 VI + . . . + xtUt . The nonsingular matrix U has an inverse, which means that we can find uniquely determined vectors

3.3.4 THE SPECTRAL TEST 93 where l4) = k(o, c, (1 + a)c, . . . , (1 + a + . + cP )c) (6) is a constant vector. The variable k1 is redundant in this representation of L, be- cause we can change (5, ki, Its, . . . , k,) to (z+kim, 0, ks –&I,. . . , kt –at- kl), reducing kl to zero without loss of generality. Therefore we obtain the compara- tively simple formula L={Vo+ylK+y2V2+~~~+ytVt I integeryl,ya,…,yt}, (7) where vl = ~(l,u,u2 ,.. , CL- ); (8) Vs=(O,l,O )…, O), vs=(o,o,1,. .) O), . ..) Vt=(O,O,O ) .) 1). (9) The points (zi, x2,. . . , xt) of L that satisfy 0 5 x3 < 1 for all j are precisely the m points of our original set (2). Note that the increment c appears only in Vo, and the effect of V. is merely to shift all elements of L without changing their relative distances; hence c does not affect the spectral test in any way, and we might as well assume that Vi = 640,. . . , 0) when we are calculating LQ. When Vo is the zero vector we have a so-called lattice of points Lo={~IK+Y~V~+.~.+Y~K I integeryl,yz,...,yt), (10) and our goal is to study the distances between adjacent (t -1)-dimensional hyperplanes, in families of parallel hyperplanes that cover all the points of Lo. A family of parallel (t -1)-dimensional hyperplanes can be defined by a nonzero vector U = (ui, . . . , Ut) that is perpendicular to all of them; and the set of points on a particular hyperplane is then {(Xl,. . . 7 xt) I 21% + . . * + xtw = 4 >, (11) where 4 is a different constant for each hyperplane in the family. In other words, each hyperplane is the set of all X for which the dot product X . U has a given value 4. In our case the hyperplanes are all separated by a fixed distance, and one of them contains (O,O, . . . ,O); hence we can adjust the magnitude of U so that the set of all integer values c~ gives all the hyperplanes in the family. Then the distance between neighboring hyperplanes is the minimum distance from (O,O,. . . , 0) to the hyperplane for o = 1, namely Cauchy s inequality (cf. exercise 1.2.3-30) tells us that (13)

92 RANDOM NUMBERS 3.3.4 in each subcube of the unit cube, when the unit cube has been divided into 64 subcubes of size $ X 4 X a; this same generator might yield completely empty subsquares of the unit square, when the unit square has been divided into 64 subsquares of size & x isl. Since we increase our expectations in lower dimensions, a separate test for each dimension is required. It is not always true that vt 5 milt, although this upper bound is valid when the points form a rectangular grid. For example, it turns out that uz = m > m in Fig. 8, because a nearly hexagonal structure brings the m points closer together than would be possible in a strictly rectangular arrrangement. In order to develop an algorithm that computes vt efficiently, we must look more deeply at the associated mathematical theory. Therefore a reader who is not mathematically inclined is advised to skip to part D of this section, where the spectral test is presented as a plug-in method accompanied by several examples. On the other hand, we shall see that the mathematics behind the spectral test requires only some elementary manipulations of vectors. Some authors have suggested using the minimum number Nt of parallel covering lines or hyperplanes as the criterion, instead of the maximum distance l/z,+ between them. However, this number does not appear to be as important as the concept of accuracy defined above, because it is biased by how nearly the slope of the lines or hyperplanes matches the coordinate axes of the cube. For example, the 20 nearly vertical lines that cover all the points of Fig. 8 are actually l/J328 units apart, and this might falsely imply an accuracy of one part in &%, or perhaps even of one part in 20. The true accuracy of only one part in $% is realized only for the larger family of 21 lines with a slope of 7/15; another family of 24 lines, with a slope of -11/13, also has a greater inter-line distance than the 20-line family, since l/m > l/m. The precise way in which families of lines act at the boundaries of the unit hypercube does not seem to be an especially clean or significant criterion; however, for those people who prefer to count hyperplanes, it is possible to compute Nt using a method quite similar to the way in which we shall calculate vt (see exercise 16). *B. Theory behind the test. In order to analyze the basic set (2), we start with the observation that ( ujx + (1 + a + *. . + aj- )c mod 1 $j(x) = m > We can get rid of the mod 1 operation by extending the set periodically, making infinitely many copies of the original t-dimensional hypercube, proceed- ing in all directions. This gives us the set L= ;+k$g+kz,…,q+kt integer x, kl, lcz, . . . , Ict {( >I I = vi+ ;+kl,;+ks,…, 5 + kt> integer :c, ICI, k2, . . . , kt , r (

3.3.4 THE SPECTRAL TEST 91 essentially good to one part in v2. Similarly, let l/z+ be the maximum distance between planes, taken over all families of parallel planes that cover all points { (4m7 s(z)lm w4Ym)); we shall call ~3 the accuracy in three dimensions. The t-dimensional accuracy vt is the reciprocal of the maximum distance between hyperplanes, taken over all families of parallel (t -1)-dimensional hyperplanes that cover all points {(s/m, s(z)/m, . . . , s - (z)/m)}. The essential difference between periodic sequences and truly random se- quences that have been truncated to multiples of l/v is that the accuracy of truly random sequences is the same in all dimensions, while the accuracy of periodic sequences decreases as t increases. Indeed, since there are only m points in the t-dimensional cube when m is the period length, we can t achieve a t-dimensional accuracy of more than about milt. When the independence of t consecutive values is considered, computer- generated random numbers will behave essentially as if we took truly random numbers and truncated them to lgv, bits, where tit decreases with increasing t. In practice, such varying accuracy is usually all we need. We don t insist that the lo-dimensional accuracy be 235, in the sense that all (235)10 possible lo-tuples (Un, &x+1,. , ) should be equally likely on a 35-bit machine; for such large . . &x+9 values of t we want only a few of the leading bits of (Un, Un+i, . . . , Un+t-i) to behave as if they were independently random. On the other hand when an application demands high resolution of the random number sequence, simple linear congruential sequences will necessarily be inadequate; a generator with larger period should be used instead, even though only a small fraction of the period will actually be generated. Squaring the period will essentially square the accuracy in higher dimensions, i.e., it will double the effective number of bits of precision. The spectral test is based on the values of z,+ for small t, say 2 2 t 2 6. Dimensions 2, 3, and 4 seem to be adequate to detect important deficiencies in a sequence, but since we are considering the entire period it seems best to be somewhat cautious and go up into another dimension or two; on the other hand the values of Vt for t > 10 seem to be of no practical significance whatever. (This is fortunate, because it appears to be rather difficult to calculate Vt when t 2 10.) Note that there is a vague relation between the spectral test and the serial test; for example, a special case of the serial test, taken over the entire period as in exercise 3.3.3-19, counts the number of boxes in each of 64 subsquares of Fig. 8(a). The main difference is that the spectral test rotates the dots so as to discover the least favorable orientation. We shall return to a consideration of the serial test later in this section. It may appear at first that we should apply the spectral test only for one suitably high value of t; if a generator passes the test in three dimensions, it seems plausible that it should also pass the 2-D test, hence we might as well omit the latter. The fallacy in this reasoning occurs because we apply more stringent conditions in lower dimensions. A similar situation occurs with the serial test: Consider a generator that (quite properly) has almost the same number of points

90 RANDOM NUMBERS 3.3.4 Fig. 8. (a) The two-dimensional grid formed by all pairs of successive points (Xn, X,+1), when Xn+l = (137X, + 187) mod 256. (b) The three-dimensional grid of triplets (Xn, X=+1, Xn+z). [Illustrations courtesy of Bruce G. Baumgart.] Perhaps the most striking thing about the pattern of boxes in Fig. 8 is that we can cover them all by a fairly small number of parallel lines; indeed, there are many different families of parallel lines that will hit all the points. For example, a set of 20 nearly vertical lines will do the job, as will a set of 21 lines that tilt upward at roughly a 300 angle. We commonly observe similar patterns when driving past farmlands that have been planted in a systematic manner. If the same generator is considered in three dimensions, we obtain 256 points in a cube, obtained by appending a height component s(s(z)) to each of the 256 points (5, S(X)) in the plane of Fig. 8(a), as shown in Fig. 8(b). Let s imagine that this 3-D crystal structure has been made into a physical model, a cube that we can turn in our hands; as we rotate it, we will notice various families of parallel planes that encompass all of the points. In the words of Wallace Givens, the random numbers stay mainly in the planes. At first glance we might think that such systematic behavior is so nonrandom as to make congruential generators quite worthless; but more careful reflection, remembering that m is quite large in practice, provides a better insight. The regular structure in Fig. 8 is essentially the grain we see when examining our random numbers under a high-power microscope. If we take truly random numbers between 0 and 1, and round or truncate them to finite accuracy so that each is an integer multiple of l/v for some given number V, then the t- dimensional points (1) we obtain will have an extremely regular character when viewed through a microscope. Let l/v2 be the maximum distance between lines, taken over all families of parallel straight lines that cover the points {(z/m, s(z)/m)} in two dimen- sions. We shall call ~2 the two-dimensional accuracy of the random number generator, since the pairs of successive numbers have a fine structure that is

3.3.4 THE SPECTRAL TEST 89 3.3.4. The Spectral Test In this section we shall study an especially important way to check the quality of linear congruential random number generators; not only do all good generators pass this test, all generators now known to be bad actually fail it. Thus it is by far the most powerful test known, and it deserves particular attention. Our discussion will also bring out some fundamental limitations on the degree of ran- domness we can expect from linear congruential sequences and their generaliza- tions. The spectral test embodies aspects of both the empirical and theoretical tests studied in previous sections: it is like the theoretical tests because it deals with properties of the full period of the sequence, and it is like the empirical tests because it requires a computer program to determine the results. A. Ideas underlying the test. The most important randomness criteria seem to rely on properties of the joint distribution of t consecutive elements of the sequence, and the spectral test deals directly with this distribution. If we have a sequence (Un) of period m, the idea is to analyze the set of all m points {v&L, &x+1,. . . , Un+t-1) > (1) in t-dimensional space. For simplicity we shall assume that we have a linear congruential sequence (X~,a,c,m) of maximum period length m (so that c # 0), or that m is prime and c = 0 and the period length is m -1. In the latter case we shall add the point (0, 0, . . . ,O) to the set (l), so that there are always m points in all; this extra point has a negligible effect when m is large, and it makes the theory much simpler. Under these assumptions, (1) can be rewritten as -I~(z,s(x),s(s(z)),. . . , s - (x)) /0 5 z < m}, where s(x) = (ax + c) mod m (3) is the successor of 5. Note that we are considering only the set of all such points in t dimensions, not the order in which those points are actually generated. But the order of generation is reflected in the dependence between components of the vectors; and the spectral test studies such dependence for various dimensions t by dealing with the totality of all points (2). For example, Fig. 8 shows a typical small case in 2 and 3 dimensions, for the generator with s(z) = (1372 + 187) mod 256. (4 Of course a generator with period length 256 will hardly be random, but 256 is small enough that we can draw the diagram and gain some understanding before we turn to the larger m s that are of practical interest.

88 RANDOM NUMBERS 3.3.3 EXERCISES-Second Set In many cases, exact computations with integers are quite difficult to carry out, but we can attempt to study the probabilities that arise when we take the average over all real values of 1: instead of restricting the calculation to integer values. Although these results are only approximate, they shed some light on the subject. It is convenient to deal with numbers U,, between zero and one; for linear con- gruential sequences, U, = Xn/m, and we have U,+l = {au,, + 0}, where 0 = c/m and {z} denotes x mod 1. For example, the formula for serial correlation now becomes 1 C= x(ax+8)dx-(~1xdx)1)/(~1~2dx-(~xdx)1) b 21. [HA&B] (R. R. Coveyou.) What is the value of C in the formula just given? b 22. [M22] Let a be an integer, and let 0 5 6 < 1. If x is a real number between 0 and 1, and if s(x) = {ax + 0}, what is the probability that s(x) < x? (This is the real number analog of Theorem P.) 23. [68] The previous exercise gives the probability that Un+l < U,. What is the probability that Un+z < Un+l < Un, assuming that Un is a random real number between zero and one? 24. [A4691 Under the assumptions of the preceding problem, except with 0 = 0, show that U, > Un+l > .. > Un++i occurs with probability What is the average length of a descending run starting at U,, assuming that U, is selected at random between zero and one? b 25. [A4951 Let QI, p, Q , /3 be real numbers with 0 5 o < p < 1, 0 2 cr < /I < 1. Under the assumptions of exercise 22, what is the probability that o 5 z < p and Q 2 s(x) < /I ? (This is the real number analog of exercise 19.) 26. [MZY] Consider a Fibonacci generator, where Un+l = {U,+U,-I}. Assuming that U1 and Uz are independently chosen at random between 0 and 1, find the prob- ability that VI < Uz < Us, Ul < U3 < UZ, UZ < Ui < UX, etc. [Hint: Divide the unit square, i.e., the points of the plane {(sly) ] 0 2 x,y < l}, into six parts, depending on the relative order of x, y, and {x $ y}, and determine the area of each part.] 27. [A&Z?] In the Fibonacci generator of the preceding exercise, let UO and UI be chosen independently in the unit square except that UO > VI. Determine the probabil- ity that Vi is the beginning of an upward run of length k, so that UO > VI < . . . < Uk > Uk+l. Compare this with the corresponding probabilities for a random sequence. 28. [M%] According to Eq. 3.2.1.3-5, a linear congruential generator with potency 2 satisfies the condition X,-l -2X,, + X n+i G (a -1)c (modulo m). Consider a generator that abstracts this situation: let U,+l = {a + 2U, -U,-I}. As in exercise 26, divide the unit square into parts that show the relative order of Vi, UZ, and Us for each pair (VI, Us). Are there any values of QI for which all six possible orders are achieved with probability &, assuming that UI and U2 are chosen at random in the unit square?

3.3.3 THEORETICAL TESTS 87 10. [M%] Show that when 0 < h < k it is possible to express a(k -h, Ic, c) and 0(h, k, -c) easily in terms of a(h, k, c). 11. [ARfO] The formulas given in the text show us how to evaluate c~(h, Ic, c) when h and Ic are relatively prime and c is an integer. For the general case, prove that a) ~(dh, dk, dc) = a(h, Ic, c), integer d > 0; b) a(h, k c + 0) = G, k c) + W cl~h integer c, real 0 < 0 < 1, when h and Ic are relatively prime and hh = 1 (modulo k). 12. [A&?,$] Show that if h is relatively prime to k and c is an integer, la(h, Ic, c)l 2 (k -l)(k -2)/k. 13. [M,z~] Generalize Eq. (26) so that it gives an expression for a(h, k, c). b 14. [MZO] The linear congruential generator that has m = 235, a = 2 + 1, c = 1, was given the serial correlation test on three batches of 1000 consecutive numbers, and the result was a very high correlation, between 0.2 and 0.3, in each case. What is the serial correlation of this generator, taken over all 235 numbers of the period? 15. [A&Z] Generalize Lemma B so that it applies to all real values of c, 0 5 c < k. 16. [MZd] Given the Euclidean tableau defined in (33), let po = 1, pl = al, and p, = a,p,-1 + p,-2 for 1 < j L t. Show that the complicated portion of the sum in Theorem D can be rewritten as follows, making it possible to avoid noninteger computations: c (-1)3+1& = & c (-q + b3(C3 + C,+1)P,-1. 3 3 ll?lt ll3lt [Hint: Prove that we have ~11j17(-1)3+1/m3m,+l = (-1)7f1p,-l/mlm,+tl for 1 5 r 2 t.1 17. [A&%?] Design an algorithm that evaluates cr(h, k, c) for integers h, k, c satisfying the hypotheses of Theorem D. Your algorithm should use only integer arithmetic (of unlimited precision), and it should produce the answer in the form A + B/k where A and B are integers. (Cf. exercise 16.) If possible, use only a finite number of variables for temporary storage, instead of maintaining arrays such as al, a2, . . . , at. b 18. [A&?31 (U. Dieter.) Given positive integers h, k, z, let Show that this sum can be evaluated in closed form, in terms of generalized Dedekind sums and the sawtooth function. [Hint: When z 5 k, the quantity lj/k] -L(j -z)/kJ equals 1 for 0 2 j < z, and it equals 0 for z 5 j < k, so we can introduce this factor and sum over 0 5 j < k. ] b 19. [MZ3] Show that the serial test can be analyzed over the full period, in terms of generalized Dedekind sums, by finding a formula for the probability that a 2 X, < /3 and a 5 Xn+l < p when CY, p, CY , p are given integers with 0 5 cy < /3 2 m, 0 5 a < p 5 m. [Hint: Consider the quantity L(z -cu)/m] -L(z -P)/m].] 20. [MZ9] (U. Dieter.) Extend Theorem P by obtaining a formula for the probability that X, > Xntl > Xn+z, in terms of generalized Dedekind sums.

86 RANDOM NUMBERS 3.3.3 not on c) have small partial quotients. In particular, the result of exercise 19 implies that the serial test on pairs will be satisfactorily passed if and only if a/m has no large partial quotients. The book Dedekind Sums by Hans Rademacher and Emil Grosswald (Math. Assoc. of America, Carus Monograph No. 16, 1972) discusses the history and properties of Dedekind sums and their generalizations. Further theoretical tests, including the serial test in higher dimensions, are discussed in Section 3.3.4. EXERCISES-First Set 1. [A4101 Express Z mod y in terms of the sawtooth and 6 functions. 2. [M20] Prove the replicative law, Eq. (10). 3. [HM%%?] What is the Fourier series expansion (in terms of sines and cosines) of the function f(s) = ((x))? b 4. [Ml91 If m = lOlo, what is the highest possible value of d (in the notation of Theorem P), given that the potency of the generator is lo? 5. [A4211 Carry out the derivation of Eq. (17). 6. [A4271 Let hh + kk = 1. (a) Show, without using Lemma B, that a@, k, c) = a(h, k, 0) + 12 for all integers c 2 0. (b) Show that if 0 < j < k, (c) Under the assumptions of Lemma B, prove Eq. (21). b 7. [A424] Give a proof of the reciprocity law (19), when c = 0, by using the general reciprocity law of exercise 1.2.4-45. b 8. [M.!?4] (L. Carlitz.) Let By generalizing the method of proof used in Lemma B, prove the following beautiful identity due to H. Rademacher: If each of p, Q, r is relatively prime to the other two, (The reciprocity law for Dedekind sums, with c = 0, is the special case r = 1.) 9. [A4401 Is there a simple proof of Rademacher s identity (exercise along the lines of the proof in exercise 7 of a special case?