f7-1

7.1UniformDeviates 267Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).As for references on this subject, the one to turn to first is Knuth [1]. Then try[2]. Only a few of the standard books on numerical methods [3-4]treat topics relating to random numbers. CITED REFERENCES AND FURTHER READING: Knuth,D.E.1981, SeminumericalAlgorithms ,2nded.,vol.2of TheArtofComputerProgramming (Reading, MA: Addison-Wesley), Chapter 3, especially §3.5. [1] Bratley, P., Fox, B.L., and Schrage, E.L. 1983, A Guide to Simulation (New York: Springer- Verlag). [2] Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), Chapter 11. [3] Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 10. [4] 7.1 Uniform Deviates Uniform deviates are just random numbers that lie within a specified range (typically0to 1),withanyonenumberintherangejust as likelyas anyother. They are, in other words, what you probably think “random numbers” are. However, we want to distinguish uniform deviates from other sorts of random numbers, for example numbers drawn from a normal (Gaussian) distribution of specified meanandstandarddeviation. Theseothersortsofdeviatesarealmostalwaysgeneratedby performing appropriate operations on one or more uniform deviates, as we will see insubsequentsections. So,areliablesourceofrandomuniformdeviates,thesubjectof this section, is an essential building block for any sort of stochastic modeling or Monte Carlo computer work. System-Supplied RandomNumberGenerators Yourcomputerverylikelyhaslurkingwithinitalibraryroutinewhichiscalled a“randomnumbergenerator.” Thatroutinetypicallyhas anunforgettablenamelike “ran,” and a calling sequence like x=ran(iseed) sets xto the next random number and updates iseed You initialize iseedto a (usually) arbitrary value before the first call to ran. Eachinitializingvaluewill typicallyreturnadifferentsubsequentrandomsequence,oratleastadifferentsubsequenceofsomeoneenormouslylongsequence. The same initializingvalueof iseedwill always returnthe samerandomsequence,however. Now our first, and perhaps most important, lesson in this chapter is: Be very, verysuspiciousofasystem-supplied ranthatresemblestheonejustdescribed. Ifall scientific papers whose results are in doubt because of bad rans were to disappear from library shelves, there would be a gap on each shelf about as big as your fist. System-supplied rans arealmost always linearcongruentialgenerators ,which 268 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).generate a sequence of integers I1,I2,I3,..., each between 0andm−1(a large number) by the recurrence relation Ij+1=aI j+c(mod m)( 7.1.1 ) Here mis called the modulus,and aandcare positiveintegers called the multiplier and theincrement , respectively. The recurrence(7.1.1)will eventuallyrepeat itself, withaperiodthatisobviouslynogreaterthan m.I fm, a,andcareproperlychosen, thentheperiodwillbeofmaximallength,i.e.,oflength m. Inthatcase,allpossible integersbetween0and m−1occuratsomepoint,soanyinitial“seed”choiceof I0 isasgoodasanyother: Thesequencejusttakesofffromthatpoint. Therealnumber between0and1whichis returnedis generally Ij+1/m,sothatit isstrictlyless than 1, but occasionally(once in mcalls) exactly equal to zero. iseedis set to Ij+1(or someencodingofit),sothatitcanbeusedonthenextcalltogenerate Ij+2,andsoon. The linear congruentialmethod has the advantageof being very fast, requiring onlyafewoperationspercall,henceitsalmostuniversaluse. Ithasthedisadvantage thatitisnotfreeofsequentialcorrelationonsuccessivecalls. If krandomnumbersat a time are used to plotpoints in kdimensionalspace (with eachcoordinatebetween 0 and 1), then the points will not tend to “fill up” the k-dimensional space, but rather will lie on (k−1)-dimensional “planes.” There will be at mostabout m1/k such planes. If the constants m,a, and care not very carefully chosen, there will bemany fewer than that. The number mis usually close to the machine’s largest representable integer, e.g., ∼232. So, for example, the number of planes on which triples of points lie in three-dimensional space is usually no greater than about the cube root of 232, about 1600. You might well be focusing attention on a physical process that occurs in a small fraction of the total volume, so that the discretenessof the planes can be very pronounced. Even worse, you might be using a ranwhose choices of m, a,andchave been botched. One infamous such routine, RANDU, with a= 65539 andm=2 31, was widespread on IBM mainframe computers for many years, and widely copied onto other systems [1]. One of us recalls producing a “random” plot with only 11 planes, and being told by his computer center’s programming consultant that he had misused the random number generator: “We guarantee that each number is randomindividually,butwedon’tguaranteethat morethanoneofthemis random.”Figure that out. Correlationin k-spaceisnottheonlyweaknessoflinearcongruentialgenerators. Such generatorsoften havetheir low-order(least significant) bits much less random than their high-order bits. If you want to generate a random integer between 1 and 10, you should always do it using high-order bits, as in j=1+int(10.*ran(iseed)) and never by anything resembling j=1+mod(int(1000000.*ran(iseed)),10) 7.1UniformDeviates 269Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).(which uses lower-order bits). Similarly you should never try to take apart a “ran” number into several supposedly random pieces. Instead use separate calls for every piece. Portable Random NumberGenerators ParkandMiller [1]havesurveyedalargenumberofrandomnumbergenerators that have been used over the last 30 years or more. Along with a good theoretical review,theypresentananecdotalsamplingofanumberofinadequategeneratorsthat havecome into widespreaduse. Thehistorical recordis nothingif notappalling. There is good evidence, both theoretical and empirical, that the simple multi- plicative congruential algorithm Ij+1=aI j(mod m)( 7.1.2 ) can be as good as any of the more general linear congruential generators that have c/negationslash=0(equation 7.1.1) — ifthe multiplier aand modulus mare chosen exquisitely carefully. Park and Miller propose a “Minimal Standard” generator based on the choices a=75= 16807 m=231−1 = 2147483647 ( 7.1.3 ) First proposed by Lewis, Goodman, and Miller in 1969, this generator has in subsequent years passed all new theoretical tests, and (perhaps more importantly) hasaccumulatedalargeamountofsuccessfuluse. ParkandMillerdonotclaimthatthe generatoris “perfect”(we will see below that it is not), but onlythat it is a good minimal standard against which other generators should be judged. It is not possible to implement equations (7.1.2) and (7.1.3) directly in a high-level language, since the product of aandm−1exceeds the maximum value for a 32-bit integer. Assembly language implementation using a 64-bit productregister is straightforward, but not portable from machine to machine. A trick due to Schrage [2,3]for multiplying two 32-bit integers modulo a 32-bit constant, withoutusinganyintermediateslargerthan32bits (includinga signbit)is thereforeextremely interesting: It allows the Minimal Standard generator to be implemented in essentially any programming language on essentially any machine. Schrage’s algorithm is based on an approximate factorization ofm, m=aq+r,i.e., q=[m/a ],r =mmod a (7.1.4 ) with square brackets denoting integer part. If ris small, specifically r<q, and 0<z<m −1, it can be shown that both a(zmod q)andr[z/q ]lie in the range 0,...,m −1, and that azmod m=/braceleftbigg a(zmod q)−r[z/q ]if itis ≥0, a(zmod q)−r[z/q ]+motherwise(7.1.5 ) The application of Schrage’s algorithm to the constants (7.1.3) uses the values q= 127773 andr= 2836. Here is an implementation of the Minimal Standard generator: 270 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).FUNCTION ran0(idum) INTEGER idum,IA,IM,IQ,IR,MASK REAL ran0,AM PARAMETER (IA=16807,IM=2147483647,AM=1./IM, * IQ=127773,IR=2836,MASK=123459876) “Minimal” random number generator of Park and Miller. Returns a uniform random deviate between 0.0 and 1.0. Set or reset idumto any integer value (except the unlikely value MASK) to initialize the sequence; idummust not be altered between calls for successive deviates in a sequence. INTEGER k idum=ieor(idum,MASK) XORing with MASKallows use of zero and other simple bit patterns for idum. k=idum/IQ idum=IA*(idum-k*IQ)-IR*k Compute idum=mod(IA*idum,IM) without overflows by Schrage’s method. if (idum.lt.0) idum=idum+IM ran0=AM*idum Convert idumto a floating result. idum=ieor(idum,MASK) Unmask before return. return END The period of ran0is231−2≈2.1×109. A peculiarity of generators of the form (7.1.2) is that the value 0must never be allowed as the initial seed — it perpetuates itself — and it never occurs for any nonzero initial seed. Experience hasshownthatusersalwaysmanagetocallrandomnumbergeneratorswiththeseed idum=0. Thatis why ran0performsits exclusive-orwith anarbitraryconstantboth onentryandexit. Ifyouarethe first userinhistorytobe proofagainsthumanerror, you can remove the two lines with the ieorfunction. Park and Miller discuss two other multipliers athat can be used with the same m=231−1. Theseare a= 48271 (with q= 44488 andr= 3399)and a= 69621 (with q= 30845 andr= 23902). These can be substituted in the routine ran0 if desired; they may be slightly superior to Lewis et al.’s longer-tested values. No values other than these should be used. The routine ran0is a Minimal Standard,satisfactory for the majority of appli- cations,butwedonotrecommenditasthefinalwordonrandomnumbergenerators. OurreasonispreciselythesimplicityoftheMinimalStandard. Itisnothardtothink ofsituationswheresuccessiverandomnumbersmightbeusedinawaythatacciden- tallyconflictswiththegenerationalgorithm. Forexample,sincesuccessivenumbers differbyamultipleofonly 1.6×104outofamodulusofmorethan 2×109,verysmall randomnumberswill tend to be followed by smaller than average values. One time in106, for example, there will be a value <10−6returned(as there should be), but this willalwaysbefollowedbya valueless thanabout 0.0168. Onecaneasily think ofapplicationsinvolvingrareeventswherethispropertywouldleadtowrongresults. There are other, more subtle, serial correlations present in ran0. For example, if successive points (Ii,Ii+1)are binned into a two-dimensional plane for i= 1,2,...,N, thenthe resulting distributionfails the χ2test when Nis greaterthan a few×107,muchlessthantheperiod m−2. Sincelow-orderserialcorrelationshave historically been such a bugaboo, and since there is a very simple way to remove them, we think that it is prudent to do so. The following routine, ran1, uses the Minimal Standard for its random value, but it shuffles the output to removelow-orderserial correlations. A random deviate derivedfromthe jthvalueinthesequence, Ij,isoutputnotonthe jthcall,butrather onarandomizedlatercall, j+3 2onaverage. TheshufflingalgorithmisduetoBays and Durham as described in Knuth [4], and is illustrated in Figure 7.1.1. 7.1UniformDeviates 271Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).FUNCTION ran1(idum) INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV REAL ran1,AM,EPS,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, * NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS) “Minimal” random number generator of Park and Miller with Bays-Durham shuffle and added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive ofthe endpoint values). Call with iduma negative integer to initialize; thereafter, do not alter idumbetween successive deviates in a sequence. RNMXshould approximate the largest floating value that is less than 1. INTEGER j,k,iv(NTAB),iySAVE iv,iyDATA iv /NTAB*0/, iy /0/ if (idum.le.0.or.iy.eq.0) then Initialize. idum=max(-idum,1) Be sure to prevent idum =0. do 11j=NTAB+8,1,-1 Load the shuffle table (after 8 warm-ups). k=idum/IQ idum=IA*(idum-k*IQ)-IR*kif (idum.lt.0) idum=idum+IMif (j.le.NTAB) iv(j)=idum enddo 11 iy=iv(1) endif k=idum/IQ Start here when not initializing. idum=IA*(idum-k*IQ)-IR*k Compute idum=mod(IA*idum,IM) without overflows by Schrage’s method. if (idum.lt.0) idum=idum+IM j=1+iy/NDIV Will be in the range 1:NTAB . iy=iv(j) Output previously stored value and refill the shuffle ta- ble. iv(j)=idum ran1=min(AM*iy,RNMX) Because users don’t expect endpoint values. return END The routine ran1passes those statistical tests that ran0is known to fail. In fact, we do not know of any statistical test that ran1fails to pass, except when the numberofcalls starts tobecomeontheorderof theperiod m, say >108≈m/20. Forsituationswhenevenlongerrandomsequencesareneeded,L’Ecuyer [6]has given a good way of combining two different sequences with different periods so as to obtain a new sequence whose period is the least common multiple of the two periods. The basic idea is simply to add the two sequences, modulothe modulusofeitherof them (call it m). A trick to avoid an intermediate value that overflows the integerwordsizeistosubtractratherthanadd,andthenaddbacktheconstant m−1 if the result is ≤0, so as to wrap aroundinto the desired interval 0,...,m −1. Notice that it is not necessary that this wrapped subtraction be able to reach all values 0,...,m −1fromeveryvalue of the first sequence. Consider the absurd extreme case where the value subtracted was only between 1 and 10: The resulting sequence would still be no less random than the first sequence by itself. As a practical matter it is only necessary that the second sequence have a range coveringsubstantially all of the range of the first. L’Ecuyer recommends the use of the two generators m 1= 2147483563 (with a1= 40014,q1= 53668,r1= 12211) and m2= 2147483399 (with a2= 40692,q2= 52774,r2= 3791). Both moduli are slightly less than 231. The periods m1−1=2 ×3×7×631×81031and m2−1=2 ×19×31×1019×1789share only the factor 2, so the period of the combined generator is ≈2.3×1018. For present computers, period exhaustion is a practical impossibility. 272 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).OUTPUT RAN1 32iy iv1 iv32 Figure 7.1.1. Shuf fling procedure used in ran1to break up sequential correlations in the Minimal Standard generator. Circled numbers indicate the sequence of events: On each call, the random numberiniyis used to choose a random element in the array iv. That element becomes the output random number, and also is the next iy. Its spot in ivis refilled from the Minimal Standard routine. Combining the two generators breaks up serial correlations to a considerable extent. We nevertheless recommend the additional shuf fle that is implemented in the following routine, ran2. We think that, within the limits of its floating-point precision, ran2providesperfectrandomnumbers;apracticalde finitionof“perfect” isthatwewillpay $1000tothe firstreaderwhoconvincesusotherwise(by findinga statistical test that ran2fails in a nontrivial way, excluding the ordinarylimitations of a machine ’sfloating-point representation). FUNCTION ran2(idum) INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV REAL ran2,AM,EPS,RNMXPARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, * IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211, * IR2=3791,NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS) Long period ( >2×10 18) random number generator of L’Ecuyer with Bays-Durham shuffle and added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint values). Call with iduma negative integer to initialize; thereafter, do not alter idumbetween successive deviates in a sequence. RNMXshould approximate the largest floating value that is less than 1. INTEGER idum2,j,k,iv(NTAB),iy SAVE iv,iy,idum2DATA idum2/123456789/, iv/NTAB*0/, iy/0/if (idum.le.0) then Initialize. idum=max(-idum,1) Be sure to prevent idum =0. idum2=idumdo 11j=NTAB+8,1,-1 Load the shuffle table (after 8 warm-ups). k=idum/IQ1 7.1UniformDeviates 273Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 if (j.le.NTAB) iv(j)=idum enddo 11 iy=iv(1) endif k=idum/IQ1 Start here when not initializing. idum=IA1*(idum-k*IQ1)-k*IR1 Compute idum=mod(IA1*idum,IM1) without over- flows by Schrage’s method. if (idum.lt.0) idum=idum+IM1 k=idum2/IQ2 idum2=IA2*(idum2-k*IQ2)-k*IR2 Compute idum2=mod(IA2*idum2,IM2) likewise. if (idum2.lt.0) idum2=idum2+IM2j=1+iy/NDIV Will be in the range 1:NTAB . iy=iv(j)-idum2 Here idumis shuffled, idumandidum2are com- bined to generate output. iv(j)=idum if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) Because users don’t expect endpoint values. returnEND L’Ecuyer[6]lists additional short generators that can be combined into longer ones, includinggeneratorsthat can be implementedin 16-bitintegerarithmetic. Finally, we give you Knuth ’s suggestion [4]for a portable routine, which we have translated to the present conventions as ran3. This is not based on the linear congruential method at all, but rather on a subtractive method (see also [5]). One might hope that its weaknesses, if any, are therefore of a highly different character from the weaknesses, if any, of ran1above. If you ever suspect trouble with one routine, it is a good idea to try the other in the same application. ran3has one nice feature: if your machine is poor on integer arithmetic (i.e., is limited to 16-bitintegers),substitutionofthethree “commented ”linesfortheonesdirectlypreceding them will render the routine entirely floating-point. FUNCTION ran3(idum) Returns a uniform random deviate between 0.0and 1.0.S e t idumto any negative value to initialize or reinitialize the sequence. INTEGER idumINTEGER MBIG,MSEED,MZ C REAL MBIG,MSEED,MZ REAL ran3,FACPARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG) C PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG) According to Knuth, any large mbig, and any smaller (but still large) mseed can be sub- stituted for the above values. INTEGER i,iff,ii,inext,inextp,k INTEGER mj,mk,ma(55) The value 55 is special and should not be modified; see Knuth. C REAL mj,mk,ma(55) SAVE iff,inext,inextp,maDATA iff /0/ if(idum.lt.0.or.iff.eq.0)then Initialization. iff=1mj=abs(MSEED-abs(idum)) Initialize ma(55) using the seed idumand the large num- bermseed. mj=mod(mj,MBIG) ma(55)=mjmk=1do 11i=1,54 Now initialize the rest of the table, ii=mod(21*i,55) in a slightly random order, ma(ii)=mk with numbers that are not especially random. mk=mj-mk if(mk.lt.MZ)mk=mk+MBIG 274 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).mj=ma(ii) enddo 11 do13k=1,4 We randomize them by “warming up the generator.” do12i=1,55 ma(i)=ma(i)-ma(1+mod(i+30,55)) if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG enddo 12 enddo 13 inext=0 Prepare indices for our first generated number. inextp=31 The constant 31 is special; see Knuth. idum=1 endifinext=inext+1 Here is where we start, except on initialization. Increment inext, wrapping around 56 to 1. if(inext.eq.56)inext=1 inextp=inextp+1 Ditto for inextp . if(inextp.eq.56)inextp=1 mj=ma(inext)-ma(inextp) Now generate a new random number subtractively. if(mj.lt.MZ)mj=mj+MBIG Be sure that it is in range. ma(inext)=mj Store it, ran3=mj*FAC and output the derived uniform deviate. return END Quick andDirty Generators Onesometimeswouldlikea “quickanddirty ”generatortoembedinaprogram,perhaps taking only one or two lines of code, just to somewhat randomize things. One might wish to process data from an experiment not always in exactly the same order, for example, so thatthefirst output is more “typical”than might otherwise be the case. For this kind of application, all we really need is a list of “good”choices for m,a, and cin equation (7.1.1). If we don ’t need a period longer than 10 4to106, say, we can keep the value of (m−1)a+csmall enough to avoid over flows that would otherwise mandate the extra complexity of Schrage ’s method (above). We can thus easily embed in our programs jran=mod(jran*ia+ic,im) ran=float(jran)/float(im) whenever we want a quick and dirty uniform deviate, or jran=mod(jran*ia+ic,im) j=jlo+((jhi-jlo+1)*jran)/im whenever we want an integer between jloandjhi, inclusive. (In both cases jranwas once initialized to any seed value between 0 and im-1.) Be sure to remember, however, that when imis small, the kth root of it, which is the number of planes in k-space, is even smaller! So a quick and dirty generator should never be used to select points in k-space with k> 1. Withthesecaveats,some “good”choicesfortheconstantsaregivenintheaccompanying table. These constants (i) give a period of maximal length im, and, more important, (ii) pass Knuth’s“spectral test ”for dimensions 2, 3, 4, 5, and 6. The increment icis a prime, close to the value (1 2−1 6√ 3)im;actually almost any value of icthatis relativelyprime to imwilldo just as well, but there is some “lore”favoring this choice (see [4], p. 84). 7.1UniformDeviates 275Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).ConstantsforQuickandDirtyRandomNumberGenerators overflow at im ia ic 6075 106 1283 220 7875 211 1663 221 7875 421 1663 222 6075 1366 12836655 936 1399 11979 430 2531 2 23 14406 967 304129282 419 617353125 171 11213 2 24 12960 1741 273114000 1541 295721870 1291 462131104 625 6571 139968 205 29573 2 25 29282 1255 6173 81000 421 17117 134456 281 28411 226overflow at im ia ic 86436 1093 18257 121500 1021 25673259200 421 54773 2 27 117128 1277 24749121500 2041 25673312500 741 66037 2 28 145800 3661 30809175000 2661 36979233280 1861 49297244944 1597 51749 2 29 139968 3877 29573214326 3613 45289714025 1366 150889 2 30 134456 8121 28411259200 7141 54773 2 31 233280 9301 49297 714025 4096 150889 232 AnEven Quicker and DirtierGenerator Many FORTRAN compilerscanbeabusedinsuchawaythattheywillmultiplytwo32-bit integersignoringany resultingoverflow . Insuchcases,onmanymachines, thevalue returned is predictably the low-order 32 bits of the true 64-bit product. ( Ccompilers, incidentally, can do this without the requirement of abuse —it is guaranteed behavior for so-called unsigned long int integers. On VMS VAXes, the necessary FORTRAN command is FORTRAN/CHECK=NOOVERFLOW .) If we now choose m=232, the“mod”in equation (7.1.1) is free, and we have simply Ij+1=aI j+c (7.1.6 ) Knuth suggests a= 1664525 as a suitable multiplier for this value of m. H.W. Lewis has conducted extensive testsofthisvalue of awith c= 1013904223 , which isaprime close to(√ 5−2)m. The resulting in-line generator (we will call it ranqd1) is simply idum=1664525*idum+1013904223 This is about as good as any 32-bit linear congruential generator, entirely adequate for many uses. And, with only a single multiply and add, it is veryfast. To check whether your compiler and machine have the desired over flow proper- ties, see if you can generate the following sequence of 32-bit values (given here inhex): 00000000, 3C6EF35F, 47502932, D1CCF6E9, AAF95334, 6252E503, 9F2EC686,57FE6C2D, A3D95FA8, 81FDBEE7, 94F0AF1A, CBF633B1. Ifyouneed floating-pointvaluesinsteadof32-bitintegers,andwanttoavoidadivideby floating-point 2 32,adirtytrickistomaskinanexponentthatmakesthevalueliebetween1and 2, then subtract 1.0. The resulting in-line generator (callit ranqd2) willlook something like 276 Chapter7. RandomNumbersSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).INTEGER idum,itemp,jflone,jflmsk REAL ftemp EQUIVALENCE (itemp,ftemp) DATA jflone /Z’3F800000’/, jflmsk /Z’007FFFFF’/... idum=1664525*idum+1013904223 itemp=ior(jflone,iand(jflmsk,idum))ran=ftemp-1.0 The hex constants 3F800000 and 007FFFFF are the appropriate ones for computers using the IEEE representation for 32-bit floating-point numbers (e.g., IBM PCs and most UNIX workstations). For DEC VAXes, the correct hex constants are, respectively, 00004080 andFFFF007F. Notice that the IEEE mask results in the floating-point number being constructed out of the 23 low-order bits of the integer, which is not ideal. Also notice that your compilermayrequireadifferentnotationforhexconstants,e.g., x’3f800000’ ,’3F800000’X ,ore ven 16#3F800000 . (Your authors have tried very hard to make almost all of the material in this bookmachineandcompilerindependent —indeed,evenprogramminglanguageindependent. This subsection is a rare aberration. Forgive us. Once in a great while the temptation tobereally dirty is just irresistible.) Relative Timingsand Recommendations Timings are inevitably machine dependent. Nevertheless the following table is indicative of the relativetimings, for typical machines, of the various uniform generatorsdiscussedinthissection,plus ran4from§7.5. Smallervaluesinthetable indicate faster generators. The generators ranqd1andranqd2refer to the “quick and dirty”generators immediately above. Generator RelativeExecutionTime ran0 ≡1.0 ran1 ≈1.3 ran2 ≈2.0 ran3 ≈0.6 ranqd1 ≈0.10 ranqd2 ≈0.25 ran4 ≈4.0 On balance, we recommend ran1for general use. It is portable, based on ParkandMiller ’sMinimalStandardgeneratorwithanadditionalshuf fle,andhas no known (to us) flaws other than period exhaustion. If you are generating more than 100,000,000 random numbers in a single calculation (that is, more than about 5% of ran1’s period), we recommend the use ofran2, with its much longer period. Knuth’ssubtractiveroutine ran3seemstobethetimingwinneramongportable routines. Unfortunately the subtractive method is not so well studied, and not astandard. Weliketokeep ran3inreservefora “secondopinion, ”substitutingitwhen wesuspectanothergeneratorofintroducingunwantedcorrelationsintoacalculation. The routine ran4generates extremely good random deviates, and has some other nice properties, but it is slow. See §7.5 for discussion. 7.2TransformationMethod: ExponentialandNormalDeviates 277Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Finally, the quick and dirty in-line generators ranqd1andranqd2are very fast, but they are machine dependent, nonportable, and at best only as good as a32-bitlinearcongruentialgeneratoreveris —in ourviewnot goodenoughinmany situations. We woulduse these onlyin veryspecial cases, wherespeed is critical. CITED REFERENCES AND FURTHER READING: Park, S.K., and Miller, K.W. 1988, Communications of the ACM , vol. 31, pp. 1192–1201. [1] Schrage, L. 1979, ACM Transactions on Mathematical Software , vol. 5, pp. 132–138. [2] Bratley, P., Fox, B.L., and Schrage, E.L. 1983, A Guide to Simulation (New York: Springer- Verlag). [3] Knuth,D.E.1981, SeminumericalAlgorithms ,2nded.,vol.2of TheArtofComputerProgramming (Reading, MA: Addison-Wesley), §§3.2–3.3. [4] Kahaner,D.,Moler,C.,andNash,S.1989, NumericalMethods andSoftware (EnglewoodCliffs, NJ: Prentice Hall), Chapter 10. [5] L’Ecuyer, P. 1988, Communications of the ACM , vol. 31, pp. 742–774. [6] Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 10. 7.2 Transformation Method: Exponential and Normal Deviates In the previous section, we learned how to generate random deviates with a uniform probability distribution, so that the probability of generating a number between xandx+dx, denoted p(x)dx, is given by p(x)dx=/braceleftBigdx 0<x< 1 0otherwise(7.2.1 ) The probability distribution p(x)is of course normalized, so that /integraldisplay∞ −∞p(x)dx=1 ( 7.2.2 ) Nowsupposethatwegenerateauniformdeviate xandthentakesomeprescribed functionofit, y(x). Theprobabilitydistributionof y,denoted p(y)dy,isdetermined by the fundamental transformation law of probabilities, which is simply |p(y)dy|=|p(x)dx| (7.2.3 ) or p(y)=p(x)/vextendsingle/vextendsingle/vextendsingle/vextendsingledx dy/vextendsingle/vextendsingle/vextendsingle/vextendsingle(7.2.4 )