37 namespace random_details
41 static ActiveDSFMT::Context thread_local_context;
42 #pragma omp threadprivate(thread_local_context) 44 static bool is_thread_local_context_initialized =
false;
45 #pragma omp threadprivate(is_thread_local_context_initialized) 49 return thread_local_context;
54 return is_thread_local_context_initialized;
59 is_thread_local_context_initialized =
true;
70 static unsigned int differ = 0;
73 unsigned char *p = (
unsigned char *) &t;
74 for(
size_t i = 0; i <
sizeof(t); ++i) {
79 p = (
unsigned char *) &c;
80 for(
size_t j = 0; j <
sizeof(c); ++j) {
84 return (h1 + differ++) ^ h2;
97 const __m128i ActiveDSFMT::sse2_param_mask = _mm_set_epi32(ActiveDSFMT::MSK32_3, ActiveDSFMT::MSK32_4, ActiveDSFMT::MSK32_1, ActiveDSFMT::MSK32_2);
106 class GlobalSeedProvider
108 static const unsigned int default_first_seed = 4257U;
112 GlobalSeedProvider(): _dsfmt(_c), _first_seed_given(false) {
113 _dsfmt.init_gen_rand(default_first_seed);
116 void set_seed(
unsigned int s) {_dsfmt.init_gen_rand(s);}
118 unsigned int get_seed() {
return _c.last_seed;}
120 void reset() {
if(_first_seed_given) _dsfmt.init_gen_rand(get_seed());}
122 void set_state(
const ivec& st) {
123 int size = (DSFMT::N + 1) * 4;
124 it_assert(st.size() ==
size + 1,
"GlobalSeedProvider::state(): " 125 "Invalid state initialization vector");
126 uint32_t *psfmt = &_c.status[0].u32[0];
127 for(
int i = 0; i <
size; ++i) {
128 psfmt[i] =
static_cast<uint32_t
>(st(i));
131 _first_seed_given =
true;
135 int size = (DSFMT::N + 1) * 4;
136 uint32_t *psfmt = &_c.status[0].u32[0];
137 ivec state(
size + 1);
138 for(
int i = 0; i <
size; ++i) {
139 state(i) =
static_cast<int>(psfmt[i]);
141 state(
size) = _c.idx;
145 void randomize() {_dsfmt.init_gen_rand(random_details::hash_time_to_seed(time(0), clock())); _first_seed_given =
true;}
147 unsigned int generate() {
148 if(_first_seed_given)
149 return _dsfmt.genrand_uint32();
154 _first_seed_given =
true;
155 return default_first_seed;
164 bool _first_seed_given;
169 GlobalSeedProvider& global_seed_provider()
171 static GlobalSeedProvider global_seed_provider_instance;
172 return global_seed_provider_instance;
180 global_seed_provider().set_seed(seed);
189 global_seed_provider().reset();
198 s = global_seed_provider().generate();
208 global_seed_provider().randomize();
217 state = global_seed_provider().get_state();
226 global_seed_provider().set_state(state);
257 dsfmt.
init_gen_rand(random_details::hash_time_to_seed(time(0), clock()));
264 int size = (random_details::ActiveDSFMT::N + 1) * 4;
266 state.set_size(
size + 1);
267 for(
int i = 0; i <
size; ++i) {
268 state(i) =
static_cast<int>(psfmt[i]);
276 int size = (random_details::ActiveDSFMT::N + 1) * 4;
278 "Invalid state initialization vector");
280 for(
int i = 0; i <
size; ++i) {
281 psfmt[i] =
static_cast<uint32_t
>(state(i));
317 for(
int i = 0; i < n; i++)
328 for(i = 0; i < h; i++)
329 for(j = 0; j < w; j++)
375 for(
int i = 0; i < n; i++)
386 for(i = 0; i < h; i++)
387 for(j = 0; j < w; j++)
396 void Gamma_RNG::init_state()
398 const static double sqrt32 = 5.656854;
400 const static double q1 = 0.04166669;
401 const static double q2 = 0.02083148;
402 const static double q3 = 0.00801191;
403 const static double q4 = 0.00144121;
404 const static double q5 = -7.388e-5;
405 const static double q6 = 2.4511e-4;
406 const static double q7 = 2.424e-4;
408 double r = 1.0 / alpha;
411 it_error_if(!std::isfinite(alpha) || !std::isfinite(_scale) || (alpha < 0.0)
412 || (_scale <= 0.0),
"Gamma_RNG::init_state() - wrong parameters");
416 _d = sqrt32 - _s * 12.0;
418 _q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
425 _b = 0.463 + _s + 0.178 * _s2;
427 _c = 0.195 / _s - 0.079 + 0.16 * _s;
429 else if(alpha <= 13.022) {
430 _b = 1.654 + 0.0076 * _s2;
431 _si = 1.68 / _s + 0.275;
432 _c = 0.062 / _s + 0.024;
444 for(
int i = 0; i < n; i++)
452 for(
int i = 0; i < r * c; i++)
463 const static double exp_m1 = 0.36787944117144232159;
465 const static double a1 = 0.3333333;
466 const static double a2 = -0.250003;
467 const static double a3 = 0.2000062;
468 const static double a4 = -0.1662921;
469 const static double a5 = 0.1423657;
470 const static double a6 = -0.1367177;
471 const static double a7 = 0.1233795;
473 double e, p, q, t, u, v, w, x, ret_val;
475 double scale = _scale;
480 e = 1.0 + exp_m1 * a;
505 return scale * ret_val;
509 if((_d * u) <= (t * t * t))
510 return scale * ret_val;
517 if(std::fabs(v) <= 0.25)
518 q = _q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v
519 + a3) * v + a2) * v + a1) * v;
521 q = _q0 - _s * t + 0.25 * t * t + (_s2 + _s2) *
log(1.0 + v);
524 if(
log(1.0 - u) <= q)
525 return scale * ret_val;
540 if(t >= -0.71874483771719) {
543 if(std::fabs(v) <= 0.25)
544 q = _q0 + 0.5 * t * t *
545 ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
548 q = _q0 - _s * t + 0.25 * t * t + (_s2 + _s2) *
log(1.0 + v);
556 if((_c * std::fabs(u)) <= (w *
std::exp(e - 0.5 * t * t)))
562 return scale * x * x;
576 const double Normal_RNG::ytab[128] = {
577 1, 0.963598623011, 0.936280813353, 0.913041104253,
578 0.892278506696, 0.873239356919, 0.855496407634, 0.838778928349,
579 0.822902083699, 0.807732738234, 0.793171045519, 0.779139726505,
580 0.765577436082, 0.752434456248, 0.739669787677, 0.727249120285,
581 0.715143377413, 0.703327646455, 0.691780377035, 0.68048276891,
582 0.669418297233, 0.65857233912, 0.647931876189, 0.637485254896,
583 0.62722199145, 0.617132611532, 0.607208517467, 0.597441877296,
584 0.587825531465, 0.578352913803, 0.569017984198, 0.559815170911,
585 0.550739320877, 0.541785656682, 0.532949739145, 0.524227434628,
586 0.515614886373, 0.507108489253, 0.498704867478, 0.490400854812,
587 0.482193476986, 0.47407993601, 0.466057596125, 0.458123971214,
588 0.450276713467, 0.442513603171, 0.434832539473, 0.427231532022,
589 0.419708693379, 0.41226223212, 0.404890446548, 0.397591718955,
590 0.390364510382, 0.383207355816, 0.376118859788, 0.369097692334,
591 0.362142585282, 0.355252328834, 0.348425768415, 0.341661801776,
592 0.334959376311, 0.328317486588, 0.321735172063, 0.31521151497,
593 0.308745638367, 0.302336704338, 0.29598391232, 0.289686497571,
594 0.283443729739, 0.27725491156, 0.271119377649, 0.265036493387,
595 0.259005653912, 0.253026283183, 0.247097833139, 0.241219782932,
596 0.235391638239, 0.229612930649, 0.223883217122, 0.218202079518,
597 0.212569124201, 0.206983981709, 0.201446306496, 0.195955776745,
598 0.190512094256, 0.185114984406, 0.179764196185, 0.174459502324,
599 0.169200699492, 0.1639876086, 0.158820075195, 0.153697969964,
600 0.148621189348, 0.143589656295, 0.138603321143, 0.133662162669,
601 0.128766189309, 0.123915440582, 0.119109988745, 0.114349940703,
602 0.10963544023, 0.104966670533, 0.100343857232, 0.0957672718266,
603 0.0912372357329, 0.0867541250127, 0.082318375932, 0.0779304915295,
604 0.0735910494266, 0.0693007111742, 0.065060233529, 0.0608704821745,
605 0.056732448584, 0.05264727098, 0.0486162607163, 0.0446409359769,
606 0.0407230655415, 0.0368647267386, 0.0330683839378, 0.0293369977411,
607 0.0256741818288, 0.0220844372634, 0.0185735200577, 0.0151490552854,
608 0.0118216532614, 0.00860719483079, 0.00553245272614, 0.00265435214565
615 const unsigned int Normal_RNG::ktab[128] = {
616 0, 12590644, 14272653, 14988939,
617 15384584, 15635009, 15807561, 15933577,
618 16029594, 16105155, 16166147, 16216399,
619 16258508, 16294295, 16325078, 16351831,
620 16375291, 16396026, 16414479, 16431002,
621 16445880, 16459343, 16471578, 16482744,
622 16492970, 16502368, 16511031, 16519039,
623 16526459, 16533352, 16539769, 16545755,
624 16551348, 16556584, 16561493, 16566101,
625 16570433, 16574511, 16578353, 16581977,
626 16585398, 16588629, 16591685, 16594575,
627 16597311, 16599901, 16602354, 16604679,
628 16606881, 16608968, 16610945, 16612818,
629 16614592, 16616272, 16617861, 16619363,
630 16620782, 16622121, 16623383, 16624570,
631 16625685, 16626730, 16627708, 16628619,
632 16629465, 16630248, 16630969, 16631628,
633 16632228, 16632768, 16633248, 16633671,
634 16634034, 16634340, 16634586, 16634774,
635 16634903, 16634972, 16634980, 16634926,
636 16634810, 16634628, 16634381, 16634066,
637 16633680, 16633222, 16632688, 16632075,
638 16631380, 16630598, 16629726, 16628757,
639 16627686, 16626507, 16625212, 16623794,
640 16622243, 16620548, 16618698, 16616679,
641 16614476, 16612071, 16609444, 16606571,
642 16603425, 16599973, 16596178, 16591995,
643 16587369, 16582237, 16576520, 16570120,
644 16562917, 16554758, 16545450, 16534739,
645 16522287, 16507638, 16490152, 16468907,
646 16442518, 16408804, 16364095, 16301683,
647 16207738, 16047994, 15704248, 15472926
651 const double Normal_RNG::wtab[128] = {
652 1.62318314817e-08, 2.16291505214e-08, 2.54246305087e-08, 2.84579525938e-08,
653 3.10340022482e-08, 3.33011726243e-08, 3.53439060345e-08, 3.72152672658e-08,
654 3.8950989572e-08, 4.05763964764e-08, 4.21101548915e-08, 4.35664624904e-08,
655 4.49563968336e-08, 4.62887864029e-08, 4.75707945735e-08, 4.88083237257e-08,
656 5.00063025384e-08, 5.11688950428e-08, 5.22996558616e-08, 5.34016475624e-08,
657 5.44775307871e-08, 5.55296344581e-08, 5.65600111659e-08, 5.75704813695e-08,
658 5.85626690412e-08, 5.95380306862e-08, 6.04978791776e-08, 6.14434034901e-08,
659 6.23756851626e-08, 6.32957121259e-08, 6.42043903937e-08, 6.51025540077e-08,
660 6.59909735447e-08, 6.68703634341e-08, 6.77413882848e-08, 6.8604668381e-08,
661 6.94607844804e-08, 7.03102820203e-08, 7.11536748229e-08, 7.1991448372e-08,
662 7.2824062723e-08, 7.36519550992e-08, 7.44755422158e-08, 7.52952223703e-08,
663 7.61113773308e-08, 7.69243740467e-08, 7.77345662086e-08, 7.85422956743e-08,
664 7.93478937793e-08, 8.01516825471e-08, 8.09539758128e-08, 8.17550802699e-08,
665 8.25552964535e-08, 8.33549196661e-08, 8.41542408569e-08, 8.49535474601e-08,
666 8.57531242006e-08, 8.65532538723e-08, 8.73542180955e-08, 8.8156298059e-08,
667 8.89597752521e-08, 8.97649321908e-08, 9.05720531451e-08, 9.138142487e-08,
668 9.21933373471e-08, 9.30080845407e-08, 9.38259651738e-08, 9.46472835298e-08,
669 9.54723502847e-08, 9.63014833769e-08, 9.71350089201e-08, 9.79732621669e-08,
670 9.88165885297e-08, 9.96653446693e-08, 1.00519899658e-07, 1.0138063623e-07,
671 1.02247952126e-07, 1.03122261554e-07, 1.04003996769e-07, 1.04893609795e-07,
672 1.05791574313e-07, 1.06698387725e-07, 1.07614573423e-07, 1.08540683296e-07,
673 1.09477300508e-07, 1.1042504257e-07, 1.11384564771e-07, 1.12356564007e-07,
674 1.13341783071e-07, 1.14341015475e-07, 1.15355110887e-07, 1.16384981291e-07,
675 1.17431607977e-07, 1.18496049514e-07, 1.19579450872e-07, 1.20683053909e-07,
676 1.21808209468e-07, 1.2295639141e-07, 1.24129212952e-07, 1.25328445797e-07,
677 1.26556042658e-07, 1.27814163916e-07, 1.29105209375e-07, 1.30431856341e-07,
678 1.31797105598e-07, 1.3320433736e-07, 1.34657379914e-07, 1.36160594606e-07,
679 1.37718982103e-07, 1.39338316679e-07, 1.41025317971e-07, 1.42787873535e-07,
680 1.44635331499e-07, 1.4657889173e-07, 1.48632138436e-07, 1.50811780719e-07,
681 1.53138707402e-07, 1.55639532047e-07, 1.58348931426e-07, 1.61313325908e-07,
682 1.64596952856e-07, 1.68292495203e-07, 1.72541128694e-07, 1.77574279496e-07,
683 1.83813550477e-07, 1.92166040885e-07, 2.05295471952e-07, 2.22600839893e-07
687 const double Normal_RNG::PARAM_R = 3.44428647676;
692 uint32_t u,
sign, i, j;
717 return sign ? x : -x;
749 for(
int i = 0; i < n; i++)
760 for(i = 0; i < h; i++)
761 for(j = 0; j < w; j++)
781 factr = -2.0 * var * (1.0 - rho * rho);
798 for(
int i = 0; i < n; i++)
809 for(i = 0; i < h; i++)
810 for(j = 0; j < w; j++)
834 mean = tgamma(1.0 + 1.0 / b) / l;
835 var = tgamma(1.0 + 2.0 / b) / (l * l) - mean;
843 for(
int i = 0; i < n; i++)
854 for(i = 0; i < h; i++)
855 for(j = 0; j < w; j++)
874 for(
int i = 0; i < n; i++)
885 for(i = 0; i < h; i++)
886 for(j = 0; j < w; j++)
905 for(
int i = 0; i < n; i++)
916 for(i = 0; i < h; i++)
917 for(j = 0; j < w; j++)
#define it_error_if(t, s)
Abort if t is true.
unsigned int GlobalRNG_get_local_seed()
Get new seed to initialize thread-local generators.
double sample()
Returns a single sample.
bool is_bigendian()
Returns true if machine endianness is BIG_ENDIAN.
Weibull_RNG(double lambda=1.0, double beta=1.0)
Constructor. Set lambda and beta.
w128_t status[N+1]
128-bit internal state array
void GlobalRNG_set_state(const ivec &state)
Resume the global seed provider state saved in memory.
void setup(double meanval, double variance)
Set mean and variance.
int idx
State array indexing.
void setup(double lambda)
Set lambda.
double sample()
Get a sample.
double genrand_open_open()
Generate uniform (0, 1) double pseudorandom number.
Laplace_RNG(double meanval=0.0, double variance=1.0)
Constructor. Set mean and variance.
double genrand_open_close()
Generate uniform (0, 1] double pseudorandom number.
double operator()()
Get one sample.
void setup(double lambda, double beta)
Set lambda, and beta.
double sample()
Get a Normal distributed (0,1) sample.
#define it_assert(t, s)
Abort if t is not true.
void get_setup(double &meanval, double &variance) const
Get mean and variance.
void GlobalRNG_reset(unsigned int seed)
Set the internal seed of the Global Seed Provider.
C++ implementation of dSFMT random number generator.
double operator()()
Get one sample.
vec log(const vec &x)
The natural logarithm of the elements.
Exponential_RNG(double lambda=1.0)
constructor. Set lambda.
Rayleigh_RNG(double sigma=1.0)
Constructor. Set sigma.
double variance(const cvec &v)
The variance of the elements in the vector. Normalized with N-1 to be unbiased.
Rice_RNG(double sigma=1.0, double v=1.0)
Constructor. Set sigma, and v (if v = 0, Rice -> Rayleigh).
T min(const Vec< T > &in)
Minimum value of vector.
void get_setup(double &meanval, double &variance) const
Get mean and variance.
double operator()()
Get one sample.
uint32_t genrand_uint32()
Generate uniform [0, UINT_MAX) integer pseudorandom number.
vec exp(const vec &x)
Exp of the elements of a vector x.
static unsigned int hash_time_to_seed(time_t t, clock_t c)
Get an unsigned int from time variables t and c.
void RNG_randomize()
Set a random seed for all Random Number Generators in the current thread.
T max(const Vec< T > &v)
Maximum value of vector.
AR1_Normal_RNG(double meanval=0.0, double variance=1.0, double rho=0.0)
Constructor. Set mean, variance, and correlation.
void lc_mark_initialized()
Function to mark thread-local context as initialized.
void RNG_reset(unsigned int seed)
Set the seed for all Random Number Generators in the current thread.
double operator()()
Get one sample.
int size(const Vec< T > &v)
Length of vector.
double operator()()
Get one sample.
void RNG_set_state(const ivec &state)
Resume Random Number generation in the current thread with previously stored context.
bool lc_is_initialized()
Function to check if thread-local context is initialized.
void get_setup(double &meanval, double &variance, double &rho) const
Get mean, variance and correlation.
IT++ compatibility types and functions.
void GlobalRNG_get_state(ivec &state)
Save current full state of global seed provider in memory.
ActiveDSFMT::Context & lc_get()
Function to access thread-local context for random numbers generation.
double operator()()
Get one sample.
void GlobalRNG_randomize()
Set a random seed for the Global Seed Provider seed.
void RNG_get_state(ivec &state)
Save Random Number generation context used in the current thread.
double sign(double x)
Signum function.
double genrand_close_open()
Generate uniform [0, 1) double pseudorandom number.
void reset()
Set memory contents to zero.
vec sqrt(const vec &x)
Square root of the elements.
void setup(double sigma, double v)
Set sigma, and v (if v = 0, Rice -> Rayleigh).
Definition of classes for random number generators.
void setup(double sigma)
Set sigma.
DSFMT_19937_RNG ActiveDSFMT
Active Generator for random (stochastic) sources.
void init_gen_rand(unsigned int seed)
Initialise the generator with a new seed.
void setup(double meanval, double variance, double rho)
Set mean, variance, and correlation.
Elementary mathematical functions - header file.
double operator()()
Get a single random sample.