#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)

#include "def.h"

float runif(long *idum){
      int j;
      long k;
      static long iy=0;
      static long iv[NTAB];
      float temp;
      if (*idum <= 0 || !iy) { //Initialize.
         if (-(*idum) < 1) *idum=1; //Be sure to prevent idum = 0.
         else *idum = -(*idum);
      
         for (j=NTAB+7;j>=0;j--) { //Load the shue table (after 8 warm-ups).
              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; //Start here when not initializing.
         *idum=IA*(*idum-k*IQ)-IR*k; //Compute idum=(IA*idum) % IM without over
         if(*idum < 0) *idum += IM; //flows by Schrage's method.
         j=iy/NDIV; //Will be in the range 0..NTAB-1.
         iy=iv[j]; //Output previously stored value and rell the
         iv[j] = *idum; //shue table.
         if ((temp=AM*iy) > RNMX) return RNMX; //Because users don't expect endpoint values.
         else return temp;      
}

float rnorm(long *idum){
      float ran1(long *idum);
      static int iset=0;
      static float gset;
      float fac,rsq,v1,v2;
      
      if (iset == 0) { //We don't have an extra deviate handy, so
      do {
         v1=2.0*runif(idum)-1.0; //pick two uniform numbers in the square extending
         v2=2.0*runif(idum)-1.0; //from -1 to +1 in each direction,
         rsq=v1*v1+v2*v2; //see if they are in the unit circle,
      } while (rsq >= 1.0 || rsq == 0.0); //and if they are not, try again.

      fac=sqrt(-2.0*log(rsq)/rsq);
      gset=v1*fac;
      iset=1; //Set flag.
      return v2*fac;
      } 
      else { //We have an extra deviate handy,
        iset=0; //so unset the flag,
        return gset; //and return it.
      }
      
}