A Discrete-Event Network Simulator
API
Loading...
Searching...
No Matches
rng-stream.cc
Go to the documentation of this file.
1//
2// Copyright (C) 2001 Pierre L'Ecuyer (lecuyer@iro.umontreal.ca)
3//
4// SPDX-License-Identifier: GPL-2.0-only
5//
6// Modified for ns-3 by:
7// - Rajib Bhattacharjea<raj.b@gatech.edu>
8// - Mathieu Lacage <mathieu.lacage@gmail.com>
9//
10
11#include "rng-stream.h"
12
13#include "fatal-error.h"
14#include "log.h"
15
16#include <cstdlib>
17#include <iostream>
18
19/**
20 * \file
21 * \ingroup rngimpl
22 * ns3::RngStream and MRG32k3a implementations.
23 */
24
25namespace ns3
26{
27
28// Note: Logging in this file is largely avoided due to the
29// number of calls that are made to these functions and the possibility
30// of causing recursions leading to stack overflow
31NS_LOG_COMPONENT_DEFINE("RngStream");
32
33} // namespace ns3
34
35/**
36 * \ingroup rngimpl
37 * @{
38 */
39/** Namespace for MRG32k3a implementation details. */
40namespace MRG32k3a
41{
42
43// clang-format off
44
45/** Type for 3x3 matrix of doubles. */
46typedef double Matrix[3][3];
47
48/** First component modulus, 2<sup>32</sup> - 209. */
49const double m1 = 4294967087.0;
50
51/** Second component modulus, 2<sup>32</sup> - 22853. */
52const double m2 = 4294944443.0;
53
54/** Normalization to obtain randoms on [0,1). */
55const double norm = 1.0 / (m1 + 1.0);
56
57/** First component multiplier of <i>n</i> - 2 value. */
58const double a12 = 1403580.0;
59
60/** First component multiplier of <i>n</i> - 3 value. */
61const double a13n = 810728.0;
62
63/** Second component multiplier of <i>n</i> - 1 value. */
64const double a21 = 527612.0;
65
66/** Second component multiplier of <i>n</i> - 3 value. */
67const double a23n = 1370589.0;
68
69/** Decomposition factor for computing a*s in less than 53 bits, 2<sup>17</sup> */
70const double two17 = 131072.0;
71
72/** IEEE-754 floating point precision, 2<sup>53</sup> */
73const double two53 = 9007199254740992.0;
74
75/** First component transition matrix. */
76const Matrix A1p0 = {
77 { 0.0, 1.0, 0.0 },
78 { 0.0, 0.0, 1.0 },
79 { -810728.0, 1403580.0, 0.0 }
80};
81
82/** Second component transition matrix. */
83const Matrix A2p0 = {
84 { 0.0, 1.0, 0.0 },
85 { 0.0, 0.0, 1.0 },
86 { -1370589.0, 0.0, 527612.0 }
87};
88
89
90//-------------------------------------------------------------------------
91/**
92 * Return (a*s + c) MOD m; a, s, c and m must be < 2^35
93 *
94 * This computes the result exactly, without exceeding the 53 bit
95 * precision of doubles.
96 *
97 * \param [in] a First multiplicative argument.
98 * \param [in] s Second multiplicative argument.
99 * \param [in] c Additive argument.
100 * \param [in] m Modulus.
101 * \returns <tt>(a*s +c) MOD m</tt>
102 */
103double MultModM (double a, double s, double c, double m)
104{
105 double v;
106 int32_t a1;
107
108 v = a * s + c;
109
110 if (v >= two53 || v <= -two53)
111 {
112 a1 = static_cast<int32_t> (a / two17);
113 a -= a1 * two17;
114 v = a1 * s;
115 a1 = static_cast<int32_t> (v / m);
116 v -= a1 * m;
117 v = v * two17 + a * s + c;
118 }
119
120 a1 = static_cast<int32_t> (v / m);
121 /* in case v < 0)*/
122 if ((v -= a1 * m) < 0.0)
123 {
124 return v += m;
125 }
126 else
127 {
128 return v;
129 }
130}
131
132
133//-------------------------------------------------------------------------
134/**
135 * Compute the vector v = A*s MOD m. Assume that -m < s[i] < m.
136 * Works also when v = s.
137 *
138 * \param [in] A Matrix argument, 3x3.
139 * \param [in] s Three component input vector.
140 * \param [out] v Three component output vector.
141 * \param [in] m Modulus.
142 */
143void MatVecModM (const Matrix A, const double s[3], double v[3],
144 double m)
145{
146 int i;
147 double x[3]; // Necessary if v = s
148
149 for (i = 0; i < 3; ++i)
150 {
151 x[i] = MultModM (A[i][0], s[0], 0.0, m);
152 x[i] = MultModM (A[i][1], s[1], x[i], m);
153 x[i] = MultModM (A[i][2], s[2], x[i], m);
154 }
155 for (i = 0; i < 3; ++i)
156 {
157 v[i] = x[i];
158 }
159}
160
161
162//-------------------------------------------------------------------------
163/**
164 * Compute the matrix C = A*B MOD m. Assume that -m < s[i] < m.
165 * Note: works also if A = C or B = C or A = B = C.
166 *
167 * \param [in] A First matrix argument.
168 * \param [in] B Second matrix argument.
169 * \param [out] C Result matrix.
170 * \param [in] m Modulus.
171 */
172void MatMatModM (const Matrix A, const Matrix B,
173 Matrix C, double m)
174{
175 int i;
176 int j;
177 double V[3];
178 Matrix W;
179
180 for (i = 0; i < 3; ++i)
181 {
182 for (j = 0; j < 3; ++j)
183 {
184 V[j] = B[j][i];
185 }
186 MatVecModM (A, V, V, m);
187 for (j = 0; j < 3; ++j)
188 {
189 W[j][i] = V[j];
190 }
191 }
192 for (i = 0; i < 3; ++i)
193 {
194 for (j = 0; j < 3; ++j)
195 {
196 C[i][j] = W[i][j];
197 }
198 }
199}
200
201
202//-------------------------------------------------------------------------
203/**
204 * Compute the matrix B = (A^(2^e) Mod m); works also if A = B.
205 *
206 * \param [in] src Matrix input argument \c A.
207 * \param [out] dst Matrix output \c B.
208 * \param [in] m Modulus.
209 * \param [in] e The exponent.
210 */
211void MatTwoPowModM (const Matrix src, Matrix dst, double m, int32_t e)
212{
213 int i;
214 int j;
215
216 /* initialize: dst = src */
217 for (i = 0; i < 3; ++i)
218 {
219 for (j = 0; j < 3; ++j)
220 {
221 dst[i][j] = src[i][j];
222 }
223 }
224 /* Compute dst = src^(2^e) mod m */
225 for (i = 0; i < e; i++)
226 {
227 MatMatModM (dst, dst, dst, m);
228 }
229}
230
231
232//-------------------------------------------------------------------------
233/**
234 * Compute the matrix B = (A^n Mod m); works even if A = B.
235 *
236 * \param [in] A Matrix input argument.
237 * \param [out] B Matrix output.
238 * \param [in] m Modulus.
239 * \param [in] n Exponent.
240 */
241void MatPowModM (const double A[3][3], double B[3][3], double m, int32_t n)
242{
243 int i;
244 int j;
245 double W[3][3];
246
247 // initialize: W = A; B = I
248 for (i = 0; i < 3; ++i)
249 {
250 for (j = 0; j < 3; ++j)
251 {
252 W[i][j] = A[i][j];
253 B[i][j] = 0.0;
254 }
255 }
256 for (j = 0; j < 3; ++j)
257 {
258 B[j][j] = 1.0;
259 }
260
261 // Compute B = A^n mod m using the binary decomposition of n
262 while (n > 0)
263 {
264 if (n % 2)
265 {
266 MatMatModM (W, B, B, m);
267 }
268 MatMatModM (W, W, W, m);
269 n /= 2;
270 }
271}
272
273/**
274 * The transition matrices of the two MRG components
275 * (in matrix form), raised to all powers of 2 from 1 to 191
276 */
278{
279 Matrix a1[190]; //!< First component transition matrix powers.
280 Matrix a2[190]; //!< Second component transition matrix powers.
281};
282
283/**
284 * Compute the transition matrices of the two MRG components
285 * raised to all powers of 2 from 1 to 191.
286 *
287 * \returns The precalculated powers of the transition matrices.
288 */
290{
291 Precalculated precalculated;
292 for (int i = 0; i < 190; i++)
293 {
294 int power = i + 1;
295 MatTwoPowModM (A1p0, precalculated.a1[i], m1, power);
296 MatTwoPowModM (A2p0, precalculated.a2[i], m2, power);
297 }
298 return precalculated;
299}
300/**
301 * Get the transition matrices raised to a power of 2.
302 *
303 * \param [in] n The power of 2.
304 * \param [out] a1p The first transition matrix power.
305 * \param [out] a2p The second transition matrix power.
306 */
307void PowerOfTwoMatrix (int n, Matrix a1p, Matrix a2p)
308{
309 static Precalculated constants = PowerOfTwoConstants ();
310 for (int i = 0; i < 3; i ++)
311 {
312 for (int j = 0; j < 3; j++)
313 {
314 a1p[i][j] = constants.a1[n-1][i][j];
315 a2p[i][j] = constants.a2[n-1][i][j];
316 }
317 }
318}
319
320} // namespace MRG32k3a
321
322// clang-format on
323
324namespace ns3
325{
326
327using namespace MRG32k3a;
328
329double
331{
332 int32_t k;
333 double p1;
334 double p2;
335 double 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 }
347 m_currentState[2] = p1;
348
349 /* Component 2 */
350 p2 = a21 * m_currentState[5] - a23n * m_currentState[3];
351 k = static_cast<int32_t>(p2 / m2);
352 p2 -= k * m2;
353 if (p2 < 0.0)
354 {
355 p2 += m2;
356 }
359 m_currentState[5] = p2;
360
361 /* Combination */
362 u = ((p1 > p2) ? (p1 - p2) * MRG32k3a::norm : (p1 - p2 + m1) * MRG32k3a::norm);
363
364 return u;
365}
366
367RngStream::RngStream(uint32_t seedNumber, uint64_t stream, uint64_t substream)
368{
369 if (seedNumber >= m1 || seedNumber >= m2 || seedNumber == 0)
370 {
371 NS_FATAL_ERROR("invalid Seed " << seedNumber);
372 }
373 for (int i = 0; i < 6; ++i)
374 {
375 m_currentState[i] = seedNumber;
376 }
377 AdvanceNthBy(stream, 127, m_currentState);
378 AdvanceNthBy(substream, 76, m_currentState);
379}
380
382{
383 for (int i = 0; i < 6; ++i)
384 {
386 }
387}
388
389void
390RngStream::AdvanceNthBy(uint64_t nth, int by, double state[6])
391{
392 Matrix matrix1;
393 Matrix matrix2;
394 for (int i = 0; i < 64; i++)
395 {
396 int nbit = 63 - i;
397 int bit = (nth >> nbit) & 0x1;
398 if (bit)
399 {
400 PowerOfTwoMatrix(by + nbit, matrix1, matrix2);
401 MatVecModM(matrix1, state, state, m1);
402 MatVecModM(matrix2, &state[3], &state[3], m2);
403 }
404 }
405}
406
407} // namespace ns3
408
409/**@}*/ // \ingroup rngimpl
Combined Multiple-Recursive Generator MRG32k3a.
Definition rng-stream.h:39
RngStream(uint32_t seed, uint64_t stream, uint64_t substream)
Construct from explicit seed, stream and substream values.
double m_currentState[6]
The RNG state vector.
Definition rng-stream.h:74
void AdvanceNthBy(uint64_t nth, int by, double state[6])
Advance state of the RNG by leaps and bounds.
double RandU01()
Generate the next random number for this stream.
NS_FATAL_x macro definitions.
#define NS_FATAL_ERROR(msg)
Report a fatal error with a message and terminate.
#define NS_LOG_COMPONENT_DEFINE(name)
Define a Log component with a specific name.
Definition log.h:191
Debug message logging.
Namespace for MRG32k3a implementation details.
Definition rng-stream.cc:41
const double a21
Second component multiplier of n - 1 value.
Definition rng-stream.cc:64
const double m1
First component modulus, 232 - 209.
Definition rng-stream.cc:49
void MatMatModM(const Matrix A, const Matrix B, Matrix C, double m)
Compute the matrix C = A*B MOD m.
const Matrix A1p0
First component transition matrix.
Definition rng-stream.cc:76
const double a23n
Second component multiplier of n - 3 value.
Definition rng-stream.cc:67
const double two17
Decomposition factor for computing a*s in less than 53 bits, 217
Definition rng-stream.cc:70
void PowerOfTwoMatrix(int n, Matrix a1p, Matrix a2p)
Get the transition matrices raised to a power of 2.
const double norm
Normalization to obtain randoms on [0,1).
Definition rng-stream.cc:55
const double a12
First component multiplier of n - 2 value.
Definition rng-stream.cc:58
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.
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.
const Matrix A2p0
Second component transition matrix.
Definition rng-stream.cc:83
double Matrix[3][3]
Type for 3x3 matrix of doubles.
Definition rng-stream.cc:46
const double m2
Second component modulus, 232 - 22853.
Definition rng-stream.cc:52
const double a13n
First component multiplier of n - 3 value.
Definition rng-stream.cc:61
const double two53
IEEE-754 floating point precision, 253
Definition rng-stream.cc:73
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.
void MatVecModM(const Matrix A, const double s[3], double v[3], double m)
Compute the vector v = A*s MOD m.
Precalculated PowerOfTwoConstants()
Compute the transition matrices of the two MRG components raised to all powers of 2 from 1 to 191.
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 ...
Matrix a1[190]
First component transition matrix powers.
Matrix a2[190]
Second component transition matrix powers.