c7-1 (779515), страница 3
Текст из файла (страница 3)
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).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 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.2817.1 Uniform Deviatesiy1RAN3OUTPUT2iv31Figure 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, .
. . , m − 1.Notice that it is not necessary that this wrapped subtraction be able to reach allvalues 0, . . . , m − 1 from every value of the first sequence. Consider the absurdextreme case where the value subtracted was only between 1 and 10: The resultingsequence would still be no less random than the first sequence by itself. As apractical 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 twogenerators m1 = 2147483563 (with a1 = 40014, q1 = 53668, r1 = 12211) andm2 = 2147483399 (with a2 = 40692, q2 = 52774, r2 = 3791). Both moduliare slightly less than 231 . The periods m1 − 1 = 2 × 3 × 7 × 631 × 81031 andm2 − 1 = 2 × 19 × 31 × 1019 × 1789 share only the factor 2, so the period ofthe combined generator is ≈ 2.3 × 1018 . For present computers, period exhaustionis a practical impossibility.Combining the two generators breaks up serial correlations to a considerableextent. We nevertheless recommend the additional shuffle that is implemented inthe following routine, ran2. We think that, within the limits of its floating-pointprecision, ran2 provides perfect random numbers; a practical definition of “perfect”is that we will pay $1000 to the first reader who convinces us otherwise (by finding astatistical test that ran2 fails in a nontrivial way, excluding the ordinary limitationsof a machine’s floating-point representation).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).iv0282Chapter 7.IM1 2147483563IM2 2147483399AM (1.0/IM1)IMM1 (IM1-1)IA1 40014IA2 40692IQ1 53668IQ2 52774IR1 12211IR2 3791NTAB 32NDIV (1+IMM1/NTAB)EPS 1.2e-7RNMX (1.0-EPS)float ran2(long *idum)Long period (> 2 × 1018) random number generator of L’Ecuyer with Bays-Durham shuffleand added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive ofthe endpoint values). Call with idum a negative integer to initialize; thereafter, do not alteridum between successive deviates in a sequence.
RNMX should approximate the largest floatingvalue that is less than 1.{int j;long k;static long idum2=123456789;static long iy=0;static long iv[NTAB];float temp;if (*idum <= 0) {Initialize.if (-(*idum) < 1) *idum=1;Be sure to prevent idum = 0.else *idum = -(*idum);idum2=(*idum);for (j=NTAB+7;j>=0;j--) {Load the shuffle table (after 8 warm-ups).k=(*idum)/IQ1;*idum=IA1*(*idum-k*IQ1)-k*IR1;if (*idum < 0) *idum += IM1;if (j < NTAB) iv[j] = *idum;}iy=iv[0];}k=(*idum)/IQ1;Start here when not initializing.*idum=IA1*(*idum-k*IQ1)-k*IR1;Compute idum=(IA1*idum) % IM1 withoutif (*idum < 0) *idum += IM1;overflows by Schrage’s method.k=idum2/IQ2;idum2=IA2*(idum2-k*IQ2)-k*IR2;Compute idum2=(IA2*idum) % IM2 likewise.if (idum2 < 0) idum2 += IM2;j=iy/NDIV;Will be in the range 0..NTAB-1.iy=iv[j]-idum2;Here idum is shuffled, idum and idum2 areiv[j] = *idum;combined to generate output.if (iy < 1) iy += IMM1;if ((temp=AM*iy) > RNMX) return RNMX; Because users don’t expect endpoint values.else return temp;}L’Ecuyer [6] lists additional short generators that can be combined into longerones, including generators that can be implemented in 16-bit integer arithmetic.Finally, we give you Knuth’s suggestion [4] for a portable routine, which wehave translated to the present conventions as ran3.
This is not based on the linearcongruential method at all, but rather on a subtractive method (see also [5]). Onemight hope that its weaknesses, if any, are therefore of a highly different characterSample 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).#define#define#define#define#define#define#define#define#define#define#define#define#define#defineRandom Numbers7.1 Uniform Deviates283#include <stdlib.h>Change to math.h in K&R C.#define MBIG 1000000000#define MSEED 161803398#define MZ 0#define FAC (1.0/MBIG)According to Knuth, any large MBIG, and any smaller (but still large) MSEED can be substitutedfor the above values.float ran3(long *idum)Returns a uniform random deviate between 0.0 and 1.0.
Set idum to any negative value toinitialize or reinitialize the sequence.{static int inext,inextp;static long ma[56];The value 56 (range ma[1..55]) is special andstatic int iff=0;should not be modified; see Knuth.long mj,mk;int i,ii,k;if (*idum < 0 || iff == 0) {Initialization.iff=1;mj=labs(MSEED-labs(*idum));Initialize ma[55] using the seed idum and themj %= MBIG;large number MSEED.ma[55]=mj;mk=1;for (i=1;i<=54;i++) {Now initialize the rest of the table,ii=(21*i) % 55;in a slightly random order,ma[ii]=mk;with numbers that are not especially random.mk=mj-mk;if (mk < MZ) mk += MBIG;mj=ma[ii];}for (k=1;k<=4;k++)We randomize them by “warming up the generfor (i=1;i<=55;i++) {ator.”ma[i] -= ma[1+(i+30) % 55];if (ma[i] < MZ) ma[i] += MBIG;}inext=0;Prepare indices for our first generated number.inextp=31;The constant 31 is special; see Knuth.*idum=1;}Here is where we start, except on initialization.if (++inext == 56) inext=1;Increment inext and inextp, wrapping aroundif (++inextp == 56) inextp=1;56 to 1.mj=ma[inext]-ma[inextp];Generate a new random number subtractively.if (mj < MZ) mj += MBIG;Be sure that it is in range.ma[inext]=mj;Store it,return mj*FAC;and output the derived uniform deviate.}Quick and Dirty GeneratorsOne sometimes would like a “quick and dirty” generator to embed in a program, perhapstaking only one or two lines of code, just to somewhat randomize things.
One might wish toSample 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).from the weaknesses, if any, of ran1 above. If you ever suspect trouble withone routine, it is a good idea to try the other in the same application.
ran3 hasone nice feature: if your machine is poor on integer arithmetic (i.e., is limitedto 16-bit integers), you can declare mj, mk, and ma[] as float, define mbig andmseed as 4000000 and 1618033, respectively, and the routine will be renderedentirely floating-point.284Chapter 7.Random Numbersprocess data from an experiment not always in exactly the same order, for example, so thatthe first 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, andc in equation (7.1.1). If we don’t need a period longer than 104 to 106 , say, we can keep thevalue of (m − 1)a + c small enough to avoid overflows that would otherwise mandate theextra complexity of Schrage’s method (above).
We can thus easily embed in our programswhenever we want a quick and dirty uniform deviate, orjran=(jran*ia+ic) % im;j=jlo+((jhi-jlo+1)*jran)/im;whenever we want an integer between jlo and jhi, inclusive. (In both cases jran was onceinitialized to any seed value between 0 and im-1.)Be sure to remember, however, that when im is small, the kth root of it, which is thenumber of planes in k-space, is even smaller! So a quick and dirty generator should neverbe used to select points in k-space with k > 1.With these caveats, some “good” choices for the constants are given in the accompanyingtable.
These constants (i) give a period of maximal length im, and, more important, (ii) passKnuth’s “spectral√test” for dimensions 2, 3, 4, 5, and 6. The increment ic is a prime, close tothe value ( 12 − 16 3)im; actually almost any value of ic that is relatively prime to im will dojust as well, but there is some “lore” favoring this choice (see [4], p. 84).An Even Quicker GeneratorIn C, if you multiply two unsigned long int integers on a machine with a 32-bit longinteger representation, the value returned is the low-order 32 bits of the true 64-bit product.
Ifwe now choose m = 232 , the “mod” in equation (7.1.1) is free, and we have simplyIj+1 = aIj + c(7.1.6)Knuth suggests a = 1664525 as a suitable multiplier for this value of m. H.W. Lewishas √conducted extensive tests of this value of a with c = 1013904223, which is a prime closeto ( 5 − 2)m.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
