equal
deleted
inserted
replaced
22 |
22 |
23 #include <cstdlib> |
23 #include <cstdlib> |
24 #include <iostream> |
24 #include <iostream> |
25 #include "rng-stream.h" |
25 #include "rng-stream.h" |
26 #include "fatal-error.h" |
26 #include "fatal-error.h" |
|
27 #include "log.h" |
|
28 |
|
29 NS_LOG_COMPONENT_DEFINE ("RngStream"); |
27 |
30 |
28 namespace |
31 namespace |
29 { |
32 { |
30 typedef double Matrix[3][3]; |
33 typedef double Matrix[3][3]; |
31 |
34 |
68 //------------------------------------------------------------------------- |
71 //------------------------------------------------------------------------- |
69 // Return (a*s + c) MOD m; a, s, c and m must be < 2^35 |
72 // Return (a*s + c) MOD m; a, s, c and m must be < 2^35 |
70 // |
73 // |
71 double MultModM (double a, double s, double c, double m) |
74 double MultModM (double a, double s, double c, double m) |
72 { |
75 { |
|
76 NS_LOG_FUNCTION (a << s << c << m); |
73 double v; |
77 double v; |
74 int32_t a1; |
78 int32_t a1; |
75 |
79 |
76 v = a * s + c; |
80 v = a * s + c; |
77 |
81 |
103 // Works also when v = s. |
107 // Works also when v = s. |
104 // |
108 // |
105 void MatVecModM (const Matrix A, const double s[3], double v[3], |
109 void MatVecModM (const Matrix A, const double s[3], double v[3], |
106 double m) |
110 double m) |
107 { |
111 { |
|
112 NS_LOG_FUNCTION (&A << s << v << m); |
108 int i; |
113 int i; |
109 double x[3]; // Necessary if v = s |
114 double x[3]; // Necessary if v = s |
110 |
115 |
111 for (i = 0; i < 3; ++i) |
116 for (i = 0; i < 3; ++i) |
112 { |
117 { |
126 // Note: works also if A = C or B = C or A = B = C. |
131 // Note: works also if A = C or B = C or A = B = C. |
127 // |
132 // |
128 void MatMatModM (const Matrix A, const Matrix B, |
133 void MatMatModM (const Matrix A, const Matrix B, |
129 Matrix C, double m) |
134 Matrix C, double m) |
130 { |
135 { |
|
136 NS_LOG_FUNCTION (&A << &B << &C << m); |
131 int i, j; |
137 int i, j; |
132 double V[3]; |
138 double V[3]; |
133 Matrix W; |
139 Matrix W; |
134 |
140 |
135 for (i = 0; i < 3; ++i) |
141 for (i = 0; i < 3; ++i) |
157 //------------------------------------------------------------------------- |
163 //------------------------------------------------------------------------- |
158 // Compute the matrix B = (A^(2^e) Mod m); works also if A = B. |
164 // Compute the matrix B = (A^(2^e) Mod m); works also if A = B. |
159 // |
165 // |
160 void MatTwoPowModM (const Matrix src, Matrix dst, double m, int32_t e) |
166 void MatTwoPowModM (const Matrix src, Matrix dst, double m, int32_t e) |
161 { |
167 { |
|
168 NS_LOG_FUNCTION (&src << &dst << m << e); |
162 int i, j; |
169 int i, j; |
163 |
170 |
164 /* initialize: dst = src */ |
171 /* initialize: dst = src */ |
165 for (i = 0; i < 3; ++i) |
172 for (i = 0; i < 3; ++i) |
166 { |
173 { |
181 // Compute the matrix B = (A^n Mod m); works even if A = B. |
188 // Compute the matrix B = (A^n Mod m); works even if A = B. |
182 // |
189 // |
183 /* |
190 /* |
184 void MatPowModM (const double A[3][3], double B[3][3], double m, int32_t n) |
191 void MatPowModM (const double A[3][3], double B[3][3], double m, int32_t n) |
185 { |
192 { |
|
193 NS_LOG_FUNCTION (A << B << m << n); |
186 int i, j; |
194 int i, j; |
187 double W[3][3]; |
195 double W[3][3]; |
188 |
196 |
189 // initialize: W = A; B = I |
197 // initialize: W = A; B = I |
190 for (i = 0; i < 3; ++i) |
198 for (i = 0; i < 3; ++i) |
231 } |
239 } |
232 return precalculated; |
240 return precalculated; |
233 } |
241 } |
234 void PowerOfTwoMatrix (int n, Matrix a1p, Matrix a2p) |
242 void PowerOfTwoMatrix (int n, Matrix a1p, Matrix a2p) |
235 { |
243 { |
|
244 NS_LOG_FUNCTION (n << &a1p << &a2p); |
236 static struct Precalculated constants = PowerOfTwoConstants (); |
245 static struct Precalculated constants = PowerOfTwoConstants (); |
237 for (int i = 0; i < 3; i ++) |
246 for (int i = 0; i < 3; i ++) |
238 { |
247 { |
239 for (int j = 0; j < 3; j++) |
248 for (int j = 0; j < 3; j++) |
240 { |
249 { |
251 //------------------------------------------------------------------------- |
260 //------------------------------------------------------------------------- |
252 // Generate the next random number. |
261 // Generate the next random number. |
253 // |
262 // |
254 double RngStream::RandU01 () |
263 double RngStream::RandU01 () |
255 { |
264 { |
|
265 NS_LOG_FUNCTION (this); |
256 int32_t k; |
266 int32_t k; |
257 double p1, p2, u; |
267 double p1, p2, u; |
258 |
268 |
259 /* Component 1 */ |
269 /* Component 1 */ |
260 p1 = a12 * m_currentState[1] - a13n * m_currentState[0]; |
270 p1 = a12 * m_currentState[1] - a13n * m_currentState[0]; |
282 return u; |
292 return u; |
283 } |
293 } |
284 |
294 |
285 RngStream::RngStream (uint32_t seedNumber, uint64_t stream, uint64_t substream) |
295 RngStream::RngStream (uint32_t seedNumber, uint64_t stream, uint64_t substream) |
286 { |
296 { |
|
297 NS_LOG_FUNCTION (this << seedNumber << stream << substream); |
287 if (seedNumber >= m1 || seedNumber >= m2 || seedNumber == 0) |
298 if (seedNumber >= m1 || seedNumber >= m2 || seedNumber == 0) |
288 { |
299 { |
289 NS_FATAL_ERROR ("invalid Seed " << seedNumber); |
300 NS_FATAL_ERROR ("invalid Seed " << seedNumber); |
290 } |
301 } |
291 for (int i = 0; i < 6; ++i) |
302 for (int i = 0; i < 6; ++i) |
296 AdvanceNthBy (substream, 76, m_currentState); |
307 AdvanceNthBy (substream, 76, m_currentState); |
297 } |
308 } |
298 |
309 |
299 RngStream::RngStream(const RngStream& r) |
310 RngStream::RngStream(const RngStream& r) |
300 { |
311 { |
|
312 NS_LOG_FUNCTION (this << &r); |
301 for (int i = 0; i < 6; ++i) |
313 for (int i = 0; i < 6; ++i) |
302 { |
314 { |
303 m_currentState[i] = r.m_currentState[i]; |
315 m_currentState[i] = r.m_currentState[i]; |
304 } |
316 } |
305 } |
317 } |
306 |
318 |
307 void |
319 void |
308 RngStream::AdvanceNthBy (uint64_t nth, int by, double state[6]) |
320 RngStream::AdvanceNthBy (uint64_t nth, int by, double state[6]) |
309 { |
321 { |
|
322 NS_LOG_FUNCTION (this << &nth << by << state); |
310 Matrix matrix1,matrix2; |
323 Matrix matrix1,matrix2; |
311 for (int i = 0; i < 64; i++) |
324 for (int i = 0; i < 64; i++) |
312 { |
325 { |
313 int nbit = 63 - i; |
326 int nbit = 63 - i; |
314 int bit = (nth >> nbit) & 0x1; |
327 int bit = (nth >> nbit) & 0x1; |