// This code is in the public domain -- castanyo@yahoo.es

#ifndef NV_MATH_RANDOM_H
#define NV_MATH_RANDOM_H

#include <nvcore/Containers.h> // nextPowerOfTwo
#include <nvmath/nvmath.h>

namespace nv
{

/// Interface of the random number generators.
class Rand
{
public:

	virtual ~Rand() {}

	enum time_e { Time };

	/// Provide a new seed.
	virtual void seed( uint s ) { /* empty */ };
	
	/// Get an integer random number.
	virtual uint get() = 0;

	/// Get a random number on [0, max] interval.
	uint getRange( uint max )
	{
		uint n;
	//	uint mask = Bitmask( max );
	//	do { n = Get() & mask; } while( n > max );		
		uint np2 = nextPowerOfTwo( max );
		do { n = get() & (np2-1); } while( n > max );
		return n;
	}

	/// Random number on [0.0, 1.0] interval.
	float getFloat()
	{
    	union
		{
			uint32 i;
			float f;
		} pun;

		pun.i = 0x3f800000UL | (get() & 0x007fffffUL);
		return pun.f - 1.0f;
	}

	/*
	/// Random number on [0.0, 1.0] interval.
	double getReal()
	{
		return double(get()) * (1.0/4294967295.0); // 2^32-1
	}
	
	/// Random number on [0.0, 1.0) interval.
	double getRealExclusive()
	{
		return double(get()) * (1.0/4294967296.0); // 2^32
	}
	*/

	/// Get the max value of the random number.
	uint max() const { return 4294967295U; }

	// Get a random seed.
	static uint randomSeed();
	
};


/// Very simple random number generator with low storage requirements.
class SimpleRand : public Rand
{
public:

	/// Constructor that uses the current time as the seed.
	SimpleRand( time_e )
	{
		seed(randomSeed());
	}
	
	/// Constructor that uses the given seed.
	SimpleRand( uint s = 0 )
	{
		seed(s);
	}
	
	/// Set the given seed.
	virtual void seed( uint s )
	{
		current = s;
	}
	
	/// Get a random number.
	virtual uint get()
	{
		return current = current * 1103515245 + 12345;
	}

private:

	uint current;
	
};


/// Mersenne twister random number generator.
class MTRand : public Rand
{
public:

	enum { N = 624 };       // length of state vector
	enum { M = 397 };

	/// Constructor that uses the current time as the seed.
	MTRand( time_e )
	{
		seed(randomSeed());
	}
	
	/// Constructor that uses the given seed.
	MTRand( uint s = 0 )
	{
		seed(s);
	}
	
	/// Constructor that uses the given seeds.
	NVMATH_API MTRand( const uint * seed_array, uint length );

	
	/// Provide a new seed.
	virtual void seed( uint s )
	{
		initialize(s);
		reload();
	}	

	/// Get a random number between 0 - 65536.
	virtual uint get()
	{
		// Pull a 32-bit integer from the generator state
		// Every other access function simply transforms the numbers extracted here
		if( left == 0 ) { 
			reload(); 
		}
		left--;
		
		uint s1;
		s1 = *next++;
		s1 ^= (s1 >> 11);
		s1 ^= (s1 <<  7) & 0x9d2c5680U;
		s1 ^= (s1 << 15) & 0xefc60000U;
		return ( s1 ^ (s1 >> 18) );		
	};

	
private:
	
	NVMATH_API void initialize( uint32 seed );
	NVMATH_API void reload();
	
	uint hiBit( uint u ) const { return u & 0x80000000U; }
	uint loBit( uint u ) const { return u & 0x00000001U; }
	uint loBits( uint u ) const { return u & 0x7fffffffU; }
	uint mixBits( uint u, uint v ) const { return hiBit(u) | loBits(v); }
	uint twist( uint m, uint s0, uint s1 ) const { return m ^ (mixBits(s0,s1)>>1) ^ ((~loBit(s1)+1) & 0x9908b0dfU); }
		
private:

	uint state[N];	// internal state
	uint * next;	// next value to get from state
	int left;		// number of values left before reload needed		

};



/** George Marsaglia's random number generator. 
 * Code based on Thatcher Ulrich public domain source code:
 * http://cvs.sourceforge.net/viewcvs.py/tu-testbed/tu-testbed/base/tu_random.cpp?rev=1.7&view=auto
 *
 * PRNG code adapted from the complimentary-multiply-with-carry
 * code in the article: George Marsaglia, "Seeds for Random Number
 * Generators", Communications of the ACM, May 2003, Vol 46 No 5,
 * pp90-93.
 * 
 * The article says:
 * 
 * "Any one of the choices for seed table size and multiplier will
 * provide a RNG that has passed extensive tests of randomness,
 * particularly those in [3], yet is simple and fast --
 * approximately 30 million random 32-bit integers per second on a
 * 850MHz PC.  The period is a*b^n, where a is the multiplier, n
 * the size of the seed table and b=2^32-1.  (a is chosen so that
 * b is a primitive root of the prime a*b^n + 1.)"
 * 
 * [3] Marsaglia, G., Zaman, A., and Tsang, W.  Toward a universal
 * random number generator.  _Statistics and Probability Letters
 * 8_ (1990), 35-39.
 */
class GMRand : public Rand
{
public:

