00001 00029 #ifndef RANDOM_DSFMT_H 00030 #define RANDOM_DSFMT_H 00031 00032 #include <itpp/base/ittypes.h> 00033 #include <itpp/base/vec.h> 00034 #include <cstring> // required for memset() 00035 #include <ctime> 00036 #include <limits> 00037 00038 #if defined(__SSE2__) 00039 # include <emmintrin.h> 00040 #endif 00041 00042 namespace itpp 00043 { 00044 00092 template <int MEXP, int POS1, int SL1, uint64_t MSK1, uint64_t MSK2, 00093 uint32_t MSK32_1, uint32_t MSK32_2, 00094 uint32_t MSK32_3, uint32_t MSK32_4, 00095 uint64_t FIX1, uint64_t FIX2, uint64_t PCV1, uint64_t PCV2> 00096 class DSFMT { 00097 public: 00099 DSFMT() { if (!initialized) init_gen_rand(4257U); } 00101 DSFMT(unsigned int seed) { init_gen_rand(seed); } 00102 00104 void randomize() { init_gen_rand(hash(time(0), clock())); } 00106 void reset() { init_gen_rand(last_seed); } 00108 void reset(unsigned int seed) { init_gen_rand(seed); } 00109 00111 double random_01() { return genrand_open_open(); } 00113 double random_01_lclosed() { return genrand_close_open(); } 00115 double random_01_rclosed() { return genrand_open_close(); } 00117 uint32_t random_int() { return genrand_uint32(); } 00118 00120 ivec get_state() const { 00121 int size = (N + 1) * 4; 00122 uint32_t *psfmt = &status[0].u32[0]; 00123 ivec state(size + 1); // size + 1 to save idx variable in the same vec 00124 for (int i = 0; i < size; ++i) { 00125 state(i) = psfmt[i]; 00126 } 00127 state(size) = idx; 00128 return state; 00129 } 00130 00132 void set_state(const ivec &state) { 00133 int size = (N + 1) * 4; 00134 it_assert(state.size() == size + 1, "Random_Generator::set_state(): " 00135 "Invalid state initialization vector"); 00136 uint32_t *psfmt = &status[0].u32[0]; 00137 for (int i = 0; i < size; ++i) { 00138 psfmt[i] = state(i); 00139 } 00140 idx = state(size); 00141 } 00142 00150 void init_gen_rand(unsigned int seed) { 00151 uint32_t *psfmt = &status[0].u32[0]; 00152 psfmt[idxof(0)] = seed; 00153 for (int i = 1; i < (N + 1) * 4; i++) { 00154 psfmt[idxof(i)] = 1812433253UL 00155 * (psfmt[idxof(i - 1)] ^ (psfmt[idxof(i - 1)] >> 30)) + i; 00156 } 00157 initial_mask(); 00158 period_certification(); 00159 idx = Nx2; 00160 #if defined(__SSE2__) 00161 sse2_param_mask = _mm_set_epi32(MSK32_3, MSK32_4, MSK32_1, MSK32_2); 00162 #endif 00163 initialized = true; 00164 last_seed = seed; 00165 } 00166 00168 static uint32_t genrand_uint32() { 00169 uint64_t *psfmt64 = &status[0].u[0]; 00170 if (idx >= Nx2) { 00171 dsfmt_gen_rand_all(); 00172 idx = 0; 00173 } 00174 return psfmt64[idx++] & 0xffffffffU; 00175 } 00176 00186 static double genrand_close1_open2() { 00187 double *psfmt64 = &status[0].d[0]; 00188 if (idx >= Nx2) { 00189 dsfmt_gen_rand_all(); 00190 idx = 0; 00191 } 00192 return psfmt64[idx++]; 00193 } 00194 00203 static double genrand_close_open() { return genrand_close1_open2() - 1.0; } 00204 00213 static double genrand_open_close() { return 2.0 - genrand_close1_open2(); } 00214 00223 static double genrand_open_open() { 00224 double *dsfmt64 = &status[0].d[0]; 00225 union { 00226 double d; 00227 uint64_t u; 00228 } r; 00229 00230 if (idx >= Nx2) { 00231 dsfmt_gen_rand_all(); 00232 idx = 0; 00233 } 00234 r.d = dsfmt64[idx++]; 00235 r.u |= 1; 00236 return r.d - 1.0; 00237 } 00238 00239 00240 private: 00241 static const int N = (MEXP - 128) / 104 + 1; 00242 static const int Nx2 = N * 2; 00243 static const uint64_t LOW_MASK = 0x000fffffffffffffULL; 00244 static const uint64_t HIGH_CONST = 0x3ff0000000000000ULL; 00245 static const unsigned int SR = 12U; 00246 #if defined(__SSE2__) 00247 static const unsigned int SSE2_SHUFF = 0x1bU; 00248 #endif // __SSE2__ 00249 00251 union W128_T { 00252 #if defined(__SSE2__) 00253 __m128i si; 00254 __m128d sd; 00255 #endif // __SSE2__ 00256 uint64_t u[2]; 00257 uint32_t u32[4]; 00258 double d[2]; 00259 }; 00261 typedef union W128_T w128_t; 00263 static w128_t status[N + 1]; 00265 static int idx; 00267 static unsigned int last_seed; 00268 00270 static bool initialized; 00272 static bool bigendian; 00273 00274 #if defined(__SSE2__) 00275 00276 static __m128i sse2_param_mask; 00277 #endif // __SSE2__ 00278 00285 static unsigned int hash(time_t t, clock_t c) 00286 { 00287 static unsigned int differ = 0; // guarantee time-based seeds will change 00288 00289 unsigned int h1 = 0; 00290 unsigned char *p = (unsigned char *) &t; 00291 for (size_t i = 0; i < sizeof(t); ++i) { 00292 h1 *= std::numeric_limits<unsigned char>::max() + 2U; 00293 h1 += p[i]; 00294 } 00295 unsigned int h2 = 0; 00296 p = (unsigned char *) &c; 00297 for (size_t j = 0; j < sizeof(c); ++j) { 00298 h2 *= std::numeric_limits<unsigned char>::max() + 2U; 00299 h2 += p[j]; 00300 } 00301 return (h1 + differ++) ^ h2; 00302 } 00303 00308 static int idxof(int i) { return (bigendian ? (i ^ 1) : i); } 00309 00314 static void initial_mask() { 00315 uint64_t *psfmt = &status[0].u[0]; 00316 for (int i = 0; i < Nx2; i++) { 00317 psfmt[i] = (psfmt[i] & LOW_MASK) | HIGH_CONST; 00318 } 00319 } 00320 00322 static void period_certification() { 00323 uint64_t pcv[2] = {PCV1, PCV2}; 00324 uint64_t tmp[2]; 00325 uint64_t inner; 00326 00327 tmp[0] = (status[N].u[0] ^ FIX1); 00328 tmp[1] = (status[N].u[1] ^ FIX2); 00329 00330 inner = tmp[0] & pcv[0]; 00331 inner ^= tmp[1] & pcv[1]; 00332 for (int i = 32; i > 0; i >>= 1) { 00333 inner ^= inner >> i; 00334 } 00335 inner &= 1; 00336 /* check OK */ 00337 if (inner == 1) { 00338 return; 00339 } 00340 /* check NG, and modification */ 00341 #if (PCV2 & 1) == 1 00342 status[N].u[1] ^= 1; 00343 #else 00344 uint64_t work; 00345 for (int i = 1; i >= 0; i--) { 00346 work = 1; 00347 for (int j = 0; j < 64; j++) { 00348 if ((work & pcv[i]) != 0) { 00349 status[N].u[i] ^= work; 00350 return; 00351 } 00352 work = work << 1; 00353 } 00354 } 00355 #endif // (PCV2 & 1) == 1 00356 return; 00357 } 00358 00366 static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *lung) { 00367 #if defined(__SSE2__) 00368 __m128i x = a->si; 00369 __m128i z = _mm_slli_epi64(x, SL1); 00370 __m128i y = _mm_shuffle_epi32(lung->si, SSE2_SHUFF); 00371 z = _mm_xor_si128(z, b->si); 00372 y = _mm_xor_si128(y, z); 00373 00374 __m128i v = _mm_srli_epi64(y, SR); 00375 __m128i w = _mm_and_si128(y, sse2_param_mask); 00376 v = _mm_xor_si128(v, x); 00377 v = _mm_xor_si128(v, w); 00378 r->si = v; 00379 lung->si = y; 00380 #else // standard C++ 00381 uint64_t t0 = a->u[0]; 00382 uint64_t t1 = a->u[1]; 00383 uint64_t L0 = lung->u[0]; 00384 uint64_t L1 = lung->u[1]; 00385 lung->u[0] = (t0 << SL1) ^ (L1 >> 32) ^ (L1 << 32) ^ b->u[0]; 00386 lung->u[1] = (t1 << SL1) ^ (L0 >> 32) ^ (L0 << 32) ^ b->u[1]; 00387 r->u[0] = (lung->u[0] >> SR) ^ (lung->u[0] & MSK1) ^ t0; 00388 r->u[1] = (lung->u[1] >> SR) ^ (lung->u[1] & MSK2) ^ t1; 00389 #endif // __SSE2__ 00390 } 00391 00396 static void dsfmt_gen_rand_all() { 00397 int i; 00398 w128_t lung = status[N]; 00399 do_recursion(&status[0], &status[0], &status[POS1], &lung); 00400 for (i = 1; i < N - POS1; i++) { 00401 do_recursion(&status[i], &status[i], &status[i + POS1], &lung); 00402 } 00403 for (; i < N; i++) { 00404 do_recursion(&status[i], &status[i], &status[i + POS1 - N], &lung); 00405 } 00406 status[N] = lung; 00407 } 00408 00409 }; 00410 00411 00412 // ---------------------------------------------------------------------- 00413 // typedefs of different RNG 00414 // ---------------------------------------------------------------------- 00415 00416 typedef DSFMT<521, 3, 25, 00417 0x000fbfefff77efffULL, 0x000ffeebfbdfbfdfULL, 00418 0x000fbfefU, 0xff77efffU, 0x000ffeebU, 0xfbdfbfdfU, 00419 0xcfb393d661638469ULL, 0xc166867883ae2adbULL, 00420 0xccaa588000000000ULL, 0x0000000000000001ULL> DSFMT_521_RNG; 00421 00422 typedef DSFMT<1279, 9, 19, 00423 0x000efff7ffddffeeULL, 0x000fbffffff77fffULL, 00424 0x000efff7U, 0xffddffeeU, 0x000fbfffU, 0xfff77fffU, 00425 0xb66627623d1a31beULL, 0x04b6c51147b6109bULL, 00426 0x7049f2da382a6aebULL, 0xde4ca84a40000001ULL> DSFMT_1279_RNG; 00427 00428 typedef DSFMT<2203, 7, 19, 00429 0x000fdffff5edbfffULL, 0x000f77fffffffbfeULL, 00430 0x000fdfffU, 0xf5edbfffU, 0x000f77ffU, 0xfffffbfeU, 00431 0xb14e907a39338485ULL, 0xf98f0735c637ef90ULL, 00432 0x8000000000000000ULL, 0x0000000000000001ULL> DSFMT_2203_RNG; 00433 00434 typedef DSFMT<4253, 19, 19, 00435 0x0007b7fffef5feffULL, 0x000ffdffeffefbfcULL, 00436 0x0007b7ffU, 0xfef5feffU, 0x000ffdffU, 0xeffefbfcU, 00437 0x80901b5fd7a11c65ULL, 0x5a63ff0e7cb0ba74ULL, 00438 0x1ad277be12000000ULL, 0x0000000000000001ULL> DSFMT_4253_RNG; 00439 00440 typedef DSFMT<11213, 37, 19, 00441 0x000ffffffdf7fffdULL, 0x000ffffffdf7fffdULL, 00442 0x000fffffU, 0xfdf7fffdU, 0x000dffffU, 0xfff6bfffU, 00443 0xd0ef7b7c75b06793ULL, 0x9c50ff4caae0a641ULL, 00444 0x8234c51207c80000ULL, 0x0000000000000001ULL> DSFMT_11213_RNG; 00445 00446 typedef DSFMT<19937, 117, 19, 00447 0x000ffafffffffb3fULL, 0x000ffdfffc90fffdULL, 00448 0x000ffaffU, 0xfffffb3fU, 0x000ffdffU, 0xfc90fffdU, 00449 0x90014964b32f4329ULL, 0x3b8d12ac548a7c7aULL, 00450 0x3d84e1ac0dc82880ULL, 0x0000000000000001ULL> DSFMT_19937_RNG; 00451 00452 00453 } // namespace itpp 00454 00455 #endif // #ifndef RANDOM_DSFMT_H
Generated on Sat Jul 9 2011 15:21:30 for IT++ by Doxygen 1.7.4