#ifndef _INCLUDED_L3_RAND_H
#define _INCLUDED_L3_RAND_H

// *Random Generator* Copyright (C) Krzysztof Bosak, 1998-03-01...2000-12-26.
// All rights reserved.
// kbosak@box43.pl
// http://www.kbosak.prv.pl

#ifndef UNDER_CE
#include <time.h>
#endif

// TO DO:
// add save state
// move benchmarking outside l3 library

// I found very quick way of generating drandoms with given density,
// let's say 'f(x)=x': you find reverse function g(x) to Integrate(f(x), dx),
// then you use 'return g(drandom()); F.ex. for generating drandoms with
// f(x)=sin(x) you use 'return acos(drandom()/(pi/2));'

// To generate uniform distribution of real numbers on unary interval:
// - [0,1] use return irandom()/irandommax() (most common case in implementation)
// - [0,1) use return irandom()/(irandommax()+1) (most demanded case by mathematicians)
//   WARNING: when irandommax()==INT_MAX (or UINT_MAX respectively) to not exceed
//   integer type range you must use return irandom()/(irandommax()+1.0) which contains
//   conversion of denominator to double, and mdivision is slower in this case. Division
//   by constant is changed by compilators to multiplication, so please don't write whole
//   denominator in reciprocal double type form.
// - (0,1] use [0,1] method and discard zeros by if() condition
// - (0,1) use [0,1) method and discard zeros by if() condition
// ...where irandom() gives random integers ranging from 0 to irandommax() INCLUSIVE.

#include <math.h>
#include <stdlib.h>
#include "l3_mt.h"

class randomgen
{
private:
	static int _randomgen_seed_modifier;
	int _seed;

	mutable bool _is_v1_unused;
	int _irandommax_half;
	double _irandommax_plus1_inv;

	mutable double _gauss_v1;
	mutable double _gauss_v2;

public:
	randomgen();
	explicit randomgen(int new_seed);
	void randomize();
	void setseed(int new_seed);
	inline int seed() const
	{
		return _seed;
	}
	inline void reheat() const
	{	// let's forget about history of seed,
		// to make sure that even first generated random will be truly random.
		for(int reheat_stage=0; reheat_stage<7; reheat_stage++)
		{
			drandom();
		}
	}
	inline int irandom() const
	{
		//return gfsr4int();// GFSR
		return static_cast<int>(mtint());// MT
	}
	inline double drandom() const
	{	// 0...1 inclusive
		//return gfsr4();// GFSR
		return mtdouble();// MT
	}
	inline int irandom(int upperlimit) const
	{	//0...upperlimit-1
		assert(upperlimit>0 && upperlimit<=irandommax());
		return static_cast<int>(upperlimit*drandom());
	}
	inline int irandom(int low, int top) const
	{	// low...top inclusive
		assert(top>=low);
		return low+irandom(top-low+1);
	}
	inline bool brandom() const
	{	// 0, 1
		return irandom()<=_irandommax_half;
	}
	inline bool brandom(double probability_of_true) const
	{	// returns true with probability_of_true (inclusive)
		return drandom()<=probability_of_true;
	}
	inline int irandommin() const
	{
		return 0;
	}
	inline int irandommax() const
	{
		//return 2147483561;// GFSR, 2^31-87
		return 2147483647;// MT
	}
	inline double drandom(double b) const
	{	// 0...b inclusive
		assert(b>0);
		return drandom()*b;
	}
	inline double anglerandom() const
	{	// 0...2Pi inclusive
		return drandom()*3.1415926535897932384626433832795*2;
	}
	inline double drandom(double a, double b) const
	{	// a...b inclusive
		assert(b>a);
		return a+drandom()*(b-a);
	}
	inline double drandomgauss() const
	{	// N(0,1) - gauss double using polar method, generates 2 numbers at once.
		// returns one and stored another one for future usage.
		if(_is_v1_unused==true)
		{
			_is_v1_unused=false;
			return _gauss_v1;
		}
		const double result=_gauss_v2;
		double s;
		do
		{
			_gauss_v1=2*drandom()-1;
			_gauss_v2=2*drandom()-1;
			s=_gauss_v1*_gauss_v1+_gauss_v2*_gauss_v2;
		}while(s>=1);
		const double mul=sqrt(-2*log(s)/s);
		_gauss_v1*=mul;
		_gauss_v2*=mul;
		_is_v1_unused=true;
		return result;
	}
	inline double drandomgauss(double avg, double sigma) const
	{	// N(avg,sigma) - gauss double using polar method
		return avg+drandomgauss()*sigma;
	}
	inline double drandomgauss(double sigma) const
	{	// N(0,sigma) - gauss double using polar method
		return drandomgauss()*sigma;
	}
};
int randomgen::_randomgen_seed_modifier=1;
randomgen RAND;

randomgen::randomgen()
	: _is_v1_unused(false),
	_irandommax_half(irandommax()>>1),
	_irandommax_plus1_inv(1/(irandommax()+1.0))// +1.0 and conversion to double is crucial
{
	randomize();
}

