A Discrete-Event Network Simulator
API
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 
36 namespace 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
42 NS_LOG_COMPONENT_DEFINE("RngStream");
43 
44 } // namespace ns3
45 
51 namespace MRG32k3a
52 {
53 
54 // clang-format off
55 
57 typedef double Matrix[3][3];
58 
60 const double m1 = 4294967087.0;
61 
63 const double m2 = 4294944443.0;
64 
66 const double norm = 1.0 / (m1 + 1.0);
67 
69 const double a12 = 1403580.0;
70 
72 const double a13n = 810728.0;
73 
75 const double a21 = 527612.0;
76 
78 const double a23n = 1370589.0;
79 
81 const double two17 = 131072.0;
82 
84 const double two53 = 9007199254740992.0;
85 
87 const 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 
94 const 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 //-------------------------------------------------------------------------
114 double 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 //-------------------------------------------------------------------------
154 void 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 //-------------------------------------------------------------------------
183 void 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 //-------------------------------------------------------------------------
222 void 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 //-------------------------------------------------------------------------
252 void 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 }
318 void 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 
335 namespace ns3
336 {
337 
338 using namespace MRG32k3a;
339 
340 double
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  }
356  m_currentState[0] = m_currentState[1];
357  m_currentState[1] = m_currentState[2];
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  }
368  m_currentState[3] = m_currentState[4];
369  m_currentState[4] = m_currentState[5];
370  m_currentState[5] = p2;
371 
372  /* Combination */
373  u = ((p1 > p2) ? (p1 - p2) * MRG32k3a::norm : (p1 - p2 + m1) * MRG32k3a::norm);
374 
375  return u;
376 }
377 
378 RngStream::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  {
396  m_currentState[i] = r.m_currentState[i];
397  }
398 }
399 
400 void
401 RngStream::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