1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
|
#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)
/*
Long period (> 2 X 10^18 ) 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.
*/
float ran2(long *idum)
{
int j;
long k;
static long idum2=123456789;
static long iy=0;
static long iv[NTAB];
float temp;
/* Initialize. */
if (*idum <= 0)
{
/* Be sure to prevent idum = 0. */
if (-(*idum) < 1) *idum=1;
else *idum = -(*idum);
idum2=(*idum);
/* Load the shuffle table (after 8 warm-ups). */
for(j=NTAB+7;j>=0;j--)
{
k=(*idum)/IQ1;
*idum=IA1*(*idum-k*IQ1)-k*IR1;
if (*idum < 0) *idum += IM1;
if (j < NTAB) iv[j] = *idum;
}
iy=iv[0];
}
/* Start here when not initializing. */
/* Compute idum=(IA1*idum) % IM1 without overflows by Schrage's method. */
k=(*idum)/IQ1;
*idum=IA1*(*idum-k*IQ1)-k*IR1;
if(*idum < 0) *idum += IM1;
/* Compute idum2=(IA2*idum) % IM2 likewise. */
k=idum2/IQ2;
idum2=IA2*(idum2-k*IQ2)-k*IR2;
if (idum2 < 0) idum2 += IM2;
/* Will be in the range 0..NTAB-1. */
j=iy/NDIV;
/* Here idum is shuffled, idum and idum2 are combined to generate output. */
iy=iv[j]-idum2;
iv[j] = *idum;
if (iy < 1) iy += IMM1;
/* Because users don't expect endpoint values. */
if((temp=AM*iy) > RNMX) return RNMX;
else return temp;
} |
Partager