SEO		[PR] 爆速!無料ブログ 無料ホームページ開設 無料ライブ放送

               

---

dKingyoRandom.h

解説を見る。

#ifndef _dKingyoRandom__
#define _dKingyoRandom__

#include "AllLoad.h"
#include "dkutilMath.h"
#include "dKingyoRandom_MT_SSE2.h"
#include "dKingyoRandom_MT_MMXASM.h"
#include <time.h>


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

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



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(){}
    virtual T Rand(){
        //m_Seed = m_Seed * 1103515245L + 12345;
        //return (unsigned)(m_Seed / 65536L) % 32768U;
        return arg_rand(&m_ulSeed);
    }
    
    virtual T Random(T domain){
        return Rand() * domain / m_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 );
        }

    }
    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;
        if(seed<0){seed = -seed;}
        m_ulSeed= m_Seed = 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)


template<class T>
class dKingyoRandom_MT : public dKingyoRandom_BasicMember<T>{//Please long or ULONG
    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:
    dKingyoRandom_MT(T seed=timeGetTime()){
        m_Max=2^32-1;
        SetSeed(seed);
        bInitialized=0;
    }
    virtual ~dKingyoRandom_MT(){}

    virtual 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();
    }
    virtual T Rand(){
        if(mti >= N){
            if(!bInitialized) SetSeed(4357);
            generateMT();
        }
        return mtr[mti++]; 
    }
    virtual T Random(T max){return Rand() * max / m_Max;}
    //min<= x <maxまでの値を返す
    virtual T RandomDomain(T min,T max){
        return (max-min)*Rand()+min;
    }
    virtual 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;
    }
};

class dKingyoRandom_ANSI_EX : public dKingyoRandom_ANSI<int>{
public: 
    dKingyoRandom_ANSI_EX(){}
};
class dKingyoRandom_MT_EX : public dKingyoRandom_MT<ULONG>{
public:
    dKingyoRandom_MT_EX(){}
};
class dKingyoRandom_Float_EX : public dKingyoRandom_Float<float>{
public:
    dKingyoRandom_Float_EX(){}
};


#ifdef USE_DKINGYO_OBJECT_TEST
inline void RandomTest(){
    ODS("BEGIN RANDOM TEST()\n");
    char s[256],s2[256],s3[256];
    
    dKingyoRandom_Float<float> f(100);
    dKingyoRandom_MT<int> m(100);
    dKingyoRandom_ANSI<short> r(100);
    int count=0;
    while(1){
        int inum=r.RandomDomain(100,200);//r.GetRandomDomain(100,200);
        float fnum=f.RandomDomain(100,200);
        int mtnum = m.Rand();
        //int i=r.Random(200);
        //int i=arg_rand(&tane) * 200 / RAND_MAX;
        //int i=(200-100)*arg_rand((int *)&tane)+100;
        sprintf(s,"%d\n",inum);
        sprintf(s2,"%f\n",fnum);
        sprintf(s3,"MT=%d\n",mtnum);
        if(100 > inum || 200 < inum){
            MB(s);
        }
        if(100 >= fnum || 200 <= fnum){
            MB(s2);
        }
        
        ODS(s);
        ODS(s2);
        ODS(s3);
        count++;
        if(count > 20)
            break;
    }
}
inline void MTRandomTest(){
    dkutil::math::srandMT(timeGetTime());
    const int loopnum=5;
    int setSeed=100;
    ULONG ulsetSeed=100;
    ODS("//**********************************************************\n");
    ODS("//メルセンヌツイスター法の乱数生成テスト\n");
    ODS("int type\n");
    {
        dKingyoRandom_MT<int> mtc(setSeed) ;
        dKingyoRandom_MT_SSE2<int> mtsse(setSeed);
        //mtsse.SetSeed(setSeed);
        //mtc.SetSeed(setSeed);
        for(int i=0;i<loopnum;i++){
            int sse,c;
            sse=mtsse.Rand();
            c=mtc.Rand();
            dkutil::dOutputDebugString("SSE2 rand num=%d ",sse);
            dkutil::dOutputDebugString("C type rand num=%d ",c);
                
        }
    }
    
    ODS("ulong type \n");
    {
        //::dkutil::math::dKingyoRandom_MT_EX mtc(ulsetSeed) ;
        //::dkutil::math::dKingyoRandom_MT_SSE2_EX mtsse(ulsetSeed);
        dKingyoRandom_MT_EX mtc ;
        dKingyoRandom_MT_SSE2_EX mtsse;
        mtsse.SetSeed(setSeed);
        mtc.SetSeed(setSeed);
        for(int i=0;i<loopnum;i++){
            ULONG sse,c;
            sse=mtsse.Rand();
            c=mtc.Rand();
            dkutil::dOutputDebugString("SSE2 rand num=%u ",sse);
            dkutil::dOutputDebugString("MMX type rand num=%d ",dkutil::math::randMT());
            dkutil::dOutputDebugString("C type rand num=%u ",c);
                        
        }
    }
    ODS("//**********************************************************\n");
    int loopcount=100000000;
    dkutil::dOutputDebugString("乱数生成速度テスト loop回数 = %d,出力単位=CPU clock",loopcount); 
    {
        dKingyoRandom_MT_EX mtc ;
        dKingyoRandom_MT_SSE2_EX mtsse;
        
        {//SSE2
            clock_t c=clock();
            for(int i=0;i<loopcount;i++){
                mtsse.Rand();           
            }
            dkutil::dOutputDebugString("SSE2 clock is %d",clock() - c);
        }
        {//SSE2 inline
            INL_dKingyoRandom_MT_SSE2<ULONG> inlmt;
            clock_t c=clock();
            for(int i=0;i<loopcount;i++){
                inlmt.Rand();       
            }
            dkutil::dOutputDebugString("SSE2 inline type clock is %d",clock() - c);
        }
        {//SSE2 original
            myMT_t *mt;
            mt = myMT_Init();
            myMT_SetSeed(mt, timeGetTime());
            clock_t c=clock();
            for(int i=0;i<loopcount;i++){
                myMT_GetInt32(mt);      
            }
            dkutil::dOutputDebugString("SSE2 original type clock is %d",clock() - c);
            myMT_Term(mt);
        }
        {//MMX
            clock_t c=clock();
            for(int i=0;i<loopcount;i++){
                dkutil::math::randMT();         
            }
            dkutil::dOutputDebugString("MMX2 clock is %d",clock() - c);
        }
        {//C Language
            clock_t c=clock();
            for(int i=0;i<loopcount;i++){
                mtc.Rand();         
            }
            dkutil::dOutputDebugString("C type clock is %d",clock() - c);
        }
    }




}


#endif //end of USE_DKINGYO_OBJECT_TEST


}//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のインクルードを忘れないように！  


*/

---