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