/*
  File: newrand.c
  CS 410 Programming assignment 3 helper file, fall 1997
  
  A portable random number generation subroutine from 
  "Numerical Recipes" by William H. Press, etc. pp.196-197
*/

#define  M1   259200
#define  IA1  7141
#define  IC1  54773
#define  RM1  1.0/M1
#define  M2   134456
#define  IA2  8121
#define  IC2  28411
#define  RM2  1.0/M2
#define  M3   243000
#define  IA3  4561
#define  IC3  51349

/*
  double newrand (int idum)
  When the input parameter 'idum' is negative (for example,
  -17), the generator is initialized.
  Else, returns a uniform # between 0 and 1.0.
*/

double newrand (int idum)
{
    static double  r[97];
    static int ix1, ix2, ix3;
    int j;
    double result;

    if (idum < 0) {

	/* seed the first routine */
	ix1 = (IC1 - idum) % M1;
	ix1 = (IA1*ix1 + IC1) % M1;

	/* use it to seed the second routine */
	ix2 = ix1 % M2;

	ix1 = (IA1*ix1+IC1) % M1;

	/* and use it to seed the third routine */
	ix3 = ix1 % M3;

	/* fill the table with sequential uniform deviates generated */
	/* by the first two routines				     */
	for (j = 0; j < 97; j++) {
	    ix1 = (IA1*ix1+IC1) % M1;
	    ix2 = (IA2*ix2+IC2) % M2;
	    r[j] = ((double)ix1 + (double)ix2*RM2)*RM1;
	};

    };

    /* except when initializing, this is where we start */
    /* generating the next number for each sequence     */
    ix1 = (IA1*ix1+IC1) % M1;
    ix2 = (IA2*ix2+IC2) % M2;
    ix3 = (IA3*ix3+IC3) % M3;

    /* use the third sequence to get an integer between 0 and 96 */
    j = 97*ix3/M3;
    if (j < 0 || j > 96) 
	printf("error: in drnum\n");

    /* return the table entry */
    result = r[j];

    /* refill the table entry */
    r[j] = ((double)ix1 + ((double)ix2)*RM2)*RM1;

    return (result);
}