

#ifndef _dKingyoRandom__
#define _dKingyoRandom__

//#include <dkutil/dkutilBuffer.h>
#include "dkutilMath.h"
#include "dKingyoRandom_MT_SSE2.h"
#include "dKingyoRandom_MT_MMXASM.h"



namespace dkutil{//begin dkutil namespace
	namespace oldmath{//begin math namespace

template<typename SEED>
inline int arg_rand(SEED *seed)
{
    (*seed) = (*seed) * 1103515245L + 12345;
    return (unsigned)((*seed) / 65536L) % 32768U;
}
inline int arg_random(ULONG *seed,ULONG Max){
	return arg_rand<ULONG>(seed) * Max / SHRT_MAX;
}



inline double toda(double p)
{
	double     retv, y, e, f, g, h, i, j, k, l, m, n, o;
   e =  0.6936233982e-12;
   f =  0.3657763036e-10;
   g = -0.3231081277e-08;
   h =  0.8360937017e-07;
   i = -0.1045274970e-05;
   j =  0.5824238515e-05;
   k =  0.6841218299e-05;
   l = -0.2250947176e-03;
   m = -0.8364353589e-03;
   n =  0.3706987906e-01;
   o =  1.570796288;
   y = -::log(4*p*(1-p));
   retv = ::sqrt(((((((((((y*e+f)*y+g)*y+h)*y+i)*y \
          +j)*y+k)*y+l)*y+m)*y+n)*y+o)*y);
   return (p < 0.5) ? -retv : retv;
}

template<typename T>
class dKingyoRandom_ANSI : public dKingyoRandom_BasicMember<T>{
public:
	dKingyoRandom_ANSI(T seed=timeGetTime())
	{		
		m_Max = 32767;//(RAMD_MAX)
		if(seed<0){seed = -seed;}
		 m_ulSeed = m_Seed =seed;

	}
	virtual ~dKingyoRandom_ANSI(){
#ifdef _DEBUG
		//ODS("しっかり呼ばれた");
#endif
	}
	virtual T Rand(){
		//m_Seed = m_Seed * 1103515245L + 12345;
		//return (unsigned)(m_Seed / 65536L) % 32768U;
		return arg_rand(&m_ulSeed);
	}
	///domainまでの範囲を返す 0 <= x <= domain @param domain[in] ランダムにしたい範囲
	virtual T Random(T domain){
		return Rand() * (++domain) / m_Max;
	}
	/// min<= x <=maxまで返す @param min[in] min値 @param max[in] max値
	virtual T RandomDomain(T min,T max){

		//return max += Random(max + min);
		//return max-=Random(max-min);
		//return (max-min)*Rand()+min;
		{//しかたない。うまくいかないからこれで代用
			//http://members.tripod.co.jp/cprogram/various1.html
			//min<= x <=maxまで返す
			float fNum, fDenom;

			fNum   = (float)Rand() * ( (float)max - (float)min + 1.0f );
			fDenom = (float)RAND_MAX + 1.0f;

			return min + (int)( fNum / fDenom );
		}

	}
	///min値とmax値を間違ってもOKなもの
	virtual T RandomDomainSafety(T min,T max){
		MINMAX_SWAP(T,min,max);
		MINMAX_ABS(min,max);
		MINMAX_SAFETY(min,max);
		return RandomDomain(min,max);
	}
	
	
	
};
template<class T>
class dKingyoRandom_Float : public dKingyoRandom_BasicMember<T>{
public:
	dKingyoRandom_Float(T seed = timeGetTime()){
		m_Max=32767;
		seed = (float)fabs(seed);
		//if(seed<0){seed = -seed;}
		m_Seed = seed;
		m_ulSeed=(ULONG) m_Seed;
	}

