54 const double m1 = 4294967087.0;
58 const double m2 = 4294944443.0;
62 const double norm = 1.0 / (m1 + 1.0);
66 const double a12 = 1403580.0;
70 const double a13n = 810728.0;
74 const double a21 = 527612.0;
78 const double a23n = 1370589.0;
86 const double two53 = 9007199254740992.0;
93 { -810728.0, 1403580.0, 0.0 }
101 { -1370589.0, 0.0, 527612.0 }
118 double MultModM (
double a,
double s,
double c,
double m)
125 if (v >= two53 || v <= -two53)
127 a1 =
static_cast<int32_t
> (a /
two17);
130 a1 =
static_cast<int32_t
> (v / m);
132 v = v * two17 + a * s + c;
135 a1 =
static_cast<int32_t
> (v / m);
137 if ((v -= a1 * m) < 0.0)
158 void MatVecModM (
const Matrix A,
const double s[3],
double v[3],
164 for (i = 0; i < 3; ++i)
166 x[i] =
MultModM (A[i][0], s[0], 0.0, m);
167 x[i] =
MultModM (A[i][1], s[1], x[i], m);
168 x[i] =
MultModM (A[i][2], s[2], x[i], m);
170 for (i = 0; i < 3; ++i)
194 for (i = 0; i < 3; ++i)
196 for (j = 0; j < 3; ++j)
201 for (j = 0; j < 3; ++j)
206 for (i = 0; i < 3; ++i)
208 for (j = 0; j < 3; ++j)
230 for (i = 0; i < 3; ++i)
232 for (j = 0; j < 3; ++j)
234 dst[i][j] = src[i][j];
238 for (i = 0; i < e; i++)
304 for (
int i = 0; i < 190; i++)
310 return precalculated;
321 for (
int i = 0; i < 3; i ++)
323 for (
int j = 0; j < 3; j++)
325 a1p[i][j] = constants.
a1[n-1][i][j];
326 a2p[i][j] = constants.
a2[n-1][i][j];
345 k =
static_cast<int32_t
> (p1 /
m1);
351 m_currentState[0] = m_currentState[1]; m_currentState[1] = m_currentState[2]; m_currentState[2] = p1;
354 p2 =
a21 * m_currentState[5] -
a23n * m_currentState[3];
355 k =
static_cast<int32_t
> (p2 /
m2);
361 m_currentState[3] = m_currentState[4]; m_currentState[4] = m_currentState[5]; m_currentState[5] = p2;
364 u = ((p1 > p2) ? (p1 - p2) *
norm : (p1 - p2 +
m1) *
norm);
371 if (seedNumber >=
m1 || seedNumber >=
m2 || seedNumber == 0)
375 for (
int i = 0; i < 6; ++i)
385 for (
int i = 0; i < 6; ++i)
395 for (
int i = 0; i < 64; i++)
398 int bit = (nth >> nbit) & 0x1;
NS_FATAL_x macro definitions.
const Matrix A1p0
First component transition matrix.
const double a13n
First component multiplier of n - 3 value.
double MultModM(double a, double s, double c, double m)
Return (a*s + c) MOD m; a, s, c and m must be < 2^35.
const double a12
First component multiplier of n - 2 value.
void MatMatModM(const Matrix A, const Matrix B, Matrix C, double m)
Compute the matrix C = A*B MOD m.
Matrix a2[190]
Second component transition matrix powers.
The transition matrices of the two MRG components (in matrix form), raised to all powers of 2 from 1 ...
#define NS_LOG_COMPONENT_DEFINE(name)
Define a Log component with a specific name.
#define NS_FATAL_ERROR(msg)
Report a fatal error with a message and terminate.
void MatVecModM(const Matrix A, const double s[3], double v[3], double m)
Compute the vector v = A*s MOD m.
Combined Multiple-Recursive Generator MRG32k3a.
Matrix a1[190]
First component transition matrix powers.
void MatTwoPowModM(const Matrix src, Matrix dst, double m, int32_t e)
Compute the matrix B = (A^(2^e) Mod m); works also if A = B.
void AdvanceNthBy(uint64_t nth, int by, double state[6])
Advance state of the RNG by leaps and bounds.
double m_currentState[6]
The RNG state vector.
const double a23n
Second component multiplier of n - 3 value.
RngStream(uint32_t seed, uint64_t stream, uint64_t substream)
Construct from explicit seed, stream and substream values.
const double norm
Normalization to obtain randoms on [0,1).
const double a21
Second component multiplier of n - 1 value.
Declaration of class RngStream.
struct Precalculated PowerOfTwoConstants(void)
Compute the transition matrices of the two MRG components.
double RandU01(void)
Generate the next random number for this stream.
void PowerOfTwoMatrix(int n, Matrix a1p, Matrix a2p)
Get the transition matrices raised to a power of 2.
const double m2
Second component modulus, 232 - 22853.
const Matrix A2p0
Second component transition matrix.
double Matrix[3][3]
Type for 3x3 matrix of doubles.
const double m1
First component modulus, 232 - 209.
const double two17
Decomposition factor for computing a*s in less than 53 bits, 217
const double two53
IEEE-754 floating point precision, 253