Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 68
Текст из файла (страница 68)
. . , 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).282Chapter 7.#define#define#define#define#define#define#define#define#define#define#define#define#define#defineRandom NumbersIM1 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 character7.1 Uniform Deviates283from 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.#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 to284Chapter 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 programsunsigned long jran,ia,ic,im;float ran;...jran=(jran*ia+ic) % im;ran=(float) jran / (float) im;whenever 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. The resulting in-line generator (we will call it ranqd1) is simplyunsigned long idum;...idum = 1664525L*idum + 1013904223L;This is about as good as any 32-bit linear congruential generator, entirely adequate for manyuses. And, with only a single multiply and add, it is very fast.To check whether your machine has the desired integer properties, see if you cangenerate the following sequence of 32-bit values (given here in hex): 00000000, 3C6EF35F,47502932, D1CCF6E9, AAF95334, 6252E503, 9F2EC686, 57FE6C2D, A3D95FA8, 81FDBEE7, 94F0AF1A, CBF633B1.If you need floating-point values instead of 32-bit integers, and want to avoid a divide byfloating-point 232 , a dirty trick is to mask in an exponent that makes the value lie between 1 and2, then subtract 1.0.
The resulting in-line generator (call it ranqd2) will look something like2857.1 Uniform DeviatesConstants for Quick and Dirty Random Number Generatorsoverflow at2202212222232imiaic6075106128378752111663787542116636075 13666655 93611979 430128313992531144062928253125overflow at967 3041419 6173171 112132272282292422522612960 1741 273114000 1541 295721870 1291 462131104 625 6571139968 205 2957329282 1255 617381000 421 17117134456 281 28411230231232iaic86436 1093121500 1021259200 421im182572567354773117128 1277121500 2041312500 741247492567366037145800175000233280244944308093697949297517493661266118611597139968 3877 29573214326 3613 45289714025 1366 150889134456 8121259200 71412841154773233280 9301 49297714025 4096 150889unsigned long idum,itemp;float rand;#ifdef vaxstatic unsigned long jflone = 0x00004080;static unsigned long jflmsk = 0xffff007f;#elsestatic unsigned long jflone = 0x3f800000;static unsigned long jflmsk = 0x007fffff;#endif...idum = 1664525L*idum + 1013904223L;itemp = jflone | (jflmsk & idum);rand = (*(float *)&itemp)-1.0;The hex constants 3F800000 and 007FFFFF are the appropriate ones for computers usingthe IEEE representation for 32-bit floating-point numbers (e.g., IBM PCs and most UNIXworkstations).
For DEC VAXes, the correct hex constants are, respectively, 00004080 andFFFF007F. Notice that the IEEE mask results in the floating-point number being constructedout of the 23 low-order bits of the integer, which is not ideal. (Your authors have triedvery hard to make almost all of the material in this book machine and compiler independent— indeed, even programming language independent. This subsection is a rare aberration.Forgive us. Once in a great while the temptation to be really dirty is just irresistible.)Relative Timings and RecommendationsTimings are inevitably machine dependent.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