	virtual T Rand(){
		return (float)arg_rand(&m_ulSeed) / ( (float)m_Max + 1.0f );
	}
	virtual T Random(T fMax){
		return ( (float)arg_rand(&m_ulSeed) * fMax ) / ( (float)RAND_MAX + 1.0f );
	}
	//min<= x <maxまでの値を返す
	virtual T RandomDomain(T min,T max){
		return (max-min)*Rand()+min;
	}
	virtual T RandomDomainSafety(T min,T max){
		MINMAX_SWAP_FLOAT(T,min,max);
		MINMAX_ABS_FLOAT(T,min,max);
		MINMAX_SAFETY(min,max);
		return RandomDomain(min,max);
	}
};




/* Period parameters */
#define N 624
#define M 397
#define MATRIX_A 0x9908b0df   /* constant vector a */
#define UPPER_MASK 0x80000000 /* most significant w-r bits */
#define LOWER_MASK 0x7fffffff /* least significant r bits */

/* Tempering parameters */
#define TEMPERING_MASK_B 0x9d2c5680
#define TEMPERING_MASK_C 0xefc60000
#define TEMPERING_SHIFT_U(y)  (y >> 11)
#define TEMPERING_SHIFT_S(y)  (y << 7)
#define TEMPERING_SHIFT_T(y)  (y << 15)
#define TEMPERING_SHIFT_L(y)  (y >> 18)

///インライン化しているので早いMT法乱数
template<class T>
class INL_dKingyoRandom_MT{
	UINT m_Max;
	unsigned long bInitialized ;
	// Static member 
	int mti;                   // index number 
	unsigned long mt[N + 1];   // the array for the state vector 
	unsigned long mtr[N];      // the array for the random number 
public:
	INL_dKingyoRandom_MT(T seed=timeGetTime()){
		m_Max=2^32-1;
		SetSeed(seed);
		bInitialized=0;
	}
	virtual ~INL_dKingyoRandom_MT(){}
	UINT GetMax(){return m_Max;}
	void SetSeed(T seed){
		int i;
		for(i = 0; i < N; i++){
			 mt[i] = seed & 0xffff0000;
			 seed = 69069 * seed + 1;
			 mt[i] |= (seed & 0xffff0000) >> 16;
			 seed = 69069 * seed + 1;
		}
		bInitialized = 1;
		generateMT();
	}
	T Rand(){
		if(mti >= N){
			if(!bInitialized) SetSeed(4357);
			generateMT();
		}
		return mtr[mti++]; 
	}
	T Random(T max){
		max++;
		return ( T )( ( ( LONGLONG )Rand() * ( LONGLONG )max ) >> 32 ) ;
	}
	//min<= x <maxまでの値を返す
	T RandomDomain(T min,T max){
		return (max-min)*Rand()+min;
	}
	T RandomDomainSafety(T min,T max){
		MINMAX_SWAP(T,min,max);
		MINMAX_ABS(min,max);
		MINMAX_SAFETY(min,max);
		return RandomDomain(min,max);
	}
	void generateMT(void)
	{
		int kk;
		unsigned long y;
		static unsigned long mag01[2] = {0x0, MATRIX_A}; /* mag01[x] = x * MATRIX_A  for x=0,1 */
    
		for(kk = 0; kk < N - M; kk++){
			y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
			mt[kk] = mt[kk + M] ^ (y >> 1) ^ mag01[y & 0x1];
		}

		mt[N] = mt[0];

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

		for(kk = 0; kk < N; kk++){
			y = mt[kk];
			y ^= TEMPERING_SHIFT_U(y);
			y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
			y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
			y ^= TEMPERING_SHIFT_L(y);
			mtr[kk] = y;
		}
		mti = 0;
	}
};

template<class T>
class dKingyoRandom_MT : public dKingyoRandom_BasicMember<T>{//Please long or ULONG
	INL_dKingyoRandom_MT<T> m_mt;
public:
	dKingyoRandom_MT(T seed=timeGetTime()){
		m_Max = m_mt.GetMax();
		m_mt.SetSeed(seed);
	}
	virtual ~dKingyoRandom_MT(){}

