c7-1 (779515), страница 2
Текст из файла (страница 2)
One infamous such routine, RANDU, with a = 65539 and m = 231 ,was widespread on IBM mainframe computers for many years, and widely copiedonto other systems [1]. One of us recalls producing a “random” plot with only 11planes, and being told by his computer center’s programming consultant that hehad misused the random number generator: “We guarantee that each number israndom individually, but we don’t guarantee that more than one of them is random.”Figure that out.278Chapter 7.Random NumbersPortable Random Number GeneratorsIj+1 = aIj(mod m)(7.1.2)can be as good as any of the more general linear congruential generators that havec 6= 0 (equation 7.1.1) — if the multiplier a and modulus m are chosen exquisitelycarefully.
Park and Miller propose a “Minimal Standard” generator based on thechoicesm = 231 − 1 = 2147483647a = 75 = 16807(7.1.3)First proposed by Lewis, Goodman, and Miller in 1969, this generator has insubsequent years passed all new theoretical tests, and (perhaps more importantly)has accumulated a large amount of successful use. Park and Miller do not claim thatthe generator is “perfect” (we will see below that it is not), but only that it is a goodminimal standard against which other generators should be judged.It is not possible to implement equations (7.1.2) and (7.1.3) directly in ahigh-level language, since the product of a and m − 1 exceeds the maximum valuefor a 32-bit integer.
Assembly language implementation using a 64-bit productregister is straightforward, but not portable from machine to machine. A trickdue to Schrage [2,3] for multiplying two 32-bit integers modulo a 32-bit constant,without using any intermediates larger than 32 bits (including a sign bit) is thereforeextremely interesting: It allows the Minimal Standard generator to be implementedin essentially any programming language on essentially any machine.Schrage’s algorithm is based on an approximate factorization of m,m = aq + r,i.e.,q = [m/a], r = m mod a(7.1.4)with square brackets denoting integer part. If r is small, specifically r < q, and0 < z < m − 1, it can be shown that both a(z mod q) and r[z/q] lie in the range0, .
. . , m − 1, and thataz mod m =a(z mod q) − r[z/q]a(z mod q) − r[z/q] + mif it is ≥ 0,otherwise(7.1.5)The application of Schrage’s algorithm to the constants (7.1.3) uses the valuesq = 127773 and r = 2836.Here is an implementation of the Minimal Standard generator:Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Park and Miller [1] have surveyed a large number of random number generatorsthat have been used over the last 30 years or more. Along with a good theoreticalreview, they present an anecdotal sampling of a number of inadequate generatorsthat have come into widespread use.
The historical record is nothing if not appalling.There is good evidence, both theoretical and empirical, that the simple multiplicative congruential algorithm7.1 Uniform Deviates#define#define#define#define#define#define279IA 16807IM 2147483647AM (1.0/IM)IQ 127773IR 2836MASK 123459876*idum ^= MASK;k=(*idum)/IQ;*idum=IA*(*idum-k*IQ)-IR*k;if (*idum < 0) *idum += IM;ans=AM*(*idum);*idum ^= MASK;return ans;XORing with MASK allows use of zero and othersimple bit patterns for idum.Compute idum=(IA*idum) % IM without overflows by Schrage’s method.Convert idum to a floating result.Unmask before return.}The period of ran0 is 231 − 2 ≈ 2.1 × 109 .
A peculiarity of generators ofthe form (7.1.2) is that the value 0 must never be allowed as the initial seed — itperpetuates itself — and it never occurs for any nonzero initial seed. Experience hasshown that users always manage to call random number generators with the seedidum=0. That is why ran0 performs its exclusive-or with an arbitrary constant bothon entry and exit. If you are the first user in history to be proof against human error,you can remove the two lines with the ∧ operation.Park and Miller discuss two other multipliers a that can be used with the samem = 231 − 1. These are a = 48271 (with q = 44488 and r = 3399) and a = 69621(with q = 30845 and r = 23902).
These can be substituted in the routine ran0if desired; they may be slightly superior to Lewis et al.’s longer-tested values. Novalues other than these should be used.The routine ran0 is a Minimal Standard, satisfactory for the majority ofapplications, but we do not recommend it as the final word on random numbergenerators. Our reason is precisely the simplicity of the Minimal Standard. It isnot hard to think of situations where successive random numbers might be usedin a way that accidentally conflicts with the generation algorithm. For example,since successive numbers differ by a multiple of only 1.6 × 104 out of a modulus ofmore than 2 × 109 , very small random numbers will tend to be followed by smallerthan average values.
One time in 106 , for example, there will be a value < 10−6returned (as there should be), but this will always be followed by a value less thanabout 0.0168. One can easily think of applications involving rare events where thisproperty would lead to wrong results.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 , then the resulting distribution fails the χ2 test when N is greater than afew ×107 , much less than the period m − 2. Since low-order serial correlations havehistorically been such a bugaboo, and since there is a very simple way to removeSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).float ran0(long *idum)“Minimal” random number generator of Park and Miller. Returns a uniform random deviatebetween 0.0 and 1.0. Set or reset idum to any integer value (except the unlikely value MASK)to initialize the sequence; idum must not be altered between calls for successive deviates ina sequence.{long k;float ans;280Chapter 7.Random Numbers#define#define#define#define#define#define#define#define#defineIA 16807IM 2147483647AM (1.0/IM)IQ 127773IR 2836NTAB 32NDIV (1+(IM-1)/NTAB)EPS 1.2e-7RNMX (1.0-EPS)float ran1(long *idum)“Minimal” random number generator of Park and Miller with Bays-Durham shuffle and addedsafeguards.
Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpointvalues). Call with idum a negative integer to initialize; thereafter, do not alter idum betweensuccessive deviates in a sequence. RNMX should approximate the largest floating value that isless than 1.{int j;long k;static long iy=0;static long iv[NTAB];float temp;if (*idum <= 0 || !iy) {if (-(*idum) < 1) *idum=1;else *idum = -(*idum);for (j=NTAB+7;j>=0;j--) {k=(*idum)/IQ;*idum=IA*(*idum-k*IQ)-IR*k;if (*idum < 0) *idum += IM;if (j < NTAB) iv[j] = *idum;}iy=iv[0];}k=(*idum)/IQ;*idum=IA*(*idum-k*IQ)-IR*k;if (*idum < 0) *idum += IM;j=iy/NDIV;iy=iv[j];iv[j] = *idum;if ((temp=AM*iy) > RNMX) return RNMX;else return temp;Initialize.Be sure to prevent idum = 0.Load the shuffle table (after 8 warm-ups).Start here when not initializing.Compute idum=(IA*idum) % IM without overflows by Schrage’s method.Will be in the range 0..NTAB-1.Output previously stored value and refill theshuffle table.Because users don’t expect endpoint values.}The routine ran1 passes those statistical tests that ran0 is known to fail.
Infact, we do not know of any statistical test that ran1 fails to pass, except when thenumber of calls starts to become on the order of the period m, say > 108 ≈ m/20.For situations when even longer random sequences are needed, L’Ecuyer [6] hasgiven a good way of combining two different sequences with different periods soas to obtain a new sequence whose period is the least common multiple of the twoperiods. The basic idea is simply to add the two sequences, modulo the modulus ofSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
