SourceForge.jp

dKingyoRandom.h

説明を見る。
00001 
00002 
00003 #ifndef _dKingyoRandom__
00004 #define _dKingyoRandom__
00005 
00006 //#include <dkutil/dkutilBuffer.h>
00007 #include "dkutilMath.h"
00008 #include "dKingyoRandom_MT_SSE2.h"
00009 #include "dKingyoRandom_MT_MMXASM.h"
00010 
00011 
00012 
00013 namespace dkutil{//begin dkutil namespace
00014     namespace oldmath{//begin math namespace
00015 
00016 template<typename SEED>
00017 inline int arg_rand(SEED *seed)
00018 {
00019     (*seed) = (*seed) * 1103515245L + 12345;
00020     return (unsigned)((*seed) / 65536L) % 32768U;
00021 }
00022 inline int arg_random(ULONG *seed,ULONG Max){
00023     return arg_rand<ULONG>(seed) * Max / SHRT_MAX;
00024 }
00025 
00026 
00027 template<typename T>
00028 class dKingyoRandom_ANSI : public dKingyoRandom_BasicMember<T>{
00029 public:
00030     dKingyoRandom_ANSI(T seed=timeGetTime())
00031     {       
00032         m_Max = 32767;//(RAMD_MAX)
00033         if(seed<0){seed = -seed;}
00034          m_ulSeed = m_Seed =seed;
00035 
00036     }
00037     virtual ~dKingyoRandom_ANSI(){
00038 #ifdef _DEBUG
00039         //ODS("しっかり呼ばれた");
00040 #endif
00041     }
00042     virtual T Rand(){
00043         //m_Seed = m_Seed * 1103515245L + 12345;
00044         //return (unsigned)(m_Seed / 65536L) % 32768U;
00045         return arg_rand(&m_ulSeed);
00046     }
00048     virtual T Random(T domain){
00049         return Rand() * (++domain) / m_Max;
00050     }
00052     virtual T RandomDomain(T min,T max){
00053 
00054         //return max += Random(max + min);
00055         //return max-=Random(max-min);
00056         //return (max-min)*Rand()+min;
00057         {//しかたない。うまくいかないからこれで代用
00058             //http://members.tripod.co.jp/cprogram/various1.html
00059             //min<= x <=maxまで返す
00060             float fNum, fDenom;
00061 
00062             fNum   = (float)Rand() * ( (float)max - (float)min + 1.0f );
00063             fDenom = (float)RAND_MAX + 1.0f;
00064 
00065             return min + (int)( fNum / fDenom );
00066         }
00067 
00068     }
00070     virtual T RandomDomainSafety(T min,T max){
00071         MINMAX_SWAP(T,min,max);
00072         MINMAX_ABS(min,max);
00073         MINMAX_SAFETY(min,max);
00074         return RandomDomain(min,max);
00075     }
00076     
00077     
00078     
00079 };
00080 template<class T>
00081 class dKingyoRandom_Float : public dKingyoRandom_BasicMember<T>{
00082 public:
00083     dKingyoRandom_Float(T seed = timeGetTime()){
00084         m_Max=32767;
00085         seed = (float)fabs(seed);
00086         //if(seed<0){seed = -seed;}
00087         m_Seed = seed;
00088         m_ulSeed=(ULONG) m_Seed;
00089     }
00090 
00091     virtual T Rand(){
00092         return (float)arg_rand(&m_ulSeed) / ( (float)m_Max + 1.0f );
00093     }
00094     virtual T Random(T fMax){
00095         return ( (float)arg_rand(&m_ulSeed) * fMax ) / ( (float)RAND_MAX + 1.0f );
00096     }
00097     //min<= x <maxまでの値を返す
00098     virtual T RandomDomain(T min,T max){
00099         return (max-min)*Rand()+min;
00100     }
00101     virtual T RandomDomainSafety(T min,T max){
00102         MINMAX_SWAP_FLOAT(T,min,max);
00103         MINMAX_ABS_FLOAT(T,min,max);
00104         MINMAX_SAFETY(min,max);
00105         return RandomDomain(min,max);
00106     }
00107 };
00108 
00109 
00110 
00111 
00112 /* Period parameters */
00113 #define N 624
00114 #define M 397
00115 #define MATRIX_A 0x9908b0df   /* constant vector a */
00116 #define UPPER_MASK 0x80000000 /* most significant w-r bits */
00117 #define LOWER_MASK 0x7fffffff /* least significant r bits */
00118 
00119 /* Tempering parameters */
00120 #define TEMPERING_MASK_B 0x9d2c5680
00121 #define TEMPERING_MASK_C 0xefc60000
00122 #define TEMPERING_SHIFT_U(y)  (y >> 11)
00123 #define TEMPERING_SHIFT_S(y)  (y << 7)
00124 #define TEMPERING_SHIFT_T(y)  (y << 15)
00125 #define TEMPERING_SHIFT_L(y)  (y >> 18)
00126 
00128 template<class T>
00129 class INL_dKingyoRandom_MT{
00130     UINT m_Max;
00131     unsigned long bInitialized ;
00132     // Static member 
00133     int mti;                   // index number 
00134     unsigned long mt[N + 1];   // the array for the state vector 
00135     unsigned long mtr[N];      // the array for the random number 
00136 public:
00137     INL_dKingyoRandom_MT(T seed=timeGetTime()){
00138         m_Max=2^32-1;
00139         SetSeed(seed);
00140         bInitialized=0;
00141     }
00142     virtual ~INL_dKingyoRandom_MT(){}
00143     UINT GetMax(){return m_Max;}
00144     void SetSeed(T seed){
00145         int i;
00146         for(i = 0; i < N; i++){
00147              mt[i] = seed & 0xffff0000;
00148              seed = 69069 * seed + 1;
00149              mt[i] |= (seed & 0xffff0000) >> 16;
00150              seed = 69069 * seed + 1;
00151         }
00152         bInitialized = 1;
00153         generateMT();
00154     }
00155     T Rand(){
00156         if(mti >= N){
00157             if(!bInitialized) SetSeed(4357);
00158             generateMT();
00159         }
00160         return mtr[mti++]; 
00161     }
00162     T Random(T max){return Rand() * max / m_Max;}
00163     //min<= x <maxまでの値を返す
00164     T RandomDomain(T min,T max){
00165         return (max-min)*Rand()+min;
00166     }
00167     T RandomDomainSafety(T min,T max){
00168         MINMAX_SWAP(T,min,max);
00169         MINMAX_ABS(min,max);
00170         MINMAX_SAFETY(min,max);
00171         return RandomDomain(min,max);
00172     }
00173     void generateMT(void)
00174     {
00175         int kk;
00176         unsigned long y;
00177         static unsigned long mag01[2] = {0x0, MATRIX_A}; /* mag01[x] = x * MATRIX_A  for x=0,1 */
00178     
00179         for(kk = 0; kk < N - M; kk++){
00180             y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
00181             mt[kk] = mt[kk + M] ^ (y >> 1) ^ mag01[y & 0x1];
00182         }
00183 
00184         mt[N] = mt[0];
00185 
00186         for(; kk < N; kk++){
00187             y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
00188             mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1];
00189         }
00190 
00191         for(kk = 0; kk < N; kk++){
00192             y = mt[kk];
00193             y ^= TEMPERING_SHIFT_U(y);
00194             y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
00195             y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
00196             y ^= TEMPERING_SHIFT_L(y);
00197             mtr[kk] = y;
00198         }
00199         mti = 0;
00200     }
00201 };
00202 
00203 template<class T>
00204 class dKingyoRandom_MT : public dKingyoRandom_BasicMember<T>{//Please long or ULONG
00205     INL_dKingyoRandom_MT<T> m_mt;
00206 public:
00207     dKingyoRandom_MT(T seed=timeGetTime()){
00208         m_Max = m_mt.GetMax();
00209         m_mt.SetSeed(seed);
00210     }
00211     virtual ~dKingyoRandom_MT(){}
00212 
00213     virtual void SetSeed(T seed){
00214         m_mt.SetSeed(seed);
00215     }
00216     virtual T Rand(){
00217         return m_mt.Rand();
00218     }
00219     virtual T Random(T max){return m_mt.Random(max);}
00220     //min<= x <maxまでの値を返す
00221     virtual T RandomDomain(T min,T max){
00222         return m_mt.RandomDomain(min,max);
00223     }
00224     virtual T RandomDomainSafety(T min,T max){
00225         return m_mt.RandomDomainSafety(min,max);
00226     }
00227 };
00228 
00230 class dKingyoRandom_ANSI_NO_SETUP : public dKingyoRandom_ANSI<int>{
00231 public: 
00232     dKingyoRandom_ANSI_NO_SETUP(){}
00233 };
00235 class dKingyoRandom_MT_NO_SETUP : public dKingyoRandom_MT<ULONG>{
00236 public:
00237     dKingyoRandom_MT_NO_SETUP(){}
00238 };
00240 class dKingyoRandom_Float_NO_SETUP : public dKingyoRandom_Float<float>{
00241 public:
00242     dKingyoRandom_Float_NO_SETUP(){}
00243 };
00244 
00245 #undef N
00246 #undef M
00247 #undef MATRIX_A
00248 #undef UPPER_MASK
00249 #undef LOWER_MASK
00250 
00251 #undef TEMPERING_MASK_B
00252 #undef TEMPERING_MASK_C
00253 #undef TEMPERING_SHIFT_U 
00254 #undef TEMPERING_SHIFT_S
00255 #undef TEMPERING_SHIFT_T
00256 #undef TEMPERING_SHIFT_L
00257 
00258 
00259 }//end of math namespace
00260 }//end of dkutil namespace
00261 
00262 
00263 #endif
00264 
00265 
00266 /* This main() outputs first 1000 generated numbers.  
00267 int main(int argc, char *argv[])
00268 { 
00269     int i;
00270 
00271     srandMT(4357);
00272     for (i=0; i<1000; i++) {
00273         printf("%10lu ", randMT());
00274         if (i%5==4) printf("\n");
00275     }
00276     return EXIT_SUCCESS;
00277 }
00278 */
00279 
00280 
00281 
00282 
00283 
00284 
00285 
00286 
00287 
00288 
00289 
00290 
00291 
00292 /*
00293 参考資料:
00294 http://www2t.biglobe.ne.jp/~bono/c/memo1.htm
00295 http://members.tripod.co.jp/cprogram/various1.html
00296 
00297 
00298 乱数の発生
00299 
00300 
00301 関数による発生 
00302 stdlib.hをインクルードして,rand()を用いる
00303 ・Cライブラリのrand関数は,直前のrand()の呼び出しで得られた乱数の値に適当な演算を施して,次の乱数を作成するようになっている.
00304 ・一番最初にrandを呼び出すとき,一般にsrand関数で乱数の初期値となる値を指定する.これを乱数の種という.
00305 ・乱数の種を適切に設定しないと,毎回同じ乱数を発生させることになる.
00306 そこで,現在の時刻を与えることで毎回違う乱数になるように設定する例
00307 
00308 srand((unsigned)time(NULL)); //乱数の初期化
00309 
00310 
00311 線形合同法による一様乱数 
00312 適当な初期値X0から始め
00313 
00314 Xn = ( A*Xn-1 + C ) mod M
00315 という式を用いて,次々に0〜Mの範囲の値を発生させる.
00316 Aは,8で割って余りが5の整数
00317 Cは,奇数
00318 Mは,2^n という条件で,0〜M−1 までの整数が周期Mで1回ずつ現れる. 
00319 
00320 
00321 実数乱数の発生 
00322 線形合同法を使って,0〜1未満の実数乱数を発生させる関数の例
00323 double rnd(void)
00324 {
00325     unsigned rndnum=13;
00326 
00327     rndnum=(rndnum*109+1021)%32768;
00328     return rndnum/32767.1;
00329 } 
00330 
00331 1〜Mまでの乱数を発生 
00332 r=(int)(M*rand()/32767.1+1);
00333 
00334 
00335 正規乱数(ボックス・ミュラー法) 
00336 平均m,標準偏差σの正規分布N(m,σ)に従う乱数の発生例
00337 
00338 void boxrnd(double m,double sig,double *x,double *y)
00339 {
00340     double r1,r2;
00341     r1=rand()/32767.1;
00342     r2=rand()/32767.1;
00343     *x=sig*sqrt(-2*log(r1))*cos(2*3.14159*r2)+m;
00344     *y=sig*sqrt(-2*log(r1))*sin(2*3.14159*r2)+m;
00345 }
00346 stdlib.hとmath.hのインクルードを忘れないように!  
00347 
00348 
00349 */

dkutil 1.02リリース前 d金魚専用マニュアルバージョンに対してSun Dec 28 21:23:07 2003に生成されました。 doxygen 1.3.5