	enum { SEED_COUNT = 8 };

//	const uint64 a = 123471786;		// for SEED_COUNT=1024
//	const uint64 a = 123554632;		// for SEED_COUNT=512
//	const uint64 a = 8001634;		// for SEED_COUNT=255
//	const uint64 a = 8007626;		// for SEED_COUNT=128
//	const uint64 a = 647535442;		// for SEED_COUNT=64
//	const uint64 a = 547416522;		// for SEED_COUNT=32
//	const uint64 a = 487198574;		// for SEED_COUNT=16
//	const uint64 a = 716514398U;	// for SEED_COUNT=8
	enum { a = 716514398U };
	

	GMRand( time_e )
	{
		seed(randomSeed());
	}
	
	GMRand(uint s = 987654321)
	{
		seed(s);
	}


	/// Provide a new seed.
	virtual void seed( uint s )
	{
		c = 362436;
		i = SEED_COUNT - 1;
		
		for(int i = 0; i < SEED_COUNT; i++) {
			s = s ^ (s << 13);
			s = s ^ (s >> 17);
			s = s ^ (s << 5);
			Q[i] = s;
		}
	}

	/// Get a random number between 0 - 65536.
	virtual uint get()
	{
		const uint32 r = 0xFFFFFFFE;		
		
		uint64 t;
		uint32 x;

		i = (i + 1) & (SEED_COUNT - 1);
		t = a * Q[i] + c;
		c = uint32(t >> 32);
		x = uint32(t + c);
		
		if( x < c ) {
			x++;
			c++;
		}
		
		uint32  val = r - x;
		Q[i] = val;
		return val;
	};


private:

	uint32 c;
	uint32 i;
	uint32 Q[8];

};


/** Random number implementation from the GNU Sci. Lib. (GSL).
 * Adapted from Nicholas Chapman version:
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 * This is the Unix rand48() generator. The generator returns the
 * upper 32 bits from each term of the sequence,
 * 
 * x_{n+1} = (a x_n + c) mod m 
 * 
 * using 48-bit unsigned arithmetic, with a = 0x5DEECE66D , c = 0xB
 * and m = 2^48. The seed specifies the upper 32 bits of the initial
 * value, x_1, with the lower 16 bits set to 0x330E.
 * 
 * The theoretical value of x_{10001} is 244131582646046.
 * 
 * The period of this generator is ? FIXME (probably around 2^48). 
 */
class Rand48 : public Rand
{
public:

	Rand48( time_e )
	{
		seed(randomSeed());
	}
	
	Rand48( uint s = 0x1234ABCD )
	{
		seed(s);
	}	


	/** Set the given seed. */
	virtual void seed( uint s ) {
		vstate.x0 = 0x330E;
		vstate.x1 = uint16(s & 0xFFFF);
		vstate.x2 = uint16((s >> 16) & 0xFFFF);
	}
	
	/** Get a random number. */
	virtual uint get() {
		
		advance();

		uint x1 = vstate.x1;
		uint x2 = vstate.x2;
		return (x2 << 16) + x1;
	}
	
	
private:
	
	void advance()
	{
		/* work with unsigned long ints throughout to get correct integer
		promotions of any unsigned short ints */
		const uint32 x0 = vstate.x0;
		const uint32 x1 = vstate.x1;
		const uint32 x2 = vstate.x2;
		
		uint32 a;
		a = a0 * x0 + c0;
		
		vstate.x0 = uint16(a & 0xFFFF);
		a >>= 16;
		
		/* although the next line may overflow we only need the top 16 bits
		in the following stage, so it does not matter */
		
		a += a0 * x1 + a1 * x0; 
		vstate.x1 = uint16(a & 0xFFFF);
		
		a >>= 16;
		a += a0 * x2 + a1 * x1 + a2 * x0;
		vstate.x2 = uint16(a & 0xFFFF);
	}

	
private:	
	NVMATH_API static const uint16 a0, a1, a2, c0;

	struct rand48_state_t { 
		uint16 x0, x1, x2; 
	} vstate;

};

} // nv namespace

#endif // NV_MATH_RANDOM_H