Web And Email Hosting, Jsp Web Hosting Blog

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

102 RANDOM NUMBERS 3.3.4 Table 1 SAMPLE RESULTS OF THE SPECTRAL TEST ht a m - 1 23 1OS + 530 530 530 530 447 2 27 + 1 235 16642 16642 16642 15602 252 3 2 8 + 1 235 34359738368 6 4 4 4 4 3141592653 235 2997222016 1026050 27822 1118 1118 5 137 256 274 30 14 6 4 6 3141592621 10 0 4577114792 1034718 62454 1776 542 7 3141592221 10 0 4293881050 276266 97450 3366 2382 8 4219755981 10 0 10721093248 2595578 49362 5868 820 9 4160984121 1010 9183801602 4615650 16686 6840 1344 10 3141592221 235 13539813818 5795090 88134 12716 2938 11 2718281829 235 12939188896 2723830 146116 10782 2914 12 513 235 i3161885770 2925242 113374 13070 2256 13 515 235 32078865098 10274746 167558 5844 2592 14 23+ 212 + 5 235 167510120 8052254 21476 16802 1630 15 23+2l3 + 5 235 168231328 5335322 21476 2008 1134 16 23 +214 + 5 235 12256151168 5733878 21476 13316 2032 17 22 + 213 + 5 235 8201443840 1830230 21476 7786 3080 18 24+ 213 + 5 235 8364058 8364058 21476 16712 1496 19 1993538883: 7 235 32300850938 705518 22270 9558 2660 20 1175245817 235 $6436418002 7362242 95306 3006 2860 21 17059465 235 ~9341117000 9476606 202796 18758 2382 22 216+3 229 536805386 118 116 116 116 23 1812433253 232 4326934538 1462856 15082 4866 906 24 1566083941 232 4659748970 2079590 44902 4652 662 25 69069 232 4243209856 2072544 52804 6990 242 26 1664525 232 4938916874 2322494 63712 4092 1038 27 314159269 2131 -1432232969 899290 36985 3427 1144 ^_ ^ _^ 28 see I:39: I (2;J -1) 1.4 X 10 643578623 12930027 837632 29 31167285 249 3.2 X 1014 4111841446 17341510 306326 59278 30 see the text 294 8.8 x lo18 6.4 x 1012 4.1 x 10 45662836 1846368 :I Thus, in six or fewer dimensions the merit is computed as follows: We might say that the multiplier a passes the spectral test if pt is 0.1 or more for 2 < t 5 6, and it passes with flying colors if pt > 1 for all these t. A low value of pt means that we have probably picked a very unfortunate multiplier, since very few lattices will have integer points so close to the origin. Conversely, a high value of /.~t means that we have found an unusually good multiplier for the given m; but it does not mean that the random numbers are necessarily very good, since m might be too small. Only the values vt truly indicate the degree of randomness.

