A Discrete-Event Network Simulator
API
Loading...
Searching...
No Matches
rng-stream.cc
Go to the documentation of this file.
1//
2// Copyright (C) 2001 Pierre L'Ecuyer (lecuyer@iro.umontreal.ca)
3//
4// This program is free software; you can redistribute it and/or modify
5// it under the terms of the GNU General Public License version 2 as
6// published by the Free Software Foundation;
7//
8// This program is distributed in the hope that it will be useful,
9// but WITHOUT ANY WARRANTY; without even the implied warranty of
10// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11// GNU General Public License for more details.
12//
13// You should have received a copy of the GNU General Public License
14// along with this program; if not, write to the Free Software
15// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16//
17// Modified for ns-3 by:
18// - Rajib Bhattacharjea<raj.b@gatech.edu>
19// - Mathieu Lacage <mathieu.lacage@gmail.com>
20//
21
22#include "rng-stream.h"
23
24#include "fatal-error.h"
25#include "log.h"
26
27#include <cstdlib>
28#include <iostream>
29
36namespace ns3
37{
38
39// Note: Logging in this file is largely avoided due to the
40// number of calls that are made to these functions and the possibility
41// of causing recursions leading to stack overflow
42NS_LOG_COMPONENT_DEFINE("RngStream");
43
44} // namespace ns3
45
51namespace MRG32k3a
52{
53
54// clang-format off
55
57typedef double Matrix[3][3];
58
60const double m1 = 4294967087.0;
61
63const double m2 = 4294944443.0;
64
66const double norm = 1.0 / (m1 + 1.0);
67
69const double a12 = 1403580.0;
70
72const double a13n = 810728.0;
73
75const double a21 = 527612.0;
76
78const double a23n = 1370589.0;
79
81const double two17 = 131072.0;
82
84const double two53 = 9007199254740992.0;
85
87const Matrix A1p0 = {
88 { 0.0, 1.0, 0.0 },
89 { 0.0, 0.0, 1.0 },
90 { -810728.0, 1403580.0, 0.0 }
91};
92
94const Matrix A2p0 = {
95 { 0.0, 1.0, 0.0 },
96 { 0.0, 0.0, 1.0 },
97 { -1370589.0, 0.0, 527612.0 }
98};
99
100
101//-------------------------------------------------------------------------
114double MultModM (double a, double s, double c, double m)
115{
116 double v;
117 int32_t a1;
118
119 v = a * s + c;
120
121 if (v >= two53 || v <= -two53)
122 {
123 a1 = static_cast<int32_t> (a / two17);
124 a -= a1 * two17;
125 v = a1 * s;
126 a1 = static_cast<int32_t> (v / m);
127 v -= a1 * m;
128 v = v * two17 + a * s + c;
129 }
130
131 a1 = static_cast<int32_t> (v / m);
132 /* in case v < 0)*/
133 if ((v -= a1 * m) < 0.0)
134 {
135 return v += m;
136 }
137 else
138 {
139 return v;
140 }
141}
142
143
144//-------------------------------------------------------------------------
154void MatVecModM (const Matrix A, const double s[3], double v[3],
155 double m)
156{
157 int i;
158 double x[3]; // Necessary if v = s
159
160 for (i = 0; i < 3; ++i)
161 {
162 x[i] = MultModM (A[i][0], s[0], 0.0, m);
163 x[i] = MultModM (A[i][1], s[1], x[i], m);
164 x[i] = MultModM (A[i][2], s[2], x[i], m);
165 }
166 for (i = 0; i < 3; ++i)
167 {
168 v[i] = x[i];
169 }
170}
171
172
173//-------------------------------------------------------------------------
183void MatMatModM (const Matrix A, const Matrix B,
184 Matrix C, double m)
185{
186 int i;
187 int j;
188 double V[3];
189 Matrix W;
190
191 for (i = 0; i < 3; ++i)
192 {
193 for (j = 0; j < 3; ++j)
194 {
195 V[j] = B[j][i];
196 }
197 MatVecModM (A, V, V, m);
198 for (j = 0; j < 3; ++j)
199 {
200 W[j][i] = V[j];
201 }
202 }
203 for (i = 0; i < 3; ++i)
204 {
205 for (j = 0; j < 3; ++j)
206 {
207 C[i][j] = W[i][j];
208 }
209 }
210}
211
212
213//-------------------------------------------------------------------------
222void MatTwoPowModM (const Matrix src, Matrix dst, double m, int32_t e)
223{
224 int i;
225 int j;
226
227 /* initialize: dst = src */
228 for (i = 0; i < 3; ++i)
229 {
230 for (j = 0; j < 3; ++j)
231 {
232 dst[i][j] = src[i][j];
233 }
234 }
235 /* Compute dst = src^(2^e) mod m */
236 for (i = 0; i < e; i++)
237 {
238 MatMatModM (dst, dst, dst, m);
239 }
240}
241
242
243//-------------------------------------------------------------------------
252void MatPowModM (const double A[3][3], double B[3][3], double m, int32_t n)
253{
254 int i;
255 int j;
256 double W[3][3];
257
258 // initialize: W = A; B = I
259 for (i = 0; i < 3; ++i)
260 {
261 for (j = 0; j < 3; ++j)
262 {
263 W[i][j] = A[i][j];
264 B[i][j] = 0.0;
265 }
266 }
267 for (j = 0; j < 3; ++j)
268 {
269 B[j][j] = 1.0;
270 }
271
272 // Compute B = A^n mod m using the binary decomposition of n
273 while (n > 0)
274 {
275 if (n % 2)
276 {
277 MatMatModM (W, B, B, m);
278 }
279 MatMatModM (W, W, W, m);
280 n /= 2;
281 }
282}
283
289{
290 Matrix a1[190];
291 Matrix a2[190];
292};
293
301{
302 Precalculated precalculated;
303 for (int i = 0; i < 190; i++)
304 {
305 int power = i + 1;
306 MatTwoPowModM (A1p0, precalculated.a1[i], m1, power);
307 MatTwoPowModM (A2p0, precalculated.a2[i], m2, power);
308 }
309 return precalculated;
310}
318void PowerOfTwoMatrix (int n, Matrix a1p, Matrix a2p)
319{
320 static Precalculated constants = PowerOfTwoConstants ();
321 for (int i = 0; i < 3; i ++)
322 {
323 for (int j = 0; j < 3; j++)
324 {
325 a1p[i][j] = constants.a1[n-1][i][j];
326 a2p[i][j] = constants.a2[n-1][i][j];
327 }
328 }
329}
330
331} // namespace MRG32k3a
332
333// clang-format on
334
335namespace ns3
336{
337
338using namespace MRG32k3a;
339
340double
342{
343 int32_t k;
344 double p1;
345 double p2;
346 double u;
347
348 /* Component 1 */
349 p1 = a12 * m_currentState[1] - a13n * m_currentState[0];
350 k = static_cast<int32_t>(p1 / m1);
351 p1 -= k * m1;
352 if (p1 < 0.0)
353 {
354 p1 += m1;
355 }
358 m_currentState[2] = p1;
359
360 /* Component 2 */
361 p2 = a21 * m_currentState[5] - a23n * m_currentState[3];
362 k = static_cast<int32_t>(p2 / m2);
363 p2 -= k * m2;
364 if (p2 < 0.0)
365 {
366 p2 += m2;
367 }
370 m_currentState[5] = p2;
371
372 /* Combination */
373 u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm);
374
375 return u;
376}
377
378RngStream::RngStream(uint32_t seedNumber, uint64_t stream, uint64_t substream)
379{
380 if (seedNumber >= m1 || seedNumber >= m2 || seedNumber == 0)
381 {
382 NS_FATAL_ERROR("invalid Seed " << seedNumber);
383 }
384 for (int i = 0; i < 6; ++i)
385 {
386 m_currentState[i] = seedNumber;
387 }
388 AdvanceNthBy(stream, 127, m_currentState);
389 AdvanceNthBy(substream, 76, m_currentState);
390}
391
393{
394 for (int i = 0; i < 6; ++i)
395 {
397 }
398}
399
400void
401RngStream::AdvanceNthBy(uint64_t nth, int by, double state[6])
402{
403 Matrix matrix1;
404 Matrix matrix2;
405 for (int i = 0; i < 64; i++)
406 {
407 int nbit = 63 - i;
408 int bit = (nth >> nbit) & 0x1;
409 if (bit)
410 {
411 PowerOfTwoMatrix(by + nbit, matrix1, matrix2);
412 MatVecModM(matrix1, state, state, m1);
413 MatVecModM(matrix2, &state[3], &state[3], m2);
414 }
415 }
416}
417
418} // namespace ns3
419 // \ingroup rngimpl
Combined Multiple-Recursive Generator MRG32k3a.
Definition: rng-stream.h:50
RngStream(uint32_t seed, uint64_t stream, uint64_t substream)
Construct from explicit seed, stream and substream values.
Definition: rng-stream.cc:378
double m_currentState[6]
The RNG state vector.
Definition: rng-stream.h:85
void AdvanceNthBy(uint64_t nth, int by, double state[6])
Advance state of the RNG by leaps and bounds.
Definition: rng-stream.cc:401
double RandU01()
Generate the next random number for this stream.
Definition: rng-stream.cc:341
NS_FATAL_x macro definitions.
#define NS_FATAL_ERROR(msg)
Report a fatal error with a message and terminate.
Definition: fatal-error.h:179
#define NS_LOG_COMPONENT_DEFINE(name)
Define a Log component with a specific name.
Definition: log.h:202
Debug message logging.
Namespace for MRG32k3a implementation details.
Definition: rng-stream.cc:52
const double a21
Second component multiplier of n - 1 value.
Definition: rng-stream.cc:75
const double m1
First component modulus, 232 - 209.
Definition: rng-stream.cc:60
void MatMatModM(const Matrix A, const Matrix B, Matrix C, double m)
Compute the matrix C = A*B MOD m.
Definition: rng-stream.cc:183
const Matrix A1p0
First component transition matrix.
Definition: rng-stream.cc:87
const double a23n
Second component multiplier of n - 3 value.
Definition: rng-stream.cc:78
const double two17
Decomposition factor for computing a*s in less than 53 bits, 217
Definition: rng-stream.cc:81
void PowerOfTwoMatrix(int n, Matrix a1p, Matrix a2p)
Get the transition matrices raised to a power of 2.
Definition: rng-stream.cc:318
const double norm
Normalization to obtain randoms on [0,1).
Definition: rng-stream.cc:66
const double a12
First component multiplier of n - 2 value.
Definition: rng-stream.cc:69
void MatPowModM(const double A[3][3], double B[3][3], double m, int32_t n)
Compute the matrix B = (A^n Mod m); works even if A = B.
Definition: rng-stream.cc:252
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.
Definition: rng-stream.cc:222
const Matrix A2p0
Second component transition matrix.
Definition: rng-stream.cc:94
double Matrix[3][3]
Type for 3x3 matrix of doubles.
Definition: rng-stream.cc:57
const double m2
Second component modulus, 232 - 22853.
Definition: rng-stream.cc:63
const double a13n
First component multiplier of n - 3 value.
Definition: rng-stream.cc:72
const double two53
IEEE-754 floating point precision, 253
Definition: rng-stream.cc:84
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.
Definition: rng-stream.cc:114
void MatVecModM(const Matrix A, const double s[3], double v[3], double m)
Compute the vector v = A*s MOD m.
Definition: rng-stream.cc:154
Precalculated PowerOfTwoConstants()
Compute the transition matrices of the two MRG components raised to all powers of 2 from 1 to 191.
Definition: rng-stream.cc:300
Every class exported by the ns3 library is enclosed in the ns3 namespace.
ns3::RngStream declaration.
The transition matrices of the two MRG components (in matrix form), raised to all powers of 2 from 1 ...
Definition: rng-stream.cc:289
Matrix a1[190]
First component transition matrix powers.
Definition: rng-stream.cc:290
Matrix a2[190]
Second component transition matrix powers.
Definition: rng-stream.cc:291