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 /**
30 * \file
31 * \ingroup rngimpl
32 * ns3::RngStream and MRG32k3a implementations.
33 */
34
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
45 /**
46 * \ingroup rngimpl
47 * @{
48 */
49 /** Namespace for MRG32k3a implementation details. */
50 namespace MRG32k3a {
51
52 // *NS_CHECK_STYLE_OFF*
53
54 /** Type for 3x3 matrix of doubles. */
55 typedef double Matrix[3][3];
56
57 /** First component modulus, 2<sup>32</sup> - 209. */
58 const double m1 = 4294967087.0;
59
60 /** Second component modulus, 2<sup>32</sup> - 22853. */
61 const double m2 = 4294944443.0;
62
63 /** Normalization to obtain randoms on [0,1). */
64 const double norm = 1.0 / (m1 + 1.0);
65
66 /** First component multiplier of <i>n</i> - 2 value. */
67 const double a12 = 1403580.0;
68
69 /** First component multiplier of <i>n</i> - 3 value. */
70 const double a13n = 810728.0;
71
72 /** Second component multiplier of <i>n</i> - 1 value. */
73 const double a21 = 527612.0;
74
75 /** Second component multiplier of <i>n</i> - 3 value. */
76 const double a23n = 1370589.0;
77
78 /** Decomposition factor for computing a*s in less than 53 bits, 2<sup>17</sup> */
79 const double two17 = 131072.0;
80
81 /** IEEE-754 floating point precision, 2<sup>53</sup> */
82 const double two53 = 9007199254740992.0;
83
84 /** First component transition matrix. */
85 const 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
91 /** Second component transition matrix. */
92 const 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 //-------------------------------------------------------------------------
100 /**
101 * Return (a*s + c) MOD m; a, s, c and m must be < 2^35
102 *
103 * This computes the result exactly, without exceeding the 53 bit
104 * precision of doubles.
105 *
106 * \param [in] a First multiplicative argument.
107 * \param [in] s Second multiplicative argument.
108 * \param [in] c Additive argument.
109 * \param [in] m Modulus.
110 * \returns <tt>(a*s +c) MOD m</tt>
111 */
MultModM(double a,double s,double c,double m)112 double 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 //-------------------------------------------------------------------------
143 /**
144 * Compute the vector v = A*s MOD m. Assume that -m < s[i] < m.
145 * Works also when v = s.
146 *
147 * \param [in] A Matrix argument, 3x3.
148 * \param [in] s Three component input vector.
149 * \param [out] v Three component output vector.
150 * \param [in] m Modulus.
151 */
MatVecModM(const Matrix A,const double s[3],double v[3],double m)152 void 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 //-------------------------------------------------------------------------
172 /**
173 * Compute the matrix C = A*B MOD m. Assume that -m < s[i] < m.
174 * Note: works also if A = C or B = C or A = B = C.
175 *
176 * \param [in] A First matrix argument.
177 * \param [in] B Second matrix argument.
178 * \param [out] C Result matrix.
179 * \param [in] m Modulus.
180 */
MatMatModM(const Matrix A,const Matrix B,Matrix C,double m)181 void 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 //-------------------------------------------------------------------------
211 /**
212 * Compute the matrix B = (A^(2^e) Mod m); works also if A = B.
213 *
214 * \param [in] src Matrix input argument \c A.
215 * \param [out] dst Matrix output \c B.
216 * \param [in] m Modulus.
217 * \param [in] e The exponent.
218 */
MatTwoPowModM(const Matrix src,Matrix dst,double m,int32_t e)219 void 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 //-------------------------------------------------------------------------
240 /**
241 * Compute the matrix B = (A^n Mod m); works even if A = B.
242 *
243 * \param [in] A Matrix input argument.
244 * \param [out] B Matrix output.
245 * \param [in] m Modulus.
246 * \param [in] n Exponent.
247 */
MatPowModM(const double A[3][3],double B[3][3],double m,int32_t n)248 void 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
279 /**
280 * The transition matrices of the two MRG components
281 * (in matrix form), raised to all powers of 2 from 1 to 191
282 */
283 struct Precalculated
284 {
285 Matrix a1[190]; //!< First component transition matrix powers.
286 Matrix a2[190]; //!< Second component transition matrix powers.
287 };
288
289 /**
290 * Compute the transition matrices of the two MRG components
291 * raised to all powers of 2 from 1 to 191.
292 *
293 * \returns The precalculated powers of the transition matrices.
294 */
PowerOfTwoConstants(void)295 struct Precalculated PowerOfTwoConstants (void)
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 }
306 /**
307 * Get the transition matrices raised to a power of 2.
308 *
309 * \param [in] n The power of 2.
310 * \param [out] a1p The first transition matrix power.
311 * \param [out] a2p The second transition matrix power.
312 */
PowerOfTwoMatrix(int n,Matrix a1p,Matrix a2p)313 void 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
331 namespace ns3 {
332
333 using namespace MRG32k3a;
334
RandU01()335 double RngStream::RandU01 ()
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 }
348 m_currentState[0] = m_currentState[1];
349 m_currentState[1] = m_currentState[2];
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 }
360 m_currentState[3] = m_currentState[4];
361 m_currentState[4] = m_currentState[5];
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
RngStream(uint32_t seedNumber,uint64_t stream,uint64_t substream)370 RngStream::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
RngStream(const RngStream & r)384 RngStream::RngStream (const RngStream& r)
385 {
386 for (int i = 0; i < 6; ++i)
387 {
388 m_currentState[i] = r.m_currentState[i];
389 }
390 }
391
392 void
AdvanceNthBy(uint64_t nth,int by,double state[6])393 RngStream::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
411 /**@}*/ // \ingroup rngimpl
412