00001 00005 #include <bpm/bpm_nr.h> 00006 00013 double nr_ran1( long *idum ) { 00014 00015 int j; 00016 long k; 00017 static long iy = 0; 00018 static long iv[RAN1_NTAB]; 00019 double temp; 00020 00021 if (*idum <= 0 || !iy) 00022 { 00023 if (-(*idum) < 1) *idum=1; 00024 else *idum = -(*idum); 00025 for (j=RAN1_NTAB+7;j>=0;j--) 00026 { 00027 k=(*idum)/RAN1_IQ; 00028 *idum=RAN1_IA*(*idum-k*RAN1_IQ)-RAN1_IR*k; 00029 if (*idum < 0) *idum += RAN1_IM; 00030 if (j < RAN1_NTAB) iv[j] = *idum; 00031 } 00032 iy=iv[0]; 00033 } 00034 k=(*idum)/RAN1_IQ; 00035 *idum=RAN1_IA*(*idum-k*RAN1_IQ)-RAN1_IR*k; 00036 if (*idum < 0) *idum += RAN1_IM; 00037 j=iy/RAN1_NDIV; 00038 iy=iv[j]; 00039 iv[j]=*idum; 00040 if((temp=RAN1_AM*iy) > RAN1_RNMX) return RAN1_RNMX; 00041 else return temp; 00042 }