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 
35 namespace 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
40 NS_LOG_COMPONENT_DEFINE ("RngStream");
41 
42 } // namespace ns3
43 
44 
50 namespace MRG32k3a
51 {
52 
54 typedef double Matrix[3][3];
55 
57 const double m1 = 4294967087.0;
58 
60 const double m2 = 4294944443.0;
61 
63 const double norm = 1.0 / (m1 + 1.0);
64 
66 const double a12 = 1403580.0;
67 
69 const double a13n = 810728.0;
70 
72 const double a21 = 527612.0;
73 
75 const double a23n = 1370589.0;
76 
78 const double two17 = 131072.0;
79 
81 const double two53 = 9007199254740992.0;
82 
84 const Matrix A1p0 = {
85  { 0.0, 1.0, 0.0 },
86  { 0.0, 0.0, 1.0 },
87  { -810728.0, 1403580.0, 0.0 }
88 };
89 
91 const Matrix A2p0 = {
92  { 0.0, 1.0, 0.0 },
93  { 0.0, 0.0, 1.0 },
94  { -1370589.0, 0.0, 527612.0 }
95 };
96 
97 
98 //-------------------------------------------------------------------------
111 double MultModM (double a, double s, double c, double m)
112 {
113  double v;
114  int32_t a1;
115 
116  v = a * s + c;
117 
118  if (v >= two53 || v <= -two53)
119  {
120  a1 = static_cast<int32_t> (a / two17);
121  a -= a1 * two17;
122  v = a1 * s;
123  a1 = static_cast<int32_t> (v / m);
124  v -= a1 * m;
125  v = v * two17 + a * s + c;
126  }
127 
128  a1 = static_cast<int32_t> (v / m);
129  /* in case v < 0)*/
130  if ((v -= a1 * m) < 0.0)
131  {
132  return v += m;
133  }
134  else
135  {
136  return v;
137  }
138 }
139 
140 
141 //-------------------------------------------------------------------------
151 void MatVecModM (const Matrix A, const double s[3], double v[3],
152  double m)
153 {
154  int i;
155  double x[3]; // Necessary if v = s
156 
157  for (i = 0; i < 3; ++i)
158  {
159  x[i] = MultModM (A[i][0], s[0], 0.0, m);
160  x[i] = MultModM (A[i][1], s[1], x[i], m);
161  x[i] = MultModM (A[i][2], s[2], x[i], m);
162  }
163  for (i = 0; i < 3; ++i)
164  {
165  v[i] = x[i];
166  }
167 }
168 
169 
170 //-------------------------------------------------------------------------
180 void MatMatModM (const Matrix A, const Matrix B,
181  Matrix C, double m)
182 {
183  int i, j;
184  double V[3];
185  Matrix W;
186 
187  for (i = 0; i < 3; ++i)
188  {
189  for (j = 0; j < 3; ++j)
190  {
191  V[j] = B[j][i];
192  }
193  MatVecModM (A, V, V, m);
194  for (j = 0; j < 3; ++j)
195  {
196  W[j][i] = V[j];
197  }
198  }
199  for (i = 0; i < 3; ++i)
200  {
201  for (j = 0; j < 3; ++j)
202  {
203  C[i][j] = W[i][j];
204  }
205  }
206 }
207 
208 
209 //-------------------------------------------------------------------------
218 void MatTwoPowModM (const Matrix src, Matrix dst, double m, int32_t e)
219 {
220  int i, j;
221 
222  /* initialize: dst = src */
223  for (i = 0; i < 3; ++i)
224  {
225  for (j = 0; j < 3; ++j)
226  {
227  dst[i][j] = src[i][j];
228  }
229  }
230  /* Compute dst = src^(2^e) mod m */
231  for (i = 0; i < e; i++)
232  {
233  MatMatModM (dst, dst, dst, m);
234  }
235 }
236 
237 
238 //-------------------------------------------------------------------------
247 void MatPowModM (const double A[3][3], double B[3][3], double m, int32_t n)
248 {
249  int i, j;
250  double W[3][3];
251 
252  // initialize: W = A; B = I
253  for (i = 0; i < 3; ++i)
254  {
255  for (j = 0; j < 3; ++j)
256  {
257  W[i][j] = A[i][j];
258  B[i][j] = 0.0;
259  }
260  }
261  for (j = 0; j < 3; ++j)
262  {
263  B[j][j] = 1.0;
264  }
265 
266  // Compute B = A^n mod m using the binary decomposition of n
267  while (n > 0)
268  {
269  if (n % 2)
270  {
271  MatMatModM (W, B, B, m);
272  }
273  MatMatModM (W, W, W, m);
274  n /= 2;
275  }
276 }
277 
283 {
284  Matrix a1[190];
285  Matrix a2[190];
286 };
287 
295 {
296  struct Precalculated precalculated;
297  for (int i = 0; i < 190; i++)
298  {
299  int power = i + 1;
300  MatTwoPowModM (A1p0, precalculated.a1[i], m1, power);
301  MatTwoPowModM (A2p0, precalculated.a2[i], m2, power);
302  }
303  return precalculated;
304 }
312 void PowerOfTwoMatrix (int n, Matrix a1p, Matrix a2p)
313 {
314  static struct Precalculated constants = PowerOfTwoConstants ();
315  for (int i = 0; i < 3; i ++)
316  {
317  for (int j = 0; j < 3; j++)
318  {
319  a1p[i][j] = constants.a1[n-1][i][j];
320  a2p[i][j] = constants.a2[n-1][i][j];
321  }
322  }
323 }
324 
325 } // namespace MRG32k3a
326 
327 
328 namespace ns3 {
329 
330 using namespace MRG32k3a;
331 
333 {
334  int32_t k;
335  double p1, p2, u;
336 
337  /* Component 1 */
338  p1 = a12 * m_currentState[1] - a13n * m_currentState[0];
339  k = static_cast<int32_t> (p1 / m1);
340  p1 -= k * m1;
341  if (p1 < 0.0)
342  {
343  p1 += m1;
344  }
345  m_currentState[0] = m_currentState[1]; m_currentState[1] = m_currentState[2]; m_currentState[2] = p1;
346 
347  /* Component 2 */
348  p2 = a21 * m_currentState[5] - a23n * m_currentState[3];
349  k = static_cast<int32_t> (p2 / m2);
350  p2 -= k * m2;
351  if (p2 < 0.0)
352  {
353  p2 += m2;
354  }
355  m_currentState[3] = m_currentState[4]; m_currentState[4] = m_currentState[5]; m_currentState[5] = p2;
356 
357  /* Combination */
358  u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm);
359 
360  return u;
361 }
362 
363 RngStream::RngStream (uint32_t seedNumber, uint64_t stream, uint64_t substream)
364 {
365  if (seedNumber >= m1 || seedNumber >= m2 || seedNumber == 0)
366  {
367  NS_FATAL_ERROR ("invalid Seed " << seedNumber);
368  }
369  for (int i = 0; i < 6; ++i)
370  {
371  m_currentState[i] = seedNumber;
372  }
373  AdvanceNthBy (stream, 127, m_currentState);
374  AdvanceNthBy (substream, 76, m_currentState);
375 }
376 
378 {
379  for (int i = 0; i < 6; ++i)
380  {
381  m_currentState[i] = r.m_currentState[i];
382  }
383 }
384 
385 void
386 RngStream::AdvanceNthBy (uint64_t nth, int by, double state[6])
387 {
388  Matrix matrix1,matrix2;
389  for (int i = 0; i < 64; i++)
390  {
391  int nbit = 63 - i;
392  int bit = (nth >> nbit) & 0x1;
393  if (bit)
394  {
395  PowerOfTwoMatrix(by + nbit, matrix1, matrix2);
396  MatVecModM (matrix1, state, state, m1);
397  MatVecModM (matrix2, &state[3], &state[3], m2);
398  }
399  }
400 }
401 
402 } // namespace ns3
403  // \ingroup rngimpl
NS_FATAL_x macro definitions.
double Matrix[3][3]
Type for 3x3 matrix of doubles.
Definition: rng-stream.cc:54
The transition matrices of the two MRG components (in matrix form), raised to all powers of 2 from 1 ...
Definition: rng-stream.cc:282
const double a23n
Second component multiplier of n - 3 value.
Definition: rng-stream.cc:75
const double a13n
First component multiplier of n - 3 value.
Definition: rng-stream.cc:69
const Matrix A2p0
Second component transition matrix.
Definition: rng-stream.cc:91
const double two53
IEEE-754 floating point precision, 253
Definition: rng-stream.cc:81
#define NS_LOG_COMPONENT_DEFINE(name)
Define a Log component with a specific name.
Definition: log.h:204
const double two17
Decomposition factor for computing a*s in less than 53 bits, 217
Definition: rng-stream.cc:78
#define NS_FATAL_ERROR(msg)
Report a fatal error with a message and terminate.
Definition: fatal-error.h:162
const double m1
First component modulus, 232 - 209.
Definition: rng-stream.cc:57
Namespace for MRG32k3a implementation details.
Definition: rng-stream.cc:50
void MatMatModM(const Matrix A, const Matrix B, Matrix C, double m)
Compute the matrix C = A*B MOD m.
Definition: rng-stream.cc:180
const Matrix A1p0
First component transition matrix.
Definition: rng-stream.cc:84
Combined Multiple-Recursive Generator MRG32k3a.
Definition: rng-stream.h:49
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:111
void AdvanceNthBy(uint64_t nth, int by, double state[6])
Advance state of the RNG by leaps and bounds.
Definition: rng-stream.cc:386
double m_currentState[6]
The RNG state vector.
Definition: rng-stream.h:85
const double a21
Second component multiplier of n - 1 value.
Definition: rng-stream.cc:72
RngStream(uint32_t seed, uint64_t stream, uint64_t substream)
Construct from explicit seed, stream and substream values.
Definition: rng-stream.cc:363
Every class exported by the ns3 library is enclosed in the ns3 namespace.
ns3::RngStream declaration.
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:294
double RandU01(void)
Generate the next random number for this stream.
Definition: rng-stream.cc:332
Matrix a2[190]
Second component transition matrix powers.
Definition: rng-stream.cc:285
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:151
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:247
const double norm
Normalization to obtain randoms on [0,1).
Definition: rng-stream.cc:63
const double m2
Second component modulus, 232 - 22853.
Definition: rng-stream.cc:60
Debug message logging.
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:218
void PowerOfTwoMatrix(int n, Matrix a1p, Matrix a2p)
Get the transition matrices raised to a power of 2.
Definition: rng-stream.cc:312
const double a12
First component multiplier of n - 2 value.
Definition: rng-stream.cc:66
Matrix a1[190]
First component transition matrix powers.
Definition: rng-stream.cc:284