	virtual void SetSeed(T seed){
		m_mt.SetSeed(seed);
	}
	virtual T Rand(){
		return m_mt.Rand();
	}
	virtual T Random(T max){return m_mt.Random(max);}
	//min<= x <maxまでの値を返す
	virtual T RandomDomain(T min,T max){
		return m_mt.RandomDomain(min,max);
	}
	virtual T RandomDomainSafety(T min,T max){
		return m_mt.RandomDomainSafety(min,max);
	}
};

///C言語のrandと同じ乱数を発生されるクラス
class dKingyoRandom_ANSI_NO_SETUP : public dKingyoRandom_ANSI<int>{
public:	
	dKingyoRandom_ANSI_NO_SETUP(){}
};
///MT法乱数で乱数を発生させるクラス
class dKingyoRandom_MT_NO_SETUP : public dKingyoRandom_MT<ULONG>{
public:
	dKingyoRandom_MT_NO_SETUP(){}
};
///浮動小数点（float)で乱数を発生させるクラス
class dKingyoRandom_Float_NO_SETUP : public dKingyoRandom_Float<float>{
public:
	dKingyoRandom_Float_NO_SETUP(){}
};

#undef N
#undef M
#undef MATRIX_A
#undef UPPER_MASK
#undef LOWER_MASK

#undef TEMPERING_MASK_B
#undef TEMPERING_MASK_C
#undef TEMPERING_SHIFT_U 
#undef TEMPERING_SHIFT_S
#undef TEMPERING_SHIFT_T
#undef TEMPERING_SHIFT_L


}//end of math namespace
}//end of dkutil namespace


#endif


/* This main() outputs first 1000 generated numbers.  
int main(int argc, char *argv[])
{ 
    int i;

    srandMT(4357);
    for (i=0; i<1000; i++) {
        printf("%10lu ", randMT());
        if (i%5==4) printf("\n");
    }
    return EXIT_SUCCESS;
}
*/













/*
参考資料：
http://www2t.biglobe.ne.jp/~bono/c/memo1.htm
http://members.tripod.co.jp/cprogram/various1.html


乱数の発生


関数による発生 
stdlib.hをインクルードして，rand()を用いる
・Ｃライブラリのrand関数は，直前のrand()の呼び出しで得られた乱数の値に適当な演算を施して，次の乱数を作成するようになっている．
・一番最初にrandを呼び出すとき，一般にsrand関数で乱数の初期値となる値を指定する．これを乱数の種という．
・乱数の種を適切に設定しないと，毎回同じ乱数を発生させることになる．
そこで，現在の時刻を与えることで毎回違う乱数になるように設定する例

srand((unsigned)time(NULL)); //乱数の初期化


線形合同法による一様乱数 
適当な初期値X0から始め

Xn = ( A*Xn-1 + C ) mod M
という式を用いて，次々に０〜Ｍの範囲の値を発生させる．
Aは，８で割って余りが５の整数
Cは，奇数
Mは，2^n　という条件で，０〜Ｍ−１　までの整数が周期Ｍで１回ずつ現れる． 


実数乱数の発生 
線形合同法を使って，０〜１未満の実数乱数を発生させる関数の例
double rnd(void)
{
    unsigned rndnum=13;

    rndnum=(rndnum*109+1021)%32768;
    return rndnum/32767.1;
} 

1〜Mまでの乱数を発生 
r=(int)(M*rand()/32767.1+1);


正規乱数（ボックス・ミュラー法） 
平均m,標準偏差σの正規分布N(m,σ)に従う乱数の発生例

void boxrnd(double m,double sig,double *x,double *y)
{
    double r1,r2;
    r1=rand()/32767.1;
    r2=rand()/32767.1;
    *x=sig*sqrt(-2*log(r1))*cos(2*3.14159*r2)+m;
    *y=sig*sqrt(-2*log(r1))*sin(2*3.14159*r2)+m;
}
stdlib.hとmath.hのインクルードを忘れないように！  


*/
