// gprand.cc
// 3 June 1994

//--------------------------------------------------------------------------
// This code is a component of Genetic Programming in C++ (Version 0.40)
// For comments, improvements, additions (or even money !?) contact:
// Adam Fraser, Postgraduate Section, Dept of Elec & Elec Eng,
// Maxwell Building, University Of Salford, Salford, M5 4WT, United Kingdom.
// Internet: a.fraser@eee.salford.ac.uk
// Tel: (UK) 061 745 5000 x3633
// Fax: (UK) 061 745 5999
//--------------------------------------------------------------------------


// This has been hacked from Tierra which is by Tom Ray and others random number generator
// A paraphrase of their comments follow
/* 9-9-92  adapted from ``Numerical Recipes in C'' by Press,
Flannery, Teukolsky and Vetterling. */


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


// global variables an absolute necessity for system.......
unsigned long RandIx1, RandIx2, RandIx3;
double TrandArray[98];

extern void ExitSystem( char* );

void gp_srand( unsigned long seed)
{   
    unsigned int  j;

// seed the first routine........................
    RandIx1 = (IC1 + seed) % M1;     
//... which is then used to seed the second.....
    RandIx1 = (IA1 * RandIx1 + IC1) % M1;
    RandIx2 = RandIx1 % M2;   
// .. and third routines....
    RandIx1 = (IA1 * RandIx1 + IC1) % M1;
    RandIx3 = RandIx1 % M3;        
// Fill the random variable table with sequential uniform values
    for (j = 1; j <= 97; j++)   
    {   
// Deviates generated by the first two routines.............................
	RandIx1 = (IA1 * RandIx1 + IC1) % M1;
	RandIx2 = (IA2 * RandIx2 + IC2) % M2;
// Low and high order pieces combined here to create value..... 
	TrandArray[j] = (RandIx1 + RandIx2 * RM2) * RM1;
    }                    
}

// returns an unsigned long
unsigned long gp_rand()
{   
    double temp;
    unsigned int j;

//  Except when initializing, this is where we start.  Generate the next number 
// for each sequence.
    RandIx1 = (IA1 * RandIx1 + IC1) % M1;
    RandIx2 = (IA2 * RandIx2 + IC2) % M2;
    RandIx3 = (IA3 * RandIx3 + IC3) % M3;  
    
// Use the third sequence to get an integer between 1 and 97
    j = 1 + ((97 * RandIx3) / M3);
    if (j > 97 || j < 1)        ExitSystem( "gplrand()" );

// Return that table entry, and refill it. 
    temp = TrandArray[j] * M3;               
    TrandArray[j] = (RandIx1 + RandIx2 * RM2) * RM1;   

   return (unsigned long)temp;
}       

#undef M1
#undef IA1
#undef IC1
#undef RM1
#undef M2
#undef IA2
#undef IC2
#undef RM2
#undef M3
#undef IA3
#undef IC3

// gprand.cc