1 /**
2   @file
3   @ingroup nr
4 */
5 #include <bpm/bpm_nr.h>
6 
7 /**
8    Random number generator as nicked from numerical recipes, c7.1, p280
9 
10    @param idum random seed, note that the global seed is set by bpm_rseed
11    @return random number between 0 and 1
12 */
nr_ran1(long * idum)13 double nr_ran1( long *idum ) {
14 
15   int j;
16   long k;
17   static long iy = 0;
18   static long iv[RAN1_NTAB];
19   double temp;
20 
21   if (*idum <= 0 || !iy)
22     {
23       if (-(*idum) < 1) *idum=1;
24       else *idum = -(*idum);
25       for (j=RAN1_NTAB+7;j>=0;j--)
26 	{
27 	  k=(*idum)/RAN1_IQ;
28 	  *idum=RAN1_IA*(*idum-k*RAN1_IQ)-RAN1_IR*k;
29 	  if (*idum < 0) *idum += RAN1_IM;
30 	  if (j < RAN1_NTAB) iv[j] = *idum;
31 	}
32       iy=iv[0];
33     }
34   k=(*idum)/RAN1_IQ;
35   *idum=RAN1_IA*(*idum-k*RAN1_IQ)-RAN1_IR*k;
36   if (*idum < 0) *idum += RAN1_IM;
37   j=iy/RAN1_NDIV;
38   iy=iv[j];
39   iv[j]=*idum;
40   if((temp=RAN1_AM*iy) > RAN1_RNMX) return RAN1_RNMX;
41   else return temp;
42 }
43