randomgen::randomgen(int new_seed)
	: _is_v1_unused(false),
	_irandommax_half(irandommax()>>1),
	_irandommax_plus1_inv(1/(irandommax()+1.0))// +1.0 and conversion to double is crucial
{
	setseed(new_seed);
}

void randomgen::randomize()
{
	#ifdef UNDER_CE
		setseed(static_cast<int>(
		(_randomgen_seed_modifier)%static_cast<unsigned>(irandommax())+1u
		));
	#else
		setseed(static_cast<int>(
		(time(NULL)+_randomgen_seed_modifier) % static_cast<unsigned>(irandommax())+1u
		));
	#endif
	_randomgen_seed_modifier++;
}

void randomgen::setseed(int new_seed)
{
	assert(new_seed>0);
	assert(new_seed<=irandommax());
	_seed=new_seed;
	//gfsr4setseeds(_seed);// GFSR
	mtsetseeds(static_cast<unsigned>(_seed));// MT
	_is_v1_unused=false;
	drandomgauss();
	reheat();
}

#ifdef MAINRAND
#include "l3_time.h"
#include <string.h>

/*
#define RANDCOND RAND.brandom()
#define RANDICE RAND.irandom(6)+1
void test_brown(int steps, int trials)
{
	int res[4]={0};
	int x, y;
	for(int t=0; t<trials; t++)
	{
		x=0;
		y=0;
		for(int s=0; s<steps; s++)
		{
			if(RANDCOND)
				x++;
			else
				x--;
			if(RANDCOND)
				y++;
			else
				y--;
		}
		if(x>0 && y>0)
			res[0]++;
		if(x>0 && y<0)
			res[1]++;
		if(x<0 && y<0)
			res[2]++;
		if(x<0 && y>0)
			res[3]++;
	}
	cout<<res[0]<<"|"<<res[1]<<"|"<<res[2]<<"|"<<res[3]<<endl;
}
void test_chi2(int n)
{// >18.31 5%		>23.21 1%
	for(int t=0; t<20; t++)
	{
		int c[13]={0};
		const double p[13]={0,0,1.0/36,1.0/18,1.0/12,1.0/9,5.0/36,1.0/6,5.0/36,1.0/9,1.0/12,1.0/18,1.0/36};
		for(int i=0; i<n; i++)
			c[RANDICE+RANDICE]++;
		double v=0;
		for(int k=2; k<=12; k++)
		{
			//cout<<c[k]<<"|";
			v+=pow((c[k]-n*p[k]),2) / (n*p[k]);
		}
		cout<<v<<endl;
	}
}
void test_init()
{
	srand(time(NULL));//rand()
	mtsetseeds(time(NULL));//mt
}
int l25t1knuth()
{
	static int seed=1234;
	return (seed=(69069 * seed) & 0xFFFFFFFF);
}
*/

//#define RANDFUNC gfsr4int()
#define RANDFUNC mtint()
//#define RANDFUNC l25t1knuth()
void test_ranges()
{
	int min=INT_MAX;
	int max=-1;
	int old[16];
	int look;
	for(look=0; look<16; look++)
		old[look]=RANDFUNC;
	look=0;
	int count=0;
	int r;
	for(;;)
	{
		r=RANDFUNC;
		if(r<min)
		{
			min=r;
			cout<<"Count="<<count<<" Min="<<min<<" Max="<<max<<endl;
		}
		if(r>max)
		{
			max=r;
			cout<<"Count="<<count<<" Min="<<min<<" Max="<<max<<endl;
		}
		count++;
		if(( count & ((1<<24)-1) )==0)
		{
			cout<<"Count="<<count<<" Min="<<min<<" Max="<<max<<endl;
		}
		if(r==old[look])
		{
			look++;
			if(look>15)
			{
				break;
			}
		}
		else
		{
			look=0;
		}
	}
	count-=16;
	cout<<"Count="<<count<<" Min="<<min<<" Max="<<max<<endl;
}

void test_speed()
{
	/*
	SPEED TESTING: CeleronA416mhz128cache83bus 3e10 trials:
	rand();//101.72s
	gfsr4int()//95.24s
	mtint()//174.99s
	tezuka();//109.03s
	RAND.irandom();//~1.8*gfsr_time! Stupid M$ compilers.
	*/
	Timing tim;
	for(int i=0; i<300000000; i++)//300000000
	{
		//rand();
		//ran2int(&seed);
		//gfsr4int();
		//rnlx();
		RAND.irandom();
		//RAND.drandomgauss();
	}
	cout<<tim<<endl;
}
/*
int main()
{
	for(int i=0; i<20; i++)
	{
		cout<<RAND.drandomgauss()<<endl;
	}
//	test_init();
//	test_brown(2000, 15000);
//	test_chi2(100000);
	test_speed();
	return 0;
	test_ranges();
	return 0;
}
*/
#endif

#endif //_INCLUDED_L3_RAND_H