// the following program (generator.cpp) includes the test functions
// f and neighbor, plus the random-variate generation routines

#include "generator.h"

//****************************************************************************

double f(double theta )
{
// this function returns a single random value from the
// quadratic test function of Nelson at point theta.
// it calls CalcNorm, a function to generate N(0,1)
// random variates via the inverse cdf method, by
// passing a U(0,1) generated by RandGen()
//
       double response;
       response = 0.5*theta*theta + 4*CalcNorm(rand(2));
       return response;

}// end fn f


int neighbor(int theta)
{
// this routine generates a candidate solution from the neighborhood
// of theta for the feasible region THETA = [-100, -99, ..., 99, 100]
// by wrapping the feasible region.
//If a bad value is entered, the function returns -9999

	 //input error check
	 if (theta>100 || theta <-100)
	 {return -9999;}

      double U;
	  int newtheta;

      newtheta = theta + 1;
      U = rand(6);
      if ( U < 0.5 ){newtheta = theta - 1;}
// wrap if necessary
      if (newtheta == 101)
			{newtheta = -100;}

      if (newtheta == -101)
			{newtheta = 100;}

      return newtheta;
}//end fn neighbor


/////////////////////////////////////////////////////////////////
double CalcNorm(double U)
{	
// Approximate inverse cdf for standard normal
// Based on vnorm in Bratley, Fox, and Schrage, p.339

	double P		= 1 - U;

	const double PLIM	= 1.0e-018; 

	const double P0		= -0.322232431088;
	const double P1		= -1.0;
	const double P2		= -0.342242088547;
	const double P3		= -0.0204231210245;
	const double P4		= -0.0000453642210148;

	const double Q0		= 0.099348462606;
	const double Q1		= 0.588581570495;
	const double Q2		= 0.531103462366;
	const double Q3		= 0.10353775285;
	const double Q4		= 0.0038560700634;

if (P > .5) {P = U;} //end if

double VTemp = 8;  //Initialize

if (P >= PLIM) //make sure P isn't tiny
{
double Y = sqrt(-log(P*P));
VTemp = Y + 
       ((((Y*P4 +P3)*Y +P2)*Y +P1)*Y +P0)/((((Y*Q4 +Q3)*Y +Q2)*Y +Q1)*Y +Q0);
}//end if

else {return 1;} //if it is tiny, return 1

if (1.0 - U < .5) {VTemp = -VTemp;}

double VNorm = VTemp;
return VNorm;
}//end CalcNorm


//////////////////////////////////////////////////////////////////////////
//Random-Number Generator form Law and Kelton, p454
//Prime modulus multiplicative linear congruential generator
//Z[i]=(630360013 *Z[i-1]) (mod(pow2,31)-1), based on Marse and
//Roberts portable Fortran random-number generator UINRAN.  Multiple
//(100) streams are supported, with seeds spaced 100,000 apart.
//throughout, input argument "stream" must be an int giving the 
//desired stream number.  The header file rand.h must be included
// in the calling program (#include "rand.h") before using these 
// functions

//Usage: (Three functions)

//1. To obtain the next U(0,1) random number from stream "stream,"
//		execute
//			u = rand(stream);
//		where rand is a float function.  The float variable u will
//		contain  the next random number.
//
//2. To set the seed for stream "stream" to a desired value zset,
//		execute
//			randst(zset, stream);
//		where randst is a void function and szet must be a long set to
//		the desired seed, a number between 1 and 2147483646 (inclusive).
//		Default seeds for all 100 streams are given in the code.
//
//3. To get the current (most recently used)integer in the sequence
//	 being generated for stream "stream" into the long variable zget,
//		execute
//		  zget = randgt(stream)
//		where randgt is a long function



/* Generate the next random number. */
double rand(int stream)
{
    long zi, lowprd, hi31;

    zi     = zrng[stream];
    lowprd = (zi & 65535) * MULT1;
    hi31   = (zi >> 16) * MULT1 + (lowprd >> 16);
    zi     = ((lowprd & 65535) - MODLUS) +
             ((hi31 & 32767) << 16) + (hi31 >> 15);
    if (zi < 0) zi += MODLUS;
    lowprd = (zi & 65535) * MULT2;
    hi31   = (zi >> 16) * MULT2 + (lowprd >> 16);
    zi     = ((lowprd & 65535) - MODLUS) +
             ((hi31 & 32767) << 16) + (hi31 >> 15);
    if (zi < 0) zi += MODLUS;
    zrng[stream] = zi;
    return (( (double) (zi >> 7 | 1) + 1 )/ 16777216.0);
}

////////////////////////////////////////////////////////////////////////////////
//Set the current zrng for stream "stream" to zset
void randst (long zset, int stream)
{
	zrng[stream] = zset;
}
/////////////////////////////////////////////////////////////////////////////////


/////////////////////////////////////////////////////////////////////////////////
//Return the current zrng for stream "stream".

long randgt (int stream)
{
	return zrng[stream];
}
//////////////////////////////////////////////////////////////////////////////////