1 | /* -*- C++ -*- |
---|
2 | * |
---|
3 | * This file is a part of LEMON, a generic C++ optimization library |
---|
4 | * |
---|
5 | * Copyright (C) 2003-2006 |
---|
6 | * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
---|
7 | * (Egervary Research Group on Combinatorial Optimization, EGRES). |
---|
8 | * |
---|
9 | * Permission to use, modify and distribute this software is granted |
---|
10 | * provided that this copyright notice appears in all copies. For |
---|
11 | * precise terms see the accompanying LICENSE file. |
---|
12 | * |
---|
13 | * This software is provided "AS IS" with no warranty of any kind, |
---|
14 | * express or implied, and with no claim as to its suitability for any |
---|
15 | * purpose. |
---|
16 | * |
---|
17 | */ |
---|
18 | |
---|
19 | #ifndef LEMON_RANDOM_H |
---|
20 | #define LEMON_RANDOM_H |
---|
21 | |
---|
22 | #include <algorithm> |
---|
23 | |
---|
24 | #include <ctime> |
---|
25 | #include <cmath> |
---|
26 | #include <cmath> |
---|
27 | |
---|
28 | ///\ingroup misc |
---|
29 | ///\file |
---|
30 | ///\brief Mersenne Twister random number generator |
---|
31 | /// |
---|
32 | ///\author Balazs Dezso |
---|
33 | |
---|
34 | namespace lemon { |
---|
35 | |
---|
36 | #if WORD_BIT == 32 |
---|
37 | |
---|
38 | /// \ingroup misc |
---|
39 | /// |
---|
40 | /// \brief Mersenne Twister random number generator |
---|
41 | /// |
---|
42 | /// The Mersenne Twister is a twisted generalized feedback |
---|
43 | /// shift-register generator of Matsumoto and Nishimura. The period of |
---|
44 | /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in |
---|
45 | /// 623 dimensions. The time performance of this generator is comparable |
---|
46 | /// to the commonly used generators. |
---|
47 | /// |
---|
48 | /// \author Balazs Dezso |
---|
49 | class Random { |
---|
50 | |
---|
51 | static const int length = 624; |
---|
52 | static const int shift = 397; |
---|
53 | |
---|
54 | public: |
---|
55 | |
---|
56 | static const unsigned long min = 0ul; |
---|
57 | static const unsigned long max = ~0ul; |
---|
58 | |
---|
59 | /// \brief Constructor |
---|
60 | /// |
---|
61 | /// Constructor with time dependent seeding. |
---|
62 | Random() { initState(std::time(0)); } |
---|
63 | |
---|
64 | /// \brief Constructor |
---|
65 | /// |
---|
66 | /// Constructor |
---|
67 | Random(unsigned long seed) { initState(seed); } |
---|
68 | |
---|
69 | /// \brief Copy constructor |
---|
70 | /// |
---|
71 | /// Copy constructor. The generated sequence will be identical to |
---|
72 | /// the other sequence. |
---|
73 | Random(const Random& other) { |
---|
74 | std::copy(other.state, other.state + length, state); |
---|
75 | current = state + (other.current - other.state); |
---|
76 | } |
---|
77 | |
---|
78 | /// \brief Assign operator |
---|
79 | /// |
---|
80 | /// Assign operator. The generated sequence will be identical to |
---|
81 | /// the other sequence. |
---|
82 | Random& operator=(const Random& other) { |
---|
83 | if (&other != this) { |
---|
84 | std::copy(other.state, other.state + length, state); |
---|
85 | current = state + (other.current - other.state); |
---|
86 | } |
---|
87 | return *this; |
---|
88 | } |
---|
89 | |
---|
90 | /// \brief Returns an unsigned random number |
---|
91 | /// |
---|
92 | /// It returns an unsigned integer random number from the range |
---|
93 | /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$. |
---|
94 | unsigned long getUnsigned() { |
---|
95 | if (current == state) fillState(); |
---|
96 | --current; |
---|
97 | unsigned long rnd = *current; |
---|
98 | rnd ^= (rnd >> 11); |
---|
99 | rnd ^= (rnd << 7) & 0x9D2C5680ul; |
---|
100 | rnd ^= (rnd << 15) & 0xEFC60000ul; |
---|
101 | rnd ^= (rnd >> 18); |
---|
102 | return rnd; |
---|
103 | } |
---|
104 | |
---|
105 | /// \brief Returns a random number |
---|
106 | /// |
---|
107 | /// It returns an integer random number from the range |
---|
108 | /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$. |
---|
109 | long getInt() { |
---|
110 | return (long)getUnsigned(); |
---|
111 | } |
---|
112 | |
---|
113 | /// \brief Returns an unsigned integer random number |
---|
114 | /// |
---|
115 | /// It returns an unsigned integer random number from the range |
---|
116 | /// \f$ \{0, 1, \dots, 2^{31} - 1\} \f$. |
---|
117 | long getNatural() { |
---|
118 | return (long)(getUnsigned() >> 1); |
---|
119 | } |
---|
120 | |
---|
121 | /// \brief Returns a random bool |
---|
122 | /// |
---|
123 | /// It returns a random bool. |
---|
124 | bool getBool() { |
---|
125 | return (bool)(getUnsigned() & 1); |
---|
126 | } |
---|
127 | |
---|
128 | /// \brief Returns a real random number |
---|
129 | /// |
---|
130 | /// It returns a real random number from the range |
---|
131 | /// \f$ [0, 1) \f$. The double will have 32 significant bits. |
---|
132 | double getReal() { |
---|
133 | return std::ldexp((double)getUnsigned(), -32); |
---|
134 | } |
---|
135 | |
---|
136 | /// \brief Returns a real random number |
---|
137 | /// |
---|
138 | /// It returns a real random number from the range |
---|
139 | /// \f$ [0, 1) \f$. The double will have 53 significant bits, |
---|
140 | /// but it is slower than the \c getReal(). |
---|
141 | double getPrecReal() { |
---|
142 | return std::ldexp((double)(getUnsigned() >> 5), -27) + |
---|
143 | std::ldexp((double)(getUnsigned() >> 6), -53); |
---|
144 | } |
---|
145 | |
---|
146 | /// \brief Returns an unsigned random number from a given range |
---|
147 | /// |
---|
148 | /// It returns an unsigned integer random number from the range |
---|
149 | /// \f$ \{0, 1, \dots, n - 1\} \f$. |
---|
150 | unsigned long getUnsigned(unsigned long n) { |
---|
151 | unsigned long mask = n - 1, rnd; |
---|
152 | mask |= mask >> 1; |
---|
153 | mask |= mask >> 2; |
---|
154 | mask |= mask >> 4; |
---|
155 | mask |= mask >> 8; |
---|
156 | mask |= mask >> 16; |
---|
157 | do rnd = getUnsigned() & mask; while (rnd >= n); |
---|
158 | return rnd; |
---|
159 | } |
---|
160 | |
---|
161 | /// \brief Returns a random number from a given range |
---|
162 | /// |
---|
163 | /// It returns an unsigned integer random number from the range |
---|
164 | /// \f$ \{0, 1, \dots, n - 1\} \f$. |
---|
165 | long getInt(long n) { |
---|
166 | long mask = n - 1, rnd; |
---|
167 | mask |= mask >> 1; |
---|
168 | mask |= mask >> 2; |
---|
169 | mask |= mask >> 4; |
---|
170 | mask |= mask >> 8; |
---|
171 | mask |= mask >> 16; |
---|
172 | do rnd = getUnsigned() & mask; while (rnd >= n); |
---|
173 | return rnd; |
---|
174 | } |
---|
175 | |
---|
176 | private: |
---|
177 | |
---|
178 | void initState(unsigned long seed) { |
---|
179 | static const unsigned long mul = 0x6c078965ul; |
---|
180 | |
---|
181 | current = state; |
---|
182 | |
---|
183 | unsigned long *curr = state + length - 1; |
---|
184 | curr[0] = seed; --curr; |
---|
185 | for (int i = 1; i < length; ++i) { |
---|
186 | curr[0] = (mul * ( curr[1] ^ (curr[1] >> 30) ) + i); |
---|
187 | --curr; |
---|
188 | } |
---|
189 | } |
---|
190 | |
---|
191 | void fillState() { |
---|
192 | static const unsigned long mask[2] = { 0x0ul, 0x9908B0DFul }; |
---|
193 | static const unsigned long loMask = (1ul << 31) - 1; |
---|
194 | static const unsigned long hiMask = ~loMask; |
---|
195 | |
---|
196 | current = state + length; |
---|
197 | |
---|
198 | register unsigned long *curr = state + length - 1; |
---|
199 | register long num; |
---|
200 | |
---|
201 | num = length - shift; |
---|
202 | while (num--) { |
---|
203 | curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ |
---|
204 | curr[- shift] ^ mask[curr[-1] & 1ul]; |
---|
205 | --curr; |
---|
206 | } |
---|
207 | num = shift - 1; |
---|
208 | while (num--) { |
---|
209 | curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ |
---|
210 | curr[length - shift] ^ mask[curr[-1] & 1ul]; |
---|
211 | --curr; |
---|
212 | } |
---|
213 | curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^ |
---|
214 | curr[length - shift] ^ mask[curr[length - 1] & 1ul]; |
---|
215 | |
---|
216 | } |
---|
217 | |
---|
218 | unsigned long *current; |
---|
219 | unsigned long state[length]; |
---|
220 | |
---|
221 | }; |
---|
222 | |
---|
223 | #elif WORD_BIT == 64 |
---|
224 | |
---|
225 | /// \brief Mersenne Twister random number generator |
---|
226 | /// |
---|
227 | /// The Mersenne Twister is a twisted generalized feedback |
---|
228 | /// shift-register generator of Matsumoto and Nishimura. The period of |
---|
229 | /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in |
---|
230 | /// 311 dimensions. The time performance of this generator is comparable |
---|
231 | /// to the commonly used generators. |
---|
232 | class Random { |
---|
233 | |
---|
234 | static const int length = 312; |
---|
235 | static const int shift = 156; |
---|
236 | |
---|
237 | public: |
---|
238 | |
---|
239 | static const unsigned long min = 0ul; |
---|
240 | static const unsigned long max = ~0ul; |
---|
241 | |
---|
242 | /// \brief Constructor |
---|
243 | /// |
---|
244 | /// Constructor with time dependent seeding. |
---|
245 | Random() { initState(std::time(0)); } |
---|
246 | |
---|
247 | /// \brief Constructor |
---|
248 | /// |
---|
249 | /// Constructor |
---|
250 | Random(unsigned long seed) { initState(seed); } |
---|
251 | |
---|
252 | /// \brief Copy constructor |
---|
253 | /// |
---|
254 | /// Copy constructor. The generated sequence will be identical to |
---|
255 | /// the other sequence. |
---|
256 | Random(const Random& other) { |
---|
257 | std::copy(other.state, other.state + length, state); |
---|
258 | current = state + (other.current - other.state); |
---|
259 | } |
---|
260 | |
---|
261 | /// \brief Assign operator |
---|
262 | /// |
---|
263 | /// Assign operator. The generated sequence will be identical to |
---|
264 | /// the other sequence. |
---|
265 | Random& operator=(const Random& other) { |
---|
266 | if (&other != this) { |
---|
267 | std::copy(other.state, other.state + length, state); |
---|
268 | current = state + (other.current - other.state); |
---|
269 | } |
---|
270 | return *this; |
---|
271 | } |
---|
272 | |
---|
273 | /// \brief Returns an unsigned random number |
---|
274 | /// |
---|
275 | /// It returns an unsigned integer random number from the range |
---|
276 | /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$. |
---|
277 | unsigned long getUnsigned() { |
---|
278 | if (current == state) fillState(); |
---|
279 | --current; |
---|
280 | unsigned long rnd = *current; |
---|
281 | rnd ^= (rnd >> 29) & 0x5555555555555555ul; |
---|
282 | rnd ^= (rnd << 17) & 0x71D67FFFEDA60000ul; |
---|
283 | rnd ^= (rnd << 37) & 0xFFF7EEE000000000ul; |
---|
284 | rnd ^= (rnd >> 43); |
---|
285 | return rnd; |
---|
286 | } |
---|
287 | |
---|
288 | /// \brief Returns a random number |
---|
289 | /// |
---|
290 | /// It returns an integer random number from the range |
---|
291 | /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$. |
---|
292 | long getInt() { |
---|
293 | return (long)getUnsigned(); |
---|
294 | } |
---|
295 | |
---|
296 | /// \brief Returns an unsigned integer random number |
---|
297 | /// |
---|
298 | /// It returns an unsigned integer random number from the range |
---|
299 | /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$. |
---|
300 | long getNatural() { |
---|
301 | return (long)(getUnsigned() >> 1); |
---|
302 | } |
---|
303 | |
---|
304 | /// \brief Returns a random bool |
---|
305 | /// |
---|
306 | /// It returns a random bool. |
---|
307 | bool getBool() { |
---|
308 | return (bool)(getUnsigned() & 1); |
---|
309 | } |
---|
310 | |
---|
311 | /// \brief Returns a real random number |
---|
312 | /// |
---|
313 | /// It returns a real random number from the range |
---|
314 | /// \f$ [0, 1) \f$. |
---|
315 | double getReal() { |
---|
316 | return std::ldexp((double)(getUnsigned() >> 11), -53); |
---|
317 | } |
---|
318 | |
---|
319 | /// \brief Returns a real random number |
---|
320 | /// |
---|
321 | /// It returns a real random number from the range |
---|
322 | /// \f$ [0, 1) \f$. This function is identical to the \c getReal(). |
---|
323 | double getPrecReal() { |
---|
324 | return getReal(); |
---|
325 | } |
---|
326 | |
---|
327 | /// \brief Returns an unsigned random number from a given range |
---|
328 | /// |
---|
329 | /// It returns an unsigned integer random number from the range |
---|
330 | /// \f$ \{0, 1, \dots, n - 1\} \f$. |
---|
331 | unsigned long getUnsigned(unsigned long n) { |
---|
332 | unsigned long mask = n - 1, rnd; |
---|
333 | mask |= mask >> 1; |
---|
334 | mask |= mask >> 2; |
---|
335 | mask |= mask >> 4; |
---|
336 | mask |= mask >> 8; |
---|
337 | mask |= mask >> 16; |
---|
338 | mask |= mask >> 32; |
---|
339 | do rnd = getUnsigned() & mask; while (rnd >= n); |
---|
340 | return rnd; |
---|
341 | } |
---|
342 | |
---|
343 | /// \brief Returns a random number from a given range |
---|
344 | /// |
---|
345 | /// It returns an unsigned integer random number from the range |
---|
346 | /// \f$ \{0, 1, \dots, n - 1\} \f$. |
---|
347 | long getInt(long n) { |
---|
348 | long mask = n - 1, rnd; |
---|
349 | mask |= mask >> 1; |
---|
350 | mask |= mask >> 2; |
---|
351 | mask |= mask >> 4; |
---|
352 | mask |= mask >> 8; |
---|
353 | mask |= mask >> 16; |
---|
354 | mask |= mask >> 32; |
---|
355 | do rnd = getUnsigned() & mask; while (rnd >= n); |
---|
356 | return rnd; |
---|
357 | } |
---|
358 | |
---|
359 | private: |
---|
360 | |
---|
361 | void initState(unsigned long seed) { |
---|
362 | |
---|
363 | static const unsigned long mul = 0x5851F42D4C957F2Dul; |
---|
364 | |
---|
365 | current = state; |
---|
366 | |
---|
367 | unsigned long *curr = state + length - 1; |
---|
368 | curr[0] = seed; --curr; |
---|
369 | for (int i = 1; i < length; ++i) { |
---|
370 | curr[0] = (mul * ( curr[1] ^ (curr[1] >> 62) ) + i); |
---|
371 | --curr; |
---|
372 | } |
---|
373 | } |
---|
374 | |
---|
375 | void fillState() { |
---|
376 | static const unsigned long mask[2] = { 0x0ul, 0xB5026F5AA96619E9ul }; |
---|
377 | static const unsigned long loMask = (1ul << 31) - 1; |
---|
378 | static const unsigned long hiMask = ~loMask; |
---|
379 | |
---|
380 | current = state + length; |
---|
381 | |
---|
382 | register unsigned long *curr = state + length - 1; |
---|
383 | register int num; |
---|
384 | |
---|
385 | num = length - shift; |
---|
386 | while (num--) { |
---|
387 | curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ |
---|
388 | curr[- shift] ^ mask[curr[-1] & 1ul]; |
---|
389 | --curr; |
---|
390 | } |
---|
391 | num = shift - 1; |
---|
392 | while (num--) { |
---|
393 | curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ |
---|
394 | curr[length - shift] ^ mask[curr[-1] & 1ul]; |
---|
395 | --curr; |
---|
396 | } |
---|
397 | curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^ |
---|
398 | curr[length - shift] ^ mask[curr[length - 1] & 1ul]; |
---|
399 | |
---|
400 | } |
---|
401 | |
---|
402 | unsigned long *current; |
---|
403 | unsigned long state[length]; |
---|
404 | |
---|
405 | }; |
---|
406 | |
---|
407 | #else |
---|
408 | |
---|
409 | /// \brief Mersenne Twister random number generator |
---|
410 | /// |
---|
411 | /// The Mersenne Twister is a twisted generalized feedback |
---|
412 | /// shift-register generator of Matsumoto and Nishimura. There is |
---|
413 | /// two different implementation in the Lemon library, one for |
---|
414 | /// 32-bit processors and one for 64-bit processors. The period of |
---|
415 | /// the generated sequence is \f$ 2^{19937} - 1 \f$, the generated |
---|
416 | /// sequence of 32-bit random numbers is equi-distributed in 623 |
---|
417 | /// dimensions. The time performance of this generator is comparable |
---|
418 | /// to the commonly used generators. |
---|
419 | class Random { |
---|
420 | public: |
---|
421 | |
---|
422 | static const unsigned long min = 0ul; |
---|
423 | static const unsigned long max = ~0ul; |
---|
424 | |
---|
425 | /// \brief Constructor |
---|
426 | /// |
---|
427 | /// Constructor with time dependent seeding. |
---|
428 | Random() {} |
---|
429 | |
---|
430 | /// \brief Constructor |
---|
431 | /// |
---|
432 | /// Constructor |
---|
433 | Random(unsigned long seed) {} |
---|
434 | |
---|
435 | /// \brief Copy constructor |
---|
436 | /// |
---|
437 | /// Copy constructor. The generated sequence will be identical to |
---|
438 | /// the other sequence. |
---|
439 | Random(const Random& other) {} |
---|
440 | |
---|
441 | /// \brief Assign operator |
---|
442 | /// |
---|
443 | /// Assign operator. The generated sequence will be identical to |
---|
444 | /// the other sequence. |
---|
445 | Random& operator=(const Random& other) { return *this; } |
---|
446 | |
---|
447 | /// \brief Returns an unsigned random number |
---|
448 | /// |
---|
449 | /// It returns an unsigned integer random number from the range |
---|
450 | /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$ for 64-bit processors and from |
---|
451 | /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$ for 32-bit processors. |
---|
452 | unsigned long getUnsigned() { return 0ul; } |
---|
453 | |
---|
454 | /// \brief Returns a random number |
---|
455 | /// |
---|
456 | /// It returns an integer random number from the range |
---|
457 | /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$ for 64-bit processors and from |
---|
458 | /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$ for 32-bit processors. |
---|
459 | long getInt() { return 0l; } |
---|
460 | |
---|
461 | /// \brief Returns an unsigned integer random number |
---|
462 | /// |
---|
463 | /// It returns an unsigned integer random number from the range |
---|
464 | /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$ for 64-bit processors and |
---|
465 | /// from \f$ \{0, 1, \dots, 2^{31} - 1\} \f$ for 32-bit processors. |
---|
466 | long getNatural() { return 0l; } |
---|
467 | |
---|
468 | /// \brief Returns a random bool |
---|
469 | /// |
---|
470 | /// It returns a random bool. |
---|
471 | bool getBool() { return false; } |
---|
472 | |
---|
473 | /// \brief Returns a real random number |
---|
474 | /// |
---|
475 | /// It returns a real random number from the range |
---|
476 | /// \f$ [0, 1) \f$. For 32-bit processors the generated random |
---|
477 | /// number will just have 32 significant bits. |
---|
478 | double getReal() { return 0.0; } |
---|
479 | |
---|
480 | /// \brief Returns a real random number |
---|
481 | /// |
---|
482 | /// It returns a real random number from the range |
---|
483 | /// \f$ [0, 1) \f$. This function returns random numbers with 53 |
---|
484 | /// significant bits for 32-bit processors. For 64-bit processors |
---|
485 | /// it is identical to the \c getReal(). |
---|
486 | double getPrecReal() { return 0.0; } |
---|
487 | |
---|
488 | /// \brief Returns an unsigned random number from a given range |
---|
489 | /// |
---|
490 | /// It returns an unsigned integer random number from the range |
---|
491 | /// \f$ \{0, 1, \dots, n - 1\} \f$. |
---|
492 | unsigned long getUnsigned(unsigned long n) { return 0; } |
---|
493 | |
---|
494 | /// \brief Returns a random number from a given range |
---|
495 | /// |
---|
496 | /// It returns an unsigned integer random number from the range |
---|
497 | /// \f$ \{0, 1, \dots, n - 1\} \f$. |
---|
498 | long getInt(long n) { return 0; } |
---|
499 | |
---|
500 | }; |
---|
501 | |
---|
502 | #endif |
---|
503 | |
---|
504 | /// \brief Global random number generator instance |
---|
505 | /// |
---|
506 | /// A global mersenne twister random number generator instance |
---|
507 | extern Random rnd; |
---|
508 | |
---|
509 | } |
---|
510 | |
---|
511 | #endif |
---|