/* Random Number Generators from Numerical Recipes in C v2.0 */ /* (C) Copr. 1986-92 Numerical Recipes Software `2\$m.1.9-153. */ #include #define IA 16807 #define IM 2147483647 #define AM (1.0/IM) #define IQ 127773 #define IR 2836 #define NTAB 32 #define NDIV (1+(IM-1)/NTAB) #define EPS 1.2e-7 #define RNMX (1.0-EPS) float ran1(long *idum) { 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; temp = AM*iy; if (temp > RNMX) return RNMX; else return temp; } #undef IA #undef IM #undef AM #undef IQ #undef IR #undef NTAB #undef NDIV #undef EPS #undef RNMX float gasdev(long *idum) { float ran1(long *idum); static int iset=0; static double gset; double fac,rsq,v1,v2; double temp; if (iset == 0) { do { v1=2.0*(double)ran1(idum)-1.0; v2=2.0*(double)ran1(idum)-1.0; rsq=v1*v1+v2*v2; } while ((rsq >= 1.0) || (rsq == 0.0)); temp = log( rsq ) / rsq; fac = sqrt( -2.0 * temp ); gset=v1*fac; iset=1; return( (float)v2*fac ); } else { iset=0; return( (float)gset ); } }