Changeset 2242:16523135943d in lemon0.x for lemon/random.h
 Timestamp:
 10/14/06 17:26:05 (14 years ago)
 Branch:
 default
 Phase:
 public
 Convert:
 svn:c9d7d8f590d60310b91f818b3a526b0e/lemon/trunk@2991
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/random.h
r2230 r2242 21 21 22 22 #include <algorithm> 23 #include <iterator> 24 #include <vector> 23 25 24 26 #include <ctime> … … 34 36 namespace lemon { 35 37 36 #if WORD_BIT == 32 38 namespace _random_bits { 39 40 template <typename _Word, int _bits = std::numeric_limits<_Word>::digits> 41 struct RandomTraits {}; 42 43 template <typename _Word> 44 struct RandomTraits<_Word, 32> { 45 46 typedef _Word Word; 47 static const int bits = 32; 48 49 static const int length = 624; 50 static const int shift = 397; 51 52 static const Word mul = 0x6c078965u; 53 static const Word arrayInit = 0x012BD6AAu; 54 static const Word arrayMul1 = 0x0019660Du; 55 static const Word arrayMul2 = 0x5D588B65u; 56 57 static const Word mask = 0x9908B0DFu; 58 static const Word loMask = (1u << 31)  1; 59 static const Word hiMask = ~loMask; 60 61 62 static Word tempering(Word rnd) { 63 rnd ^= (rnd >> 11); 64 rnd ^= (rnd << 7) & 0x9D2C5680u; 65 rnd ^= (rnd << 15) & 0xEFC60000u; 66 rnd ^= (rnd >> 18); 67 return rnd; 68 } 69 70 }; 71 72 template <typename _Word> 73 struct RandomTraits<_Word, 64> { 74 75 typedef _Word Word; 76 static const int bits = 64; 77 78 static const int length = 312; 79 static const int shift = 156; 80 81 static const Word mul = (Word)0x5851F42Du << 32  (Word)0x4C957F2Du; 82 static const Word arrayInit = (Word)0x00000000u << 32 (Word)0x012BD6AAu; 83 static const Word arrayMul1 = (Word)0x369DEA0Fu << 32 (Word)0x31A53F85u; 84 static const Word arrayMul2 = (Word)0x27BB2EE6u << 32 (Word)0x87B0B0FDu; 85 86 static const Word mask = (Word)0xB5026F5Au << 32  (Word)0xA96619E9u; 87 static const Word loMask = ((Word)1u << 31)  1; 88 static const Word hiMask = ~loMask; 89 90 static Word tempering(Word rnd) { 91 rnd ^= (rnd >> 29) & ((Word)0x55555555u << 32  (Word)0x55555555u); 92 rnd ^= (rnd << 17) & ((Word)0x71D67FFFu << 32  (Word)0xEDA60000u); 93 rnd ^= (rnd << 37) & ((Word)0xFFF7EEE0u << 32  (Word)0x00000000u); 94 rnd ^= (rnd >> 43); 95 return rnd; 96 } 97 98 }; 99 100 template <typename _Word> 101 class RandomCore { 102 public: 103 104 typedef _Word Word; 105 106 private: 107 108 static const int bits = RandomTraits<Word>::bits; 109 110 static const int length = RandomTraits<Word>::length; 111 static const int shift = RandomTraits<Word>::shift; 112 113 public: 114 115 void initState() { 116 static const Word seedArray[4] = { 117 0x12345u, 0x23456u, 0x34567u, 0x45678u 118 }; 119 120 initState(seedArray, seedArray + 4); 121 } 122 123 void initState(Word seed) { 124 125 static const Word mul = RandomTraits<Word>::mul; 126 127 current = state; 128 129 Word *curr = state + length  1; 130 curr[0] = seed; curr; 131 for (int i = 1; i < length; ++i) { 132 curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits  2)) ) + i); 133 curr; 134 } 135 } 136 137 template <typename Iterator> 138 void initState(Iterator begin, Iterator end) { 139 140 static const Word init = RandomTraits<Word>::arrayInit; 141 static const Word mul1 = RandomTraits<Word>::arrayMul1; 142 static const Word mul2 = RandomTraits<Word>::arrayMul2; 143 144 145 Word *curr = state + length  1; curr; 146 Iterator it = begin; int cnt = 0; 147 int num; 148 149 initState(init); 150 151 num = length > end  begin ? length : end  begin; 152 while (num) { 153 curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits  2))) * mul1)) 154 + *it + cnt; 155 ++it; ++cnt; 156 if (it == end) { 157 it = begin; cnt = 0; 158 } 159 if (curr == state) { 160 curr = state + length  1; curr[0] = state[0]; 161 } 162 curr; 163 } 164 165 num = length  1; cnt = length  (curr  state)  1; 166 while (num) { 167 curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits  2))) * mul2)) 168  cnt; 169 curr; ++cnt; 170 if (curr == state) { 171 curr = state + length  1; curr[0] = state[0]; curr; 172 cnt = 1; 173 } 174 } 175 176 state[length  1] = (Word)1 << (bits  1); 177 } 178 179 void copyState(const RandomCore& other) { 180 std::copy(other.state, other.state + length, state); 181 current = state + (other.current  other.state); 182 } 183 184 Word operator()() { 185 if (current == state) fillState(); 186 current; 187 Word rnd = *current; 188 return RandomTraits<Word>::tempering(rnd); 189 } 190 191 private: 192 193 194 void fillState() { 195 static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask }; 196 static const Word loMask = RandomTraits<Word>::loMask; 197 static const Word hiMask = RandomTraits<Word>::hiMask; 198 199 current = state + length; 200 201 register Word *curr = state + length  1; 202 register long num; 203 204 num = length  shift; 205 while (num) { 206 curr[0] = (((curr[0] & hiMask)  (curr[1] & loMask)) >> 1) ^ 207 curr[ shift] ^ mask[curr[1] & 1ul]; 208 curr; 209 } 210 num = shift  1; 211 while (num) { 212 curr[0] = (((curr[0] & hiMask)  (curr[1] & loMask)) >> 1) ^ 213 curr[length  shift] ^ mask[curr[1] & 1ul]; 214 curr; 215 } 216 curr[0] = (((curr[0] & hiMask)  (curr[length  1] & loMask)) >> 1) ^ 217 curr[length  shift] ^ mask[curr[length  1] & 1ul]; 218 219 } 220 221 222 Word *current; 223 Word state[length]; 224 225 }; 226 227 228 template <typename Result, 229 int shift = (std::numeric_limits<Result>::digits + 1) / 2> 230 struct Masker { 231 static Result mask(const Result& result) { 232 return Masker<Result, (shift + 1) / 2>:: 233 mask((Result)(result  (result >> shift))); 234 } 235 }; 236 237 template <typename Result> 238 struct Masker<Result, 1> { 239 static Result mask(const Result& result) { 240 return (Result)(result  (result >> 1)); 241 } 242 }; 243 244 template <typename Result, typename Word, 245 int rest = std::numeric_limits<Result>::digits, int shift = 0, 246 bool last = rest <= std::numeric_limits<Word>::digits> 247 struct IntConversion { 248 static const int bits = std::numeric_limits<Word>::digits; 249 250 static Result convert(RandomCore<Word>& rnd) { 251 return (Result)(rnd() >> (bits  rest)) << shift; 252 } 253 254 }; 255 256 template <typename Result, typename Word, int rest, int shift> 257 struct IntConversion<Result, Word, rest, shift, false> { 258 static const int bits = std::numeric_limits<Word>::digits; 259 260 static Result convert(RandomCore<Word>& rnd) { 261 return ((Result)rnd() << shift)  262 IntConversion<Result, Word, rest  bits, shift + bits>::convert(rnd); 263 } 264 }; 265 266 267 template <typename Result, typename Word, 268 bool one_word = std::numeric_limits<Word>::digits < 269 std::numeric_limits<Result>::digits> 270 struct Mapping { 271 static Result map(RandomCore<Word>& rnd, const Result& bound) { 272 Word max = (Word)(bound  1); 273 Result mask = Masker<Result>::mask(bound  1); 274 Result num; 275 do { 276 num = IntConversion<Result, Word>::convert(rnd) & mask; 277 } while (num > max); 278 return num; 279 } 280 }; 281 282 template <typename Result, typename Word> 283 struct Mapping<Result, Word, false> { 284 static Result map(RandomCore<Word>& rnd, const Result& bound) { 285 Word max = (Word)(bound  1); 286 Word mask = Masker<Word, (std::numeric_limits<Result>::digits + 1) / 2> 287 ::mask(max); 288 Word num; 289 do { 290 num = rnd() & mask; 291 } while (num > max); 292 return num; 293 } 294 }; 295 296 template <typename Result, int exp, bool pos = (exp >= 0)> 297 struct ShiftMultiplier { 298 static const Result multiplier() { 299 Result res = ShiftMultiplier<Result, exp / 2>::multiplier(); 300 res *= res; 301 if ((exp & 1) == 1) res *= (Result)2.0; 302 return res; 303 } 304 }; 305 306 template <typename Result, int exp> 307 struct ShiftMultiplier<Result, exp, false> { 308 static const Result multiplier() { 309 Result res = ShiftMultiplier<Result, exp / 2>::multiplier(); 310 res *= res; 311 if ((exp & 1) == 1) res *= (Result)0.5; 312 return res; 313 } 314 }; 315 316 template <typename Result> 317 struct ShiftMultiplier<Result, 0, true> { 318 static const Result multiplier() { 319 return (Result)1.0; 320 } 321 }; 322 323 template <typename Result> 324 struct ShiftMultiplier<Result, 20, true> { 325 static const Result multiplier() { 326 return (Result)(1.0/1048576.0); 327 } 328 }; 329 330 template <typename Result> 331 struct ShiftMultiplier<Result, 32, true> { 332 static const Result multiplier() { 333 return (Result)(1.0/424967296.0); 334 } 335 }; 336 337 template <typename Result> 338 struct ShiftMultiplier<Result, 53, true> { 339 static const Result multiplier() { 340 return (Result)(1.0/9007199254740992.0); 341 } 342 }; 343 344 template <typename Result> 345 struct ShiftMultiplier<Result, 64, true> { 346 static const Result multiplier() { 347 return (Result)(1.0/18446744073709551616.0); 348 } 349 }; 350 351 template <typename Result, int exp> 352 struct Shifting { 353 static Result shift(const Result& result) { 354 return result * ShiftMultiplier<Result, exp>::multiplier(); 355 } 356 }; 357 358 template <typename Result, typename Word, 359 int rest = std::numeric_limits<Result>::digits, int shift = 0, 360 bool last = rest <= std::numeric_limits<Word>::digits> 361 struct RealConversion{ 362 static const int bits = std::numeric_limits<Word>::digits; 363 364 static Result convert(RandomCore<Word>& rnd) { 365 return Shifting<Result,  shift  rest>:: 366 shift((Result)(rnd() >> (bits  rest))); 367 } 368 }; 369 370 template <typename Result, typename Word, int rest, int shift> 371 struct RealConversion<Result, Word, rest, shift, false> { 372 static const int bits = std::numeric_limits<Word>::digits; 373 374 static Result convert(RandomCore<Word>& rnd) { 375 return Shifting<Result,  shift  bits>::shift((Result)rnd()) + 376 RealConversion<Result, Word, restbits, shift + bits>::convert(rnd); 377 } 378 }; 379 380 template <typename Result, typename Word> 381 struct Initializer { 382 383 template <typename Iterator> 384 static void init(RandomCore<Word>& rnd, Iterator begin, Iterator end) { 385 std::vector<Word> ws; 386 for (Iterator it = begin; it != end; ++it) { 387 ws.push_back((Word)*it); 388 } 389 rnd.initState(ws.begin(), ws.end()); 390 } 391 392 template <typename Iterator> 393 static void init(RandomCore<Word>& rnd, Result seed) { 394 rnd.initState(seed); 395 } 396 }; 397 398 template <typename Word> 399 struct BoolConversion { 400 static bool convert(RandomCore<Word>& rnd) { 401 return (rnd() & 1) == 1; 402 } 403 }; 404 405 } 37 406 38 407 /// \ingroup misc … … 41 410 /// 42 411 /// The Mersenne Twister is a twisted generalized feedback 43 /// shiftregister generator of Matsumoto and Nishimura. The period of 44 /// this generator is \f$ 2^{19937}  1 \f$ and it is equidistributed in 45 /// 623 dimensions. The time performance of this generator is comparable 46 /// to the commonly used generators. 412 /// shiftregister generator of Matsumoto and Nishimura. The period 413 /// of this generator is \f$ 2^{19937}  1 \f$ and it is 414 /// equidistributed in 623 dimensions for 32bit numbers. The time 415 /// performance of this generator is comparable to the commonly used 416 /// generators. 417 /// 418 /// This implementation is specialized for both 32bit and 64bit 419 /// architectures. The generators differ sligthly in the 420 /// initialization and generation phase so they produce two 421 /// completly different sequences. 422 /// 423 /// The generator gives back random numbers of serveral types. To 424 /// get a random number from a range of a floating point type you 425 /// can use one form of the \c operator(). If you want to get random 426 /// number from a the {0, 1, ..., n1} integer range use the \c 427 /// operator[]. And to get random number from the whole range of an 428 /// integer type you can use the \c integer() or \c uinteger() 429 /// functions. After all you can get random bool with equal chance 430 /// or given probability with the \c boolean() member function. 431 /// 432 ///\code 433 /// int a = rnd[100000]; // 0..99999 434 /// int b = rnd.uinteger<int>(); // 0 .. 2^31  1 435 /// int c = rnd.integer<int>(); //  2^31 .. 2^31  1 436 /// double d = rnd(); // [0.0, 1.0) 437 /// double e = rnd(100.0); // [0.0, 100.0) 438 /// double f = rnd(1.0, 2.0); // [1.0, 2.0) 439 /// bool g = rnd.boolean(); // P(g = true) = 0.5 440 /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 441 ///\endcode 442 /// 443 /// The lemon provides a global instance of the random number generator 444 /// which name is \c rnd. Usually it is a good programming convenience 445 /// to use this global generator to get random numbers. 47 446 /// 48 447 /// \author Balazs Dezso 49 448 class Random { 50 51 static const int length = 624; 52 static const int shift = 397; 449 private: 450 451 #if WORD_BIT == 32 452 typedef unsigned int Word; 453 #elif WORD_BIT == 64 454 typedef unsigned long Word; 455 #endif 456 457 _random_bits::RandomCore<Word> core; 53 458 54 459 public: 55 460 56 static const unsigned long min = 0ul;57 static const unsigned long max = ~0ul;58 59 461 /// \brief Constructor 60 462 /// 61 /// Constructor with time dependent seeding.62 Random() { initState(std::time(0)); }463 /// Constructor with constant seeding. 464 Random() { core.initState(); } 63 465 64 466 /// \brief Constructor 65 467 /// 66 /// Constructor 67 Random(unsigned long seed) { initState(seed); } 468 /// Constructor with seed. The current number type will be converted 469 /// to the architecture word type. 470 template <typename Number> 471 Random(Number seed) { 472 _random_bits::Initializer<Number, Word>::init(core, seed); 473 } 474 475 /// \brief Constructor 476 /// 477 /// Constructor with array seeding. The given range should contain 478 /// any number type and the numbers will be converted to the 479 /// architecture word type. 480 template <typename Iterator> 481 Random(Iterator begin, Iterator end) { 482 typedef typename std::iterator_traits<Iterator>::value_type Number; 483 _random_bits::Initializer<Number, Word>::init(core, begin, end); 484 } 68 485 69 486 /// \brief Copy constructor … … 71 488 /// Copy constructor. The generated sequence will be identical to 72 489 /// the other sequence. 73 Random(const Random& other) { 74 std::copy(other.state, other.state + length, state); 75 current = state + (other.current  other.state); 490 Random(const Random& other) { 491 core.copyState(other.core); 76 492 } 77 493 … … 82 498 Random& operator=(const Random& other) { 83 499 if (&other != this) { 84 std::copy(other.state, other.state + length, state); 85 current = state + (other.current  other.state); 500 core.copyState(other.core); 86 501 } 87 502 return *this; 88 503 } 89 504 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 505 /// \brief Returns a random real number 506 /// 507 /// It returns a random real number from the range [0, 1). The default 508 /// result type of this function is double. 509 template <typename Number> 510 Number operator()() { 511 return _random_bits::RealConversion<Number, Word>::convert(core); 512 } 513 514 double operator()() { 515 return operator()<double>(); 516 } 517 518 /// \brief Returns a random real number 519 /// 520 /// It returns a random real number from the range [0, b). 521 template <typename Number> 522 Number operator()(Number b) { 523 return operator()<Number>() * b; 524 } 525 526 /// \brief Returns a random real number 527 /// 528 /// It returns a random real number from the range [a, b). 529 template <typename Number> 530 Number operator()(Number a, Number b) { 531 return operator()<Number>() * (b  a) + a; 532 } 533 534 /// \brief Returns a random integer from a range 535 /// 536 /// It returns a random integer from the range {0, 1, ..., bound  1}. 537 template <typename Number> 538 Number operator[](const Number& bound) { 539 return _random_bits::Mapping<Number, Word>::map(core, bound); 540 } 541 542 /// \brief Returns a random nonnegative integer 543 /// 544 /// It returns a random nonnegative integer uniformly from the 545 /// whole range of the current \c Number type. The default result 546 /// type of this function is unsigned int. 547 template <typename Number> 548 Number uinteger() { 549 return _random_bits::IntConversion<Number, Word>::convert(core); 550 } 551 552 unsigned int uinteger() { 553 return uinteger<unsigned int>(); 554 } 555 556 /// \brief Returns a random integer 557 /// 558 /// It returns a random integer uniformly from the whole range of 559 /// the current \c Number type. The default result type of this 560 /// function is int. 561 template <typename Number> 562 Number integer() { 563 static const int nb = std::numeric_limits<Number>::digits + 564 (std::numeric_limits<Number>::is_signed ? 1 : 0); 565 return _random_bits::IntConversion<Number, Word, nb>::convert(core); 566 } 567 568 int integer() { 569 return integer<int>(); 570 } 571 121 572 /// \brief Returns a random bool 122 573 /// 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]; 574 /// It returns a random bool 575 bool boolean() { 576 return _random_bits::BoolConversion<Word>::convert(core); 577 } 578 579 /// \brief Returns a random bool 580 /// 581 /// It returns a random bool with given probability of true result 582 bool boolean(double p) { 583 return operator()() < p; 584 } 220 585 221 586 }; 222 587 223 #elif WORD_BIT == 64224 225 /// \ingroup misc226 ///227 /// \brief Mersenne Twister random number generator228 ///229 /// The Mersenne Twister is a twisted generalized feedback230 /// shiftregister generator of Matsumoto and Nishimura. The period of231 /// this generator is \f$ 2^{19937}  1 \f$ and it is equidistributed in232 /// 311 dimensions. The time performance of this generator is comparable233 /// to the commonly used generators.234 class Random {235 236 static const int length = 312;237 static const int shift = 156;238 239 public:240 241 static const unsigned long min = 0ul;242 static const unsigned long max = ~0ul;243 244 /// \brief Constructor245 ///246 /// Constructor with time dependent seeding.247 Random() { initState(std::time(0)); }248 249 /// \brief Constructor250 ///251 /// Constructor252 Random(unsigned long seed) { initState(seed); }253 254 /// \brief Copy constructor255 ///256 /// Copy constructor. The generated sequence will be identical to257 /// the other sequence.258 Random(const Random& other) {259 std::copy(other.state, other.state + length, state);260 current = state + (other.current  other.state);261 }262 263 /// \brief Assign operator264 ///265 /// Assign operator. The generated sequence will be identical to266 /// the other sequence.267 Random& operator=(const Random& other) {268 if (&other != this) {269 std::copy(other.state, other.state + length, state);270 current = state + (other.current  other.state);271 }272 return *this;273 }274 275 /// \brief Returns an unsigned random number276 ///277 /// It returns an unsigned integer random number from the range278 /// \f$ \{0, 1, \dots, 2^{64}  1\} \f$.279 unsigned long getUnsigned() {280 if (current == state) fillState();281 current;282 unsigned long rnd = *current;283 rnd ^= (rnd >> 29) & 0x5555555555555555ul;284 rnd ^= (rnd << 17) & 0x71D67FFFEDA60000ul;285 rnd ^= (rnd << 37) & 0xFFF7EEE000000000ul;286 rnd ^= (rnd >> 43);287 return rnd;288 }289 290 /// \brief Returns a random number291 ///292 /// It returns an integer random number from the range293 /// \f$ \{2^{63}, \dots, 2^{63}  1\} \f$.294 long getInt() {295 return (long)getUnsigned();296 }297 298 /// \brief Returns an unsigned integer random number299 ///300 /// It returns an unsigned integer random number from the range301 /// \f$ \{0, 1, \dots, 2^{63}  1\} \f$.302 long getNatural() {303 return (long)(getUnsigned() >> 1);304 }305 306 /// \brief Returns a random bool307 ///308 /// It returns a random bool.309 bool getBool() {310 return (bool)(getUnsigned() & 1);311 }312 313 /// \brief Returns a real random number314 ///315 /// It returns a real random number from the range316 /// \f$ [0, 1) \f$.317 double getReal() {318 return std::ldexp((double)(getUnsigned() >> 11), 53);319 }320 321 /// \brief Returns a real random number322 ///323 /// It returns a real random number from the range324 /// \f$ [0, 1) \f$. This function is identical to the \c getReal().325 double getPrecReal() {326 return getReal();327 }328 329 /// \brief Returns an unsigned random number from a given range330 ///331 /// It returns an unsigned integer random number from the range332 /// \f$ \{0, 1, \dots, n  1\} \f$.333 unsigned long getUnsigned(unsigned long n) {334 unsigned long mask = n  1, rnd;335 mask = mask >> 1;336 mask = mask >> 2;337 mask = mask >> 4;338 mask = mask >> 8;339 mask = mask >> 16;340 mask = mask >> 32;341 do rnd = getUnsigned() & mask; while (rnd >= n);342 return rnd;343 }344 345 /// \brief Returns a random number from a given range346 ///347 /// It returns an unsigned integer random number from the range348 /// \f$ \{0, 1, \dots, n  1\} \f$.349 long getInt(long n) {350 long mask = n  1, rnd;351 mask = mask >> 1;352 mask = mask >> 2;353 mask = mask >> 4;354 mask = mask >> 8;355 mask = mask >> 16;356 mask = mask >> 32;357 do rnd = getUnsigned() & mask; while (rnd >= n);358 return rnd;359 }360 361 private:362 363 void initState(unsigned long seed) {364 365 static const unsigned long mul = 0x5851F42D4C957F2Dul;366 367 current = state;368 369 unsigned long *curr = state + length  1;370 curr[0] = seed; curr;371 for (int i = 1; i < length; ++i) {372 curr[0] = (mul * ( curr[1] ^ (curr[1] >> 62) ) + i);373 curr;374 }375 }376 377 void fillState() {378 static const unsigned long mask[2] = { 0x0ul, 0xB5026F5AA96619E9ul };379 static const unsigned long loMask = (1ul << 31)  1;380 static const unsigned long hiMask = ~loMask;381 382 current = state + length;383 384 register unsigned long *curr = state + length  1;385 register int num;386 387 num = length  shift;388 while (num) {389 curr[0] = (((curr[0] & hiMask)  (curr[1] & loMask)) >> 1) ^390 curr[ shift] ^ mask[curr[1] & 1ul];391 curr;392 }393 num = shift  1;394 while (num) {395 curr[0] = (((curr[0] & hiMask)  (curr[1] & loMask)) >> 1) ^396 curr[length  shift] ^ mask[curr[1] & 1ul];397 curr;398 }399 curr[0] = (((curr[0] & hiMask)  (curr[length  1] & loMask)) >> 1) ^400 curr[length  shift] ^ mask[curr[length  1] & 1ul];401 402 }403 404 unsigned long *current;405 unsigned long state[length];406 407 };408 409 #else410 411 /// \ingroup misc412 ///413 /// \brief Mersenne Twister random number generator414 ///415 /// The Mersenne Twister is a twisted generalized feedback416 /// shiftregister generator of Matsumoto and Nishimura. There is417 /// two different implementation in the Lemon library, one for418 /// 32bit processors and one for 64bit processors. The period of419 /// the generated sequence is \f$ 2^{19937}  1 \f$, the generated420 /// sequence of 32bit random numbers is equidistributed in 623421 /// dimensions. The time performance of this generator is comparable422 /// to the commonly used generators.423 class Random {424 public:425 426 static const unsigned long min = 0ul;427 static const unsigned long max = ~0ul;428 429 /// \brief Constructor430 ///431 /// Constructor with time dependent seeding.432 Random() {}433 434 /// \brief Constructor435 ///436 /// Constructor437 Random(unsigned long seed) {}438 439 /// \brief Copy constructor440 ///441 /// Copy constructor. The generated sequence will be identical to442 /// the other sequence.443 Random(const Random& other) {}444 445 /// \brief Assign operator446 ///447 /// Assign operator. The generated sequence will be identical to448 /// the other sequence.449 Random& operator=(const Random& other) { return *this; }450 451 /// \brief Returns an unsigned random number452 ///453 /// It returns an unsigned integer random number from the range454 /// \f$ \{0, 1, \dots, 2^{64}  1\} \f$ for 64bit processors and from455 /// \f$ \{0, 1, \dots, 2^{32}  1\} \f$ for 32bit processors.456 unsigned long getUnsigned() { return 0ul; }457 458 /// \brief Returns a random number459 ///460 /// It returns an integer random number from the range461 /// \f$ \{2^{63}, \dots, 2^{63}  1\} \f$ for 64bit processors and from462 /// \f$ \{2^{31}, \dots, 2^{31}  1\} \f$ for 32bit processors.463 long getInt() { return 0l; }464 465 /// \brief Returns an unsigned integer random number466 ///467 /// It returns an unsigned integer random number from the range468 /// \f$ \{0, 1, \dots, 2^{63}  1\} \f$ for 64bit processors and469 /// from \f$ \{0, 1, \dots, 2^{31}  1\} \f$ for 32bit processors.470 long getNatural() { return 0l; }471 472 /// \brief Returns a random bool473 ///474 /// It returns a random bool.475 bool getBool() { return false; }476 477 /// \brief Returns a real random number478 ///479 /// It returns a real random number from the range480 /// \f$ [0, 1) \f$. For 32bit processors the generated random481 /// number will just have 32 significant bits.482 double getReal() { return 0.0; }483 484 /// \brief Returns a real random number485 ///486 /// It returns a real random number from the range487 /// \f$ [0, 1) \f$. This function returns random numbers with 53488 /// significant bits for 32bit processors. For 64bit processors489 /// it is identical to the \c getReal().490 double getPrecReal() { return 0.0; }491 492 /// \brief Returns an unsigned random number from a given range493 ///494 /// It returns an unsigned integer random number from the range495 /// \f$ \{0, 1, \dots, n  1\} \f$.496 unsigned long getUnsigned(unsigned long n) { return 0; }497 498 /// \brief Returns a random number from a given range499 ///500 /// It returns an unsigned integer random number from the range501 /// \f$ \{0, 1, \dots, n  1\} \f$.502 long getInt(long n) { return 0; }503 504 };505 506 #endif507 588 508 589 /// \brief Global random number generator instance
Note: See TracChangeset
for help on using the changeset viewer.