25 #include "rng-stream.h"
26 #include "fatal-error.h"
37 typedef double Matrix[3][3];
39 const double m1 = 4294967087.0;
40 const double m2 = 4294944443.0;
41 const double norm = 1.0 / (m1 + 1.0);
42 const double a12 = 1403580.0;
43 const double a13n = 810728.0;
44 const double a21 = 527612.0;
45 const double a23n = 1370589.0;
46 const double two17 = 131072.0;
47 const double two53 = 9007199254740992.0;
48 const double fact = 5.9604644775390625e-8;
50 const Matrix InvA1 = {
51 { 184888585.0, 0.0, 1945170933.0 },
56 const Matrix InvA2 = {
57 { 0.0, 360363334.0, 4225571728.0 },
65 { -810728.0, 1403580.0, 0.0 }
71 { -1370589.0, 0.0, 527612.0 }
78 double MultModM (
double a,
double s,
double c,
double m)
85 if (v >= two53 || v <= -two53)
87 a1 =
static_cast<int32_t
> (a / two17);
90 a1 =
static_cast<int32_t
> (v / m);
92 v = v * two17 + a * s + c;
95 a1 =
static_cast<int32_t
> (v / m);
97 if ((v -= a1 * m) < 0.0)
112 void MatVecModM (
const Matrix
A,
const double s[3],
double v[3],
118 for (i = 0; i < 3; ++i)
120 x[i] = MultModM (A[i][0], s[0], 0.0, m);
121 x[i] = MultModM (A[i][1], s[1], x[i], m);
122 x[i] = MultModM (A[i][2], s[2], x[i], m);
124 for (i = 0; i < 3; ++i)
135 void MatMatModM (
const Matrix A,
const Matrix B,
142 for (i = 0; i < 3; ++i)
144 for (j = 0; j < 3; ++j)
148 MatVecModM (A, V, V, m);
149 for (j = 0; j < 3; ++j)
154 for (i = 0; i < 3; ++i)
156 for (j = 0; j < 3; ++j)
167 void MatTwoPowModM (
const Matrix src, Matrix dst,
double m, int32_t e)
172 for (i = 0; i < 3; ++i)
174 for (j = 0; j < 3; ++j)
176 dst[i][j] = src[i][j];
180 for (i = 0; i < e; i++)
182 MatMatModM (dst, dst, dst, m);
234 for (
int i = 0; i < 190; i++)
237 MatTwoPowModM (A1p0, precalculated.a1[i], m1, power);
238 MatTwoPowModM (A2p0, precalculated.a2[i], m2, power);
240 return precalculated;
242 void PowerOfTwoMatrix (
int n, Matrix a1p, Matrix a2p)
244 static struct Precalculated constants = PowerOfTwoConstants ();
245 for (
int i = 0; i < 3; i ++)
247 for (
int j = 0; j < 3; j++)
249 a1p[i][j] = constants.a1[n-1][i][j];
250 a2p[i][j] = constants.a2[n-1][i][j];
268 p1 = a12 * m_currentState[1] - a13n * m_currentState[0];
269 k =
static_cast<int32_t
> (p1 / m1);
275 m_currentState[0] = m_currentState[1]; m_currentState[1] = m_currentState[2]; m_currentState[2] = p1;
278 p2 = a21 * m_currentState[5] - a23n * m_currentState[3];
279 k =
static_cast<int32_t
> (p2 / m2);
285 m_currentState[3] = m_currentState[4]; m_currentState[4] = m_currentState[5]; m_currentState[5] = p2;
288 u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm);
293 RngStream::RngStream (uint32_t seedNumber, uint64_t stream, uint64_t substream)
295 if (seedNumber >= m1 || seedNumber >= m2 || seedNumber == 0)
299 for (
int i = 0; i < 6; ++i)
301 m_currentState[i] = seedNumber;
303 AdvanceNthBy (stream, 127, m_currentState);
304 AdvanceNthBy (substream, 76, m_currentState);
307 RngStream::RngStream(
const RngStream& r)
309 for (
int i = 0; i < 6; ++i)
311 m_currentState[i] = r.m_currentState[i];
316 RngStream::AdvanceNthBy (uint64_t nth,
int by,
double state[6])
318 Matrix matrix1,matrix2;
319 for (
int i = 0; i < 64; i++)
322 int bit = (nth >> nbit) & 0x1;
325 PowerOfTwoMatrix(by + nbit, matrix1, matrix2);
326 MatVecModM (matrix1, state, state, m1);
327 MatVecModM (matrix2, &state[3], &state[3], m2);
#define NS_LOG_COMPONENT_DEFINE(name)
#define NS_FATAL_ERROR(msg)
fatal error handling