/***************************************************************
*    "Minimal" random number generator of Park and Miller 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.    */


#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;
  if ((temp=AM*iy) > RNMX) return RNMX;
  else return temp;
}
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX
/* After Press et al. (1992), Numerical Recipes for C */


/**************************************************************
 * This program performs a simple histogram analysis of ten million calls
 * to RAN1.  If RAN1 is correctly yielding a uniformly distributed value
 * between 0.0 and slightly less than 1.0, then the counts in each bin
 * should be approximately 100,000, plus or minus a few tenths of a
 * percent.  For the last bin on the last line of output, the count
 * should be zero.  If you see anything wrong with your output, please 
 * let me know and I'll try to help figure out why this doesn't work on 
 * your system.
 *
 * Once you have confirmed that RAN1 appears to be working correctly on
 * your system, then you can comment out or delete PROGRAM TEST and copy
 * RAN1 directly into the source code module for your Monte Carlo program. */

main()
{
  int icount[101] ;
  
  float x;
  long idum ;
  int i,j,k,ibin ;
  
  /* clear histogram */
  for (i=0; i<101; i++) icount[i] = 0 ;
  
  /* initialize with negative number */
  idum = -1  ;
  x = ran1(&idum) ;
  
  /* tabulate values returned by many calls */
  for (i=0; i<10000000; i++) {
    x = ran1(&idum) ;
    ibin = x*100 ;
    icount[ibin]++ ;
  }
  
  k = 0 ;
  for (i=0; i<10; i++) {
    for (j=0; j<10; j++) {
      printf("%8d",icount[k]) ;
      k++ ;
    }
    printf("\n") ;
  }
  printf("%8d\n",icount[100]) ;
  
}