Web And Email Hosting, Jsp Web Hosting Blog

20 RANDOM NUMBERS 3.2.1.2 Theorem C. The number a is a primitive element module p if and only if i) pe = 2, a is odd; or pe = 4, amod = 3; orpe = 8, amod = 3,5,7; orp= 2, e 2 4, amod = 3,5; Or ii) p is odd, e = 1, a $ 0 (modulo p), and a(P- )lq $ :l (modulo p) for any prime divisor q of p -1; Or iii) p is odd, e > 1, a satisfies (ii), and ap- JZ 1 (modulo p2). 1 Conditions (ii) and (iii) of this theorem are readily tested on a computer for large values of p, by using the efficient methods for evaluating powers discussed in Section 4.6.3. Theorem C applies to powers of primes only; if we are given values a3 that are primitive modulo p;, it is possible to find a single value a such that a E a3 (modulo p;,~), for 1 5 j 5 t, using the Chinese remainder algorithm discussed in Se&on 4.3.2, and this number a will be a primitive element modulo pi . . . p, . Hence there is a reasonably efficient way to construct multipliers satisfying the condition of Theorem B, for any desired value of m, although the calculations can be somewhat lengthy in the general case. In the common case m = 2e, with e 2 4, the conditions above simplify to the single requirement that a G 3 or 5 (modulo . In this case, one-fourth of all possible multipliers give the maximum period. The second most common case is when m = 10e. Using Lemmas P and Q, it is not difficult to obtain necessary and sufficient conditions for the achievement of the maximum period in the case of a decimal computer (cf. exercise 18): Theorem D. If m = 10e, e 2 5, c = 0, and X0 is not a multiple of 2 or 5, the period of the linear congruential sequence is 5 X 10ee2 if and only if a mod 200 equals one of the following 32 values: 3, 11, 13, 19, 21, 27, 29, 37, 53, 59, 61, 67, 69, 77, 83, 91, 109, 117, 123, 131, 133, 139, 141, 147, 163, 171, 173, 179, 181, 187, 189, 197. 1 (10) EXERCISES 1. [IO] What is the length of the period of the linear congruential sequence with X0 = 5772156648, a=3141592621, c= 2718281829,and m= 10000000000? 2. [IO] Are the following two conditions sufficient to guarantee the maximum length period, when m is a power of 2? (i) c is odd; (ii) a mod 4 = 1. 3. [IS] Suppose that m = 10e, where e 2 2, and suppose further that c is odd and not a multiple of 5. Show that the linear congruential sequence will have the maximum length period if and only if a mod 20 = 1.

3.2.1.2 CHOICE OF MULTIPLIER 19 It may be possible to achieve an acceptably long period even if we stipulate that c = 0. Let us now try to find conditions on the multiplier so that the period is as long as possible in this special case. According to Lemma Q, the period of the sequence depends entirely on the periods of the sequences when m = pe, so let us consider that situation. We have X, = anXo modpe, and it is clear that the period will be of length 1 if a is a multiple of p, so we take a to be relatively prime to p. Then the period is the smallest integer X such that X0 = aXXo modpe. If the greatest common divisor of X0 and pe is pf, this condition is equivalent to ux = 1 (modulo p -f). (8) By Euler s theorem (exercise 1.2.4-28), u~(P -~) = 1 (modulo p -f); hence X is a divisor of p(p -f) = pe-f - (p -1). When a is relatively prime to m, the smallest integer X for which ax = 1 (modulo m) is conventionally called the order of a module m. Any such value of a that has the maximum possible order modulo m is called a primitive element modulo m. Let x(m) denote the order of a primitive element, i.e., the maximum possible order, modulo m. The remarks above show that X(p ) is a divisor of peP1(p- 1); with a little care (see exercises 11 through 16 below) we can give the precise value of x(m) in all cases as follows: x(2) = 1, X(4) = 2, x(27 = 2e-2 if e 2 3. X(P ) = Pe-YP -11, if p>2. (9) X(pi . . .pz ) = lcm(X(pil), . . . , X(pft)). Our remarks may be summarized in the following theorem: Theorem B. [R. D. Carmichael, Bull. Amer. Math. Sot. 16 (1910), 232-238.1 The maximum period possible when c = 0 is x(m), where h(m) is defined in (9). This period is achieved if i) X0 is relatively prime to m; ii) u is a primitive element module m. 1 Note that we can obtain a period of length m -1 if m is prime; this is just one less than the maximum length, so for all practical purposes such a period is as long as we want. The question now is, how can we find primitive elements modulo m? The exercises at the close of this section tell us that there is a fairly simple answer when m is prime or a power of a prime, namely the results stated in our next theorem.

RANDOM NUMBERS 3.2.1.2 Lemma R. Assume that 1 < a < pe, where p is prime. If X is the smallest positive integer for which (aX - l)/(a -1) s 0 (modulo p ), then a s 1 (modulo p) when p > 2, x = pe if and only if a z 1 (modulo 4) when p= 2. Proof. Assume that X = pe. If a $ 1 (modulo p), then (an -l)/(a -1) G 0 (modulo p ) if and only if an -1 E 0 (modulo p ). The condition ape - 1 E 0 (modulo p ) then implies that ape E 1 (modulo p); but by Theorem 1.2.4F we have ape s a (modulo p), hence a $ 1 (modulo p) leads to a contradiction. And if p = 2 and a zz 3 (modulo 4), we have (aze- -l)/(a -1) E 0 (modulo 2e) by exercise 8. These arguments show that it is necessary in general to have a = 1 + qpf , where pf > 2 and q is not a multiple of p, whenever X = pe. It remains to be shown that this condition is sufficient to make X = pe. By repeated application of Lemma P, we find that apg E 1 (modulo P~+~), apg f 1 (modulo P~+~+ ), for all g 2 0, and therefore (up9 -l)/(a -1) E 0 (modulo pg) , (apg -l)/(a -1) $ 0 (modulo pgfl). (6) In particular, (a pe - l)/(a-1) 3 0 (modulo p ). Now the congruential sequence (0, a, 1, p ) has X, = (a -l)/(a-1) modp ; therefore it has a period of length X, that is, X, = 0 if and only if n is a multiple of X. Hence pe is a multiple of X. This can happen only if X = pg for some g, and the relations in (6) imply that X = pe, completing the proof. I The proof of Theorem A is now complete. 4 We will conclude this section by considering the special case of pure mul- tiplicative generators, when c = 0. Although the random number generation process is slightly faster in this case, Theorem A shows us that the maximum period length cannot be achieved. In fact, this is quite obvious, since the sequence now satisfies the relation X n+l = ax, mod m, (7) and the value X, = 0 should never appear, lest the sequence degenerate to zero. In general, if d is any divisor of m and if X, is a multiple of d, all succeeding elements X,+1, Xn+2, . . . of the multiplicative sequence will be multiples of d. So when c = 0, we will want X, to be relatively prime to m for all n, and this limits the length of the period to at most p(m), the number of integers between 0 and m that are relatively prime to m.

3.2.1.2 CHOICE OF MULTIPLIER 17 The length X of the period of the linear congruential sequence determined by (X0, a, c, m) is the least common multiple of the lengths Xj of the periods of the linear congruential sequences (X0 mod p$ , a mod p;3f c mod pT3, p; ), 1 < j 5 t. Proof. By induction on t, it suffices to prove that if ml and m2 are relatively prime, the length X of the linear congruential sequence determined by the param- eters (X0, a, c, mlm2) is the least common multiple of the lengths X1 and X2 of the periods of the sequences determined by (X0 mod ml, a mod ml, c mod ml, ml) and (X0 mod m2, a mod m2, c mod m2, m2). We observed in the previous section, Eq. (4), that if the elements of these three sequences are respectively denoted by X,, Yn, and 2,) we will have Y, = X, mod ml and 2, =Xnmodm2, for all n > 0. Therefore, by Law D of Section 1.2.4, we find that x, =xk if and only if Y,, = Y, and 2, = 2,. (4 Let X be the least common multiple of X1 and X2; we wish to prove that A = X. Since X, = Xn+x for all suitably large n, we have Y, = Y,+x (hence X is a multiple of X1) and 2, = &+A (hence X is a multiple of X2), so we must have X 2 X . Furthermore, we know that Y, = Yn+x, and 2, = &+A, for all suitably large n; therefore, by (4), X, = X,+X,. This proves X 5 X . 1 Now we are ready to prove Theorem A. Because of Lemma Q, it suffices to prove the theorem when m is a power of a prime number. For p l . . . pi = X = lcm(X1,. . . ,X,) 5 X1 . . . Xt 5 p: . . . pz can be true if and only if X, = p; for 1 5 j 5 t. Therefore, assume that m = pe, where p is prime and e is a positive integer. The theorem is obviously true when a = 1, so we may take a > 1. The period can be of length m if and only if each possible integer 0 5 x < m occurs in the period, since no value occurs in the period more than once. Therefore the period is of length m if and only if the period of the sequence with X0 = 0 is of length m, and we are justified in supposing that X0 = 0. By formula 3.2.1-6 we have c mod m. If c is not relatively prime to m, this value X, could never be equal to 1, so condition (i) of the theorem is necessary. The period has length m if and only if the smallest positive value of n for which X, = X0 = 0 is n = m. By (5) and condition (i), our theorem now reduces to proving the following fact:

16 RANDOM NUMBERS 3.2.1.2 Theorem A. The linear congruential sequence defined by m, a, c, and X0 has period length m if and only if i) c is relatively prime to m; ii) b = a - 1 is a multiple of p, for every prime p dividing m; iii) b is a multiple of 4, if m is a multiple of 4. The ideas used in the proof of this theorem go back at least a hundred years. The first proof of the theorem in this particular form was given by M. Greenberger in the special case m = 2e [see JACM 8 (1961), 383-3891, and the sufficiency of conditions (i), (ii), and (iii) in the general case was shown by Hull and Dobell [see SIAM Review 4 (1962), 230-2541. To prove the theorem we will first consider some auxiliary number-theoretic results that are of interest in themselves. Lemma P. Let p be a prime number, and Jet e be a positive integer, where pe > 2. If z E 1 (modulo p ), z $ 1 (modulo pe+ ), (1) then z* E 1 (modulo pefl), x* $ 1 (modulo pefa). (2) Proof. We have x = 1-t qpe for some integer q that is not a multiple of p. By the binomial formula x*=1+ I 4Pe + *** + p q*-lp(*-l)e + q*p*e 0 ( P–l > =1+qpe+l ( 1 +;(;)qp~+~~)q~p~~+…+~(;)q*-~pb-l)~). The quantity in parentheses is an integer, and, in fact, every term inside the parentheses is a multiple of p except the first term. For if 1 < k < p, the binomial coefficient (3 is divisible by p (cf. exercise 1.2.6-lo), hence is divisible by P( - )~; and the last term is q*-lp(*-l)e- , which is divisible by p since (p - 1)e > 1 when pe > 2. So zp E 1 + qpe+ (modulo pef2), and this completes the proof. (Note: A generalization of this result appears in exercise 3.2.2-11(a).) 1 Lemma Q. Let the decomposition of m into prime factors be m = pt . . . p, . (3)

3.2.1.2 CHOICE OF MULTIPLIER 15 3. [,%o] How can the constant (w -1) be specified in general in the MIX assembly language, regardless of the value of the byte size? b 4. [Z] Discuss the calculation of linear congruential sequences with m = 232 on two s-complement machines such as the System/370 series. 5. [Zoo] Given that m is less than the word size, and that 2, y are nonnegative integers less than m, show that the difference (Z -y) mod m may be computed in just four MIX instructions, without requiring any division. What is the best code for the sum (x + y) mod m? b 6. [,%%I The previous exercise suggests that subtraction mod m is easier to perform than addition mod m. Discuss sequences generated by the rule Xn+l = (ax, -c) mod m. Are these sequences essentially different from linear congruential sequences as defined in the text? Are they more suited to efficient computer calculation? 7. [MzL$] What patterns can you spot in Table l? b 8. [.20] Write a MIX program analogous to (2) that computes (ax) mod (w -1). The values 0 and w -1 are to be treated as equivalent in the input and output of your program. 9. [~?3] Write a MIX program analogous to the one in exercise 8, but it should compute (ax) mod (w -2) . 3.2.1.2. Choice of multiplier. In this section we shall show how to choose the multiplier a so as to give the period of maximum length.. A long period is essential for any sequence that is to be used as a source of random numbers; indeed, we would hope that the period contains considerably more numbers than will ever be used in a single application. Therefore we shall concern ourselves in this section with the question of period length. The reader should keep in mind, however, that a long period is only one desirable criterion for the randomness of our sequence. For example, when a = c = 1, the sequence is simply X,+, = (X, + 1) mod m, and this obviously has a period of length m, yet it is anything but random. Other considerations affecting the choice of a multiplier will be given later in this chapter. Since only m different values are possible, the period surely cannot be longer than m. Can we achieve the maximum length, m? The example above shows that it is always possible, although the choice a = c = 1 does not yield a desirable sequence. Let us investigate all possible choices of a, c, and X0 that give a period of length m. It turns out that all such values of the parameters can be characterized very simply; when m is the product of distinct primes, only a = 1 will produce the full period, but when m is divisible by a high power of some prime there is considerable latitude in the choice of a. The following theorem makes it easy to tell if the maximum period is achieved.

14 RANDOMNUMBERS 3.2.1.1 the numbers Y, = X, mod24. The gist of Eq. (4) is that the low-order four bits of (Xn) form a congruential sequence that has a period of length 16 or less. Similarly, the low-order five bits are periodic with a period of at most 32; and the least significant bit of X, is either constant or strictly alternating. This situation does not occur when m = w f 1; in such a case, the low-order bits of X, will behave just as randomly as the high-order bits do. If, for example, w = 235 and m = 235 - 1, the numbers of the sequence will not be very random if we consider only their remainders mod 31, 71, 127, or 122921 (cf. Table 1); but the low-order bit, which represents the numbers of the sequence taken mod 2, would be satisfactorily random. Another alternative is to let m be the largest prime number less than w. This prime may be found by using the techniques of Section 4.5.4, and a table of suitably large primes appears in that section. In most applications, the low-order bits are insignificant, and the choice m = w is quite satisfactory-provided that the programmer using the random numbers does so wisely. Our discussion so far has been based on a signed magnitude computer like MIX. Similar ideas apply to machines that use complement notations, although there are some instructive variations. For example, a DEC 20 computer has 36 bits with two s complement arithmetic; when it computes the product of two nonnegative integers, the lower half contains the least significant 35 bits with a plus sign. On this machine we should therefore take w = 235, not 236. The 32-bit two s complement arithmetic on IBM System/370 computers is different: the lower half of a product contains a full 32 bits. Some programmers have felt that this is a disadvantage, since the lower half can be negative when the operands are positive, and it is a nuisance to correct this; but actually it is a distinct advantage from the standpoint of random number generation, since we can take m = 232 instead of 231 (see exercise 4). EXERCISES 1. [MLZ] In exercise 3.2.1-3 we concluded that the best congruential generators will have the multiplier a relatively prime to m. Show that when m = w in this case it is possible to compute (ax + c)modw in just three MIX instructions, rather than the four in (1) with the result appearing in register X. 2. [16] Write a MIX subroutine having the following characteristics: Calling sequence: JMP RANDM Entry conditions: Location XRAND contains an integer x. Exit conditions: X + rA + (ax + c) mod 20, rX t 0, overflow off. (Thus a call on this subroutine will produce the next random number of a linear congruential sequence.)

3.2.1.1 Table 1 PRIME FACTORIZATIONS - 2 -1 7 31 151 3 5.17 257 131071 33 7 .19 73 524287 3 52 11 31 .41 7= 127. 337 3 23 89.683 47.178481 32 . 5.7 13.17. 241 31 ,601 1801 3.2731. 8191 7. 73.262657 3.5.29.43.113.127 233.110312089 32 7.11 31 .151 .331 2147483647 3.5.17 ,257 .65537 7 23.89 599479 3.43691.131071 31 71 .127.122921 33 5.7.13 19 * 37 73 109 223 616318177 3.174763 524287 7 79.8191 121369 3 52 11 17 31 41 61681 13367.164511 353 32 7 .43 127 ,337 .5419 431 .9719.2099863 3.5.23.89.397.683.2113 7.31 73. 151.631 23311 3.47 178481 .2796203 2351 4513.13264529 32 5 7. 13. 17.97. 241. 257.673 179951.3203431780337 3= 5= .7.11 . 13. 31. 41 61. 151 331 1321 7 73.127.337.92737.649657 3 5.17 257 641 65537.6700417 lae - 1 33 7.11 13 37 32 .239 4649 32 .ll 73.101.137 34. 37.333667 32 . 11 41 271.9091 3= .21649.513239 33.7.11.13 37.101.9901 32 11 .17.73.101.137.5882353 - CHOICE OF MODULUS 13 OF e 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 59 w f 1 2? + 1 32 11.331 65537 3 43691 5 13 37 109 3.174763 17 61681 32 43 5419 5 397.2113 3.2796203 97.257 673 3.11.251.4051 5.53.157.1613 34 19.87211 17.15790321 3 59.3033169 5 13.41 61 3.715827883 641.6700417 1321 60 17.241 61681 4562284561 63 33. 19 4315419.77158673929 64 274177.67280421310721 32 67 683.20857 5. 137.953 26317 3. 11 43.281 86171 17.241 433. 38737 3. 1777 25781083 5 229.457 525313 3 2731 22366891 257.4278255361 3 83.8831418697 5 13 29 113 1429 14449 3 2932031007403 17 353 2931542417 33. 11 19.331. 18837001 5.277.1013.1657.30269 3 ,283 165768537521 193.65537 22253377 3 .2833.37171 1824726041 e 10 + 1 6 101 9901 7 11*909091 f 8 17.5882353 9 7.11 13.19 10 -101 .3541.27961 52579 11 112 .23.4093.8779 12 73 137 I 99990001 16 353.449.641 .1409.69857 1

12 RANDOM NUMBERS 3.2.1.1 The following program does this: 01 LDAN X rA + -X. 02 MUL A rAX t @A). a. 03 STX TEMP 04 SUB TEMP rA+rA-rX. ( 4 05 JANN *+3 Exit if rA > 0. 06 INCA 2 rA+rA+2. 07 ADD =w -1= rA + rA f w -1. (Cf. exercise 3.) 1 Register A now contains the value (ax) mod (W + 1). Of course, this value might lie anywhere between 0 and w, inclusive, so the reader may legitimately wonder how we can represent so many values in the A-register! (The register obviously cannot hold a number larger than w -1.) The answer is that overflow will be on after the above program if and only if the result equals w, assuming that overflow was initially off. It is convenient simply to reject the value w if it appears in the congruential sequence modulo w + 1; this will happen if lines 05 and 06 of (2) are replaced by JANN ++4; INCA 2; JAP *-5 . To prove that code (2) actually does determine (a~) mod (w + l), note that in line 04 we are subtracting the lower half of the product from the upper half. No overflow can occur at this step; and if aX = qw + r, with 0 2 r < w, we will have the quantity r -q in register A after line 04. Now ax = q(w + 1) + (r - 91, and since q < w, we have -w < r -q < w; hence (ax) mod(w + 1) equals either r -q or r -q + (w + l), depending on whether r -q 2 0 or r -q < 0. A similar technique can be used to get the product of two numbers modulo (w -1); see exercise 8. In later sections we shall require a knowledge of the prime factors of m in order to choose the multiplier a correctly. Table 1 lists the complete factorization of w f 1 into primes for nearly every known computer word size; the methods of Section 4.5.4 can be used to extend this table if desired. The reader may well ask why we bother to consider using m = w f 1, when the choice m = w is so manifestly convenient. The reason is t,hat when m = w, the right-hand digits of X, are much less random than the left-hand digits. If d is a divisor of m. and if Y, = X, mod d, (3) we can easily show that Yn+l = (aY, + c) mod d. (4 (For, Xn+l = ax, + c - qm for some integer q, and taking both sides mod d causes the quantity qm to drop out when d is a factor of m.) To illustrate the significance of Eq. (4), let us suppose, for example, that we have a binary computer. If m = 20 = 2e, the low-order four bits of X, are

3.2.1.1 CHOICE OF MODULUS 11 b 2. [M.zo] Show that if a and m are relatively prime, the number X0 will always appear in the period. 3. [Ml01 If a and m are not relatively prime, explain why the sequence will be somewhat handicapped and probably not very random; hence we will generally want the multiplier a to be relatively prime to the modulus m. 4. [II] Prove Eq. (6). 5. [M?O] Equation (6) holds for k > 0. If possible, give a formula that expresses Xn+k in terms of X, for negative values of k. 3.2.1.1. Choice of modulus. Our current goal is to find good values for the parameters that define a linear congruential sequence. Let us first consider the proper choice of the number m. We want m to be rather large, since the period cannot have more than m elements. (Even if a person wants to generate only random zeros and ones, he should not take m = 2, for then the sequence would at best have the form. . . , 0, 1, 0, 1, 0, 1, . . . ! Methods for modifying random numbers to get random zeros and ones are discussed in Section 3.4. ) Another factor that influences our choice of m is speed of generation: We want to pick a value so that the computation of (ax, + c) mod m is fast. Consider MIX as an example. We can compute ymod m by putting y in registers A and X and dividing by m; assuming that y and m are positive, we see that y mod m will then appear in register X. But division is a comparatively slow operation, and it can be avoided if we take m to be a value that is especially s convenient, such as the word size of our computer. Let w be the computer s word size, namely, 2e on an e-bit binary computer or 10e on an e-digit decimal machine. (In this book we shall often use the letter e to denote an arbitrary integer exponent, instead of the base of natural logarithms, hoping that the context will make our notation unambiguous. Physicists have a similar problem when they use e for the charge on an electron.) The result of an addition operation is usually given modulo w, except on ones -complement machines; and multiplication mod w is also quite simple, since the desired result is the lower half of the product. Thus, the following program computes the quantity (ax + c) mod w efficiently: LDA A rA + a. MUL X rAX + (rA).X. SLAX 5 rA t rAXmodw. (1) ADD C rAt(rA+c)modw. m The result appears in register A. The overflow toggle might be on at the conclu- sion of the above sequence of instructions, and if this is undesirable, the code should be followed by, e.g., JOV *+l to turn it off. A clever technique that is less commonly known can be used to perform computations modulo (w + 1). For reasons to be explained later, we will generally want c = 0 when m = w + 1, so we merely need to compute (ax) mod (w + 1).