3.3.4 THE SPECTRAL TEST 101 (Note that the first row Vi has actually gotten longer in this transformation, although eventually the rows of U should get shorter.) The next fourteen iterations of step S5 have (j, Q, qs, qs) = (2, -2, *, 0), (37% 3, *), (1, f, -10, -11, (2, -1, *, –6), (3, –1,O, $1, (1, *t 0,2), (2,0, *, -1): (3,3,4, $1, (1, *, O,O), (2, -5, *I 01, (3,1,0, *It (1, *, -3, -11, (2,0, *, 01, (3,&O, *I. Now the transformation process is stuck, but the rows of the matrices have become significantly shorter: -1479 -888874 601246 -2994234 u = -3022 -2809871 438109 (34) -227 -983 -854296 -9749816 -1707736 The search limits (zi, 22, zs) in step S8 turn out to be (O,O, l), so Us is the shortest solution to (33); we have v3 = J2272 + 9832 + 1302 zz 1017.21089. Note that only a few iterations were needed to find this value, although condition (33) looks quite difficult to deal with at first glance. All points (Un, Un+l, Un+2) produced by this random number generator lie on a family of parallel planes about 0.001 units apart. E. Ratings for various generators. So far we haven t really given a criterion that tells us whether or not a particular random number generator passes or flunks the spectral test. In fact, this depends on the application, since some applications demand higher resolution than others. It appears that vt 2 230/t for 2 2 t 5 6 will be quite adequate in most applications (although the author must admit choosing this criterion partly because 30 is conveniently divisible by 2, 3, 5, and 6). For some purposes we would like a criterion that is relatively independent of m, so we can say that a particular multiplier is good or bad with respect to the set of all other multipliers for the given m, without examining any others. A reasonable figure of merit for rating the goodness of a particular multiplier seems to be the volume of the ellipsoid in t-space defined by the relation (zim -22~ - * * . - xtut-1)2 + x; + . . . + xt < u;, since this volume tends to indicate how likely it is that nonzero integer points (zi, . . . , ~-corresponding to solutions of (15)-are in the ellipsoid. We therefore propose to calculate this volume, namely +I2 u; pt = (t/2)! m (35) as an indication of the effectiveness of the multiplier a for the given m. In this formula. (a)! = (a)(; -l)…(&G, for t odd. (36)

100 RANDOM NUMBERS 3.3.4 especially in higher dimensions. During the exhaustive search, the vector Y will always be equal to zrUr + . . . + xJJt, so that f(xl, . . . , xt) = Ye Y. Since f(-x1,. . . , -xt) = f(xi,. . . , xt), we shall examine only vectors whose first nonzero component is positive. The method is essentially that of counting in steps of one, regarding (xi, . . . , xt) as the digits in a balanced number system with mixed radices (221 + 1,. . . ,2.z2 + 1); cf. Section 4.1.) S9. [Advance xk.1 If z k = Zk, go to Sll. Otherwise increase xk by 1 and set Y+Y+Uk. SlO. [Advance k.] Set k + k+l. Then if k 5 t, set xk t -Zk, Y c Y-2zkuk, and repeat step SlO. But if k > t, set s t min(s,Y . Y). Sll. [Decrease k.] Set k c k -1. If k 2 1, return to S9. Otherwise output vt = fi (the exhaustive search is completed) and return to S4. m In practice Algorithm S is applied for T = 5 or 6, say; it usually works reasonably well when T = 7 or 8, but it can be terribly slow when T 2 9 since the exhaustive search tends to make the running time grow as 3T. (If the minimum value vt occurs at many different points, the exhaustive search will hit them all; hence we typically find that all zk = 1 for large t. As remarked above, the values of vt are generally irrelevant for practical purposes when t is large.) An example will help to make Algorithm S clear. Consider the linear congruential sequence defined by m = lOlo, a = 3141592621, c = 1, x0 = 0. (32) Six cycles of the Euclidean algorithm in steps S2 and S3 suffice to prove that the minimum nonzero value of x1 + xi with xi + 3141592621×2 = 0 (modulo lOlo) occurs for x1 = 67654, 22 = 226; hence the two-dimensional accuracy of this generator is v2 = d676542 + 2262 e 67654.37748. Passing to three dimensions, we seek the minimum nonzero value of xy +xz +x: such that x1 + 3141592621×2 + 31415926212xs = 0 (modulo lOlo). (33) Step S4 sets up the matrices -67654 -226 0 -191 -44190611 2564918569 u = -44190611 191 0 , v = -226 67654 1307181134 . ( 5793866 33 11 ( 0 0 10000000000 1 The first iteration of step S5, with q = 1 for i = 2 and q = 4 for i = 3, changes them to -21082801 97 4 -191 -44190611 2564918569 u = -44190611 191 0 -35 44258265 -1257737435 . 5793866 33 1 764 176762444 -259674276

3.3.4 THE SPECTRAL TEST 99 S3. [Compute vs.] Set u c u -h, v t v -p; and if u2 + v2 < s, set s +- u2 + v2, h +- u, p + v. Then output fi = v2. (The validity of this calculation for the two-dimensional case is proved in exercise 5. Now we will set up the U and V matrices satisfying (28) and (29), in preparation for calculations in higher dimensions.) Set where the -sign is chosen for V if and only if p > 0. S4. [Advance t.] If t = T, the algorithm terminates. (Otherwise we want to increase t by 1. At this point U and V are t X t matrices satisfying (28) and (29) and we must enlarge them by adding an appropriate new row and column.) Set t t t + 1 and r e (or) modm. Set Ut to the new row (–r, 0, 0, . . . ,O, 1) of t elements, and set uit t-0 for 1 5 i < t. Set V, to the new row (O,O, 0,. . . ,O, m). Finally, for 1 5 i < t, set g c round(v,l r/m), vzt + vile -gm, and Ut t Ut + qUi. (Here round(z) denotes the nearest integer to z, e.g., 1~ + l/2]. We are essentially setting vit + viir and immediately applying transformation (23) with j = t, since the numbers ]vilr] are so large they ought to be reduced at once.) Finally set s c min(s, U,. Ut), k t t, and j + 1. (In the following steps, j denotes the current row index for transformation (23), and Ic denotes the last such index where the transformation shortened at least one of the Vi.) S5. [Transform.] For 1 5 i 5 t, do the following operations: If i # j and 2lVi . Vjl > Vy . Vy, set q + round(Vi . Vj / I+ . Vj), Vi +-Vi -qV7, Uy + Vi + qUi, and k + j. (The fact that we omit this transformation, when 21Vi . Vjl exactly equals V, . V,, prevents the algorithm from looping endlessly; see exercise 19.) S6. [Examine new bound.] If k = j (i.e., if the transformation in S5 has just done something useful), set s +- min(s, Uj . Uj). S7. [Advance j.] If j = t, set j c 1; otherwise set j c j + 1. Now if j # k, return to step S5. (If j = k, we have gone through t -1 consecutive cycles of no transformation, so the transformation process is stuck.) S8. [Prepare for search.] (Now the absolute minimum will be determined, using an exhaustive search over all (Q, . . . , Q) satisfying condition (21) of Lemma A.) Set X t Y t (0, . . . ,O), set k t t, and set (We will examine all X = (51, . . . ,zt) with ]zrj] 2 zj for 1 5 j 5 t. In hundreds of applications of this algorithm, no zJ has yet turned out to be greater than 1, nor has the exhaustive search in the following steps ever reduced s; however, such phenomena are probably possible in weird cases,

98 RANDOM NUMBERS 3.3.4 be quite small in most cases. Occasionally these bounds (21) will be poor, and another type of transformation will usually get the algorithm unstuck again and reduce the bounds (see exercise 18). However, transformation (23) by itself has proved to be quite adequate for the spectral test; in fact, it has proved to be amazingly powerful when the computations are arranged as in the algorithm discussed below. *D. How to perform the spectral test. Here now is an efficient computational procedure that follows from our considerations. R. W. Gosper and U. Dieter have observed that it is possible to use the results of lower dimensions to make the spectral test significantly faster in higher dimensions. This refinement has been incorporated into the following algorithm, together with a significant simp- lification of the two-dimensional case. Algorithm S (The spectral test). This algorithm determines the value of v,=rnin{fi+…+xz ).~+.~~+…+at- s~O(modulorn) } (27) for 2 2 t < T, given a, m, and T, where 0 < a < m and a is relatively prime to m. (The number vt measures the t-dimensional accuracy of random number generators, as discussed in the text above.) All arithmetic within this algorithm is done on integers whose magnitudes rarely if ever exceed m2, except in step S8; in fact, nearly all of the integer variables will be less than m in absolute value during the computation. When vt is being calculated for t 2 3, the algorithm works with two t x t matrices U and V, whose row vectors are denoted by Vi = (uir, . . . , uit) and vi = (Wil,…, vit) for 1 5 i 5 t. These vectors satisfy the conditions uil + uuis + 1.. + ut-ruit = 0 (modulo m), 1

3.3.4 THE SPECTRAL TEST 97 diagram makes this plain. (25) Turning to question (b), we want to choose the q2 so that U, +&+ qiiIJi has minimum length; geometrically, we want to start with Uj and add some vector in the (t -1)-dimensional hyperplane whose points are the sums of multiples of { Ui 1 i # j }. Again the best solution is to choose things so that Uj’ is perpendicular to the hyperplane, i.e., so that Uj . Uk = 0 for all k # j, i.e., U, . Uk + C qi(Ui ‘ uk) O? 1

96 RANDOM NUMBERS 3.3.4 has the same answer. For example, Euclid s algorithm has this form; if we don t know the gcd of the input numbers, we change them into smaller numbers having the same gtd. (In fact, a slightly more general approach probably underlies the discovery of nearly all algorithms: If you can t solve a problem directly, change it into one or more simpler problems, from whose solution you can solve the original one. ) In our case, a simpler problem is one that requires less searching because the right-hand side of (22) is smaller. The key idea we shall use is that it is possible to change one quadratic form into another one that is equivalent for all practical purposes. Let j be any fixed subscript, 1 5 j 5 t; let (qi, . . . , CJ~-…1, . . . , qt) oj+l, be any sequence of t - 1 integers; and consider the following transformation of the vectors: Vi’ = vi -qivj, xi = xi -qixj, Vi’ = vi, for i # j; Vj = Vj, X[1 = Xj, Uj = Uj + Cifj q%Ui. (23) It is easy to see that the new vectors VI , . . . , U, define a quadratic form f for which f/(x$, . . . ,xi) = f(xl,. . . , xt); furthermore the basic orthogonality condition (19) remains valid, because it is easy to check that Vi . Vj = Sij. As (Xl,…, xt) runs through all nonzero integer vectors, so does (xi,. . . , xi); hence the new form j has the same minimum as f. Our goal is to use transformation (23), replacing Vi by Vi and Vi by Vi for all i, in order to make the right-hand side of (22) small; and the right-hand side of (22) will be small when both Uj . Uj and Vk . vk are small. Therefore it is natural to ask the following two questions about the transformation (23): a) What choice of qi makes Vi . Vi as small as possible? b) What choice of 41, . . . , qj-1, qj+l, . . . , qt makes Uj . Uj as small as possible? It is easiest to solve these questions first for real values of the qi. Question (a) is quite simple, since and the minimum occurs when qi=Vi Vj/Vj Vj. (24) Geometrically, we are asking what multiple of Vj should be subtracted from Vi so that the resulting vector Vi has minimum length, and the answer is to choose qi so that Vi is perpendicular to Vj (i.e., so that Vi . Vj = 0); the following