module std.random;

/*
pseudo random number generator with a period of (2^19937 - 1)
using Takuji Nishimura and Makoto Matsumoto's "Mersenne Twister 19937" algorithm

written by Wang Zhen (nehzgnaw@gmail.com)
based on Takuji Nishimura and Makoto Matsumoto's mt19937ar.c (2002-version)
freely available at http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
*/

private import std.math;

version (Win32)
{
	extern(Windows) int QueryPerformanceCounter(ulong *count);
}

version (linux)
{
	private import std.c.linux.linux;
}


class Random
{
private:
	const int N = 624, M = 397;
	const uint UPPER_MASK = 0x80000000U, LOWER_MASK = 0x7fffffffU;
	uint mt[N];
	int mti;
	bit next_gaussian_available;
	double next_gaussian_value;
	uint generate_seed()
	// returns a seed value based on time
	{
		// using Walter Bright's code in this method
		ulong s;
		version(Win32)
		{
			QueryPerformanceCounter(&s);
		}
		version(linux)
		{
			// time.h
			// sys/time.h

			timeval tv;

			if (gettimeofday(&tv, null))
			{   // Some error happened - try time() instead
				s = time(null);
			}
			else
			{
				s = tv.tv_usec;	// changed from : s = (cast(long)tv.tv_sec << 32) + tv.tv_usec;
			}
		}
		return cast(uint)s;
	}
public:
	this()
	{
		set_seed();
	}
	this(uint seed)
	{
		set_seed(seed);
	}
	void set_seed()
	// initializes the random sequence with a seed based on time
	{
		set_seed(generate_seed());
	}
	void set_seed(uint seed)
	// initializes the random sequence with the given seed value
	{
		mt[0] = seed;
		for(mti = 1; mti < N; ++mti)
			mt[mti] = (1812433253U * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
	}
	bit next_bit()
	// returns the next uniformly distributed random bit
	{
		return (next_uint() & 1) == 1;
	}
	ubyte next_ubyte()
	// returns the next uniformly distributed random ubyte
	{
		return cast(ubyte)next_uint();
	}
	ushort next_ushort()
	// returns the next uniformly distributed random ushort
	{
		return cast(ushort)next_uint();
	}
	uint next_uint()
	// returns the next uniformly distributed random uint
	{
		uint result;
		static uint mag01[2] = [0, 0x9908b0dfU];
		if(mti == N)
		{
			uint i, y;
			for(i = 0; i < N-M; i++)
			{
				y = (mt[i] & UPPER_MASK) | (mt[i+1] & LOWER_MASK);
				mt[i] = mt[i+M] ^ (y >> 1) ^ mag01[y & 1];
			}
			for(; i < N-1; i++)
			{
				y = (mt[i] & UPPER_MASK) | (mt[i+1] & LOWER_MASK);
				mt[i] = mt[i+(M-N)] ^ (y >> 1) ^ mag01[y & 1];
			}
			y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
			mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 1];
			mti = 0;
		}

		result = mt[mti++];

		result ^= (result >> 11);
		result ^= (result << 7) & 0x9d2c5680U;
		result ^= (result << 15) & 0xefc60000U;
		result ^= (result >> 18);

		return result;
	}
	ulong next_ulong()
	// returns the next uniformly distributed random ulong
	{
		return (cast(ulong)next_uint() << 32) + next_uint();
	}
	/* reserved for future use
	ucent next_ucent()
	// returns the next uniformly distributed random ucent
	{
		return (cast(ucent)next_ulong() << 64) + next_ulong();
	}
	*/
	float next_float()
	// returns the next uniformly distributed random float between 0.0f (inclusive) and 1.0f (exclusive)
	{
		return (next_uint() & 0x00ffffffU) / (cast(float)(1 << 24));
	}
	double next_double()
	// returns the next uniformly distributed random double between 0.0 (inclusive) and 1.0 (exclusive)
	{
		return (((cast(ulong)next_uint() & 0x03ffffffUL) << 27) + (next_uint() & 0x07ffffffUL))
				/ cast(double)(1L << 53);
	}
	double next_gaussian()
	// returns the next normally distributed random double with mean 0.0 and standard deviation 1.0
	{		
		if(next_gaussian_available)
		{
			next_gaussian_available = false;
			return next_gaussian_value;
		}
		else
		{
			next_gaussian_available = true;
			double d0, d1, m;
			do
			{ 
				d0 = 2.0 * next_double() - 1.0;
				d1 = 2.0 * next_double() - 1.0;
				m = d0 * d0 + d1 * d1;
			}while(m >= 1.0 || m == 0.0);
			m = sqrt(-2.0 * log(m) / m);
			next_gaussian_value = d1 * m;
			return d0 * m;
		}
	}
}

unittest
{
	static uint testcases[][] =
	[
		[0x00000000, 0x0000, 0x8C7F0AAC],
		[0x00000000, 0x06B0, 0x2D504023],
		[0x00000001, 0x0000, 0x6AC1F425],
		[0x00000001, 0x0001, 0xFF4780EB],
		[0x00000001, 0x000A, 0x17A38090],
		[0x0CFFF866, 0x0000, 0x18209ACC],
		[0x411524A5, 0x1F41, 0xC1012C18],
		[0xDFF9D2A2, 0x0400, 0xCAF33CC3],
		[0xF552D63C, 0x0001, 0x0BE03D6D],
		[0xFFFF0000, 0x03E7, 0x95F36C9C],
		[0xFFFFFFFF, 0x0100, 0x131574D3],
	];

	Random rng = new Random();
	foreach(uint[] tc; testcases)
	{
		rng.set_seed(tc[0]);
		for(uint i = 0; i < tc[1]; ++i)
			rng.next_uint();
		assert(rng.next_uint() == tc[2]);
	}

}