00001 
00059 #ifndef _MTRAND_H_
00060 #define _MTRAND_H_
00061 
00062 #include <climits>
00063 #include <ctime>
00064 
00065 #pragma warning (disable : 4146 )
00066 
00069 
00070 class MTRand {
00071 
00072   public:
00074     inline MTRand();
00075     inline MTRand(unsigned int oneSeed);
00076     inline MTRand(unsigned int bigSeed[], unsigned int seedLength = N);
00077 
00079     inline unsigned int randInt(); 
00080     inline unsigned int randInt(unsigned int n); 
00081     inline unsigned int operator()(unsigned int n); 
00082 
00084     inline float rand(); 
00085     inline float randExc(); 
00086     inline float randDblExc(); 
00087 
00089     inline double rand53(); 
00090 
00092     inline void seed();
00093     inline void seed(unsigned int oneSeed);
00094     inline void seed(unsigned int bigSeed[], unsigned int seedLength = N);
00095 
00096     const static int N = 624; 
00097 
00098   protected:
00099     const static int M = 397; 
00100 
00101     unsigned int state[N]; 
00102     unsigned int *pNext; 
00103     int left; 
00104 
00105     inline void initialize(unsigned int oneSeed);
00106     inline void reload();
00107     inline static unsigned int hash( time_t t, clock_t c );
00108 
00109     inline unsigned int hiBit(unsigned int u) const {
00110         return u & 0x80000000U;
00111     }
00112 
00113     inline unsigned int loBit(unsigned int u) const {
00114         return u & 0x00000001U;
00115     }
00116 
00117     inline unsigned int loBits(unsigned int u) const {
00118         return u & 0x7fffffffU;
00119     }
00120 
00121     inline unsigned int mixBits(unsigned int u, unsigned int v) const {
00122         return hiBit(u) | loBits(v);
00123     }
00124 
00125     inline unsigned int twist(unsigned int m, unsigned int s0, unsigned int s1) const {
00126         return m ^ (mixBits(s0, s1) >> 1) ^ (-loBit(s1) & 0x9908b0dfU);
00127     }
00128 
00129 };
00130 
00131 
00135 inline MTRand::MTRand() {
00136     seed();
00137 }
00138 
00139 inline MTRand::MTRand(unsigned int oneSeed) {
00140     seed(oneSeed);
00141 }
00142 
00143 inline MTRand::MTRand(unsigned int bigSeed[], unsigned int seedLength) {
00144     seed(bigSeed, seedLength);
00145 }
00146 
00147 
00156 inline unsigned int MTRand::randInt() {
00157     if (left == 0)
00158         reload();
00159     --left;
00160 
00161     register unsigned int s1;
00162     s1 = *pNext++;
00163     s1 ^= (s1 >> 11);
00164     s1 ^= (s1 << 7) & 0x9d2c5680U;
00165     s1 ^= (s1 << 15) & 0xefc60000U;
00166     return (s1 ^ (s1 >> 18));
00167 }
00168 
00172 inline unsigned int MTRand::randInt(unsigned int n) {
00173     
00174     unsigned int used = n;
00175     used |= used >> 1;
00176     used |= used >> 2;
00177     used |= used >> 4;
00178     used |= used >> 8;
00179     used |= used >> 16;
00180 
00181     
00182     unsigned int i;
00183     do {
00184         i = randInt() & used; 
00185     } while (i >= n);
00186     return i;
00187 }
00188 
00189 inline unsigned int MTRand::operator()(unsigned int n) {
00190     return randInt(n);
00191 }
00192 
00193 
00198 inline float MTRand::rand() {
00199     return static_cast<float>(randInt()) / 4294967295.0f;
00200 }
00201 
00202 inline float MTRand::randExc() {
00203     return static_cast<float>(randInt()) / 4294967296.0f;
00204 }
00205 
00206 inline float MTRand::randDblExc() {
00207     return (static_cast<float>(randInt()) + 0.5f) / 4294967296.0f;
00208 }
00209 
00213 inline double MTRand::rand53() {
00214     unsigned int a(randInt() >> 5), b(randInt() >> 6);
00215     return (a * 67108864.0 + b) / 9007199254740992.0;
00216 }
00217 
00218 
00223 inline void MTRand::seed() {
00224     seed(hash(time(NULL), clock()));
00225 }
00226 
00227 inline void MTRand::seed(unsigned int oneSeed) {
00228     initialize(oneSeed);
00229     reload();
00230 }
00231 
00238 inline void MTRand::seed(unsigned int bigSeed[], unsigned int seedLength) {
00239     initialize(19650218UL);
00240     register int i(1);
00241     register unsigned int j(0);
00242     register int k((N > static_cast<int>(seedLength)) ? N : static_cast<int>(seedLength));
00243 
00244     for ( ; k; --k) {
00245         state[i] = state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525U);
00246         state[i] += (bigSeed[j] & 0xffffffffU) + j;
00247         state[i] &= 0xffffffffU;
00248         ++i;
00249         ++j;
00250         if (i >= N) {
00251             state[0] = state[N - 1];
00252             i = 1;
00253         }
00254         if (j >= seedLength)
00255             j = 0;
00256     }
00257     for (k = (N - 1); k; --k) {
00258         state[i] = state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941U);
00259         state[i] -= i;
00260         state[i] &= 0xffffffffU;
00261         ++i;
00262         if (i >= N) {
00263             state[0] = state[N - 1];
00264             i = 1;
00265         }
00266     }
00267     state[0] = 0x80000000UL; 
00268     reload();
00269 }
00270 
00271 
00282 inline void MTRand::initialize(unsigned int seed) {
00283     register unsigned int *s(state);
00284     register unsigned int *r(state);
00285     register int i(1);
00286     *s++ = seed & 0xffffffffU;
00287     for ( ; i < N; ++i) {
00288         *s++ = (1812433253U * (*r ^ (*r >> 30)) + i) & 0xffffffffU;
00289         r++;
00290     }
00291 }
00292 
00297 inline void MTRand::reload() {
00298     register unsigned int *p(state);
00299     register int i;
00300     for (i = N - M; i--; ++p)
00301         *p = twist(p[M], p[0], p[1]);
00302     for (i = M; --i; ++p)
00303         *p = twist(p[M - N], p[0], p[1]);
00304     *p = twist(p[M - N], p[0], state[0] );
00305 
00306     left = N;
00307     pNext = state;
00308 }
00309 
00310 
00311 inline unsigned int MTRand::hash(time_t t, clock_t c) {
00312     
00313     static unsigned int differ(0);
00314 
00315     unsigned int h1(0);
00316     unsigned char *p(reinterpret_cast<unsigned char *>(&t));
00317     for (size_t i = 0; i < sizeof(t); ++i) {
00318         h1 *= UCHAR_MAX + 2U;
00319         h1 += p[i];
00320     }
00321 
00322     unsigned int h2(0);
00323     p = reinterpret_cast<unsigned char *>(&c);
00324     for (size_t j = 0; j < sizeof(c); ++j) {
00325         h2 *= UCHAR_MAX + 2U;
00326         h2 += p[j];
00327     }
00328 
00329     return (h1 + differ++) ^ h2;
00330 }
00331 
00332 #endif