IT++ Logo
random_dsfmt.h
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4