42 bool Random_Generator::initialized =
false;
46 Random_Generator::w128_t Random_Generator::status[N + 1] = { };
48 int Random_Generator::idx = 0;
50 unsigned int Random_Generator::last_seed = 0U;
53 __m128i Random_Generator::sse2_param_mask = _mm_set_epi32(0, 0, 0, 0);
123 for (
int i=0; i < n; i++)
134 for (i = 0; i < h; i++)
135 for (j = 0; j < w; j++)
181 for (
int i=0; i < n; i++)
192 for (i = 0; i < h; i++)
193 for (j = 0; j < w; j++)
206 for (
int i = 0; i < n; i++)
214 for (
int i = 0; i < r * c; i++)
225 const static double sqrt32 = 5.656854;
226 const static double exp_m1 = 0.36787944117144232159;
228 const static double q1 = 0.04166669;
229 const static double q2 = 0.02083148;
230 const static double q3 = 0.00801191;
231 const static double q4 = 0.00144121;
232 const static double q5 = -7.388e-5;
233 const static double q6 = 2.4511e-4;
234 const static double q7 = 2.424e-4;
236 const static double a1 = 0.3333333;
237 const static double a2 = -0.250003;
238 const static double a3 = 0.2000062;
239 const static double a4 = -0.1662921;
240 const static double a5 = 0.1423657;
241 const static double a6 = -0.1367177;
242 const static double a7 = 0.1233795;
245 static double aa = 0.;
246 static double aaa = 0.;
247 static double s, s2, d;
248 static double q0, b, si, c;
250 double e, p, q, r, t, u, v, w, x, ret_val;
252 double scale = 1.0 / beta;
254 it_error_if(!std::isfinite(a) || !std::isfinite(scale) || (a < 0.0)
255 || (scale <= 0.0),
"Gamma_RNG wrong parameters");
260 e = 1.0 + exp_m1 * a;
284 d = sqrt32 - s * 12.0;
292 return scale * ret_val;
296 if ((d * u) <= (t * t * t))
297 return scale * ret_val;
303 q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
310 b = 0.463 + s + 0.178 * s2;
312 c = 0.195 / s - 0.079 + 0.16 * s;
314 else if (a <= 13.022) {
315 b = 1.654 + 0.0076 * s2;
316 si = 1.68 / s + 0.275;
317 c = 0.062 / s + 0.024;
330 if (std::fabs(v) <= 0.25)
331 q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v
332 + a3) * v + a2) * v + a1) * v;
334 q = q0 - s * t + 0.25 * t * t + (s2 + s2) *
log(1.0 + v);
337 if (
log(1.0 - u) <= q)
338 return scale * ret_val;
353 if (t >= -0.71874483771719) {
356 if (std::fabs(v) <= 0.25)
357 q = q0 + 0.5 * t * t *
358 ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
361 q = q0 - s * t + 0.25 * t * t + (s2 + s2) *
log(1.0 + v);
369 if ((c * std::fabs(u)) <= (w *
std::exp(e - 0.5 * t * t)))
375 return scale * x * x;
385 variance = sigma * sigma;
389 const double Normal_RNG::ytab[128] = {
390 1, 0.963598623011, 0.936280813353, 0.913041104253,
391 0.892278506696, 0.873239356919, 0.855496407634, 0.838778928349,
392 0.822902083699, 0.807732738234, 0.793171045519, 0.779139726505,
393 0.765577436082, 0.752434456248, 0.739669787677, 0.727249120285,
394 0.715143377413, 0.703327646455, 0.691780377035, 0.68048276891,
395 0.669418297233, 0.65857233912, 0.647931876189, 0.637485254896,
396 0.62722199145, 0.617132611532, 0.607208517467, 0.597441877296,
397 0.587825531465, 0.578352913803, 0.569017984198, 0.559815170911,
398 0.550739320877, 0.541785656682, 0.532949739145, 0.524227434628,
399 0.515614886373, 0.507108489253, 0.498704867478, 0.490400854812,
400 0.482193476986, 0.47407993601, 0.466057596125, 0.458123971214,
401 0.450276713467, 0.442513603171, 0.434832539473, 0.427231532022,
402 0.419708693379, 0.41226223212, 0.404890446548, 0.397591718955,
403 0.390364510382, 0.383207355816, 0.376118859788, 0.369097692334,
404 0.362142585282, 0.355252328834, 0.348425768415, 0.341661801776,
405 0.334959376311, 0.328317486588, 0.321735172063, 0.31521151497,
406 0.308745638367, 0.302336704338, 0.29598391232, 0.289686497571,
407 0.283443729739, 0.27725491156, 0.271119377649, 0.265036493387,
408 0.259005653912, 0.253026283183, 0.247097833139, 0.241219782932,
409 0.235391638239, 0.229612930649, 0.223883217122, 0.218202079518,
410 0.212569124201, 0.206983981709, 0.201446306496, 0.195955776745,
411 0.190512094256, 0.185114984406, 0.179764196185, 0.174459502324,
412 0.169200699492, 0.1639876086, 0.158820075195, 0.153697969964,
413 0.148621189348, 0.143589656295, 0.138603321143, 0.133662162669,
414 0.128766189309, 0.123915440582, 0.119109988745, 0.114349940703,
415 0.10963544023, 0.104966670533, 0.100343857232, 0.0957672718266,
416 0.0912372357329, 0.0867541250127, 0.082318375932, 0.0779304915295,
417 0.0735910494266, 0.0693007111742, 0.065060233529, 0.0608704821745,
418 0.056732448584, 0.05264727098, 0.0486162607163, 0.0446409359769,
419 0.0407230655415, 0.0368647267386, 0.0330683839378, 0.0293369977411,
420 0.0256741818288, 0.0220844372634, 0.0185735200577, 0.0151490552854,
421 0.0118216532614, 0.00860719483079, 0.00553245272614, 0.00265435214565
428 const unsigned int Normal_RNG::ktab[128] = {
429 0, 12590644, 14272653, 14988939,
430 15384584, 15635009, 15807561, 15933577,
431 16029594, 16105155, 16166147, 16216399,
432 16258508, 16294295, 16325078, 16351831,
433 16375291, 16396026, 16414479, 16431002,
434 16445880, 16459343, 16471578, 16482744,
435 16492970, 16502368, 16511031, 16519039,
436 16526459, 16533352, 16539769, 16545755,
437 16551348, 16556584, 16561493, 16566101,
438 16570433, 16574511, 16578353, 16581977,
439 16585398, 16588629, 16591685, 16594575,
440 16597311, 16599901, 16602354, 16604679,
441 16606881, 16608968, 16610945, 16612818,
442 16614592, 16616272, 16617861, 16619363,
443 16620782, 16622121, 16623383, 16624570,
444 16625685, 16626730, 16627708, 16628619,
445 16629465, 16630248, 16630969, 16631628,
446 16632228, 16632768, 16633248, 16633671,
447 16634034, 16634340, 16634586, 16634774,
448 16634903, 16634972, 16634980, 16634926,
449 16634810, 16634628, 16634381, 16634066,
450 16633680, 16633222, 16632688, 16632075,
451 16631380, 16630598, 16629726, 16628757,
452 16627686, 16626507, 16625212, 16623794,
453 16622243, 16620548, 16618698, 16616679,
454 16614476, 16612071, 16609444, 16606571,
455 16603425, 16599973, 16596178, 16591995,
456 16587369, 16582237, 16576520, 16570120,
457 16562917, 16554758, 16545450, 16534739,
458 16522287, 16507638, 16490152, 16468907,
459 16442518, 16408804, 16364095, 16301683,
460 16207738, 16047994, 15704248, 15472926
464 const double Normal_RNG::wtab[128] = {
465 1.62318314817e-08, 2.16291505214e-08, 2.54246305087e-08, 2.84579525938e-08,
466 3.10340022482e-08, 3.33011726243e-08, 3.53439060345e-08, 3.72152672658e-08,
467 3.8950989572e-08, 4.05763964764e-08, 4.21101548915e-08, 4.35664624904e-08,
468 4.49563968336e-08, 4.62887864029e-08, 4.75707945735e-08, 4.88083237257e-08,
469 5.00063025384e-08, 5.11688950428e-08, 5.22996558616e-08, 5.34016475624e-08,
470 5.44775307871e-08, 5.55296344581e-08, 5.65600111659e-08, 5.75704813695e-08,
471 5.85626690412e-08, 5.95380306862e-08, 6.04978791776e-08, 6.14434034901e-08,
472 6.23756851626e-08, 6.32957121259e-08, 6.42043903937e-08, 6.51025540077e-08,
473 6.59909735447e-08, 6.68703634341e-08, 6.77413882848e-08, 6.8604668381e-08,
474 6.94607844804e-08, 7.03102820203e-08, 7.11536748229e-08, 7.1991448372e-08,
475 7.2824062723e-08, 7.36519550992e-08, 7.44755422158e-08, 7.52952223703e-08,
476 7.61113773308e-08, 7.69243740467e-08, 7.77345662086e-08, 7.85422956743e-08,
477 7.93478937793e-08, 8.01516825471e-08, 8.09539758128e-08, 8.17550802699e-08,
478 8.25552964535e-08, 8.33549196661e-08, 8.41542408569e-08, 8.49535474601e-08,
479 8.57531242006e-08, 8.65532538723e-08, 8.73542180955e-08, 8.8156298059e-08,
480 8.89597752521e-08, 8.97649321908e-08, 9.05720531451e-08, 9.138142487e-08,
481 9.21933373471e-08, 9.30080845407e-08, 9.38259651738e-08, 9.46472835298e-08,
482 9.54723502847e-08, 9.63014833769e-08, 9.71350089201e-08, 9.79732621669e-08,
483 9.88165885297e-08, 9.96653446693e-08, 1.00519899658e-07, 1.0138063623e-07,
484 1.02247952126e-07, 1.03122261554e-07, 1.04003996769e-07, 1.04893609795e-07,
485 1.05791574313e-07, 1.06698387725e-07, 1.07614573423e-07, 1.08540683296e-07,
486 1.09477300508e-07, 1.1042504257e-07, 1.11384564771e-07, 1.12356564007e-07,
487 1.13341783071e-07, 1.14341015475e-07, 1.15355110887e-07, 1.16384981291e-07,
488 1.17431607977e-07, 1.18496049514e-07, 1.19579450872e-07, 1.20683053909e-07,
489 1.21808209468e-07, 1.2295639141e-07, 1.24129212952e-07, 1.25328445797e-07,
490 1.26556042658e-07, 1.27814163916e-07, 1.29105209375e-07, 1.30431856341e-07,
491 1.31797105598e-07, 1.3320433736e-07, 1.34657379914e-07, 1.36160594606e-07,
492 1.37718982103e-07, 1.39338316679e-07, 1.41025317971e-07, 1.42787873535e-07,
493 1.44635331499e-07, 1.4657889173e-07, 1.48632138436e-07, 1.50811780719e-07,
494 1.53138707402e-07, 1.55639532047e-07, 1.58348931426e-07, 1.61313325908e-07,
495 1.64596952856e-07, 1.68292495203e-07, 1.72541128694e-07, 1.77574279496e-07,
496 1.83813550477e-07, 1.92166040885e-07, 2.05295471952e-07, 2.22600839893e-07
500 const double Normal_RNG::PARAM_R = 3.44428647676;
505 uint32_t u,
sign, i, j;
530 return sign ? x : -x;
540 setup(meanval, variance);
562 for (
int i=0; i < n; i++)
573 for (i = 0; i < h; i++)
574 for (j = 0; j < w; j++)
586 setup(meanval, variance, rho);
594 factr = -2.0 * var * (1.0 - rho * rho);
611 for (
int i=0; i < n; i++)
622 for (i = 0; i < h; i++)
623 for (j = 0; j < w; j++)
647 mean = tgamma(1.0 + 1.0 / b) / l;
648 var = tgamma(1.0 + 2.0 / b) / (l * l) - mean;
656 for (
int i=0; i < n; i++)
667 for (i = 0; i < h; i++)
668 for (j = 0; j < w; j++)
687 for (
int i=0; i < n; i++)
698 for (i = 0; i < h; i++)
699 for (j = 0; j < w; j++)
718 for (
int i=0; i < n; i++)
729 for (i = 0; i < h; i++)
730 for (j = 0; j < w; j++)