282 Chapter 7.Random Numbers #define IM12147483563 #def1neIN22147483399 #define AM (1.0/IM1) #define IMM1 (IM1-1) #define IA1 40014 #define IA2 40692 #define IQ1 53668 #define IQ2 52774 #define IR1 12211 #define IR2 3791 #define NTAB 32 #define NDIV (1+IMM1/NTAB) http://www.nr. #define EPS 1.2e-7 #define RNMX (1.0-EPS) 83 float ran2(long *idum) granted for 19881992 Long period (>2x 1018)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 1800 the endpoint values).Call with idum a negative integer to initialize:thereafter,do not alter idum between successive deviates in a sequence.RNMX should approximate the largest floating value that is less than 1. from NUMERICAL RECIPESI int ji long ki static long idum2=123456789; (Nort static long iy=0; static long iv[NTAB]; America server computer, to make e University Press. THE float temp; one paper ART if (*idum <=0){ Initialize. if (-(*idum)<1)*idum=1; Be sure to prevent idum =0. 9 else *idum =-(*idum); idum2=(*idum); for (j=NTAB+7;j>=0;j--){ Load the shuffle table (after 8 warm-ups). SCIENTIFIC k=(*idum)/IQ1; *idum=IA1*(*idum-k*IQ1)-k*IR1; to dir if (*idum 0)*idum +IM1; if (j NTAB)iv[j]tidum; 188 iy=iv[o]; 1920 COMPUTING(ISBN k=(*idum)/IQ1; Start here when not initializing. *idum=IA1*(*idum-k*IQ1)-k*IR1; Compute idum=(IA1*idum)%IM1 without 0621 if (*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. Numerical Recipes -43108 if (idum2 0)idum2 +IM2; i=iy/NDIV; Will be in the range 0..NTAB-1. iy=iv[j]-idum2; Here idum is shuffled,idum and idum2 are (outside iv[j]*idum; combined to generate output. if (iy 1)iy +IMM1; North Software. if ((temp=AM*iy)>RNMX)return RNMX; Because users don't expect endpoint values. else return temp; Ame visit website machine L'Ecuyer [6]lists additional short generators that can be combined into longer ones,including generators that can be implemented in 16-bit integer arithmetic. 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 character282 Chapter 7. Random Numbers Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #define IM1 2147483563 #define IM2 2147483399 #define AM (1.0/IM1) #define IMM1 (IM1-1) #define IA1 40014 #define IA2 40692 #define IQ1 53668 #define IQ2 52774 #define IR1 12211 #define IR2 3791 #define NTAB 32 #define NDIV (1+IMM1/NTAB) #define EPS 1.2e-7 #define RNMX (1.0-EPS) float ran2(long *idum) Long period (> 2 × 1018) 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 idum a negative integer to initialize; thereafter, do not alter idum between successive deviates in a sequence. RNMX should approximate the largest floating value 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 without if (*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 are iv[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 longer ones, including generators that can be implemented in 16-bit integer arithmetic. 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