ifii txHa m ■ If! mm Brail MM ill m Jfnjfl I mi L I B R.ARY OF THE "• -i UNIVLR.5ITY Of ILLINOIS . 510. 84 no.324-33d cop. 2 * * *> The person charging this material is re- sponsible for its return on or before the Latest Date stamped below. Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. University of Illinois Library m: SEP 3 1977 MAY 2 9 RET |M* 1 4 1978 OCT 2 \ 1970 NOV 1 4 1970 NOV j'iak NOV 1 2 #74 1 i JIN < fS7f OEC 1 9 REft 004 1371 FEB MAR 28 197 APR 2 Nftl 1 JUL 1 19?S MAR 8 1976 MAR 8 r£d ■ to 1377 *f*4 .APR A-WECH L161— O-1096 Digitized by the Internet Archive in 2013 http://archive.org/details/randomnumbergene329rich J3^ c ? tw5 ±wro No. 329 April , 1969 RANDOM NUMBER GENERATION ON THE I.B.M. 360 Beth Carson Richardson JIM ! 6 f§;g RANDOM NUMBER GENERATION ON THE I.B.M. 3^0 BY BETH CARSON RICHARDSON B.S., University of Illinois, 1966 THESIS Submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science in the Graduate College of the University of Illinois, 19&9 Urbana, Illinois > a •P I «i w d o m o •H a o cq .d -P CD pq •H 5 a H 4*- X tf I t M co • CQ OJ + CQ OJ CO M H II % OJ o II CO OJ W H o JL, o o H OJ OJ X O -r-D H OJ M ^i H O co ^H CQ t>0 o OJ <+H 11 X + •H 01 <^™ ' V CQ > a III OJ O H OO s — ' X II M d •H ^J OJ o II CO OJ X! o o o H OJ OJ t*i o o H H II OJ M d •H 0\ H OJ H H OJ VO OJ CD bO t— 00 O H -cj- H ir\ 0- H r O CO o K H (-3 OJ X I OJ OJ H II •H X CD H o3 H CD O CO LTN CO ~H" ir\ LT\| CO II OJ V/ H i M + •H o5 A + •H A + •H "ST +3 CO O a s o3 ■P •H ■s o3 O ft o5 o 0) CO ft VO H rH ON on o\ on ... ONMD LTN ON ON ON ON ^ OJ mo ON ON ON ON ,,. OOONOJ ON ON ON ON ^ OJ OJ o CO IB ■ OJ OJ X M CD H c3 H CD O CO LTN CO LTN CO O LTN OJ V H i + •H 05 A Pi + •H A + •H "ft" -P •H •p •s ro -H Q fl 1 0} •\ £ Tf o a ft 1 U O w ft MD H rH ON ON ON ON ONVO Ln ON ON ON ON OJ 00 O ON ON ON ON CO ONOJ ON ON ON ON OJ OJ OJ •H iH a, OJ CO H OJ OJ -4" OJ ON OJ OJ m OJ kn. OJ VO OJ OJ OJ I 00 LfN no -=1- 21 o H i-3 EH i O H AV ** H 52! a a • |oj a £ CO CD H fl ft a3 S 3 ft o P S ■P •s cd a ■H d° ,d •H •p -d 1) S 1 ft d Ti a CD ■H H ,d f> ^ -p -d £ a H 0) o 1 n II & o i j I O d o ■P & CQ •N o H (U CD • iaq K !> CQ 2 5°> 000 > \ +2 = 2 5°° Other drawbacks are that the method is slow and the period is short. For these reasons, the mid-square method was superseded by the congruential method. The congruential method is based on the relation x. , = ax.+c(mod m) which means that ax.+c is to be divided by m and x. n X+X X X X+X n set equal to the remainder. Lehmer first suggested this method in 19^-9 with c=0 in which case it is called the multiplicative congruential method. Q The mixed congruential method with c/^0 was later introduced by Rotenberg. The modulus m, the starting value x n , the multiplier a and the constant c are chosen so as to provide a long period and good statistical behavior. The problem of choosing constants is discussed in section two. Since x. , , i=0,l, ... is equal to the remainder of ax.+c upon division by the modulus m, each generated number lies somewhere in the range (0,1, . . . ,m-l) so the period of a generator is less than or equal to m. In this regard the mixed method has an apparent advantage over the multiplicative method because for a suitable choice of constants the full period m can be achieved whereas the full period can never be achieved for the multiplicative 9 10 method. Unfortunately, extensive testing ' has shown that the multiplicative method generally behaves better statistically than the mixed method. The poor statistical behavior of the mixed method led Hull and Dobell to conclude that multiplicative methods were to be preferred. In an effort to speed up congruent ial methods even more, the multiplication operation was replaced by an addition so that the re- currence relation became x. , = x.+x. (mod m) . This relation is known l+l i i-n N as the additive method and when n = 1, x = and x, = 1 one has the Fibonacci generator. Although the Fibonacci generator is faster and has a longer period than the multiplicative method, testing has shown that it gives poorer statistical results. When n was increased the generator produced numbers with better statistical properties but the program required indexing so the speed advantage of the additive method was 12 lost, 13 MacLaren and Marsaglia found mixed congruential generators to have poor statistical qualities when used in Monte Carlo calculations and so proposed two improved congruential generators. One of the genera- tors uses a stored table of random numbers. When a number is used it is replaced by a scrambled version of itself provided by the relation x. = ax.+c(mod m) . This method eliminates the main problem of stored tables of numbers, that of exhausting the list, but requires some inconvenient tape -handling. The second generator is designed to improve the distribution of n-tuples of random numbers x,,...,x . A list of 128 numbers u_,...,u p n is produced by the relation u . _ = a 1 u_+c..(mo& m) , A second set of numbers v, ,...,v where v. , = a v.+c (mod m) is also 1' ' n l+l 2 l 2. needed. The first seven bits of v are used to select the k'th random number, x , from the list of the u's. This position is then refilled with v. , . Although the numbers produced by this generator were shown J£+_L to have "better statistical properties than numbers produced "by the mixed or multiplicative methods, this generator has other problems. Its theoretical properties are unknown and the generation time is about twice that of the mixed or multiplicative methods because two numbers must he produced to generate every one random number. In conclusion, it appears that there is no generator that provides random numbers with the fastest generation time, the longest period, and the best statistical properties. The generator which appears best to combine these three properties, though, is the multi- plicative congruential generator. FOOTNOTES 1. R. D. Richtmyer, "Monte Carlo Methods," Nuclear Reactor Theory Am . Math . Soc , Proc . Symposium Appl. Math. 11 (1961), pp. 190-205. 2. Student, "The Probable Error of a Mean," Biometrika 6 (1908), PP. 1-25. 3. L. H. G. Tippett, Random Sampling Numbers , Tracts for Computers, no. 15 (Cambridge, ~1927y. h. M. G. Kendall and B. Babington-Smith, Tables of Random Sampling Numbers, Tracts for Computers, no. 2k ( Cambridge , 1939). 5. Rand Corporation, A Million Random Digits with 100,000 Normal Deviates , (Free Press, Glencoe, 1955) • 6. B. Jannson, Random Number Generators , (Victor Pettersons Bokindustri Aktiebolag, Stockholm, 1966), p. 31. 7. D. H. Lehmer, "Mathematical Models in Large-Scale Computing Units," Annals Comp. Laboratory , Harvard Univ . 26 (I95l)> PP« lhl-lk6. 8. A. Rotenberg, "A New Pseudo-Random Number Generator," J. Assoc . Comp . Mach . 7 (i960), pp. 75-77- 9. T. E. Hull and A. R. Dobell, "Mixed Congruential Random Number Generators for Binary Machines," J. Assoc . Comp . Mach . 11 (1964), pp. 31-40. 10. M. D. MacLaren and G. Marsaglia, "Uniform Random Number Generators," J. Assoc . Comp . Mach . 12 (1965), pp. 83-88. 11. Hull and Dobell, J. Assoc . Comp . Mach . 11. 12. B. F. Green, Jr., J. E. Smith, and L. Klem, "Empirical Tests of an Additive Random Number Generator," J. Assoc . Comp . Mach . 6 (1959) PP. 527-536. 13. MacLaren and Marsaglia, J. Assoc. Comp. Mach. 12. THEORETICAL PROPERTIES OF THE MULTIPLICATIVE CONGRUENTIAL METHOD The congruence relation for the multiplicative method, x. .. = ax. (mod m), requires the programmer to define three parameters: the modulus m, the starting value x n , and the multiplier a. If chosen wisely, these parameters will provide speed, a long period and good statistical properties. The following paragraphs describe some rules by which the parameters should be chosen if these goals are to be attained. r The modulus m is usually chosen to be 2 where r is a positive integer less than or equal to the number of bits in a computer word. This choice provides a speed advantage for suppose that ax. = b 2 r+n +...+b 2 r +b n 2 r " 1 +...+b n 2 1 +b^ , b = or 1. l r+n r r-1 1 ' Now x. , = ax. (mod m) means that ax. is to be divided by m and x. n set l+l i x ' l J l+l x equal to the remainder, but if m = 2 then the division is eliminated because reduction modulo m can be achieved by merely setting to zero all but the low order r bits of ax. • l The starting value x n should be chosen relatively prime to 2 the modulus m or the period may be reduced. The proof of this statement requires the following fact from number theory. Let (d,m) = g represent the fact that the greatest common divisor of d and m is g. [Now if (d,m) = 1, that is d and m are relatively prime, and dx = dy(mod m) then x a y(mod m) if (d,m) = g 1 and dx = dy(mod m) then x = y(mod m/g)]. Suppose (x n ,m) = 1 and the period is n. 8 x = ax (mod m) = x p = ax., (mod m) = a x (mod m) • x = a x_(mod m) n Since x is the n+l'st number in the sequence it is equal to x . So a x Q = x (mod m) or a = l(mod m) . Now suppose (x„,m) = g / 1 and the period is n x = a x_(mod m) n a x Q = x Q (mod m) a = l(mod m/g) The period may be shortened because the modulus has been reduced by the factor g. The proof of this statement hinges on this fact: If n is the period modulo m then it is the period modulo m/g, but the reverse is not always true. Suppose a = l(mod m), then a =X • m+1 where X is an integer or a = X.'g»m/g+l so a = l(mod m/g). Now suppose a = l(mod m/g), then a = \,»m/g+l and a = l(mod m) if and only if x/g is an integer. This fact proves that the period when (x n ,m) 4 1 will be the same or less than the period when (x„,m) = 1. If the modulus m is chosen to be 2 , the requirement on x reduces to choosing x such that it is any odd positive integer. The multiplier a should be chosen so that it provides random numbers with the longest period and the best statistical properties. Greenberger showed that the longest period is obtained when a = 8k-3 where k is any positive integer. The following results from number k theory are needed to explain this rule . (l) The Euler Phi function cp(m) is defined as the number of positive integers less than m and relatively prime to m. (2) Euler's theorem states that if (a,m) = 1 then q)(m) is the largest integer such that a ^ ' a l(mod m) . (3) If the least positive integer n such that a = l(mod m) is such that n = cp(m) then a is called a primitive root of m. (k) If n < cp(m) then n is a divisor of cp(m). As was shown in the preceding paragraph the period is the smallest posi- tive integer n such that a = l(mod m) where we will assume that the modulus m is 2 . We would like the smallest n to he as large as possible. ■v* ^» Now if a is a primitive root of 2 then n = cp(2 ) by (3) and if in addition a is odd then this n will he maximum by (2). However, the only moduli which have primitive roots are h, 2»p , and p for p and odd prime so since 2 cannot have a primitive root n < cp(2 ). Now cp(2 ) = 2 by (l) and since a is not a primitive root n must be of the form 2 for s <; r-2 by (k) . The theory requires that a be an odd positive integer and all odd positive integers can be written in the form 8k-l or 8k-3 for k a non-negative integer. We would like to know which form of a will give the maximum period. Suppose that a is 8k-l. The period is 2 for s ^ r-2 and to determine the period we want to find the minimum s satisfying (8k^l) 2 a l(mod 2 r ) or (8k-l) -1 a O(mod 2 r ) This means that (8k-l) -1 is evenly divisible by 2 r . Expanding by the binomial theorem, (tl + 8k) 2S -l = l2 S .8k + 2S ( 2S -jH 8k ) 2 ± ... = (2 s .8k) -1 +(2 S -1) 8k + . . . 2 10 All the terms in square "brackets that include an 8k will 'be even as long as one of the denominators does not cancel out a numerator. The fact that this can never happen is seen by setting n = 2 and rewriting the j+lst term. We have ( (n-l)(n-2) t "- (n-J+l) 1 . (8k) J The part in curly brackets is the binomial coefficient I n-lj and every binomial coefficient is an integer. Therefore every term in square brackets except the first term will be even so the contents of the square r brackets is odd. Since 2 cannot evenly divide an odd number it must be true that 2 S *8k = 0(mod 2 r ) or 2 S+3 'k = 0(mod 2 r ). The smallest s satisfying this relation depends on the choice of k and is such that s ^ r-3« So the maximum period for a multiplier of the form 8k-l is 2 r "3. Now suppose instead that a is 8k-3« We want to find the smallest 2 s s such that (8k-3) = l(mod 2 ). Again expand by the binomial theorem (-3+8k) 2 = 3 2S + 2 S '8k.3 2S-1 + 2 s (2 s -l)(8k) 2 .3 2S - 2 + ... 2! 2 s s =3 + 2 s 8k 1 3 2S - 1 + (2 s -l)(8k) 3 2S " 2 Using the same reasoning as for the 8k-l case, all terms in square brackets except the first will be even so the contents of the square brackets is odd. Then (-3+8k) 2 = 2 S+3 . k -[odd] + 3' 2 S 2 2 s But 3 = (-1+2 ) and expanding this by the binomial theorem 11 3 2S = 1-2 S .2 2 + 2 S (2 S -D(2 2 ) 2 21 1-2 S+2 1- (2 S -1)(2 2 ) . The contents of the square "brackets is again odd. = l-2 S+2 [odd] + 2 S r Now (8k-3) -Is O(mod 2 ) g But (8k-3) 2 - l=-2 S+2 [odd] + 2 S+3 -k.[odd] =2 S+2 (-[odd] + 2k- [odd]) =2 S+2 [odd] r s+2 Since 2 cannot evenly divide an odd term, it must divide the 2 " term. The smallest s satisfying this is s = r-2. Therefore the period for a + r-2 multiplier of the form 8k-3 is 2 . This proves that the longest period for the multiplicative congruential method is obtained -when the multiplier a is of the form 8k-3. Unfortunately, one of the properties of the multiplicative method is that only the most significant bit has the maximum period and the length of period decreases with the significance of the bit. This is ex- r 4- plained in the following way. Let x = ax (mod 2 ). If a = 8k-3 then r-2 the period of the x's is of length 2 . Suppose ax = ...b n 2 r+1 +b 2 r +b n 2 r " 1 +...+b 2 1 + h2° n r+1 r r-1 1 then ax (mod 2 r ) = b 1 2 r ' 1 +...+b n 2 1 +h2° n v ' r-1 1 12 orx , = "b ..2 + y where y = ax (mod 2 ). n+1 r-1 J n J n n v ' r-3 r-2 Now the sequence of y's has the period 2 , hut y = h ? 2 + z r-2\ T-k- where z = ax (mod 2 ) so the sequence of z's has the period 2 r-i-1 Clearly then, "b . has period 2 . Finally, since the multiplier and the starting value are odd, each generated number will be odd and b will always be 1. Although a multiplier of the form 8k-3 will provide random numbers with the longest period, these numbers may not have good statistical properties. Unfortunately, there are no simple rules to ensure the latter and more work needs to be done in this area. From a theoretical standpoint the multiplier should not be close to a simple 5 rational multiple of the modulus m or serial correlation will result. If a = m/k for k an integer then x. , = (m/k)x.(mod m) so x. .. = x./k. i Contrary to the recommendation of some authors the multiplier should not be close to the square root of the modulus m either, because this choice results in strong correlation of lag 2. If a = vm/k for k an integer then x. , = (vm/k)x. (mod m) x. 2 = (vm/k) x.(mod m) = (m/k )x.(mod m) so x. „ = x./k i+2 i' Of course no final selection of the multiplier can be approved until a careful check has been made on the generated numbers. For this purpose a series of appropriate statistical tests is discussed in the next section. 13 FOOTNOTES 1. T. E. Hull and A. R. Dobell, "Random Number Generators/' SIAM Review, k (1962), pp. 230-25^. 2. International Business Machines Corporation, Random Number Genera - tion and Testing, Reference manual C20-8011 (New York, 1959), p. k. 3. M. Greenberger, "Notes on a New Pseudo-Random Number Generator," J. Assoc . Comp. Mach . 8 (1961), pp. 163-166. k. 0. Ore, Number Theory and Its History, (McGraw-Hill, New York, 19^8), pp. 272-302. 5. R. R. Coveyou and R. D. MacPherson, "Fourier Analysis of Uniform Random Number Generators," J. Assoc . Comp . Mach. 1^ (1967), pp. 100-119. 6. International Business Machines Corporation, p. 5« Ik STATISTICAL TESTS The methods discussed so far are designed to produce uniformly- distributed random integers. In addition, a set of uniformly distributed random floating point numbers on the unit interval (0,1) can be produced by interpreting each binary integer as a binary fraction and making the appropriate conversion. Once a supply of uniformly distributed random numbers on the unit interval is available random numbers "with other statistical distributions can he easily obtained. The hypothesis to be tested, then, is that the generated numbers are independent and uniformly distributed. For each statistical test the expected results if this hypothesis is true are compared with the actual results by a Chi-Square test. In order to use the Chi-Square statistic the generated numbers x,,...,x must be separated into distinct classes where the kinds of classes vary with each statistical test. Suppose there are k cells r,,...,r in which the generated numbers can fall and let x.er. represent the fact that x . is contained in cell r. . Let p. be the expected proba- bility that x.er. where j = 1, ...,N$ i = 1, ...k and p. is independent of j. Let f . be the actual frequency of generated numbers in cell i. Then X^ = Z (f . -N«p.)/N«p. has a Chi-Square distribution with k-1 degrees of freedom if the hypothesis of uniform randomness is true. These results can be extended to more than one dimension. Let x.,...,x. .. be an m sequence of generated numbers. As before, let there y j+m-l be k cells r n ,...,r, in which an individual number can fall. There are 1 k 15 then k cells s..,...,s , in -which an m sequence can fall where (x.,...,x. -,)es. = (r ,r , ...,r ) represents the fact that the first v y ' j+m-1' 1 v x y' z member of the sequence fell in cell x, the second member fell in cell y and the id 'th member fell in cell z. Let p. be the expected probability that x.,...,x. -.es.,(j = 1, ...,N-m+l; i = 1, ...,k ) when the hypothesis is true. Let f . be the actual number of m sequences of generated numbers which are in class i. When the m sequences are taken consecutively, i.e., (x 1 ,...,x m ), ( x m+1 > ' - > x 2x) > ( x 2m+l' * * * ^ > " ' then m 2 2 X =Z (f . -[N/m]«p.) /[N/m] «p. has a Chi-Square distribution. However, .111 -i m | J- _1_ _L when the m sequences are taken successively, i.e., (x, ,...,x ), x 2' ' ' ' > x -i ) > ( x q^ • • • )t • • • "then the probability that x . , . . . ,x . _ €s . k m 2 2/ is not independent of j so X = E (f . -[N-m+l].p.)7[N-m+l]#p. does not i=l 1 X 2 have a Chi-Square distribution. For this case , Good shows that X -"* . has asymptotically a Chi-Square distribution with k m -k m ~ 2 degrees of freedom. By defining X = 0, the case for m = 1 holds as before. 2 Once the Chi-Square statistic X has been calculated the next r> oo step is to evaluate / p f(x)dx where J X n/2-1 -x/2 f(x) = ¥-rr+n — is 'the Chi-Square probability distribution for n 2 n/d T(n/2) degrees of freedom. This integral denotes the probability that a 2 Chi-Square random variable will be greater than X , or equivalently, 2 denotes the area under the Chi-Square curve to the right of X . For degrees of freedom n ^ 30 the probability can be found in a table of the 16 Chi -Square distribution. For degrees of freedom n > 30 the quantity v2X is approximately normally distributed about mean ^2n-l with unit 3 variance. A significance level of .05 or .01 is usually selected and if the probability is less than this significance level then the hy- pothesis of uniform randomness is rejected. Various statistical procedures that employ such a test have been suggested for testing the randomness of the generated numbers. In 1938 Kendall and Babington-Smith proposed four widely used tests : the frequency test, serial correlation test, poker test, and gap test. Since then a vast number of other tests have also been proposed. Many authors feel that exhaustive testing of a generator is impractical, that it is better to select a generator that satisfies the more straight- forward tests like the frequency and serial tests. However, it is the opinion of this author that the greater the number of tests that are passed, the more confidence can be placed in the generated numbers. One of the most important statistical tests, the frequency test shows how close the generated numbers are to being uniformly distributed. In this test the unit interval is divided into 100 equal cells and the frequency of numbers in each cell f., i = 1,...,100 is calculated. If the numbers are uniformly distributed then each digit should occur an equal number of times so the probability that a generated number is in any cell is l/lOO. For a sample of N numbers the expected frequency of generated numbers in each cell is N/100. The hypothesis of uniform distribution can be tested by comparing the calculated fre- quencies with the expected frequencies by a Chi-Square test. The sta- 2 10 ° / 2 tistic X , =• (IOO/N) Z (f.-N/lOO) has a Chi-Square distribution with 1 i=l 1 17 99 degrees of freedom. Note that this procedure tests the uniform distribution of the first two digits of every generated number. It is a much stricter test than the usual procedure of dividing the unit interval into 10 subintervals and thereby only testing the frequency of the first digit of every number. A frequency test can also be performed on the octal digits of each random integer. Suppose there are m possible octal digit positions in a random integer and let f . . be the frequency of the i'th octal digit in the j'th digit position, i = 0, ...,7j J = 1, ...,m. If the integers are uniformly distributed then each octal digit should occur an equal number of times so the expected frequency of each octal digit in each digit position is N/8. A Chi-Square test can be used to compare the expected frequencies with the calculated frequencies. The statistic tions j = 1, ...,m has a Chi-Square distribution with 7 degrees of free- X . =(8/n)E (f . .-N/8) calculated for each of the m octal digit posi- .J \ ' -;_ n 1.1 dom. The serial correlation test determines whether any digit tends to be followed by any other digit. This test is defined by dividing the unit interval into 10 equal subintervals and letting f . . be the fre- ij quency of generated numbers in the i'th interval which are followed by a generated number in the j'th interval. Suppose a sample of N generated numbers is tested. If the numbers are truly random then no digit will tend to be followed by any other digit and the frequency of random num- bers in each cell will be the same N/lOO. The expected frequencies can be compared with the calculated frequencies by a Chi-Square test. Let 2 10 9 o I \ 10 , 2 X P =.100/N Z (f, ,-N/lOO) and t n = lO/N'Z (f -N/10) vh ere X 2 , ^ ' i,j=l 1J 1 i,j=l 1 18 corresponds to a frequency test on the random numbers -with 10 divisions 2 of the unit interval. Then X -X has asymptotically a Chi-Square distribution with 90 degrees of freedom. A visual test for serial correlation is derived by noting that the above procedure is equivalent to dividing the unit square into 100 equal cells. If successive pairs of numbers are plotted as points in the unit square then there should be an even scattering of points in all areas of the square. Any preponderance of points in an area represents serial correlation. The mutual independence of the octal digits in each random integer can be checked by a poker test. To perform this test, the first five octal digits of each integer are treated as a poker hand without re- gard to suit. Each hand is tallied as either a bust (abcde), one pair (aabcd), two pair (aabbc), three of a kind (aaabc), full house (aaabb), four of a kind ( aaaab ) } or five of a kind ( aaaaa ) . For instance , if the first five octal digits are 62521 this hand is tallied as one pair. The expected frequencies of poker hands assuming uniform distribution and mutual independence can be computed with the use of the multinomial dis- tribution. The multinomial distribution is defined as follows. If there are n independent trials with k outcomes per trial and 0.,i=l, ...,k is the probability of the i'th outcome and a.,i=l, ...,k is the frequency of the i'th outcome in n trials then the joint distribution of a , ...,a is: J. K. f(a n ,...,a, ) = n! ft. 1 . . . fl k <"W : *i 1 k T~* T^ 1 k 1 k For the poker test n = 5, k = 8, 6. = l/8 i=l, . . . ,8 p(5 of a kind) = 8 y. (1/8) 5 [(l/8)° ••• (l/8)°] (aaaaa) 51(0! • • «0!) 19 p(4 of a kind) (aaaab) p(Full House) : (aaabb) p(3 of a kind) (aaabc) p(2 pair) ( aabbc ) p(l pair) (aabcd) p(bust) (abede) 8-7' 8.7- 8 A 8 ft 5! (l/S^d/S) 1 TUT 51 (l/8) 3 (l/8)< 312! 5\ (l/8) 3 (l/8) 1 (l/8) 1 3!1!1I _5I_ (l/S^Cl/S^d/S) 1 2I2I1. ^ (l/8) 2 (l/8) 1 (l/8) 1 (l/8) 1 2111111! 51 (1/8) 5 1! For a sample of random integers of size N the expected frequencies of the 7 poker hands are N»p. ,i=l, . . . ,7* The actual frequencies f. can be compared with the expected frequencies e. by a Chi-Square test. The 7 ■f.-e. N 2 statistic X = Z ( i i ) has a Chi-Square distribution with 6 degrees i=l e~. l of freedom. In the interpretation of this test a positive correlation in the octal digits 'would result in too many good hands while a negative correlation would result in too many bad hands. A run test shows how the generated numbers oscillate among themselves. To describe this test suppose we replace the sample of N random numbers x , . . . ,x by an N-l sequence S=s , . . . , s where s.=0 if x. x., n . Runs can be ins ide runs or end i i ^ l+l i i l+l runs depending on whether they occur inside or at the two ends of the N-sequence. A subsequence of k zeros bounded by ones at each end forms an inside run up of length k. A subsequence of k ones bounded by zeros at each end forms an inside run down of length k. End runs have to be 20 "bounded only on one side. For example, the sequence 001110 begins with a run up of length 2 and is followed by a run down of length 3* The run test involves counting the number of runs of different lengths and com- paring them by a Chi-Square test with expected results. In the interpre- tation of this test a high frequency of long runs indicates non-randomness of the generated numbers . The expected number of runs of length k is : 2N. k +3k+l -2, k 3 +3k 2 -k-*i- (k+3)! (k+3): This is derived in the following manner. Let r be the number of runs of length k. Let r, . =/ 1 if there is a run of length k starting in the i'th position of the S sequence otherwise Then the expected number of runs of length k is: N N E(r ) = E(Z r ki ). Now E(r^) = 0, i > N-k so E(Z r ) = i=l i=l N-k-1 E(r n , ) + Z E(r, . ) + E(r. , T , ). For 2 < i < N-k-1 then x kl , p v ki' k,N-k' v x E(r, .) = PClO 1 ^) + P(01 k 0) = 2-P(lO k l) = 2-P(x._ 1 > x. x. +k+1 ) 1 Since < x . < 1, then P(x . ) = / dx . J J Jri J P(X. , < X.) a / J dX. _ X 21 r 1 P(x. .. > x.) = / dx. , J" 1 J J x J" 1 j So E(r. , . ) = 2 / " ' f' I ' Xi+k . • • / Xl+1 f ** ± ^3x ± ~ '^i+k-l^i+k^i+k+l i+k+1 x = 2- k 2 +3k+l (k+3)! Nov E(r ) = P(O k l) + P(l k O) = 2-P(0 k l) = 2*2{x ± ^ x 2 ^ . . . ^ x k < x k+1 > x k+2 ) = 2 / / / / . . . / dx_, . . . dx, , dx, dx, , dx. 1 k-1 k k+1 k+2 x, . n „ k+2 = 2. k+1 Xk+2T: But since E(r, , ) = E ( r v w t^ "then N E(r k ) = E(Z_ r k .) = S.E^) + (N-k-2)E(r kj _) = 2N k 2 +3k+l - 2. k 3 +3k 2 -k-+ (k+3)! (k+3)! Another type of run test counts the number of runs above and below the mean. To describe this test suppose we replace the sample of N random numbers x, , . . . t x«, by an N sequence S = s, , . . . ,s. T where s . = 1' 'TI IN l if x. ^ l/2 and s. = 1 if x. > l/2. The number of runs up and down of different lengths are counted as described in the previous paragraphs and compared with expected results by a Chi-Square test. The expected -k-1 number of runs of length k is (N-k+3)2 . This is calculated in the 22 following manner. Let r "be the number of runs of length k. Let r.. =(l if there is a run of length k starting in the i'th position of the S-sequence. otherwise Then the expected number of runs of length k is : N N E(r, ) = E(Z r ). Now E(r . ) = 0, i > N-k+1 so E(Z r, . ) = K 1=1 1CL Kx i=l K± N-k E(r kl } + &>*(**) + E(r k,N-k + l ) ' *or 2 1/2, x.. 1/2) But P(x. ^ 1/2) = 1/2, P(x. > 1/2) = 1/2 J J So E(r k± ) = 2-P(2~ k " 2 ) = 2' k_1 Now ECr^) = F(o\) + P(l k 0) = 2-P(0 k l) = 2»P(x 1 ^ 1/2, ...,x k 1/2) = 2.(2- k_1 ) = 2" k But since ECr^) = E(r k N _ k+1 ) "then N E(r k ) = E(Z rfcL ) = S.ECr^) + (N-k-l)E(r ]d ) = (N-k + 3)2" k_1 A test which is similar to the run test is the gap test. If a maximum in a sequence of generated numbers is a number that is surrounded by two smaller numbers and a minimum is a number surrounded "by two larger numbers then a gap is defined as the distance between two successive 23 maxima or minima. In determining the length of a gap the two maxima or minima are included as -well as the numbers "between them so the smallest gap is of length 3. The gap test is performed by counting the number of gaps of different lengths and comparing them with expected results by a Chi-Square test. For a sample of generated numbers of size N the expected number of gaps of length r is : r-1 3«N-2 - r-2 . This formula is derived very clearly by Kermack- rl r+2 McKendrick and it will not be discussed here. 7 MacLaren and Marsaglia were dissatisfied with the afore- mentioned statistical tests because they thought the tests were irrelevant to actual uses of random numbers. One of the most common ways of using a sequence of random numbers is to sample consecutive n-tuples and so MacLaren and Marsaglia devised frequency tests on the range of n-tuples and on the range of maximums and minimums of n-tuples. If ¥. = max(x. _,..., x. ) i = 0, . ...N-n where the x. . are independent and 1 v l+l' ' i+n 7 ' ' l's uniformly distributed on the unit interval (0,1) then W. is not uniformly distributed because : P(W. < a) = P(x. . < a)'-- P(x. < a) v 1 ' v l+l ' v 1+n ' pa pa = / dx. . . . . / dx. J i+l J Q i+n n = a , whereas uniformity would require P(W. < a) = a. However, P(W. < a) = P(W. < a ) = a so ¥. is theoretically uniformly dis- 11 l tributed on (0,1). In a sample of N generated numbers there are N/n n-tuples and the observed W. , i = 1, ...,N/n can be tested for a uniform distribution by the frequency test as described earlier. In a 2k similar manner the minimum of n-tuples can be tested. In this case, if W. = min(x. .,,..., x. ) then 1 x i+l ' i+n' P(W. < a) = 1-P(W. > a) =1- [p(x. +1 >a) P(x i+n > aij = !■ dx a i+l J ■1 dx. a l+n ,n n . is = l-(l-a) Now P(W j , < a) = P l-fl-W^ 11 < l-(l-a) n = l-(l-a) n so l-(l-W ± ) theoretically uniformly distributed. A frequency test can again be employed to determine how close the actual results are to being uniformly distributed. Unfortunately, a random number generator which has passed all of these standard statistical tests may still fail when a particular application is sensitive to a property of the numbers that was not tested. In the end, the most relevant test of a random number generator for a particular application is to use the generator on a similar problem for which the answer is known. 25 FOOTNOTES 1. B. Jannson, pp. 170-191. 2. I. J. Good, "The Serial Test for Sampling Numbers and Other Tests y u for Randomness," Proc . Camb . Phil . Soc . % (1953), pp. 276-28^. 3. M. G. Kendall and A. Stuart, The Advanced Theory of Statistics , (Griffin, London, 1958), I, p.~~500. ,. x 4. M. G. Kendall and B. Babingt on- Smith, Randomness and Random Sampling Numbers," J. Roy. Stat . Soc . 101 (1938), pp. IV7-I66. 5. H. Levene and J. Wolfowitz, "The Covariance Matrix of Runs Up and ^ l Down," Annals Math . Stat . 15 (19^), pp. 58-69. 6. ¥. 0. Kermack and A. G. McKendrick, "Some Distributions Associated with a Randomly Arranged Set of Numbers," Proc . Roy . Soc . Edinburgh 57 (1937), PP. 332-376. 7. MacLaren and Marsaglia, J. Assoc . Comp. Mach . 12. 26 EXPERIMENTAL RESULTS The multiplicative method of random number generation has been generally accepted as a good one because its theoretical properties are known, it is fast and easy to program, it has a long period, and> most 1 2 important, it produces numbers with good statistical properties. ' Unfortunately, published results of the application of the multiplicative method so far have been for computers with word lengths of 36 bits or larger. There is some speculation whether the multiplicative method will be an adequate method of random number generation for the IBM 360. A prob- lem inherent within the method concerns the selection of a multiplier which will produce numbers with good statistical properties. From a theoretical standpoint it is known that the multiplier 1) should be of the form 8k-3 for k a non-negative integer in order to ensure the maximum length of period 2) should not be close to a simple rational multiple of the modulus or serial correlation will result 3) should not be close to a simple rational multiple of the square root of the modulus or correlation of lag 2 will result but, aside from these rules, little is known about why some multipliers produce numbers with better statistical properties than others. This prob- lem is magnified by the 3^0 's small word length of 32 bits and it is believed that random number generation on the 3^0 by the multiplicative method will be extremely sensitive to choice of a multiplier. In an attempt to learn more about random number generation on the IBM 360 the multiplicative method was programmed for the 360. A total 27 of 1,500 different multipliers was generated randomly and used to generate sequences of pseudorandom numbers. Each sequence was statis- tically tested and the test results were compared. Integer arithmetic was used to generate fixed point numbers and floating point numbers on the unit interval were obtained by setting to zero the first eight bits of each integer number and inserting the appropriate exponent. In order that the bit patterns for the integer random numbers and floating point random numbers be identical for the statistical tests the integers were stored as twenty four bit results with the first eight bits of each 2k word set to zero. The modulus of the method thus was reduced to 2 and since the starting value used was odd and all multipliers were of the form 8k-3 for k a non-negative integer the period of all the generators 22 tested was the maximum possible of 2 . To test the 1,500 random multipliers the decimal number 2173 was arbitrarily selected as the starting value and for each multiplier a sequence of 5>000 or 10,000 numbers was generated. The frequency test and serial correlation test were performed on each sequence of generated numbers and for both tests a Chi-Square statistic and the associated Chi-Square probability (i.e., the probability that if the numbers are random a Chi-Square variable will exceed this test statistic) were calculated. Of the 1,500 multipliers tested the ten which had the largest Chi-Square probabilities were selected for further testing. These ten multipliers were compared in three different ways with the nine multi- pliers which gave the lowest Chi-Square probabilities and for which the hypothesis of uniform randomness was rejected. First, since a multiplier 2k could be as large as 2 -1 it could have as many as eight decimal digits and it was thought that larger decimal numbers were probably better 28 multipliers than smaller ones. However, some of the five -digit numbers were "better multipliers than the eight -digit numbers so this idea had 3 to be rejected. Coveyou stated that the multiplier should not have a small number of l's and it was conjectured by this author that the "best multipliers probably had an equal number of O's and l's in the binary representation. After tallying O's and l's it appeared that this idea also had to be rejected. It was then suggested that "by looking at the octal digits of the multipliers it might be observed that certain octal digits were present or were lacking in particular digit positions of the ten acceptable multipliers. This hypothesis again proved false and so the search for finding rules to predict when a multiplier would be a good one was concluded. Next it was desired to learn more about the kinds of generators that were produced by using the ten selected multipliers as defined in the previous paragraph. Passing tests like the frequency and serial correlation tests is necessary but certainly not sufficient for the acceptance of a generator so each multiplier was subjected to an exten- sive set of tests. In all, sixteen tests were performed on the first 10,000 numbers generated with starting value 2173* The tests are described as follows : TEST 1 - Frequency test on the generated numbers with 100 divisions of the unit interval. TEST 2 - Serial correlation test. Successive 2-tuples of numbers were taken as points in the unit square. The unit square was divided into 100 equal cells and the frequency of 9>999 pairs was tallied. TEST 3 - Successive 3 -tuples of numbers were taken as points in the unit 29 TEST 7 eube. The unit cube was divided into 1,000 equal cells and the frequency of 9>998 triples was tallied. TEST k - Consecutive, non-overlapping, 2-tuples of numbers were taken as points in the unit square. The unit square was divided into 100 equal cells and the frequency of 5>000 pairs was tallied. \V TEST 5 - Consecutive, non-overlapping, 3-tuples of numbers were taken as points in the unit cube. The unit cube was divided into 1,000 equal cells and the frequency of 3>333 triples was tallied. TEST 6 - Consecutive, non- overlapping, 4-tuples of numbers were taken as points in the unit cube in four space. The unit cube in four space was divided into 10,000 equal cells and the fre- quency of 2,500 4-tuples was tallied. Frequency test on the eight octal digit positions of each generated integer number. Poker test. Test of runs up and down. Test of runs above and below the mean. The maximum number in consecutive, non- overlapping, 2-tuples of numbers was calculated. The unit interval was divided into 100 equal cells and the frequency of the 5>000 maximums was tallied. TEST 12 - The maximum number in consecutive, non -overlapping, 3-tuples of numbers was calculated. The unit interval was divided into 100 equal cells and the frequency of the 3>333 maximums was tallied. TEST 13 - The maximum number in consecutive, non- overlapping, ^-tuples of numbers was calculated. The unit interval was divided into 100 equal cells; frequency of the 2,500 minimums was tallied. < ] ,vV> "; TEST 8 TEST 10 TEST 11 30 For each test the calculated results for the generated sequence of numbers were compared with results that would have "been expected from a truly random sequence of numbers "by a Chi-Square test. The results of the Chi-Square test are given in TABLE 1. Each column of the table represents a particular multiplier, shown in decimal at the head of the column, and each row represents a particular test. The numbers listed represent the Chi-Square probability or percentage probability of exceed- ing the calculated Chi-Square statistic. The appearance of e in the table means that the Chi-Square probability was less than .05 and that the hypothesis of uniform randomness was rejected. The eight probabilities listed for TEST 7 are the results of the frequency test performed separately on each of the eight octal digit positions of the 2k bit random integers. The first number is for bits 22-2^ (counting beginning from the right end of the word), the second number bits 19-21, and the eighth number bits 1-3 • It is obvious that the Chi-Square probability for the eighth digit position, bits 1-3, will be extremely low since the least significant bit must always be a 1 because the generated numbers are odd. The two tests on the digits, tests 7 and 8, are not as important as the other tests because most applications require numbers rather than digits. In this regard, the manner in which TABLE 1 is viewed depends largely on the particular application for which a generator is needed. If an application requires randomness of consecutive 2 -tuples then the results of tests k, 11, and 1^ indicate that multiplier 9 should not be chosen. If, instead, good randomness of consecutive 3-tuples is needed then tests 5, 12, and 15 indicate that any of the multipliers is adequate. Tests 6, 13, and l6 show that for randomness of consecutive ^-tuples multiplier 10 31 TEST 1 TEST 2 TEST 3 TEST k TEST 5 TEST 6 TEST 7 TEST 8 TEST 9 TEST 10 TEST 11 TEST 12 TEST 13 TEST Ik TEST 15 TEST 16 TABLE 1 -1- 13651723 -2- 12316023 -3- 50515 -k- 29355 -5- 5^27 .9* .<* .98 • 97 • 99 .90 .92 .98 •97 .89 .93 • 57 • 98 Al .76 -99 .68 .87 • 99 .57 • 79 • 55 .kQ • 53 .k2 .80 .27 M • 33 .83 • 95 .1£ • 51 • 99 • 99 • 99 • 99 .85 .76 .16 •97 • 99 •99 •99 .76 • 55 .k2 • 99 •99 • 99 -99 .81 • 77 • 98 •99 • 99 -99 -99 .69 .50 • 58 -99 -99 • 99 • 99 6 G e e € .08 .68 .07 -19 .27 e .07 •35 .56 .11 .26 .ko •93 .76 .90 •97 •96 •97 .83 AO .9^ • 57 .62 • 75 -95 .65 .86 .96 .29 .18 • 53 • 05 .98 • 37 .92 •27 • 30 -9^> .76 • 99 .90 .13 .88 .72 .7* 32 TABLE 1 (Continued) TEST 1 TEST 2 TEST 3 TEST k TEST 5 TEST 6 TEST 7 TEST 8 TEST 9 TEST 10 TEST 11 TEST 12 TEST 13 TEST Ik TEST 15 TEST 16 -6- 223227 -7- 220915 -8- 220615 -9- 3739387 -10- 1^59251 • 95 • 98 .<* .96 •96 • 93 .89 • 93 .88 .98 • 52 • 39 .88 .81 .60 .kk .89 • 51 .72 •91 .kd .66 • 93 .92 .3^ .76 .90 • 97 .62 •95 • 59 • 55 .65 • 99 • 99 • 99 • 99 .83 • 55 .89 • 93 • 99 • 99 • 99 .72 • 3^ .66 • 99 • 99 • 99 -99 .84 .63 .09 • 99 • 99 • 99 • 99 .5^ •93 • 38 • 99 -99 -99 -99 e € £ € e •15 .kk Al G .92 .85 •32 .8lt- e .20 .68 .Ik • 27 .28 •89 • 52 .71 • 70 e .ko .60 .kz -9^ .22 • 57 .54 .k6 .66 .06 .kl • 19 • 98 .06 .85 •15 • 70 .84 .73 .kd .21 .1+8 .70 .35 .17 e 33 should not "be chosen. Some applications require uniform numbers that oscillate randomly among themselves. This property of the generated num- bers is tested for by tests 9 and 10 and the results of these tests show that multipliers 1 and 9 should not be used for this application. Although the ten carefully chosen multipliers of TABLE 1 gave good results it must be remembered that the 10,000 numbers tested in each case comprise a very small sample of the ^,19^,30^- possible numbers in the period of the generator. Consequently, further tests were performed using different sample sizes and different starting values. TABLE 2 shows the effect of a change in starting value. Each row repre- sents a particular decimal starting value shown at the left of each row and each column represents a particular decimal multiplier shown at the head of the column. The numbers listed in the table are the Chi-Square probabilities for the frequency test performed on a sample of 10,000 numbers generated with each particular starting value . TABLE 3 shows the effect of a change in sample size. Each row represents a particular size sample shown at the left of the row and each column represents a particu- lar multiplier shown in decimal at the head of the column. The numbers listed in the table are the Chi-Square probabilities for the frequency test performed on particular size samples generated with a starting value of 2173. As the statistics in TABLES 2 and 3 show, when one repeats the same test on different samples one should get a range of Chi-Square values with the corresponding Chi-Square probabilities distributed uni- formly. Since the interpretation of the statistics does not vary drastically with a change in sample size or starting value there is no reason to believe that the original tests on samples of 10,000 numbers with starting value 2173 were biased. OJ g H co a> -d- OJ H CO -=t" ON 1 LP* OJ LPv J=f -* VO ON H O OJ • • • • • • H On 1 ltn -4" H t— CO ON t- o\ ON w OJ -d- 1 CO CO CO ON LTN LTN co 00 ON CO • • • • • • • i ON co t- - co H o W w H OJ -eh co H LTN co LTN ON ir\ H VO 1 VO • • • • • • CO O 1 OJ OJ LTN H H r£ t— LTN H H CO H CO VO VO CO LTN VO VO c- i ON • • • • • • • f-O I OJ OJ t— co H VO ON ON VO CO -d; ,A °J ON OJ LTN OJ -=f co [- CO VO OJ • • • • • • • i en OJ OJ r- Lf\ OJ OJ CO OJ LTN (D ON 1 OJ ON LTN CO 3F VO H H lf\Jf • • • i 3- LT\ LT\ VO OJ -=*■ ON CO 3 co ON 1 LTN -=h LTN VO CO OJ ON C— -d" CO • • • • • 1 ON OJ lf\ o\ CO a -d- vo LTN w ITN i H VO CO t— VO -=f [*- CO LTN • • • • • 1 O LT\ CO OJ o co OJ Ol ON co VO .9 -4- VO VO ON LTN CO • o OJVQ I H CO -4- co ?1 co VO ON u; C— CO OJ ON J* £vi t- CO ltn c— ON CO CO r- ■ • • • • • i H H LTN 1 VO CO ri H H OJ H CM CO OJ CO J- H LTN VO H & -4- LTN VO p- H LTN VO t- CO oo H t- VO t— CO o LTv UA D— 1 UA ON ON CO c*- H CJ H O 3& • • * • • • • * 1 UA -* H t- VO VO OO CO UA ON ON CO co ON ON CO ON ON ON ON CO 1 oo • • • • • • • • ON f- 1 OO t— oo H c- -=h LT\ UA c— VO D— -d" ua ON ON CO c— VO VO _=}- ON 1 VO • • • • • • • GO O i OJ CM UA oo CO LTN c— VO H if\ O H CO ON ON CO ON ON ON CO I ON • • • • • • • • t>- O I OJ OJ t- -=»■ LTN LTN OO -4- o o t— OJ 3- ON CO UA UA -3- oo OO 1 OJ • • • • • • • • VQ OO 1 OJ OJ t- o ON VO oo c— VO t- d 1 OJ ON ON c— VO OO OJ H UA,d- • • • • • • • 1 -4- ua UA OJ c— ON -d- oo o UA -d- 1 UA VO ON t- H H OJ UA OJ ^t OO • • • • • 1 ON OJ UA ON CO t- OJ ON jdr OO VO i H ON ON VO -4- CVJ VO t- CO OO LTN * • • • • • i o LTN OO ON -=h ON OJ D— r— oo -* ■ SI -d- ON t— UA t- UA CO CO • • • • • • • • OJ vo 1 H OO SI OO H -3- OO J" OJ OO UA s SM ON ON ON 00 r- VO f- I c— • • • • • • H H 1 LTN vo OO H UA OO O O O o o o o o o O O o o o o o o O o o o o o o o UA o UA o UA o UA H H OJ OJ oo OO -d- 36 This "battery of tests was next used to learn more about the IBM random number generator for the 360, RANDU. There had "been some hesitation to use RANDU because its statistical properties were not known. The same sixteen tests were applied to the first 10,000 numbers generated by RANDU with starting value 2173 and the Chi-Square probabilities are listed in TABLE 4 under multiplier 11. The results were theoretically predictable. The generator RANDU uses the multiplicative method with the 31 modulus 2 and the multiplier 65539* Note, however, that 65539 is l6 2 +3* IBM made a mistake in choosing the multiplier so close to the square root of the modulus and this accounts for the extremely high correlation among 3-tuples as shown by the results of TEST 3» One of the most recent developments in the random number genera- 4 tion field is a paper by Marsaglia in which he reports that all multi- plicative generators have a defect in that when successive n-tuples of numbers (x_ , . . . ,x ), (x p , . . . ,x ), . . . are plotted in the unit n-cube all the points will fall into a small number of parallel hyperplanes. This defect is not too serious since most users will not select random numbers in this manner anyway. If an application does require randomness of successive n-tuples, however, some improved methods of generation have been suggested. One of these methods, suggested earlier by MacLaren and Marsaglia, involves the mixing of two separate multiplicative genera- tors. This method was programmed for the 360 and produced random numbers x.,...,x in the following manner. A list of 128 numbers u_, ...,11.-0 vas 24 produced by the relation u. 1 = a u. (mod 2 ). Another relation 24 v. = a v. (mod 2 ) was used to generate a second set of numbers v, ,...,v . Bits 18-24 of vl were used to select the k'th random number, 37 TEST 1 TEST 2 TEST 3 TEST k TEST 5 TEST 6 TEST 7 TEST 8 TEST 9 TEST 10 TEST 11 TEST 12 TEST 13 TEST Ik TEST 15 TEST 16 TABLE 4 /2#i>J>4 MWi -11- 65539 -12- 220651 13651723 .6? .9k .38 .58 \ .86 .85 • 90 • 31 .k6 .07 • 91 • 32 • 93 .20 -99 ■ 99 • 99 • 99 .96 .61 • 51 .99 -99 -99 -99 • Ik .09 .21 .20 .26 .50 .64 • 70 • 71 •97 .81 .21 • 97 •3^ • 71 .6k • 37 .68 38 x , from the list of the u's and this position was then refilled with v . The two multipliers selected for a and a were multipliers 8 and 1 respectively. The sixteen tests were performed on the first 10,000 numbers generated with starting value u = 2173 and v = 2173 and the Chi-Square probabilities are listed in TABLE k- under the heading multi- plier 12. Though this method eliminates the inherent non- randomness of successive n-tuples in the multiplicative method it requires 128 storage locations and is slower than the multiplicative method. The time for generation of one random number for multipliers 1-10 and RANDU was .0022 hundredths of a second, but the generation time for the Marsaglia mixing process was .0031 hundredths of a second. With the speed and storage capacities of modern-day computers, however, these two problems should not be serious for most users . In conclusion, a battery of programs was written to provide stringent statistical tests for a random number generator. A group of random number generators was devised, each one meeting the needs of a different type of application. It is hoped that the user will consult TABLES 1-k and then select the generator which best satisfies his requirements . In the likely event of a further breakthrough in the field of random number generation the battery of statistical programs will undoubtedly be of assistance in testing future random number generators. 39 FOOTNOTES 1. Coveyou and MacPherson, J. Assoc . Comp . Mach. Ik. 2. Hull and Dobell, J. Assoc . Comp . Mach . 11. 3« Coveyou and MacPherson, lk. k. Marsaglia, Proc. N. A. S. 6l. d ****•»* ■ 5. MacLaren and Marsaglia, J. Assoc. Comp. Mach. 12. ko BIBLIOGRAPHY 1. Coveyou, R. R. and MacPherson, R. D., "Fourier Analysis of Uniform Random Number Generators," J. Assoc . Comp . Mach . Ik (1967), pp. 100-119. 2. Good, I. J., "The Serial Test for Sampling Numbers and Other Tests for Randomness," Proc . Camb . Phil . Soc . k-9 (1953), PP- 276-284. 3. Green, Jr., B. F., Smith, J. E., KLem, L., "Empirical Tests of an Additive Random Number Generator," J. Assoc . Comp . Mach . 6 (1959), PP. 527-536. k. Greenberger, M. "Notes on a New Pseudo-Random Number Generator," J. Assoc . Comp . Mach . 8 (1961), pp. 163-166. 5- Hull, T. E. and Dobell, A. R., "Random Number Generators," SIAM Review h (1962), pp. 230-25^. 6. Hull, T. E. and Dobell, A. R., "Mixed Congruent ial Random Number Generators for Binary Machines," J. Assoc . Comp . Mach . 11 (1964) pp. 31-lK). 7. International Business Machines Corporation, Random Number Generation and Testing, Reference manual C20-8011 (New York, 1959) . 8. Jannson, B., Random Number Generators, (Victor Pettersons Bokindustri Aktiebolag, Stockholm, 1966). 9- Kendall, M. G. and Babington-Smith, B., "Randomness and Random Sampling Numbers," J. Roy . Stat . Soc . 101 (1938), pp. lVf-l66. 10. Kendall, M. G. and Babington-Smith, B., Tables of Random Sampling Numbers , Tracts for Computers, no. 2k (Cambridge, 1939)' 11. Kendall, M. G., and Stuart, A., The Advanced Theory of Statistics, (Griffin, London, 19 58), I. 12. Kermack, W. 0. and McKendrick, A. G., "Some Distributions Associated with a Randomly Arranged Set of Numbers," Proc . Roy . Soc . Edinburgh 57 (1937), PP- 332-376. 13- Lehmer, D. H., "Mathematical Models in Large-Scale Computing Units," Annals Comp . Laboratory , Harvard Univ . 26 (1951), pp. lkl-lk6. Ik. Levene H. and Wolfowitz, J., "The Covariance Matrix of Runs Up and Down," Annals Math. Stat . 15 (1 N. RN - The one dimensional array in which the N floating point random numbers are stored by RANB. Its contents may be arbitrary upon entry. IRN - The one dimensional array in which the N integer random numbers are stored by RANB. Its contents may be arbitrary upon entry. Note : When a particular multiplier is desired, the card labelled MULT in the program should be changed to the correct multiplier. kk W AMI IN 7>ku C H /\ K 7 inn M|l|. | "I I— ! P L L L L I L l> LP I. S« i\i SI 1 1 SI l> AJ- STh L. ST L LA Hf. I SI Ll-flVH DC UC i»C DC US t- iv I ) MM. I K f 7) *.M 1 ) !i .0(3) ft,H( 1 ) M , ] ? ( 1 ) ?,()( ] ) 5.0(?) 10,1 <>,H I 1 , V, V «t. 1 1 h, /i ill I 4 . ( . I U\ K !j, I H-'iH I ) . I h i" P , / t- K I I ( > , 'I I- l* k •j, I PnP ■> , i ) ( ft , 7 ) !j • i)I * » 7 J U, I. ' II IH U),0(?) h'O.'H X • Pn()()l X »00|-H-H>h • X * » )f^( 1 -\ •-> f- »- -i * 1 h KHTS AOUKPSS HI- N PUTS N IN KM,, -i (it- IS AuUKhSS ul- WWII) In Pbf;. ft (;HTS ADDPPSS UP IKim(1) IM RPG. « (iPTS AUDPPSS (it- SlAKllMI, VALUP PUTS STARTING, WALHH In hK;, S SAVhS 1HI- STARTING VALUt- 7HKMS lllll THt- t-jKSl H HITS HI- THH I'MTHGhK SI i IK I- S TMh iNfhKPR Ii\i Ihi- AKKAY IKN CIinVPR'I I li PLUA'f I>\'fi H C 1 1 ivi I f;«- I S IHf- KAiVUtiH iMUmMfcR I'M A hLTINU PI RK; NliKf- AL 1 Zl-S THt- KAMDilW MHKiHhK GETS THfc NWF'tritn* rtA(.K I'M GbNL. KPG. b PUTS KAlMDOivi IMUmHpK Jim I HP AKKAY K IM VII IS THp I'M'IbOK I l\i KH;. 4 KhSlUWhS THt- SIAKlI'v(; \/A|.U(- ChAkACTHK IS1 1L ^5 SUBROUTINE *** UNINUM IDENTIFICATION *** Frequency test for random numbers PURPOSE "**■* Determines how close a set of random numbers are to being uniformly distributed by dividing the unit interval into 100 subinter- vals, counting the frequency of random numbers in each subinterval and comparing with expected results by a Chi-Square test. OTHER SUBROUTINES USED *** PROLIM USAGE *** REAlM RN (NX) CALL UNINUM (N,RN,PR0B) Where the parameters in the calling sequence are: N - The number of random numbers to be tested, NX ^. N. RN - The one dimensional array containing the N random numbers to be tested. PROB - The outputed value of the Chi-Square probability. It is a number between and 1 which indicates the probability that a Chi-Square variable will exceed the calculated test statistic. k6 S l • •-» W I H II IWl- NHlPVlIN, ( |\. t KM, PWIIH ) IhSI '»»- O'vlHIKM'fY hF NUmH^RS WITH 100 DlVlSIIMS UF 1 HK UNIT INTERVAL H(=4L*4 HM(w| f HIV|lfil) ?iM| R;^k*<. (>LLI 100) niv(i )«o.() dii ^ T=2.mi oi w( I )=mv( T-i ) + .oi on h l=l,) oo Cf-LL(J)=n III I H [ = 1 , I>1 DII 7 .1=1,100 1 1- (w>'( I ) .<;i-.i>]\/(.|) .ANI),KN( I ) .LT.DIVJ J + l ) )f;n TH « L.'jNl IMIh CHLL ( J ) =CHLL ( .1 ) + 1 |-XPMlS = M/]00. XrtJllHsM.O mi w i = i , loo XW(ji»isXNt)M+(CHLL ( I )-t-XM0IS)**2 iiim(. S"=xi\i(ii"i/eXh , ii! s CaLL i j kiiL !""< ww, I hvf.Su, PWIIR) MHTUKim (-Ml) SUBROUTINE NAME *** FRESOC IDENTIFICATION *** Frequency test on successive or consecutive n-tuples of random numbers . PURPOSE *** Determines how close a set of n-tuples of random numbers are to being uniformly distributed. The n-tuples can be taken either successively or consecutively and n must be 1, 2, 3> or 4. The unit 2 (interval, square, cube or cube in 4-space) is divided into (10, 10 , 3 h- 10 , or 10 ) subintervals respectively depending on the value of n. The frequency of random numbers in each sub interval is tallied and compared "with expected results by a Chi-Square test. USAGE *** REAlM RN(NX) GALL FRESOC (%RN?NTUPLE?=IS0C=CHISQ^PR0B) Where the parameters in the calling sequence are : N - The number of random numbers to be used in the test, NX ^N. RN - The one dimensional array containing the N random numbers to be tested. NTUPLE - The size of the tuple, NTUPLE = 1,2,3,4. ISOC - If ISOC = 1 then the N random numbers are to be divided into N/NTUPLE-eonsecutive, non- overlapping, tuples of numbers. CHISQ - On exit from FRESOC CHISQ will contain the calculated Chi- Square statistic for the test. If ISOC = 1 the input value for CHISQ must be the CHISQ value obtained from a previous call to FRESOC with the NTUPLE value one less than the NTUPLE value for this call. If the NTUPLE value for this call is 1 then the input value for CHISQ is 0. If ISOC = 2 the input value for CHISQ may be arbitrary. PROB - The outputed value of the Chi Square probability. It is a number between and 1 which indicates the probability that a Chi -Square variable will exceed the calculated test statistic CHISQ. 1*9 7 H 9 10 20 ?S 30 40 50 53 SOBkOOl I .Mi- |-wt-SOC(N,RN,l\!TOPLF, I SoC,OK'.Su, PwilB ) FRFOUHNCY I>S1 UN SUCCESSIVE OR CONSFCUTIVF N TUPLES KFftL*4 HN(W) , UI V( 11 ) Iivjl FGF«*4 CFLL ( 10, 10,10,10) ,LPM{4) l)TV( 1 )=o.o l)D 5 1 = /, 1 1 I ) T w ( T ) = i ) 1 v I I - 1 ) + . 1 J = l _ K=l L = l mi y i=i,io no r j = i , l o 1)11 7 K=l , 10 Dll ft L=l,10 CFLLI I ,.), K,l. )=0 IFJNTUPLE.FO. 1 )GU TO 9 IF(NTUPLF.FO.2)G0 fo fi IF(n"IUPLE.F0.3)GU TO 7 COM IIiMtIF COMT IiMllfc CON "I IlMllh CUM T IIMOH 00 10 1=1,4 LPM( I )=l IF( ISOC.fcu.2 )(,n" TO To" INC«=1 1 k M A X = m - im IUPLb+I Gil TO ?>i MOi»iTOP = lM/ivTOPLt:" IKMAX = imiimiip#nT0PLF-imTUPLF+_1 INCR=NTUPLE DO 50 I«=l , IkMAX, INCR 00 40 MT=1 . fivTltPLF IRP= IH+MJ-I DO 30 J=l , 10 IRKM IRP) ,GE.l>TV( J) ,AMD.RN( IRP) .LT.DI V( J+l ) )G(.I TO 40 CONTINUE „LPM(NT)=J CFLULPMI 1 ) ,LPM( 2) ,LPM( 3) , LPM(4) )=CtLL(LPM(l) ,LPM(2) ,LPM( 3) ,LPM(4) 1 ) + l 1 F ( IS0C.bi-).2)G0 TU~~bi~ FXPDIS=lRMAX/( 10.**IMTUPLE) GU TO 54 FXPDIS=M0mT0P/ ( 10.*#NTUPLE) 50 s^ S7 NUI'I = 1 = 1 = 1 n *6 II b n 5 N 1 1 M F ( M F ( N F ( M [}i\iT = 0.0 1=1,10 4 -1 = 1,10 H K = l , 10 7 L=l,10 =XNUM+(CELL( It J»K,L)-EX PETS') "I UPLE.EQ. 1 )0(i Til 60 riJPLE.E0.2)r,f) TU b9 rUPLE.E0.3)Gf] TM 68 1 M I I h **2 5H COMTIMUE ^V CUNT INI J E 60 CHN'I I IMliE IE( IS!)C.EQ.2)GD In yo TEwiCSO=XMUM/EXPI)I S MmCS"= I thCSO-iiiMC S'i i-i i ) I- = 1 * * N T U P L E- 1 0** ( N f UPL E- 1 GO T(l HO HMCSn=XNl.iM/EXPI)TS Mf)F=] 0**n riJPLfc CALL ^41 iL J M ( NDH , UNCStl , PROFS ) R F T UR M FlMI'l 70 HO SUBROUTINE NAME *** SERPLT 51 IDENTIFICATION *** Plotter routine for random number generator PURPOSE *** Successive 2 -tuples of random numbers are plotted as points in the unit square . Any preponderance of points in an area of the plot indicates serial correlation. Note that the numbers printed along each axis should be divided by 10 to obtain their actual values . OTHER SUBROUTINES USED *** CCP1PL,CCP5AX USAGE *** REAlM RN(NX) GALL SERPLT (N,RN) Where the parameters in the calling sequence are : N - The number of random numbers to be used in the test, NX ^. N. RN - The one dimensional array containing the N random numbers to be tested. C". 40 50 C 52 S lift ROUTINE SERPLT ( N,RN) REAL*A RN(N) ? T(2) .EDGE [mTEGER*4 UP, DOWN, LAST ( 16,16) LGRID=16 _. _ C THINK BASK LGRID C LGRIO MUST.BE _POWER OF TWO_ IMJ 10 1=1, LGRID DO 10 J = l, LGRIO 10 LAST(T,J)=0 Mnls=N-i C PUT POINTS IN BOXES on so 1=1, NM1 LAB=LGR Il)*RN~( I ) + 1 LORD=LGRID*RN( 1+1 )+l LPRFV=LAST< LAB,LOWU) fh(LPKEV.E0.0)G0 TO 40 LDIFF=I-LPREV R N ( I ) = R N ( I ) + L D I F F STORED BACK POINTER IN RN( I j LAS f ( LAB, LORD) =1 CI ImT INUE HAVF POINTS IN BOXES, NOW DECODE AND PLOT. I ( ] )=0.0 T(2)=1.2b CALL CCPlPL(0.b,0.t>,-3) CALL CCPbAX(0.0,0.0f 21HSEC0ND MEMBER OF P A I R , 2 1 , 8 . , 90 . , T ) C/U.L CCPbAX ( 0.0,0.0, 20HFIRST MEMBER OF PA I R , -20 , B . 0, . , "I ) 0P = 3 i ii IWN=2 FDGF=0.01 DO 100 LAB=1,LGRID OH 99 L0RD=1,LGR ID I=LAST (LAB, LORD) 60 1 F ( 1 .Hi.O )G0 TM 99 JX=RM( I ) JY = RN'( 1 + 1 ) RN( 1 )=RN{ I )-JX BY=RN( 1+1 )-JY XC=B.O*RN( I ) YC=B.O*FY XL=XC-EDGE Xm=XC+EDGF YL = YC-EDGF Y*«=YC + EDGE CALL CCP1PL ( XL, YL, DP ) (.ALL CCP1 PL ( XM, YW» DOWN) CALL CCP1PL ( XM, YL,OP ) CALL CCP1PL ( XL , Ym f DOWN) In! PP=JX IF( I DIFF. h(i.O )GI I H> 99 I = I-|J)IFF GO id 60 99 l".i i. • I T •Of- 1 00 CI I I 1 1 < I IE RETURN f (v n 53 SUBROUTINE NAME *** POKDIG IDENTIFICATION *** Poker test and frequency test on the octal digits of a set of random integers. PURPOSE *** The 2k bit random integers are divided into 8 octal digits each. In the poker test the first 5 octal digits of each integer are treated as a poker hand without regard to suit. The frequencies of certain poker hands are tallied and compared with expected results by a Chi-Square test. In the frequency test the frequency of each octal digit in each of the 8 digit positions is tallied and compared with expected results by a Chi-Square test. OTHER SUBROUTINES USED *** FREPOK, PROLIM USAGE *** REAlM PR0Bl(8) INTEGER *K IRN (NX) • CALL POKDIG(N, IRN, IFOP, PR0B1,PR0B2) Where the parameters in the calling sequence are : N - The number of random integers to be used in the test, NX > N. IRN - The one dimensional array containing the N random integers to be tested. IFOP - If IFOP-1 then the frequency test is to be performed. If IFOP = 2 then the poker test is to be performed. PR0B1 - If IFOP = 1 then PROBl(l), I = 1, ...,8 contains the Chi-Square probability of the frequency test for the i'th digit position where the 1st digit position is bits 22-24 (counting from the right end of the word), the 2nd digit position is bits 19-21,..., the 8'th digit position is bits 1-3. The Chi-Square probability is a number between and 1 which indicates the 5^ probability that a Chi -Square variable "will exceed the calculated test statistic. If IFOP = 2., then PROKL(l), I = 1,...,8 will be arbitrary upon exit from POKDIG. PR0B2 - If IFOP = 1 then PR0B2 will be arbitrary upon exit from POKDIG. If IPOP = 2 then PR0B2 contains the Chi-Square probability for the poker test. It is a number between and 1 which indicates the probability that a Chi-Square variable will exceed the calcu- lated test statistic. 55 4 14 lb 20 26 1000 SUBROU1 IMP POKDIG (N, IRN, IF0P,PRr>Bl,PK(JB2 ) FREQUENCY TEST AMI) P OKER TEST UN T HE 01 GITS REAL*4 PXPPOK ( 7 ) , PRUB1 ( 8 ) INTEGERS MCT ( H, 8) ,POK(7) , IRN(N) i) 1 = 1, H 5 Dfl 3 J = l,8 ()CT( I ,J)=0 Dfl 4 1 = 1,7 PMK ( T )=() CALL PR F Pf. IK ( N , I RN.flLQ J_ t RJJKjLl FOP-l: I F ( IFUP.E0.2 )GU TO 20 FXPFRE=)M/8, 00 lb J=1,R XNUM=0 , . DU 14 1=1, H XNUM=XNUM+(UC I ( I , J )-PXPPRE)**2 CHI FRE = Xi\iUM/EXPERE CALL PR OL I M ( 7 , C H I PR P , P R U B 1 ( J ) ) GO IN 1000 ,?0i>07R*N i>12695*N lb3809#l\! 102b39*N 01 7090#N 008b45*N 000244*IM EXPPUK ( 1 )=0 PXPPOK ( 2 ) =0 PXPPUK ( 3 )=0, PXPPMK ( 4) =0 PXPPUK ( b)=0. PXPPMK (6)=0 PXPPUK ( 7 )=0, CHI PMK = 0.() Ull 2 b 1 = 1,7 CHIP(IK = CHIPl)K+(PHK( I )-PXPP()K( I ) ) CALL PRULIto{6,CHIP0K,PR0B2) RETURN FNO 2/PXPPl)K( I ) 56 SUBROUTINE NAME *** FREPOK IDENTIFICATION *** A subroutine called by POKDIG PURPOSE **"* The 2k bit random integers are divided into 8 octal digits each. In the poker test the first 5 octal digits of each integer are treated as a poker hand without regard to suit. The frequencies of cer- tain poker hands are tallied by FREPOK. In the frequency test the fre- quency of each octal digit in each of the 8 digit positions is tallied by FREPOK. OTHER SUBROUTINES USED *** None USAGE *** INTEGER*^- IRN(NX),0CT(8,8),P0K(7) • CALL FREP0K(N,IRN,0CT,P0K,IF0P) Where the parameters in the calling sequence are : N - The number of random integers to be used in the test, NX ^ N. IRN - The one dimensional array containing the N random integers to be tested. OCT - If IFOP = 1 this two dimensional array is set to zero on entry to FREPOK. On exit from FREPOK OCT(l,j) ; I,J = 1, ...,8 contains the frequency of the octal digit 1-1 in the J'th digit position where bits 22.-2h are the 1st digit position,..., bits 1-3 are the 8th digit position. If IFOP = 2 the contents of this array are arbitrary on entry and exit from FREPOK. POK - If IFOP = 2 then POK(l), I = 1, . . . ,7 is set to zero on entry to FREPOK. On exit from FREPOK POK contains the frequencies of 7 poker hands where : POK(l) = frequency of bust P0K(2) = frequency of one pair 57 P0K(3) = frequency of two pair POK^) = frequency of 3 of a kind P0K(5) = frequency of a full house P0K(6) = frequency of k of a kind P0K(7) = frequency of 5 of a kind If IFOP = 1 then the contents of this array are arbitrary on entry and exit from FKEPOK. IFOP - If IFOP = 1 then the frequency of digits is to be tallied. If IFOP 2 then the frequency of poker hands is to be tallied. 58 FR hPHK i\i (- U K iv i T N I I t-^HSP I MI-WI IC I" ii Ik Hhkt- ERFTST Pf IK IS! r E G I w LA CR RC SR SR S I LA CR kC I. )L SI. SR SR SR SLDL SR x _>•••> x , < x , , . is a run down of length k, n-1 n n+1 n+k n+k+1 ' and x -. > x ^ x n <«««^x n > x , - is a run up of length k. n-1 n ^ n+1 ** ^ n+k n+k+1 ° This test involves tallying runs up and down of different lengths and comparing with expected results by a Chi-Square test. OTHER SUBROUTINES USED *** PROLIM USAGE *** REAlM RN(NX),EXPRUN(100) INTEGER^ RUNLTH(IOO) CALL RNSUAD(N,RN,EXPRUN,RUNLTH,PR0B) Where the parameters in the calling sequence are : N - The number of random numbers to be tested, NX ^ N. RN - The one dimensional array containing the N random numbers to be tested. EXPRUN - Internal array used by RNSUAD, its contents are arbitrary on entry to RNSUAD. On exit from RNSUAD it contains the expected frequencies of runs of lengths 1, . . .m where m is the maximum run length found in the N random numbers. RUNLTH - Internal array used by RNSUAD, its contents are arbitrary on entry to RNSUAD. On exit from RNSUAD it contains the actual frequencies of runs of lengths 1, . . .m where m is the maximum run length found in the N random numbers. 62 PROB - The outputed value of the Chi -Square probability. It is a number between and 1 which indicates the probability that a Chi -Square variable will exceed the calculated test statistic. 63 10 13 14 ] S ?3 ?4 ?5 SO 41 M) I Id 1 A kn'( i+T) )gi TO 10 SUBROUTINE RNSUAD(N, RN, EX PR UN , RUNLTH , PROB ) TEST DF RUNS U P AND DOWN REAL#4 RN(N) , EXPRUN( 1 ) INTEJ2ER.*^_RJJNLTHJ 1) M4XLTH=6 DO 2 1 = 1 , 100 RUNLTH ( I )=0 LAS "I 7=0 . LAST 1 = 1 = 1 TF(RM( I ) ,GT.RN( T + l ) )G0 TO 10 1 = 1 + 1 IF( I ,EO.N)G(J Tl) 13 IF(RN( I ) .LE.RN( 1+1) )G0 T b K= I -LAST 1-1 L A S T 7 = I - 1 MSWICH=1 GO Tm IS 1=1 + 1 I F ( I ,EO.I\l)G( IF{RN{ I ) .GT K= I -LA ST 7.-1 L A S T 1 = I - 1 NSWICH=0 GO III IS K=M-LAST1-1 GO TO IS K=m-LAST7-1 RUNLTH(K )=RUi\iLTH( K ) + 1 IF ( K.G1 .MAXLTH) mAXLTH=K IF( I .Em. m) GU TO ?3 IF(NSWlCH.EO._0JGg Id 5 GO "in 10 U(i ^s I = 1,MAXLTH IP3=I+3 EACT=1.0 DO ?A J=l, IP3 FACT=FACT*J ISU*T**2 EXPRON( I ) = ?.*(N = I M AX = i"AXLIH Ull Sf) 1 = 1, I MAX I 1= lMAX-I+1 IE (ROiMLTHl I I ) .GE.SlGO TO SI RUNLTH ( I 1-1 ) =RUNLTH ( I I — 1 ) +RUNLTH ( I I ) EXPRUN( I I -1 ) = EXPRUN( I I -1 ) +EXPRUN( I I ) MAXLTH=MAXLT H-l Chi sn = f).o on (So I = i,MAXLTH CHI SO = CHI S0+( RUNLTH ( I ) -EX PR UN ( I ) )_**2/ EXPRUNI I ) NUF=i"iAXLTH-l CALL PRUL IM( MDFtCHISQf PROB) R E T R N E M I ) ( TSU+3.*I + 1. )-( I SO* I + 3.':= I SO- I -4. ) ) /FACT 6k SUBROUTINE NAME *** RNSABM IDENTIFICATION *** Test of runs above and below the mean on a set of random numbers. PURPOSE *** If x, , . . .x„ are the set of random numbers then x _ l/2, x p > l/2,...,x > l/2 is a run above the mean of length k. This test involves tallying runs above and below the mean of different lengths and comparing with expected results by a Chi-Square test. OTHER SUBROUTINES USED *** PROLIM USAGE *** REAlM RN(N), EXPRUN(IOO) INTEGER *k RUNLTH(IOO) • GALL RNSABM(N,RN,EXPRUN,RUNLTH,PR0B) Where the parameters in the calling sequence are: N - The number of random numbers to be tested, NX > N. RN - The one dimensional array containing the N random numbers to be tested. EXPRUN - Internal array used by RNSABM, its contents are arbitrary on entry to RNSABM. On exit from RNSABM it contains the expected frequencies of runs of lengths 1, . . .m where m is the maximum run length found in the N random numbers . RUNLTH - Internal array used by RNSABM, its contents may be arbitrary on entry to RNSABM. On exit from RNSABM it contains the actual frequencies of runs of lengths 1, ...,m where m is the maximum run length found in the N random numbers. 65 PROB - The outputed -value of the Chi -Square probability. It is a number between and 1 which indicates the probability that a Chi-Square variable will exceed the calculated test statistic. 66 si 10 13 14 ] i 33 3 b bO bl 60 Rl F M c I) c c IB PS HA Ml AX I I DM AS AS = 1 PI F( = 1 F( F( = 1 AS s w II = 1 F< F( = 1 AS sw 1 1 = M = l\l I I'M F( F( F( I I 1 1 XP M A I I 1 = F( IM XP AX Ml WIUITINF RNSARM( N v RN , HXPRUN, RUNLTH, PRUR) T (IF RUMS ABOVh AMI) BF LUW THE MEAN L*4 RM(M), FXPRI)M( 1 ) FGFR#4 RUMLTH( 1 ) LTH=0 2 1=1,100 L T H ( I ) = TZ = T1 = HI OF AL F I Ml) = M+1 RM( I ) + 1 I .En. R IV ( I ) -LAST 17=1- ICH=1 Til lb + 1 I . F . R N ( I ) -LAS T II = 1- ICH = Tfl lb - L A S T Tfl lb -LAST L I H(K K . G T . I .to. MS WIG III 10 3b 1 = RDM ( I X = MAX 50 1 = I M A X - RUN LI LTH( I RUN( I LTH=M sn=0. (SO i = Sn = LH = MAXL L FPU URN ,GT.0.5)GU 1U 10. MF1 )Gll TU 13 .LE.0.b)GU 111 b 1-1 1 mfi )c;u IN 14 .GT.0..b)G() IN 10 Z-l 1 1 7 )=RHMLTH( K )+l MAXLTH)MAXLTH=K NPDGll Tfl 33 H.EO.0)Gll TO b 1 (MAXLTH )=(n-1+3. )#2.**(-I-l ) L I H 1 , I M A X 1 + 1 H( II ) .Gt.bJGH Tf.) bl 1-1 )=FIIML I M ( I I — 1 )+RHMLTH( I I ) 1-1 )=EXPRUN( I 1-1 )+FXFRUN( I I ) AXLTH-1 1 , M A X L 1 HI I SO+ ( RUWLTHI I ) -FXFRUM( I ) }**2/EXPRUN< I ) IH-1 L I M { NO F f C H I SO t P R fl R ) 67 £> A P' SUBROUTINE NAME *** maxmin IDENTIFICATION *** Frequency test on maximums or minimums of n-tuples of a set of random numbers. If W is the maximum of the n numbers in an n-tuple then X = w is uniformly distributed^ if W is the minimum of the n numbers in an n-tuple then Y = l-(l-W) is uniformly dis- tributed. In this test the unit interval is divided into 100 sub- intervals and the frequency of X's and Y's is tallied and compared against expected results by a Chi-Square test. OTHER SUBROUTINES USED *** UNINUM USAGE *** REAlM RN(NX),MX0MN(MX) • GALL MAXMIN(N,RN,]OTUPI£,MXOMN,IXON,PROB) Where the parameters in the calling sequence are : N - The number of random numbers to be tested, NX ;> N. RN - The one dimensional array containing the N random numbers to be tested. NTUPLE - The number of random numbers in a tuple. Note that the tuples are taken consecutively from the N random numbers. MXOMN - An internal array used by MAXMIN, MX ^.n/NTUPLE. Its contents are arbitrary on entry to MAXMIN; on exit from MAXMIN it con- tains the N/NTUPLE X's or n/NTUPLE Y's. IXON - If IXON = 1 then the maximum of n-tuples is to be tested. If LXON =s 2 then the minimum of n-tuples is to be tested. PROB - The outputed value of the Chi-Square probability. A number between and 1 which indicates the probability that a Chi-Square variable will exceed the calculated test statistic. 68 id ?() 30 40 SUBROUTINE MAXHIN( N f RN, N TUPLE , MXOMN, IXON, PROB) REAL*4 RN(N) , MXOMIM( 1) _. I M = N / 1\| i I J P L E I k iv, a X = IN*NTyPLE-NTUPLE+l i\lTMl=NTUPLE-I 1 = 1 l)(") 40 J = l , TRMAX,NTUPLE MXUMM( I ) =«IM( J ) ______ ._ ,'__ I E( IX0N.E0.2 )GU TO 20 DO 10 K=l t MThl IF(RM( J + K) . GT.MXOM!M( 1 ) ) MXOMN ( I ) =RN ( J + K ) MXUMNt I )=MXOMN( 1 )**N.TUPLE. GO TO 40 D\\ 30 K = l , iMTMl IF(RM( J+K) .LT.MXOMNI I ) )MXOMN( I )=RN( J + K ) MXOMN( T ) =1 .-( 1 . -MXOMN( I ) ) **NTOPLE 1 = 1 + 1 CALL UN I IMUM ( lN f MXI IMM , PKUH ) R E T U R N FND 6 9 SUBROUTINE NAME *** PROLIM IDENTIFICATION *** Chi -Square probability PURPOSE *** Given a Chi-Square statistic X and the number of degrees of freedom n, this subroutine calculates the probability that a Chi- 2 Square variable will exceed X . For n ^ 30 this is done by evaluating /DO II f(x)dx where f(x) = x e is the Chi-Square probability x 2 2n/2r(n/2) distribution. For n > 30 the quantity n/2X - \/2n-l is normally distributed f °° - 2 /2 so / p f (x)dx where f(x) = 1 e ' ' X n/2^ is evaluated. OTHER SUBROUTINES USED *** None USAGE *** CALL PROLIM(NDF CRTSQ,PROB) Where the parameters in the calling sequence are : NDF - The number of degrees of freedom, an input parameter. 2 CHISQ - The Chi-Square statistic X , an input parameter. PROB - The Chi-Square probability, an output parameter. It is a number between and 1 which indicates the probability that a Chi-Square variable will exceed the statistic CRTSQ. 70 SI muni JUNE PRML Im( NUF,UNCSQ, PROB) K( L AL*4 X(96') f W(%) INC1 ( I ) = (?**< (NUH-2)/2. ) )*EXP(-T/2. ) IMC2 ( T )=bXP( (-1**2) /2. ) 01) = 00. 0162767 A 02) = -0,01627674 03) = 00.04Rfi].29R 04) = -0.04HR129H Oh) = 00. OR 129749 06) = -0.08129749 07 ) = 00. 1 1 369585 08) = -0.11369585 09) = 00.1459737] 10) = -0.14597371 11) = 00. 17R0968R 12) = -0.178096R8 13) = 00.21003131 14) = -0.21003131 l b ) = 00.241 74315 16) = -0.24174315 17) = 00.27319881 in) = -0.27319RR1 19) = 00.30436494 20) = -0.30436494 21 ) = 00. 3352085 2 22) = -0.33520852 23) = 00.36569686 24) = -0.36569686 25) = 00.39579764 26) = -0.39579764 27) = 00.4254789 8 ?H) = -0.42547898 29) = 00.45470942 30) = -0.45470942 3] ) = 00.48345797 32) = -0.48345797 33) = 00.51169417 34) = -0. hi 1694 17 35) = 00.53938810 36) = -0.53938810 37) = 00.36651041 38) = -0.56651041 3S ) = 00.59 30 32 36 X(40) = -0.59303236 71 X( 41 ) X ( 42) X ( 43) X ( '-4) X ( 4 3 ) X ( 46) X( 47) X( 48) X ( z.9) X ( SO) X( SI ) X ( 52) x ( 5 3 ) x ( 54) x ( 55 ) X (56) x (57) X ( SK ) X ( 59 ) X ( 60) X (61 ) X ( 62 ) X ( 6 3 ) X ( 64) X ( 6b ) X ( 66 ) X (67) x [ 6 8 ) X ( 6 9 ) X 70) X [71) X 7?) X 73 ) X I 74) X 75 ) X ( 76 ) X 1 7 7 ) X ( 78 ) x ( 79 ) X ( MO ) X ( 81 ) x ( 82 ) X ( 8 3 ) X ( 84) x ( 83 ) X ( 8 6 ) x ( 8 7 ) X ( « 8 ) x ( 8 9 ) X ( 9(1) > ( 91 ) = 00.6 = -0.6 = 00.6 = -0.6 = 00.6 = -0.6 = 00.6 = -0.6 = 00.7 = -0.7 = 00.7 = -0.7 = 00.7 1892584 1H92SH4 4416340 4416 340 6 8 718 3] 6M71M31 9256433 9256463 1567681 1367681 3803064 3 M 3 6 4 3 4 6 2 34 -Q. 75960234 00. 78036904 -0. 780369)04 00. 800308 74 -0. 80030874 00.81940031 -0. 8144003] 00. 837623&1 -0.83762 331 00. 85495903 -0. 8 34-95903 .87138850 -0. 871 38830 00. 8868 9451 -0.R8689451 00. 90146063 -0.90146063 00.9 13071 42 -0.9] 307142 00.92771248 -0.92771245 00.9 3937033 -0.9 34 37 33 00.93003271 -0.9 300 32 71 00.93968829 -0.95968 829 00.96*32682 -0.96832682 00.9739 391 7 -0.97393917 00.982.51726 -0.48231726 00.98805412 -0.988034] 2 0.442^4340 -0.99 2 34390 00. 44348 1 P4 72 X 92) X (93) X 94) X [9b ) X 96) w (01 ) l»J 02 ) Id [03) w 04) w [ OS ) w 06) I.I [07 ) W OH ) w [ 09 ) u;' 10) 1,1 ( 11 ) Ul 1 2) N Hi) lr ] 4) -0.9 0.9 -0 . 9 00.9 -0.9 0.0 3 03 03 03 03 3 3 03 03 03 03 3 3 03 9b9R184 98 36437 9836437 9968950 99689b6 255061 255061 251611 251611 2 44716 244716 ? 343 H? 23438 2 2 2 062 220620 20344b 20344b 182 87b 18287b '-'(IS) W (16) ta ( 17 ) W (18) i-i ( l (-0 L-J ( 20 ) W (2 1) w(22) 11 ( 2 3 ) '-(24) W (25) W (26) w( 27 ) w (28) H (29) W( 30) W (31) i- 1 ( 32 ) W ( 3 3 ) W ( 34 ) w (3b) H ( 36 ) '" (37) W ( 3H) i,i ( 39 ) W ( 40 ) ■ ' ( 4 1 ) t" ( 4 2 ) H ( 4 3 ) '•' ( LyU ) 11 (4b ) |,| ( Uh ) 3 1 b 8 9 3 3 3 1 b 8 9 3 3 03] 3] 642 31316 4 2 310 10 3 3 3 1 1 3 3 306 7137 0306713 f 03029941 30 2 9991 02989634 02989634 02946108 = 0.024461.08 = 0.02899461 = 0.02899461 = 0.0284974] = 0.02849741 = (J. 02 79 7000 = 0.02797000 = 0.02741296 = U. 02 74 12 96 = 0.02682686 = 0.02682686 = 0.02621234 = 0.02621234 02 bb 700 3 ,02557003 02490063 02490063 02420484 = 0. 02420484 73 (47) (48) (49) ( SO) (51) ( 52 ) (63) ( 64) ( 6 t 5 ) ( 56) (57) ( 58 ) (59) ( 60) (61 ) (62) ( 63 ) ( 64) (65 ) ( 66) (67 ) ( 6R) (69) (70) (71 ) (72) (73) ( 74) (75 ) (76) (77) (7R) (79) ( R0) (5] ) ( R2 ) (83) ( 84) ( 85 ) ( 86) ( h 7 ) ( HM ) (89) ( 90 ) (9] ) ( 92 ) (93) ( 94) ( 95 ) (96) 0.02 0.02 0.02 O.J) 2 0.02 0.0 2 0.02 0.0 2 0.02 ,02 01 ,01 01 01 01 1 01 01 01 01 01 01 01 01 01 348339 348 3 39 273706 2 7 37 06 196664 196664 1.17293 117293 03 5679 3 6679 961908 9 51908 866067 R 6 6 6 7 778260 778250 688847 6RR547 59 70 56 697056 503R72 503872 409094 409094 3 1 2 R 2 2 0.01312822 0.01215160 0.012 1 5160 0.01 116210 0.011162 10 0.01016077 0.010 160 7 7 0.00914H67 .00914867 0.00812687 0.00R126R7 0. 00709647 0.00709647 0. 00608854 0.00605854 0.00 50 142 0.00501 42 0.00396455 0.00396455 0.00291073 0.0029 1.07 3 0. 001 R5396 0.00185396 0.00079679 0.0007^679 / PKI 18 = 0.0 7^ IF ( NDF.G1 ,30)G0 TU 20 B = UNCSO_/2_. S = . IT] 10 J = 1,96 __ Y=B*X( J )+B 10 S=S+HJNC1(Y) *W( J) _ PROB=B*S+PROB l.)EN= ( ?.** ( NDF/2. ) ) *G_AMMAJ_N_DF_/_2._i Pkf)B=l .-PROB/'DEN R H T U R N 20 >il)p = B=SOR I ( 2.*UNCSO)-SOR7 ( 2.*NDF-1. ) IF(R.GF.0.0)GU TU 2 5 B="-B '»U)P=1 25 B=B/2. S = 0.0 DfJ 30 J = 1,96 Y = B * X ( J ) + B 40 S = S + FI.INC2 ( Y) *w{ J ) PR()B = H*S + PRUB PROB=. 398 9422*PRnR I F ( MDP.EO. 1 )GU TO 40 PR()B=. 5-PRGB n>F ] URi\j 40 PR| IB = . ^+pkl)H RFTURN FMO