Chapter 9Random NumbersThis chapter describes algorithms for the generation of pseudorandom numbers withboth uniform and normal distributions.9.1 Pseudorandom NumbersHere is an interesting number:0.95012928514718This is the first number produced by the Matlab random number generator withits default settings. Start up a fresh Matlab, set format long, type rand, and it’sthe number you get.If all Matlab users, all around the world, all on different computers, keepgetting this same number, is it really “random”? No, it isn’t. Computers are (inprinciple) deterministic machines and should not exhibit random behavior. If yourcomputer doesn’t access some external device, like a gamma ray counter or a clock,then it must really be computing pseudorandom numbers. Our favorite definitionwas given in 1951 by Berkeley professor D. H. Lehmer, a pioneer in computing and,especially, computational number theory:A random sequence is a vague notion . . . in which each term is unpre-dictable to the uninitiated and whose digits pass a certain number oftests traditional with statisticians . . .9.2 Uniform DistributionLehmer also invented the multiplicative congruential algorithm, which is the basis formany of the random number generators in use today. Lehmer’s generators involvethree integer parameters, a, c, and m, and an initial value, x0, called the seed. ADecember 26, 200512 Chapter 9. Random Numberssequence of integers is defined byxk+1= axk+ c mod m.The operation “mod m” means take the remainder after division by m. For example,with a = 13, c = 0, m = 31, and x0= 1, the sequence b egins with1, 13, 14, 27, 10, 6, 16, 22, 7, 29, 5, 3, . . . .What’s the next value? Well, it looks pretty unpredictable, but you’ve been ini-tiated. So you can compute (13 · 3) mod 31, which is 8. The first 30 terms inthe sequence are a permutation of the integers from 1 to 30 and then the sequencerepeats itself. It has a period equal to m − 1.If a pseudorandom integer sequence with values between 0 and m is scaledby dividing by m, the result is floating-point numbers uniformly distributed in theinterval [0, 1]. Our simple example begins with0.0323, 0.4194, 0.4516, 0.8710, 0.3226, 0.1935, 0.5161, . . . .There is only a finite number of values, 30 in this case. The smallest value is 1/31;the largest is 30/31. Each one is equally probable in a long run of the sequence.In the 1960s, the Scientific Subroutine Package (SSP) on IBM mainframecomputers included a random number generator named RND or RANDU. It was amultiplicative congruential with parameters a = 65539, c = 0, and m = 231. Witha 32-bit integer word size, arithmetic mod 231can be done quickly. Furthermore,because a = 216+3, the multiplication by a can be done with a shift and an addition.Such considerations were important on the computers of that era, but they gavethe resulting sequence a very undesirable property. The following relations are alltaken mod 231:xk+2= (216+ 3)xk+1= (216+ 3)2xk= (232+ 6 · 216+ 9)xk= [6 · (216+ 3) − 9]xk.Hencexk+2= 6xk+1− 9xkfor all k.As a result, there is an extremely high correlation among three successive randomintegers of the sequence generated by RANDU.We have implemented this defective generator in an M-file randssp. A demon-stration program randgui tries to compute π by generating random points in a cubeand counting the fraction that actually lie within the inscribed sphere. With theseM-files on your path, the statementrandgui randsspwill show the consequences of the correlation of three successive terms. The resultingpattern is far from random, but it can still be used to compute π from the ratio ofthe volumes of the cube and sphere.9.2. Uniform Distribution 3For many years, the Matlab uniform random number function, rand, wasalso a multiplicative congruential generator. The parameters werea = 75= 16807,c = 0,m = 231− 1 = 2147483647.These values are recommended in a 1988 paper by Park and Miller [10].This old Matlab multiplicative congruential generator is available in the M-file randmcg. The statementrandgui randmcgshows that the points do not suffer the correlation of the SSP generator. Theygenerate a much better “random” cloud within the cube.Like our toy generator, randmcg and the old version of the Matlab functionrand generate all real numbers of the form k/m for k = 1, . . . , m − 1. The smallestand largest are 0.00000000046566 and 0.99999999953434. The sequence repeatsitself after m − 1 values, which is a little over 2 billion numb ers. A few years ago,that was regarded as plenty. But today, an 800 MHz Pentium laptop can exhaustthe period in less than half an hour. Of course, to do anything useful with 2 billionnumbers takes more time, but we would still like to have a longer perio d.In 1995, version 5 of Matlab introduced a completely different kind of randomnumber generator. The algorithm is based on work of George Marsaglia, a professorat Florida State University and author of the classic analysis of random numbergenerators, “Random numb ers fall mainly in the planes” [6].Marsaglia’s generator [9] does not use Lehmer’s congruential algorithm. Infact, there are no multiplications or divisions at all. It is specifically designed toproduce floating-point values. The results are not just scaled integers. In place of asingle seed, the new generator has 35 words of internal memory or state. Thirty-twoof these words form a cache of floating-point numbers, z, between 0 and 1. Theremaining three words contain an integer index i, which varies between 0 and 31, asingle random integer j, and a “borrow” flag b. This entire state vector is built upa bit at a time during an initialization phase. Different values of j yield differentinitial states.The generation of the ith floating-point number in the sequence involves a“subtract-with-borrow” step, where one number in the cache is replaced with thedifference of two others:zi= zi+20− zi+5− b.The three indices, i, i+ 20, and i+5, are all interpreted mod 32 (by using just theirlast five bits). The quantity b is left over from the previous step; it is either zero ora small positive value. If the computed ziis positive, b is set to zero for the nextstep. But if the computed ziis negative, it is made positive by adding 1.0 before itis saved and b is set to 2−53for the next step. The quantity 2−53, which is half ofthe Matlab constant eps, is called one ulp because it is one unit in the last placefor floating-point numbers