/////////////////////////////////////////////////
//
// File:
//  trandom.cpp	
//
// Description:
//	Random number generator based on program
//      in Numerical Recipes in C.
//
// Last modified:
//	18 August 1997
// (c) 1997, Thomas F. El-Maraghi
//
/////////////////////////////////////////////////


#include <math.h>
#include <time.h>

#ifndef GNUGPP
 #ifndef BOOL_DEFN
  #include "bool.h"
 #endif
#endif

#define _trandom_cpp
#include "trandom.h"
#undef _trandom_cpp


/////////////////////////////////////////////////
// Global random number generator

cRandom Random;

// 
// To avoid possible correlations you should use only this one
// global object Random.  In particular, do not use several different
// cRandom objects for parallel streams, since there could be
// unwanted correlations between the different streams.
/////////////////////////////////////////////////

/////////////////////////////////////////////////
// Random number class implementation


// The following definitions are identical to those
// in the original implementation as found in
// "Numerical Recipies"
#define IA 16807
#define IM 2147483647
#define IQ 127773
#define IR 2836
#define AM (1.0/IM)
#define NDIV (1+(IM-1)/SEED_TABLE_SIZE)
#define EPS 1.2e-7
#define RandomMax (1.0-EPS)


void cRandom::Seed( long number )
{
	// The Seed() function accepts a seed from the user
	// and uses it to fill in the mSeedTable shuffle table.  Values
	// in the mSeedTable shuffle table are then used by the
	// get function to generate the pseudo-random numbers.
	mSeed = (number<0) ? number : -number;
	if( -mSeed < 1)
		mSeed = 1;
	else
		mSeed = -mSeed;
	for( int i = SEED_TABLE_SIZE + 7 ; i >= 0 ; i-- ) {
		long k = mSeed / IQ;
		mSeed = IA*(mSeed - k * IQ) - IR*k;
		if( mSeed < 0 )
			mSeed += IM;
		if( i < SEED_TABLE_SIZE )
			mSeedTable[i] = mSeed;
	}
	mNextSeed = mSeedTable[0];
}


void cRandom::Seed()
{
	// A useful variant of the Seed() function which uses the system
	// clock to generate the initial seed.
	Seed( (long)time( 0 ) );
}


float cRandom::Uniform()
{
	// Teturns a pseudo-random number evenly distributed in the range (0,1)
	if( mSeed <= 0 || !mNextSeed )
		Seed();
	long k = mSeed / IQ;
	mSeed = IA*(mSeed - k * IQ) - IR*k;
	if( mSeed < 0 )
		mSeed += IM;
	int Pick = (int)( mNextSeed / NDIV );
	mNextSeed = mSeedTable[Pick];
	mSeedTable[Pick] = mSeed;
	float RetVal = (float)AM * mNextSeed;
	if( RetVal > RandomMax )
		return (float)RandomMax;
	else
		return RetVal;
}


float cRandom::Exponential( float Lambda )
{
	// The Exponential() funciton returns the waiting time to the
	// next event of a Poisson distribution which has
	// Lambda events per second.
	float Dummy;
	do
	   Dummy = Uniform();
	while( Dummy == 0.0 );
	return (float)(-log(Dummy)/Lambda);
}


float cRandom::Gaussian()
{
	// Returns a normally distributed random number with zero
	// mean and unit variance
	if( fExtraGaussian )
	{
		// if we have an extra deviate available, then
		// use it
		fExtraGaussian = false;
		return ExtraGaussian;
	}
	
	// Get two uniform deviates inside the unit circle
	float uniform1, uniform2, radius_square;
	do {
		uniform1 = (float)(2.0*Uniform()-1.0);
		uniform2 = (float)(2.0*Uniform()-1.0);
		radius_square = uniform1*uniform1 + uniform2*uniform2;
	} while( radius_square >= 1.0 || radius_square == 0.0 );

	// Compute 2 gaussian deviates, saving one for later
	float factor = (float)sqrt(-2.0*log(radius_square)/radius_square);
	ExtraGaussian = uniform1*factor;
	fExtraGaussian = true;
	return uniform2*factor;
}

