29 #ifndef RANDOM_DSFMT_H 30 #define RANDOM_DSFMT_H 37 #include <itpp/itexports.h> 40 # include <emmintrin.h> 46 namespace random_details
96 template <
int MEXP,
int POS1,
int SL1, uint64_t MSK1, uint64_t MSK2,
97 uint64_t FIX1_V, uint64_t FIX2_V, uint64_t PCV1_V, uint64_t PCV2_V >
103 static const int N = (MEXP - 128) / 104 + 1;
104 static const uint64_t
FIX1 = FIX1_V;
105 static const uint64_t FIX2 = FIX2_V;
106 static const uint64_t PCV1 = PCV1_V;
107 static const uint64_t PCV2 = PCV2_V;
109 #if defined(__SSE2__) 110 static const uint32_t MSK32_1 =
static_cast<uint32_t
>((MSK1 >> 32) & (0xffffffffULL));
111 static const uint32_t MSK32_2 =
static_cast<uint32_t
>(MSK1 & (0xffffffffULL));
112 static const uint32_t MSK32_3 =
static_cast<uint32_t
>((MSK2 >> 32) & (0xffffffffULL));
113 static const uint32_t MSK32_4 =
static_cast<uint32_t
>(MSK2 & (0xffffffffULL));
125 #if defined(__SSE2__) 155 uint32_t *psfmt = &_context.status[0].u32[0];
156 psfmt[idxof(0)] = seed;
157 for(
int i = 1; i < (N + 1) * 4; i++) {
158 psfmt[idxof(i)] = 1812433253UL
159 * (psfmt[idxof(i - 1)] ^ (psfmt[idxof(i - 1)] >> 30)) + i;
162 period_certification();
164 _context.last_seed = seed;
169 uint64_t *psfmt64 = &_context.status[0].u[0];
170 if(_context.idx >= Nx2) {
171 dsfmt_gen_rand_all();
174 return (uint32_t)(psfmt64[_context.idx++] & 0xffffffffU);
187 double *psfmt64 = &_context.status[0].d[0];
188 if(_context.idx >= Nx2) {
189 dsfmt_gen_rand_all();
192 return psfmt64[_context.idx++];
204 double *dsfmt64 = &_context.status[0].d[0];
210 if(_context.idx >= Nx2) {
211 dsfmt_gen_rand_all();
214 r.d = dsfmt64[_context.idx++];
220 static const int Nx2 = N * 2;
221 static const unsigned int SR = 12U;
223 static const bool bigendian;
225 #if defined(__SSE2__) 226 static const __m128i sse2_param_mask;
238 static int idxof(
int i) {
return (bigendian ? (i ^ 1) : i); }
244 void initial_mask() {
246 const uint64_t LOW_MASK = 0x000fffffffffffffULL;
247 const uint64_t HIGH_CONST = 0x3ff0000000000000ULL;
249 uint64_t *psfmt = &_context.
status[0].u[0];
250 for(
int i = 0; i < Nx2; i++) {
251 psfmt[i] = (psfmt[i] & LOW_MASK) | HIGH_CONST;
256 void period_certification() {
257 uint64_t pcv[2] = {PCV1, PCV2};
262 tmp[1] = (_context.
status[N].u[1] ^ FIX2);
264 inner = tmp[0] & pcv[0];
265 inner ^= tmp[1] & pcv[1];
266 for(
int i = 32; i > 0; i >>= 1) {
276 _context.
status[N].u[1] ^= 1;
279 for(
int i = 1; i >= 0; i--) {
281 for(
int j = 0; j < 64; j++) {
282 if((work & pcv[i]) != 0) {
283 _context.
status[N].u[i] ^= work;
289 #endif // (PCV2 & 1) == 1 301 #if defined(__SSE2__) 302 #define SSE2_SHUFF 0x1bU 305 __m128i z = _mm_slli_epi64(x, SL1);
306 __m128i y = _mm_shuffle_epi32(lung->si, SSE2_SHUFF);
307 z = _mm_xor_si128(z, b->si);
308 y = _mm_xor_si128(y, z);
310 __m128i v = _mm_srli_epi64(y, SR);
311 __m128i w = _mm_and_si128(y, sse2_param_mask);
312 v = _mm_xor_si128(v, x);
313 v = _mm_xor_si128(v, w);
316 #else // standard C++ 317 uint64_t t0 = a->u[0];
318 uint64_t t1 = a->u[1];
319 uint64_t L0 = lung->u[0];
320 uint64_t L1 = lung->u[1];
321 lung->u[0] = (t0 << SL1) ^ (L1 >> 32) ^ (L1 << 32) ^ b->u[0];
322 lung->u[1] = (t1 << SL1) ^ (L0 >> 32) ^ (L0 << 32) ^ b->u[1];
323 r->u[0] = (lung->u[0] >> SR) ^ (lung->u[0] & MSK1) ^ t0;
324 r->u[1] = (lung->u[1] >> SR) ^ (lung->u[1] & MSK2) ^ t1;
332 void dsfmt_gen_rand_all() {
336 do_recursion(&status[0], &status[0], &status[POS1], &lung);
337 for(i = 1; i < N - POS1; i++) {
338 do_recursion(&status[i], &status[i], &status[i + POS1], &lung);
341 do_recursion(&status[i], &status[i], &status[i + POS1 - N], &lung);
353 typedef DSFMT < 521, 3, 25,
354 0x000fbfefff77efffULL, 0x000ffeebfbdfbfdfULL,
355 0xcfb393d661638469ULL, 0xc166867883ae2adbULL,
356 0xccaa588000000000ULL, 0x0000000000000001ULL >
DSFMT_521_RNG;
358 typedef DSFMT < 1279, 9, 19,
359 0x000efff7ffddffeeULL, 0x000fbffffff77fffULL,
360 0xb66627623d1a31beULL, 0x04b6c51147b6109bULL,
363 typedef DSFMT < 2203, 7, 19,
364 0x000fdffff5edbfffULL, 0x000f77fffffffbfeULL,
365 0xb14e907a39338485ULL, 0xf98f0735c637ef90ULL,
368 typedef DSFMT < 4253, 19, 19,
369 0x0007b7fffef5feffULL, 0x000ffdffeffefbfcULL,
370 0x80901b5fd7a11c65ULL, 0x5a63ff0e7cb0ba74ULL,
373 typedef DSFMT < 11213, 37, 19,
374 0x000ffffffdf7fffdULL, 0x000dfffffff6bfffULL,
375 0xd0ef7b7c75b06793ULL, 0x9c50ff4caae0a641ULL,
378 typedef DSFMT < 19937, 117, 19,
379 0x000ffafffffffb3fULL, 0x000ffdfffc90fffdULL,
380 0x90014964b32f4329ULL, 0x3b8d12ac548a7c7aULL,
422 #endif // #ifndef RANDOM_DSFMT_H w128_t status[N+1]
128-bit internal state array
int idx
State array indexing.
ITPP_EXPORT double sd(const vec &In1, const vec &In2)
Spectral distortion between two vectors, in dB.
Data structure to hold 128-bit values.
const Fix_Factory FIX1(1, TC, WRAP)
Fix_Factories for signed Fix/CFix with wrap-around (FIX1, FIX2, ..., FIX64)
C++ implementation of dSFMT random number generator.
void lc_mark_initialized()
Function to mark thread-local context as initialized.
double genrand_open_open()
Generate uniform (0, 1) double pseudorandom number.
bool lc_is_initialized()
Function to check if thread-local context is initialized.
DSFMT(Context &c)
Constructor using a certain context.
ActiveDSFMT::Context & lc_get()
Function to access thread-local context for random numbers generation.
unsigned int last_seed
Last known seed used to initialize context.
double genrand_close1_open2()
Generate uniform [1, 2) double pseudorandom number.
DSFMT_19937_RNG ActiveDSFMT
Active Generator for random (stochastic) sources.
void init_gen_rand(unsigned int seed)
Initialise the generator with a new seed.
uint32_t genrand_uint32()
Generate uniform [0, UINT_MAX) integer pseudorandom number.
Templated Vector Class Definitions.