32 const double m1 = 4294967087.0;
33 const double m2 = 4294944443.0;
34 const double norm = 1.0 / (m1 + 1.0);
35 const double a12 = 1403580.0;
36 const double a13n = 810728.0;
37 const double a21 = 527612.0;
38 const double a23n = 1370589.0;
39 const double two17 = 131072.0;
40 const double two53 = 9007199254740992.0;
41 const double fact = 5.9604644775390625e-8;
46 const double InvA1[3][3] = {
47 { 184888585.0, 0.0, 1945170933.0 },
52 const double InvA2[3][3] = {
53 { 0.0, 360363334.0, 4225571728.0 },
58 const double A1p0[3][3] = {
61 { -810728.0, 1403580.0, 0.0 }
64 const double A2p0[3][3] = {
67 { -1370589.0, 0.0, 527612.0 }
70 const double A1p76[3][3] = {
71 { 82758667.0, 1871391091.0, 4127413238.0 },
72 { 3672831523.0, 69195019.0, 1871391091.0 },
73 { 3672091415.0, 3528743235.0, 69195019.0 }
76 const double A2p76[3][3] = {
77 { 1511326704.0, 3759209742.0, 1610795712.0 },
78 { 4292754251.0, 1511326704.0, 3889917532.0 },
79 { 3859662829.0, 4292754251.0, 3708466080.0 }
82 const double A1p127[3][3] = {
83 { 2427906178.0, 3580155704.0, 949770784.0 },
84 { 226153695.0, 1230515664.0, 3580155704.0 },
85 { 1988835001.0, 986791581.0, 1230515664.0 }
88 const double A2p127[3][3] = {
89 { 1464411153.0, 277697599.0, 1610723613.0 },
90 { 32183930.0, 1464411153.0, 1022607788.0 },
91 { 2824425944.0, 32183930.0, 2093834863.0 }
94 const double InvA1p76[3][3] = {
95 { 2585822061.0, 2346541846.0, 600781890.0},
96 { 42385315.0, 4257896290.0, 2346541846.0},
97 { 1248824805.0, 2390631828.0, 4257896290.0}
100 const double InvA2p76[3][3] = {
101 { 855407695.0, 4134906251.0, 112088500.0},
102 { 2897599610.0, 855407695.0, 1987588141.0},
103 { 854109890.0, 2897599610.0, 1099731892.0}
106 const double InvA1p127[3][3] = {
107 { 1737602145.0, 4223429791.0, 3623540388.0},
108 { 296220492.0, 2329649265.0, 4223429791.0},
109 { 2348335727.0, 4013827785.0, 2329649265.0}
112 const double InvA2p127[3][3] = {
113 { 3244453311.0, 924928292.0, 95182267.0},
114 { 3169310862.0, 3244453311.0, 3740757140.0},
115 { 2096686917.0, 3169310862.0, 2895877872.0}
122 double MultModM (
double a,
double s,
double c,
double m)
129 if (v >= two53 || v <= -two53) {
130 a1 =
static_cast<int32_t
> (a / two17); a -= a1 * two17;
132 a1 =
static_cast<int32_t
> (v / m); v -= a1 * m;
133 v = v * two17 + a * s + c;
136 a1 =
static_cast<int32_t
> (v / m);
138 if ((v -= a1 * m) < 0.0)
return v += m;
else return v;
146 void MatVecModM (
const double A[3][3],
const double s[3],
double v[3],
152 for (i = 0; i < 3; ++i) {
153 x[i] = MultModM (A[i][0], s[0], 0.0, m);
154 x[i] = MultModM (A[i][1], s[1], x[i], m);
155 x[i] = MultModM (A[i][2], s[2], x[i], m);
157 for (i = 0; i < 3; ++i)
166 void MatMatModM (
const double A[3][3],
const double B[3][3],
167 double C[3][3],
double m)
170 double V[3], W[3][3];
172 for (i = 0; i < 3; ++i) {
173 for (j = 0; j < 3; ++j)
175 MatVecModM (A, V, V, m);
176 for (j = 0; j < 3; ++j)
179 for (i = 0; i < 3; ++i)
180 for (j = 0; j < 3; ++j)
188 void MatTwoPowModM (
const double A[3][3],
double B[3][3],
double m, int32_t e)
194 for (i = 0; i < 3; ++i)
195 for (j = 0; j < 3; ++j)
199 for (i = 0; i < e; i++)
200 MatMatModM (B, B, B, m);
207 void MatPowModM (
const double A[3][3],
double B[3][3],
double m, int32_t n)
213 for (i = 0; i < 3; ++i)
214 for (j = 0; j < 3; ++j) {
218 for (j = 0; j < 3; ++j)
223 if (n % 2) MatMatModM (W, B, B, m);
224 MatMatModM (W, W, W, m);
234 int CheckSeed (
const double seed[6])
238 for (i = 0; i < 3; ++i) {
246 for (i = 3; i < 6; ++i) {
254 if (seed[0] == 0 && seed[1] == 0 && seed[2] == 0) {
260 if (seed[3] == 0 && seed[4] == 0 && seed[5] == 0) {
282 p1 = a12 *
Cg[1] - a13n *
Cg[0];
283 k =
static_cast<int32_t
> (p1 / m1);
285 if (p1 < 0.0) p1 += m1;
289 p2 = a21 *
Cg[5] - a23n *
Cg[3];
290 k =
static_cast<int32_t
> (p2 / m2);
292 if (p2 < 0.0) p2 += m2;
296 u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm);
298 return (
anti ==
false) ? u : (1 - u);
311 u += (
U01() - 1.0) * fact;
312 return (u < 0.0) ? u + 1.0 : u;
315 return (u < 1.0) ? u : (u - 1.0);
330 12345.0, 12345.0, 12345.0, 12345.0, 12345.0, 12345.0
349 for (
int i = 0; i < 6; ++i) {
363 for (
int i = 0; i < 6; ++i)
373 for (
int i = 0; i < 6; ++i)
383 MatVecModM(A1p76,
Bg,
Bg, m1);
384 MatVecModM(A2p76, &
Bg[3], &
Bg[3], m2);
385 for (
int i = 0; i < 6; ++i)
393 if (CheckSeed (seed))
395 for (
int i = 0; i < 6; ++i)
404 if (CheckSeed (seed))
406 for (
int i = 0; i < 6; ++i)
407 Cg[i] =
Bg[i] =
Ig[i] = seed[i];
419 const double A1[3][3],
const double A2[3][3],
420 const double InvA1[3][3],
const double InvA2[3][3])
422 double B1[3][3], C1[3][3], B2[3][3], C2[3][3];
425 MatTwoPowModM (A1, B1, m1, e);
426 MatTwoPowModM (A2, B2, m2, e);
428 MatTwoPowModM (InvA1, B1, m1, -e);
429 MatTwoPowModM (InvA2, B2, m2, -e);
433 MatPowModM (A1, C1, m1, c);
434 MatPowModM (A2, C2, m2, c);
436 MatPowModM (InvA1, C1, m1, -c);
437 MatPowModM (InvA2, C2, m2, -c);
441 MatMatModM (B1, C1, C1, m1);
442 MatMatModM (B2, C2, C2, m2);
445 MatVecModM (C1,
Cg,
Cg, m1);
446 MatVecModM (C2, &
Cg[3], &
Cg[3], m2);
456 for (
int i = 0; i < 6; ++i)
462 for (
int i = 0; i < 6; ++i)
469 double B1[3][3], B2[3][3];
472 MatTwoPowModM (A1p0, B1, m1, e);
473 MatTwoPowModM (A2p0, B2, m2, e);
476 MatTwoPowModM (InvA1, B1, m1, -e);
477 MatTwoPowModM (InvA2, B2, m2, -e);
481 MatPowModM (A1p0, C1, m1, c);
482 MatPowModM (A2p0, C2, m2, c);
485 MatPowModM (InvA1, C1, m1, -c);
486 MatPowModM (InvA2, C2, m2, -c);
490 MatMatModM (B1, C1, C1, m1);
491 MatMatModM (B2, C2, C2, m2);
498 for (
int i = 0; i < 6; ++i)
580 return low +
static_cast<int> ((high - low + 1.0) *
RandU01 ());