A Discrete-Event Network Simulator
API
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
rng-stream.cc
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 
29 // Note: Logging in this file is largely avoided due to the
30 // number of calls that are made to these functions and the possibility
31 // of causing recursions leading to stack overflow
32 
33 NS_LOG_COMPONENT_DEFINE ("RngStream");
34 
35 namespace
36 {
37 typedef double Matrix[3][3];
38 
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; /* 1 / 2^24 */
49 
50 const Matrix InvA1 = { // Inverse of A1p0
51  { 184888585.0, 0.0, 1945170933.0 },
52  { 1.0, 0.0, 0.0 },
53  { 0.0, 1.0, 0.0 }
54 };
55 
56 const Matrix InvA2 = { // Inverse of A2p0
57  { 0.0, 360363334.0, 4225571728.0 },
58  { 1.0, 0.0, 0.0 },
59  { 0.0, 1.0, 0.0 }
60 };
61 
62 const Matrix A1p0 = {
63  { 0.0, 1.0, 0.0 },
64  { 0.0, 0.0, 1.0 },
65  { -810728.0, 1403580.0, 0.0 }
66 };
67 
68 const Matrix A2p0 = {
69  { 0.0, 1.0, 0.0 },
70  { 0.0, 0.0, 1.0 },
71  { -1370589.0, 0.0, 527612.0 }
72 };
73 
74 
75 //-------------------------------------------------------------------------
76 // Return (a*s + c) MOD m; a, s, c and m must be < 2^35
77 //
78 double MultModM (double a, double s, double c, double m)
79 {
80  double v;
81  int32_t a1;
82 
83  v = a * s + c;
84 
85  if (v >= two53 || v <= -two53)
86  {
87  a1 = static_cast<int32_t> (a / two17);
88  a -= a1 * two17;
89  v = a1 * s;
90  a1 = static_cast<int32_t> (v / m);
91  v -= a1 * m;
92  v = v * two17 + a * s + c;
93  }
94 
95  a1 = static_cast<int32_t> (v / m);
96  /* in case v < 0)*/
97  if ((v -= a1 * m) < 0.0)
98  {
99  return v += m;
100  }
101  else
102  {
103  return v;
104  }
105 }
106 
107 
108 //-------------------------------------------------------------------------
109 // Compute the vector v = A*s MOD m. Assume that -m < s[i] < m.
110 // Works also when v = s.
111 //
112 void MatVecModM (const Matrix A, const double s[3], double v[3],
113  double m)
114 {
115  int i;
116  double x[3]; // Necessary if v = s
117 
118  for (i = 0; i < 3; ++i)
119  {
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);
123  }
124  for (i = 0; i < 3; ++i)
125  {
126  v[i] = x[i];
127  }
128 }
129 
130 
131 //-------------------------------------------------------------------------
132 // Compute the matrix C = A*B MOD m. Assume that -m < s[i] < m.
133 // Note: works also if A = C or B = C or A = B = C.
134 //
135 void MatMatModM (const Matrix A, const Matrix B,
136  Matrix C, double m)
137 {
138  int i, j;
139  double V[3];
140  Matrix W;
141 
142  for (i = 0; i < 3; ++i)
143  {
144  for (j = 0; j < 3; ++j)
145  {
146  V[j] = B[j][i];
147  }
148  MatVecModM (A, V, V, m);
149  for (j = 0; j < 3; ++j)
150  {
151  W[j][i] = V[j];
152  }
153  }
154  for (i = 0; i < 3; ++i)
155  {
156  for (j = 0; j < 3; ++j)
157  {
158  C[i][j] = W[i][j];
159  }
160  }
161 }
162 
163 
164 //-------------------------------------------------------------------------
165 // Compute the matrix B = (A^(2^e) Mod m); works also if A = B.
166 //
167 void MatTwoPowModM (const Matrix src, Matrix dst, double m, int32_t e)
168 {
169  int i, j;
170 
171  /* initialize: dst = src */
172  for (i = 0; i < 3; ++i)
173  {
174  for (j = 0; j < 3; ++j)
175  {
176  dst[i][j] = src[i][j];
177  }
178  }
179  /* Compute dst = src^(2^e) mod m */
180  for (i = 0; i < e; i++)
181  {
182  MatMatModM (dst, dst, dst, m);
183  }
184 }
185 
186 
187 //-------------------------------------------------------------------------
188 // Compute the matrix B = (A^n Mod m); works even if A = B.
189 //
190 /*
191 void MatPowModM (const double A[3][3], double B[3][3], double m, int32_t n)
192 {
193  NS_LOG_FUNCTION (A << B << m << n);
194  int i, j;
195  double W[3][3];
196 
197  // initialize: W = A; B = I
198  for (i = 0; i < 3; ++i)
199  {
200  for (j = 0; j < 3; ++j)
201  {
202  W[i][j] = A[i][j];
203  B[i][j] = 0.0;
204  }
205  }
206  for (j = 0; j < 3; ++j)
207  {
208  B[j][j] = 1.0;
209  }
210 
211  // Compute B = A^n mod m using the binary decomposition of n
212  while (n > 0)
213  {
214  if (n % 2)
215  {
216  MatMatModM (W, B, B, m);
217  }
218  MatMatModM (W, W, W, m);
219  n /= 2;
220  }
221 }
222 */
223 
224 // The following are the transition matrices of the two MRG components
225 // (in matrix form), raised to all powers of 2 from 1 to 191
227 {
228  Matrix a1[190];
229  Matrix a2[190];
230 };
231 struct Precalculated PowerOfTwoConstants (void)
232 {
233  struct Precalculated precalculated;
234  for (int i = 0; i < 190; i++)
235  {
236  int power = i + 1;
237  MatTwoPowModM (A1p0, precalculated.a1[i], m1, power);
238  MatTwoPowModM (A2p0, precalculated.a2[i], m2, power);
239  }
240  return precalculated;
241 }
242 void PowerOfTwoMatrix (int n, Matrix a1p, Matrix a2p)
243 {
244  static struct Precalculated constants = PowerOfTwoConstants ();
245  for (int i = 0; i < 3; i ++)
246  {
247  for (int j = 0; j < 3; j++)
248  {
249  a1p[i][j] = constants.a1[n-1][i][j];
250  a2p[i][j] = constants.a2[n-1][i][j];
251  }
252  }
253 }
254 
255 } // end of anonymous namespace
256 
257 
258 namespace ns3 {
259 //-------------------------------------------------------------------------
260 // Generate the next random number.
261 //
263 {
264  int32_t k;
265  double p1, p2, u;
266 
267  /* Component 1 */
268  p1 = a12 * m_currentState[1] - a13n * m_currentState[0];
269  k = static_cast<int32_t> (p1 / m1);
270  p1 -= k * m1;
271  if (p1 < 0.0)
272  {
273  p1 += m1;
274  }
275  m_currentState[0] = m_currentState[1]; m_currentState[1] = m_currentState[2]; m_currentState[2] = p1;
276 
277  /* Component 2 */
278  p2 = a21 * m_currentState[5] - a23n * m_currentState[3];
279  k = static_cast<int32_t> (p2 / m2);
280  p2 -= k * m2;
281  if (p2 < 0.0)
282  {
283  p2 += m2;
284  }
285  m_currentState[3] = m_currentState[4]; m_currentState[4] = m_currentState[5]; m_currentState[5] = p2;
286 
287  /* Combination */
288  u = ((p1 > p2) ? (p1 - p2) * norm : (p1 - p2 + m1) * norm);
289 
290  return u;
291 }
292 
293 RngStream::RngStream (uint32_t seedNumber, uint64_t stream, uint64_t substream)
294 {
295  if (seedNumber >= m1 || seedNumber >= m2 || seedNumber == 0)
296  {
297  NS_FATAL_ERROR ("invalid Seed " << seedNumber);
298  }
299  for (int i = 0; i < 6; ++i)
300  {
301  m_currentState[i] = seedNumber;
302  }
303  AdvanceNthBy (stream, 127, m_currentState);
304  AdvanceNthBy (substream, 76, m_currentState);
305 }
306 
307 RngStream::RngStream(const RngStream& r)
308 {
309  for (int i = 0; i < 6; ++i)
310  {
311  m_currentState[i] = r.m_currentState[i];
312  }
313 }
314 
315 void
316 RngStream::AdvanceNthBy (uint64_t nth, int by, double state[6])
317 {
318  Matrix matrix1,matrix2;
319  for (int i = 0; i < 64; i++)
320  {
321  int nbit = 63 - i;
322  int bit = (nth >> nbit) & 0x1;
323  if (bit)
324  {
325  PowerOfTwoMatrix(by + nbit, matrix1, matrix2);
326  MatVecModM (matrix1, state, state, m1);
327  MatVecModM (matrix2, &state[3], &state[3], m2);
328  }
329  }
330 }
331 
332 } // namespace ns3
#define NS_LOG_COMPONENT_DEFINE(name)
Definition: log.h:122
#define NS_FATAL_ERROR(msg)
fatal error handling
Definition: fatal-error.h:72
double RandU01(void)
Definition: rng-stream.cc:262