Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 67
Текст из файла (страница 67)
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.Correlation in k-space is not the only weakness of linear congruential generators.Such generators often have their low-order (least significant) bits much less randomthan their high-order bits.
If you want to generate a random integer between 1 and10, you should always do it using high-order bits, as inj=1+(int) (10.0*rand()/(RAND_MAX+1.0));and never by anything resemblingj=1+(rand() % 10);(which uses lower-order bits). Similarly you should never try to take apart a“rand()” number into several supposedly random pieces. Instead use separatecalls for every piece.278Chapter 7.Random NumbersPortable Random Number GeneratorsPark 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 algorithmIj+1 = aIj(mod m)(7.1.2)can be as good as any of the more general linear congruential generators that havec = 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 that+az 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:7.1 Uniform Deviates#define#define#define#define#define#define279IA 16807IM 2147483647AM (1.0/IM)IQ 127773IR 2836MASK 123459876float 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;*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 remove280Chapter 7.Random Numbersthem, 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 remove low-order serial correlations.
A random deviatederived from the jth value in the sequence, Ij , is output not on the jth call, but ratheron a randomized later call, j + 32 on average. The shuffling algorithm is due to Baysand Durham as described in Knuth [4], and is illustrated in Figure 7.1.1.#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 of2817.1 Uniform Deviatesiy1iv0RAN3OUTPUT2iv31Figure 7.1.1. Shuffling procedure used in ran1 to break up sequential correlations in the MinimalStandard generator.
Circled numbers indicate the sequence of events: On each call, the random numberin iy is used to choose a random element in the array iv. That element becomes the output randomnumber, and also is the next iy. Its spot in iv is refilled from the Minimal Standard routine.either of them (call it m). A trick to avoid an intermediate value that overflows theinteger wordsize is to subtract rather than add, and then add back the constant m − 1if the result is ≤ 0, so as to wrap around into the desired interval 0, .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
