module y4d_math.rand;

private import std.date; // 乱数の種にgetUTCtime()を用いる
private import ytl.exception;

/**
	Mersenne Twister法による高精度の乱数発生ルーチンです。
	2^19937-1という天文学的な大きさの周期性を持った乱数を高速に生成します。

	コンストラクタは、２種類あり、パラメータ無しのほうは、
	乱数の初期化を行ないません。パラメータ有りのほうは、
	乱数の種を引数として取り、それに基づいた乱数を生成します。
	（内部的にsetSeedを行なうということです。
	よって、setSeedをこのあと行なう必要はありません）
	前回と同じ乱数系列を再現したい場合などにこれを使います。
*/
///	高速で長い周期の乱数発生クラス。
class Rand {

	///	乱数の取得。0～n-1までの乱数を取得。
	/**
		n==0の場合は y4d_error がthrowされる
	*/
	uint	get(uint n) {			//	0～n-1の乱数の取得
		if (n==0) { 
			throw new y4d_error(this,cast(char[]) "getで引数が0");
		}
		return get() % n;
	}

	///	乱数の取得。uint全域。
	uint	get()					///	32bit精度の乱数の取得
	{
		uint y;
		static uint mag01[2] = [ 0x0, MATRIX_A ];
		/* mag01[x] = x * MATRIX_A	for x=0,1 */

		if (m_nMti >= N) { /* generate N words at one time */
			int kk;
			if (m_nMti == N+1)	 /* if sgenrand() has not been called, */
				setSeed(4357);	 /* a default initial seed is used	 */

			for (kk=0;kk<N-M;kk++) {
				y = (m_dwMt[kk]&UPPER_MASK)|(m_dwMt[kk+1]&LOWER_MASK);
				m_dwMt[kk] = m_dwMt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
			}
			for (;kk<N-1;kk++) {
				y = (m_dwMt[kk]&UPPER_MASK)|(m_dwMt[kk+1]&LOWER_MASK);
				m_dwMt[kk] = m_dwMt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
			}
			y = (m_dwMt[N-1]&UPPER_MASK)|(m_dwMt[0]&LOWER_MASK);
			m_dwMt[N-1] = m_dwMt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
			m_nMti = 0;
		}
		y = m_dwMt[m_nMti++];

		y ^= (y >> 11);
		y ^= (y << 7) & TEMPERING_MASK_B;
		y ^= (y << 15) & TEMPERING_MASK_C;
		y ^= (y >> 18);
		return y; 
	}

	///	乱数の発生。[0,1] 範囲で乱数生成 ←0,1を含むの意味
	double	get1()
	{
		return cast(double)get() * (cast(real)1.0/4294967295.0);	
		/* divided by 2^32-1 */ 
	}

	///	乱数の発生。[0,1) 範囲で乱数生成 ←0は含む,1は含まないの意味
	double	get2()
	{
		return cast(double)get() * (cast(real)1.0/4294967296.0); 
		/* divided by 2^32 */
	}

	/// 乱数の発生。(0,1] 範囲で乱数生成 ←0は含まない,1は含むの意味
	double get3()
	{
		return (cast(double)get()+1) * (cast(real)1.0/4294967296.0); 
		/* divided by 2^32 */
	}

	///	乱数の発生。(0,1) 範囲で乱数生成 ←0,1は含まないの意味
	double get4()
	{
		return (cast(double)get() + 0.5) * (cast(double)1.0/4294967296.0);	
		/* divided by 2^32 */
	}

	///　乱数の種を設定します。
	void	setSeed(uint dwSeed)
	/**
		設定された種に従って、乱数は生成されていきます。
		必ず一度は呼び出す必要があります。
		呼び出さないときは、
			SetSeed(4357);
		が、一番最初の乱数取得のときに実行されます。
	*/
	{
		//	乱数の種の設定。必ず一度呼び出す必要がある
		const int N = 624;
		for (int i=0;i<N;i++) {
			m_dwMt[i] = dwSeed & 0xffff0000;
			dwSeed = 69069 * dwSeed + 1;
			m_dwMt[i] |= (dwSeed & 0xffff0000) >> 16;
			dwSeed = 69069 * dwSeed + 1;
		}
		m_nMti = N;
	}

	///	乱数の種として、現在時刻を与えます。
	void	randomize()
	/**
		要するに、再現性の無い乱数が出来ます。
		setSeed(タイマーの値)とやっているので、
		setSeedを呼び出す必要はありません。
	*/
	{
		setSeed(cast(uint)(getUTCtime()));
	}

	/**
		コンストラクタは、２種類あり、パラメータ無しのほうは、
		乱数の初期化を行ないません。パラメータ有りのほうは、
		乱数の種を引数として取り、それに基づいた乱数を生成します。
		（内部的にSetSeedを行なうということです。
		よって、SetSeedをこのあと行なう必要はありません）
		前回と同じ乱数系列を再現したい場合などにこれを使います。
	*/
	this() { m_nMti = 624+1; } /* means m_dwMt is not initialized */
	this(uint dwSeed) { setSeed(dwSeed); }

private:
	const uint N = 624;
	const uint M = 397;

	uint	m_dwMt[N];	// the array for the state vector
	int		m_nMti;			// initialization counter

	const uint MATRIX_A = 0x9908b0df;		/* constant vector a */
	const uint UPPER_MASK = 0x80000000;		/* most significant w-r bits */
	const uint LOWER_MASK = 0x7fffffff;		/* least significant r bits */

	/* Tempering parameters */	 
	const uint TEMPERING_MASK_B = 0x9d2c5680;
	const uint TEMPERING_MASK_C = 0xefc60000;
}
