1.1 --- a/lemon/random.h Fri Oct 13 15:10:50 2006 +0000
1.2 +++ b/lemon/random.h Sat Oct 14 15:26:05 2006 +0000
1.3 @@ -20,6 +20,8 @@
1.4 #define LEMON_RANDOM_H
1.5
1.6 #include <algorithm>
1.7 +#include <iterator>
1.8 +#include <vector>
1.9
1.10 #include <ctime>
1.11 #include <cmath>
1.12 @@ -33,46 +35,460 @@
1.13
1.14 namespace lemon {
1.15
1.16 -#if WORD_BIT == 32
1.17 + namespace _random_bits {
1.18 +
1.19 + template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
1.20 + struct RandomTraits {};
1.21 +
1.22 + template <typename _Word>
1.23 + struct RandomTraits<_Word, 32> {
1.24 +
1.25 + typedef _Word Word;
1.26 + static const int bits = 32;
1.27 +
1.28 + static const int length = 624;
1.29 + static const int shift = 397;
1.30 +
1.31 + static const Word mul = 0x6c078965u;
1.32 + static const Word arrayInit = 0x012BD6AAu;
1.33 + static const Word arrayMul1 = 0x0019660Du;
1.34 + static const Word arrayMul2 = 0x5D588B65u;
1.35 +
1.36 + static const Word mask = 0x9908B0DFu;
1.37 + static const Word loMask = (1u << 31) - 1;
1.38 + static const Word hiMask = ~loMask;
1.39 +
1.40 +
1.41 + static Word tempering(Word rnd) {
1.42 + rnd ^= (rnd >> 11);
1.43 + rnd ^= (rnd << 7) & 0x9D2C5680u;
1.44 + rnd ^= (rnd << 15) & 0xEFC60000u;
1.45 + rnd ^= (rnd >> 18);
1.46 + return rnd;
1.47 + }
1.48 +
1.49 + };
1.50 +
1.51 + template <typename _Word>
1.52 + struct RandomTraits<_Word, 64> {
1.53 +
1.54 + typedef _Word Word;
1.55 + static const int bits = 64;
1.56 +
1.57 + static const int length = 312;
1.58 + static const int shift = 156;
1.59 +
1.60 + static const Word mul = (Word)0x5851F42Du << 32 | (Word)0x4C957F2Du;
1.61 + static const Word arrayInit = (Word)0x00000000u << 32 |(Word)0x012BD6AAu;
1.62 + static const Word arrayMul1 = (Word)0x369DEA0Fu << 32 |(Word)0x31A53F85u;
1.63 + static const Word arrayMul2 = (Word)0x27BB2EE6u << 32 |(Word)0x87B0B0FDu;
1.64 +
1.65 + static const Word mask = (Word)0xB5026F5Au << 32 | (Word)0xA96619E9u;
1.66 + static const Word loMask = ((Word)1u << 31) - 1;
1.67 + static const Word hiMask = ~loMask;
1.68 +
1.69 + static Word tempering(Word rnd) {
1.70 + rnd ^= (rnd >> 29) & ((Word)0x55555555u << 32 | (Word)0x55555555u);
1.71 + rnd ^= (rnd << 17) & ((Word)0x71D67FFFu << 32 | (Word)0xEDA60000u);
1.72 + rnd ^= (rnd << 37) & ((Word)0xFFF7EEE0u << 32 | (Word)0x00000000u);
1.73 + rnd ^= (rnd >> 43);
1.74 + return rnd;
1.75 + }
1.76 +
1.77 + };
1.78 +
1.79 + template <typename _Word>
1.80 + class RandomCore {
1.81 + public:
1.82 +
1.83 + typedef _Word Word;
1.84 +
1.85 + private:
1.86 +
1.87 + static const int bits = RandomTraits<Word>::bits;
1.88 +
1.89 + static const int length = RandomTraits<Word>::length;
1.90 + static const int shift = RandomTraits<Word>::shift;
1.91 +
1.92 + public:
1.93 +
1.94 + void initState() {
1.95 + static const Word seedArray[4] = {
1.96 + 0x12345u, 0x23456u, 0x34567u, 0x45678u
1.97 + };
1.98 +
1.99 + initState(seedArray, seedArray + 4);
1.100 + }
1.101 +
1.102 + void initState(Word seed) {
1.103 +
1.104 + static const Word mul = RandomTraits<Word>::mul;
1.105 +
1.106 + current = state;
1.107 +
1.108 + Word *curr = state + length - 1;
1.109 + curr[0] = seed; --curr;
1.110 + for (int i = 1; i < length; ++i) {
1.111 + curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i);
1.112 + --curr;
1.113 + }
1.114 + }
1.115 +
1.116 + template <typename Iterator>
1.117 + void initState(Iterator begin, Iterator end) {
1.118 +
1.119 + static const Word init = RandomTraits<Word>::arrayInit;
1.120 + static const Word mul1 = RandomTraits<Word>::arrayMul1;
1.121 + static const Word mul2 = RandomTraits<Word>::arrayMul2;
1.122 +
1.123 +
1.124 + Word *curr = state + length - 1; --curr;
1.125 + Iterator it = begin; int cnt = 0;
1.126 + int num;
1.127 +
1.128 + initState(init);
1.129 +
1.130 + num = length > end - begin ? length : end - begin;
1.131 + while (num--) {
1.132 + curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1))
1.133 + + *it + cnt;
1.134 + ++it; ++cnt;
1.135 + if (it == end) {
1.136 + it = begin; cnt = 0;
1.137 + }
1.138 + if (curr == state) {
1.139 + curr = state + length - 1; curr[0] = state[0];
1.140 + }
1.141 + --curr;
1.142 + }
1.143 +
1.144 + num = length - 1; cnt = length - (curr - state) - 1;
1.145 + while (num--) {
1.146 + curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2))
1.147 + - cnt;
1.148 + --curr; ++cnt;
1.149 + if (curr == state) {
1.150 + curr = state + length - 1; curr[0] = state[0]; --curr;
1.151 + cnt = 1;
1.152 + }
1.153 + }
1.154 +
1.155 + state[length - 1] = (Word)1 << (bits - 1);
1.156 + }
1.157 +
1.158 + void copyState(const RandomCore& other) {
1.159 + std::copy(other.state, other.state + length, state);
1.160 + current = state + (other.current - other.state);
1.161 + }
1.162 +
1.163 + Word operator()() {
1.164 + if (current == state) fillState();
1.165 + --current;
1.166 + Word rnd = *current;
1.167 + return RandomTraits<Word>::tempering(rnd);
1.168 + }
1.169 +
1.170 + private:
1.171 +
1.172 +
1.173 + void fillState() {
1.174 + static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
1.175 + static const Word loMask = RandomTraits<Word>::loMask;
1.176 + static const Word hiMask = RandomTraits<Word>::hiMask;
1.177 +
1.178 + current = state + length;
1.179 +
1.180 + register Word *curr = state + length - 1;
1.181 + register long num;
1.182 +
1.183 + num = length - shift;
1.184 + while (num--) {
1.185 + curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
1.186 + curr[- shift] ^ mask[curr[-1] & 1ul];
1.187 + --curr;
1.188 + }
1.189 + num = shift - 1;
1.190 + while (num--) {
1.191 + curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
1.192 + curr[length - shift] ^ mask[curr[-1] & 1ul];
1.193 + --curr;
1.194 + }
1.195 + curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
1.196 + curr[length - shift] ^ mask[curr[length - 1] & 1ul];
1.197 +
1.198 + }
1.199 +
1.200 +
1.201 + Word *current;
1.202 + Word state[length];
1.203 +
1.204 + };
1.205 +
1.206 +
1.207 + template <typename Result,
1.208 + int shift = (std::numeric_limits<Result>::digits + 1) / 2>
1.209 + struct Masker {
1.210 + static Result mask(const Result& result) {
1.211 + return Masker<Result, (shift + 1) / 2>::
1.212 + mask((Result)(result | (result >> shift)));
1.213 + }
1.214 + };
1.215 +
1.216 + template <typename Result>
1.217 + struct Masker<Result, 1> {
1.218 + static Result mask(const Result& result) {
1.219 + return (Result)(result | (result >> 1));
1.220 + }
1.221 + };
1.222 +
1.223 + template <typename Result, typename Word,
1.224 + int rest = std::numeric_limits<Result>::digits, int shift = 0,
1.225 + bool last = rest <= std::numeric_limits<Word>::digits>
1.226 + struct IntConversion {
1.227 + static const int bits = std::numeric_limits<Word>::digits;
1.228 +
1.229 + static Result convert(RandomCore<Word>& rnd) {
1.230 + return (Result)(rnd() >> (bits - rest)) << shift;
1.231 + }
1.232 +
1.233 + };
1.234 +
1.235 + template <typename Result, typename Word, int rest, int shift>
1.236 + struct IntConversion<Result, Word, rest, shift, false> {
1.237 + static const int bits = std::numeric_limits<Word>::digits;
1.238 +
1.239 + static Result convert(RandomCore<Word>& rnd) {
1.240 + return ((Result)rnd() << shift) |
1.241 + IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);
1.242 + }
1.243 + };
1.244 +
1.245 +
1.246 + template <typename Result, typename Word,
1.247 + bool one_word = std::numeric_limits<Word>::digits <
1.248 + std::numeric_limits<Result>::digits>
1.249 + struct Mapping {
1.250 + static Result map(RandomCore<Word>& rnd, const Result& bound) {
1.251 + Word max = (Word)(bound - 1);
1.252 + Result mask = Masker<Result>::mask(bound - 1);
1.253 + Result num;
1.254 + do {
1.255 + num = IntConversion<Result, Word>::convert(rnd) & mask;
1.256 + } while (num > max);
1.257 + return num;
1.258 + }
1.259 + };
1.260 +
1.261 + template <typename Result, typename Word>
1.262 + struct Mapping<Result, Word, false> {
1.263 + static Result map(RandomCore<Word>& rnd, const Result& bound) {
1.264 + Word max = (Word)(bound - 1);
1.265 + Word mask = Masker<Word, (std::numeric_limits<Result>::digits + 1) / 2>
1.266 + ::mask(max);
1.267 + Word num;
1.268 + do {
1.269 + num = rnd() & mask;
1.270 + } while (num > max);
1.271 + return num;
1.272 + }
1.273 + };
1.274 +
1.275 + template <typename Result, int exp, bool pos = (exp >= 0)>
1.276 + struct ShiftMultiplier {
1.277 + static const Result multiplier() {
1.278 + Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
1.279 + res *= res;
1.280 + if ((exp & 1) == 1) res *= (Result)2.0;
1.281 + return res;
1.282 + }
1.283 + };
1.284 +
1.285 + template <typename Result, int exp>
1.286 + struct ShiftMultiplier<Result, exp, false> {
1.287 + static const Result multiplier() {
1.288 + Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
1.289 + res *= res;
1.290 + if ((exp & 1) == 1) res *= (Result)0.5;
1.291 + return res;
1.292 + }
1.293 + };
1.294 +
1.295 + template <typename Result>
1.296 + struct ShiftMultiplier<Result, 0, true> {
1.297 + static const Result multiplier() {
1.298 + return (Result)1.0;
1.299 + }
1.300 + };
1.301 +
1.302 + template <typename Result>
1.303 + struct ShiftMultiplier<Result, -20, true> {
1.304 + static const Result multiplier() {
1.305 + return (Result)(1.0/1048576.0);
1.306 + }
1.307 + };
1.308 +
1.309 + template <typename Result>
1.310 + struct ShiftMultiplier<Result, -32, true> {
1.311 + static const Result multiplier() {
1.312 + return (Result)(1.0/424967296.0);
1.313 + }
1.314 + };
1.315 +
1.316 + template <typename Result>
1.317 + struct ShiftMultiplier<Result, -53, true> {
1.318 + static const Result multiplier() {
1.319 + return (Result)(1.0/9007199254740992.0);
1.320 + }
1.321 + };
1.322 +
1.323 + template <typename Result>
1.324 + struct ShiftMultiplier<Result, -64, true> {
1.325 + static const Result multiplier() {
1.326 + return (Result)(1.0/18446744073709551616.0);
1.327 + }
1.328 + };
1.329 +
1.330 + template <typename Result, int exp>
1.331 + struct Shifting {
1.332 + static Result shift(const Result& result) {
1.333 + return result * ShiftMultiplier<Result, exp>::multiplier();
1.334 + }
1.335 + };
1.336 +
1.337 + template <typename Result, typename Word,
1.338 + int rest = std::numeric_limits<Result>::digits, int shift = 0,
1.339 + bool last = rest <= std::numeric_limits<Word>::digits>
1.340 + struct RealConversion{
1.341 + static const int bits = std::numeric_limits<Word>::digits;
1.342 +
1.343 + static Result convert(RandomCore<Word>& rnd) {
1.344 + return Shifting<Result, - shift - rest>::
1.345 + shift((Result)(rnd() >> (bits - rest)));
1.346 + }
1.347 + };
1.348 +
1.349 + template <typename Result, typename Word, int rest, int shift>
1.350 + struct RealConversion<Result, Word, rest, shift, false> {
1.351 + static const int bits = std::numeric_limits<Word>::digits;
1.352 +
1.353 + static Result convert(RandomCore<Word>& rnd) {
1.354 + return Shifting<Result, - shift - bits>::shift((Result)rnd()) +
1.355 + RealConversion<Result, Word, rest-bits, shift + bits>::convert(rnd);
1.356 + }
1.357 + };
1.358 +
1.359 + template <typename Result, typename Word>
1.360 + struct Initializer {
1.361 +
1.362 + template <typename Iterator>
1.363 + static void init(RandomCore<Word>& rnd, Iterator begin, Iterator end) {
1.364 + std::vector<Word> ws;
1.365 + for (Iterator it = begin; it != end; ++it) {
1.366 + ws.push_back((Word)*it);
1.367 + }
1.368 + rnd.initState(ws.begin(), ws.end());
1.369 + }
1.370 +
1.371 + template <typename Iterator>
1.372 + static void init(RandomCore<Word>& rnd, Result seed) {
1.373 + rnd.initState(seed);
1.374 + }
1.375 + };
1.376 +
1.377 + template <typename Word>
1.378 + struct BoolConversion {
1.379 + static bool convert(RandomCore<Word>& rnd) {
1.380 + return (rnd() & 1) == 1;
1.381 + }
1.382 + };
1.383 +
1.384 + }
1.385
1.386 /// \ingroup misc
1.387 ///
1.388 /// \brief Mersenne Twister random number generator
1.389 ///
1.390 /// The Mersenne Twister is a twisted generalized feedback
1.391 - /// shift-register generator of Matsumoto and Nishimura. The period of
1.392 - /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in
1.393 - /// 623 dimensions. The time performance of this generator is comparable
1.394 - /// to the commonly used generators.
1.395 + /// shift-register generator of Matsumoto and Nishimura. The period
1.396 + /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
1.397 + /// equi-distributed in 623 dimensions for 32-bit numbers. The time
1.398 + /// performance of this generator is comparable to the commonly used
1.399 + /// generators.
1.400 + ///
1.401 + /// This implementation is specialized for both 32-bit and 64-bit
1.402 + /// architectures. The generators differ sligthly in the
1.403 + /// initialization and generation phase so they produce two
1.404 + /// completly different sequences.
1.405 + ///
1.406 + /// The generator gives back random numbers of serveral types. To
1.407 + /// get a random number from a range of a floating point type you
1.408 + /// can use one form of the \c operator(). If you want to get random
1.409 + /// number from a the {0, 1, ..., n-1} integer range use the \c
1.410 + /// operator[]. And to get random number from the whole range of an
1.411 + /// integer type you can use the \c integer() or \c uinteger()
1.412 + /// functions. After all you can get random bool with equal chance
1.413 + /// or given probability with the \c boolean() member function.
1.414 + ///
1.415 + ///\code
1.416 + /// int a = rnd[100000]; // 0..99999
1.417 + /// int b = rnd.uinteger<int>(); // 0 .. 2^31 - 1
1.418 + /// int c = rnd.integer<int>(); // - 2^31 .. 2^31 - 1
1.419 + /// double d = rnd(); // [0.0, 1.0)
1.420 + /// double e = rnd(100.0); // [0.0, 100.0)
1.421 + /// double f = rnd(1.0, 2.0); // [1.0, 2.0)
1.422 + /// bool g = rnd.boolean(); // P(g = true) = 0.5
1.423 + /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8
1.424 + ///\endcode
1.425 + ///
1.426 + /// The lemon provides a global instance of the random number generator
1.427 + /// which name is \c rnd. Usually it is a good programming convenience
1.428 + /// to use this global generator to get random numbers.
1.429 ///
1.430 /// \author Balazs Dezso
1.431 class Random {
1.432 + private:
1.433
1.434 - static const int length = 624;
1.435 - static const int shift = 397;
1.436 +#if WORD_BIT == 32
1.437 + typedef unsigned int Word;
1.438 +#elif WORD_BIT == 64
1.439 + typedef unsigned long Word;
1.440 +#endif
1.441 +
1.442 + _random_bits::RandomCore<Word> core;
1.443
1.444 public:
1.445
1.446 - static const unsigned long min = 0ul;
1.447 - static const unsigned long max = ~0ul;
1.448 -
1.449 /// \brief Constructor
1.450 ///
1.451 - /// Constructor with time dependent seeding.
1.452 - Random() { initState(std::time(0)); }
1.453 + /// Constructor with constant seeding.
1.454 + Random() { core.initState(); }
1.455
1.456 /// \brief Constructor
1.457 ///
1.458 - /// Constructor
1.459 - Random(unsigned long seed) { initState(seed); }
1.460 + /// Constructor with seed. The current number type will be converted
1.461 + /// to the architecture word type.
1.462 + template <typename Number>
1.463 + Random(Number seed) {
1.464 + _random_bits::Initializer<Number, Word>::init(core, seed);
1.465 + }
1.466 +
1.467 + /// \brief Constructor
1.468 + ///
1.469 + /// Constructor with array seeding. The given range should contain
1.470 + /// any number type and the numbers will be converted to the
1.471 + /// architecture word type.
1.472 + template <typename Iterator>
1.473 + Random(Iterator begin, Iterator end) {
1.474 + typedef typename std::iterator_traits<Iterator>::value_type Number;
1.475 + _random_bits::Initializer<Number, Word>::init(core, begin, end);
1.476 + }
1.477
1.478 /// \brief Copy constructor
1.479 ///
1.480 /// Copy constructor. The generated sequence will be identical to
1.481 /// the other sequence.
1.482 - Random(const Random& other) {
1.483 - std::copy(other.state, other.state + length, state);
1.484 - current = state + (other.current - other.state);
1.485 + Random(const Random& other) {
1.486 + core.copyState(other.core);
1.487 }
1.488
1.489 /// \brief Assign operator
1.490 @@ -81,429 +497,94 @@
1.491 /// the other sequence.
1.492 Random& operator=(const Random& other) {
1.493 if (&other != this) {
1.494 - std::copy(other.state, other.state + length, state);
1.495 - current = state + (other.current - other.state);
1.496 + core.copyState(other.core);
1.497 }
1.498 return *this;
1.499 }
1.500
1.501 - /// \brief Returns an unsigned random number
1.502 + /// \brief Returns a random real number
1.503 ///
1.504 - /// It returns an unsigned integer random number from the range
1.505 - /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$.
1.506 - unsigned long getUnsigned() {
1.507 - if (current == state) fillState();
1.508 - --current;
1.509 - unsigned long rnd = *current;
1.510 - rnd ^= (rnd >> 11);
1.511 - rnd ^= (rnd << 7) & 0x9D2C5680ul;
1.512 - rnd ^= (rnd << 15) & 0xEFC60000ul;
1.513 - rnd ^= (rnd >> 18);
1.514 - return rnd;
1.515 + /// It returns a random real number from the range [0, 1). The default
1.516 + /// result type of this function is double.
1.517 + template <typename Number>
1.518 + Number operator()() {
1.519 + return _random_bits::RealConversion<Number, Word>::convert(core);
1.520 }
1.521
1.522 - /// \brief Returns a random number
1.523 + double operator()() {
1.524 + return operator()<double>();
1.525 + }
1.526 +
1.527 + /// \brief Returns a random real number
1.528 ///
1.529 - /// It returns an integer random number from the range
1.530 - /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$.
1.531 - long getInt() {
1.532 - return (long)getUnsigned();
1.533 + /// It returns a random real number from the range [0, b).
1.534 + template <typename Number>
1.535 + Number operator()(Number b) {
1.536 + return operator()<Number>() * b;
1.537 + }
1.538 +
1.539 + /// \brief Returns a random real number
1.540 + ///
1.541 + /// It returns a random real number from the range [a, b).
1.542 + template <typename Number>
1.543 + Number operator()(Number a, Number b) {
1.544 + return operator()<Number>() * (b - a) + a;
1.545 + }
1.546 +
1.547 + /// \brief Returns a random integer from a range
1.548 + ///
1.549 + /// It returns a random integer from the range {0, 1, ..., bound - 1}.
1.550 + template <typename Number>
1.551 + Number operator[](const Number& bound) {
1.552 + return _random_bits::Mapping<Number, Word>::map(core, bound);
1.553 + }
1.554 +
1.555 + /// \brief Returns a random non-negative integer
1.556 + ///
1.557 + /// It returns a random non-negative integer uniformly from the
1.558 + /// whole range of the current \c Number type. The default result
1.559 + /// type of this function is unsigned int.
1.560 + template <typename Number>
1.561 + Number uinteger() {
1.562 + return _random_bits::IntConversion<Number, Word>::convert(core);
1.563 + }
1.564 +
1.565 + unsigned int uinteger() {
1.566 + return uinteger<unsigned int>();
1.567 + }
1.568 +
1.569 + /// \brief Returns a random integer
1.570 + ///
1.571 + /// It returns a random integer uniformly from the whole range of
1.572 + /// the current \c Number type. The default result type of this
1.573 + /// function is int.
1.574 + template <typename Number>
1.575 + Number integer() {
1.576 + static const int nb = std::numeric_limits<Number>::digits +
1.577 + (std::numeric_limits<Number>::is_signed ? 1 : 0);
1.578 + return _random_bits::IntConversion<Number, Word, nb>::convert(core);
1.579 + }
1.580 +
1.581 + int integer() {
1.582 + return integer<int>();
1.583 }
1.584
1.585 - /// \brief Returns an unsigned integer random number
1.586 + /// \brief Returns a random bool
1.587 ///
1.588 - /// It returns an unsigned integer random number from the range
1.589 - /// \f$ \{0, 1, \dots, 2^{31} - 1\} \f$.
1.590 - long getNatural() {
1.591 - return (long)(getUnsigned() >> 1);
1.592 + /// It returns a random bool
1.593 + bool boolean() {
1.594 + return _random_bits::BoolConversion<Word>::convert(core);
1.595 }
1.596
1.597 /// \brief Returns a random bool
1.598 ///
1.599 - /// It returns a random bool.
1.600 - bool getBool() {
1.601 - return (bool)(getUnsigned() & 1);
1.602 + /// It returns a random bool with given probability of true result
1.603 + bool boolean(double p) {
1.604 + return operator()() < p;
1.605 }
1.606 -
1.607 - /// \brief Returns a real random number
1.608 - ///
1.609 - /// It returns a real random number from the range
1.610 - /// \f$ [0, 1) \f$. The double will have 32 significant bits.
1.611 - double getReal() {
1.612 - return std::ldexp((double)getUnsigned(), -32);
1.613 - }
1.614 -
1.615 - /// \brief Returns a real random number
1.616 - ///
1.617 - /// It returns a real random number from the range
1.618 - /// \f$ [0, 1) \f$. The double will have 53 significant bits,
1.619 - /// but it is slower than the \c getReal().
1.620 - double getPrecReal() {
1.621 - return std::ldexp((double)(getUnsigned() >> 5), -27) +
1.622 - std::ldexp((double)(getUnsigned() >> 6), -53);
1.623 - }
1.624 -
1.625 - /// \brief Returns an unsigned random number from a given range
1.626 - ///
1.627 - /// It returns an unsigned integer random number from the range
1.628 - /// \f$ \{0, 1, \dots, n - 1\} \f$.
1.629 - unsigned long getUnsigned(unsigned long n) {
1.630 - unsigned long mask = n - 1, rnd;
1.631 - mask |= mask >> 1;
1.632 - mask |= mask >> 2;
1.633 - mask |= mask >> 4;
1.634 - mask |= mask >> 8;
1.635 - mask |= mask >> 16;
1.636 - do rnd = getUnsigned() & mask; while (rnd >= n);
1.637 - return rnd;
1.638 - }
1.639 -
1.640 - /// \brief Returns a random number from a given range
1.641 - ///
1.642 - /// It returns an unsigned integer random number from the range
1.643 - /// \f$ \{0, 1, \dots, n - 1\} \f$.
1.644 - long getInt(long n) {
1.645 - long mask = n - 1, rnd;
1.646 - mask |= mask >> 1;
1.647 - mask |= mask >> 2;
1.648 - mask |= mask >> 4;
1.649 - mask |= mask >> 8;
1.650 - mask |= mask >> 16;
1.651 - do rnd = getUnsigned() & mask; while (rnd >= n);
1.652 - return rnd;
1.653 - }
1.654 -
1.655 - private:
1.656 -
1.657 - void initState(unsigned long seed) {
1.658 - static const unsigned long mul = 0x6c078965ul;
1.659 -
1.660 - current = state;
1.661 -
1.662 - unsigned long *curr = state + length - 1;
1.663 - curr[0] = seed; --curr;
1.664 - for (int i = 1; i < length; ++i) {
1.665 - curr[0] = (mul * ( curr[1] ^ (curr[1] >> 30) ) + i);
1.666 - --curr;
1.667 - }
1.668 - }
1.669 -
1.670 - void fillState() {
1.671 - static const unsigned long mask[2] = { 0x0ul, 0x9908B0DFul };
1.672 - static const unsigned long loMask = (1ul << 31) - 1;
1.673 - static const unsigned long hiMask = ~loMask;
1.674 -
1.675 - current = state + length;
1.676 -
1.677 - register unsigned long *curr = state + length - 1;
1.678 - register long num;
1.679 -
1.680 - num = length - shift;
1.681 - while (num--) {
1.682 - curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
1.683 - curr[- shift] ^ mask[curr[-1] & 1ul];
1.684 - --curr;
1.685 - }
1.686 - num = shift - 1;
1.687 - while (num--) {
1.688 - curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
1.689 - curr[length - shift] ^ mask[curr[-1] & 1ul];
1.690 - --curr;
1.691 - }
1.692 - curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
1.693 - curr[length - shift] ^ mask[curr[length - 1] & 1ul];
1.694 -
1.695 - }
1.696 -
1.697 - unsigned long *current;
1.698 - unsigned long state[length];
1.699
1.700 };
1.701
1.702 -#elif WORD_BIT == 64
1.703 -
1.704 - /// \ingroup misc
1.705 - ///
1.706 - /// \brief Mersenne Twister random number generator
1.707 - ///
1.708 - /// The Mersenne Twister is a twisted generalized feedback
1.709 - /// shift-register generator of Matsumoto and Nishimura. The period of
1.710 - /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in
1.711 - /// 311 dimensions. The time performance of this generator is comparable
1.712 - /// to the commonly used generators.
1.713 - class Random {
1.714 -
1.715 - static const int length = 312;
1.716 - static const int shift = 156;
1.717 -
1.718 - public:
1.719 -
1.720 - static const unsigned long min = 0ul;
1.721 - static const unsigned long max = ~0ul;
1.722 -
1.723 - /// \brief Constructor
1.724 - ///
1.725 - /// Constructor with time dependent seeding.
1.726 - Random() { initState(std::time(0)); }
1.727 -
1.728 - /// \brief Constructor
1.729 - ///
1.730 - /// Constructor
1.731 - Random(unsigned long seed) { initState(seed); }
1.732 -
1.733 - /// \brief Copy constructor
1.734 - ///
1.735 - /// Copy constructor. The generated sequence will be identical to
1.736 - /// the other sequence.
1.737 - Random(const Random& other) {
1.738 - std::copy(other.state, other.state + length, state);
1.739 - current = state + (other.current - other.state);
1.740 - }
1.741 -
1.742 - /// \brief Assign operator
1.743 - ///
1.744 - /// Assign operator. The generated sequence will be identical to
1.745 - /// the other sequence.
1.746 - Random& operator=(const Random& other) {
1.747 - if (&other != this) {
1.748 - std::copy(other.state, other.state + length, state);
1.749 - current = state + (other.current - other.state);
1.750 - }
1.751 - return *this;
1.752 - }
1.753 -
1.754 - /// \brief Returns an unsigned random number
1.755 - ///
1.756 - /// It returns an unsigned integer random number from the range
1.757 - /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$.
1.758 - unsigned long getUnsigned() {
1.759 - if (current == state) fillState();
1.760 - --current;
1.761 - unsigned long rnd = *current;
1.762 - rnd ^= (rnd >> 29) & 0x5555555555555555ul;
1.763 - rnd ^= (rnd << 17) & 0x71D67FFFEDA60000ul;
1.764 - rnd ^= (rnd << 37) & 0xFFF7EEE000000000ul;
1.765 - rnd ^= (rnd >> 43);
1.766 - return rnd;
1.767 - }
1.768 -
1.769 - /// \brief Returns a random number
1.770 - ///
1.771 - /// It returns an integer random number from the range
1.772 - /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$.
1.773 - long getInt() {
1.774 - return (long)getUnsigned();
1.775 - }
1.776 -
1.777 - /// \brief Returns an unsigned integer random number
1.778 - ///
1.779 - /// It returns an unsigned integer random number from the range
1.780 - /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$.
1.781 - long getNatural() {
1.782 - return (long)(getUnsigned() >> 1);
1.783 - }
1.784 -
1.785 - /// \brief Returns a random bool
1.786 - ///
1.787 - /// It returns a random bool.
1.788 - bool getBool() {
1.789 - return (bool)(getUnsigned() & 1);
1.790 - }
1.791 -
1.792 - /// \brief Returns a real random number
1.793 - ///
1.794 - /// It returns a real random number from the range
1.795 - /// \f$ [0, 1) \f$.
1.796 - double getReal() {
1.797 - return std::ldexp((double)(getUnsigned() >> 11), -53);
1.798 - }
1.799 -
1.800 - /// \brief Returns a real random number
1.801 - ///
1.802 - /// It returns a real random number from the range
1.803 - /// \f$ [0, 1) \f$. This function is identical to the \c getReal().
1.804 - double getPrecReal() {
1.805 - return getReal();
1.806 - }
1.807 -
1.808 - /// \brief Returns an unsigned random number from a given range
1.809 - ///
1.810 - /// It returns an unsigned integer random number from the range
1.811 - /// \f$ \{0, 1, \dots, n - 1\} \f$.
1.812 - unsigned long getUnsigned(unsigned long n) {
1.813 - unsigned long mask = n - 1, rnd;
1.814 - mask |= mask >> 1;
1.815 - mask |= mask >> 2;
1.816 - mask |= mask >> 4;
1.817 - mask |= mask >> 8;
1.818 - mask |= mask >> 16;
1.819 - mask |= mask >> 32;
1.820 - do rnd = getUnsigned() & mask; while (rnd >= n);
1.821 - return rnd;
1.822 - }
1.823 -
1.824 - /// \brief Returns a random number from a given range
1.825 - ///
1.826 - /// It returns an unsigned integer random number from the range
1.827 - /// \f$ \{0, 1, \dots, n - 1\} \f$.
1.828 - long getInt(long n) {
1.829 - long mask = n - 1, rnd;
1.830 - mask |= mask >> 1;
1.831 - mask |= mask >> 2;
1.832 - mask |= mask >> 4;
1.833 - mask |= mask >> 8;
1.834 - mask |= mask >> 16;
1.835 - mask |= mask >> 32;
1.836 - do rnd = getUnsigned() & mask; while (rnd >= n);
1.837 - return rnd;
1.838 - }
1.839 -
1.840 - private:
1.841 -
1.842 - void initState(unsigned long seed) {
1.843 -
1.844 - static const unsigned long mul = 0x5851F42D4C957F2Dul;
1.845 -
1.846 - current = state;
1.847 -
1.848 - unsigned long *curr = state + length - 1;
1.849 - curr[0] = seed; --curr;
1.850 - for (int i = 1; i < length; ++i) {
1.851 - curr[0] = (mul * ( curr[1] ^ (curr[1] >> 62) ) + i);
1.852 - --curr;
1.853 - }
1.854 - }
1.855 -
1.856 - void fillState() {
1.857 - static const unsigned long mask[2] = { 0x0ul, 0xB5026F5AA96619E9ul };
1.858 - static const unsigned long loMask = (1ul << 31) - 1;
1.859 - static const unsigned long hiMask = ~loMask;
1.860 -
1.861 - current = state + length;
1.862 -
1.863 - register unsigned long *curr = state + length - 1;
1.864 - register int num;
1.865 -
1.866 - num = length - shift;
1.867 - while (num--) {
1.868 - curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
1.869 - curr[- shift] ^ mask[curr[-1] & 1ul];
1.870 - --curr;
1.871 - }
1.872 - num = shift - 1;
1.873 - while (num--) {
1.874 - curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
1.875 - curr[length - shift] ^ mask[curr[-1] & 1ul];
1.876 - --curr;
1.877 - }
1.878 - curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
1.879 - curr[length - shift] ^ mask[curr[length - 1] & 1ul];
1.880 -
1.881 - }
1.882 -
1.883 - unsigned long *current;
1.884 - unsigned long state[length];
1.885 -
1.886 - };
1.887 -
1.888 -#else
1.889 -
1.890 - /// \ingroup misc
1.891 - ///
1.892 - /// \brief Mersenne Twister random number generator
1.893 - ///
1.894 - /// The Mersenne Twister is a twisted generalized feedback
1.895 - /// shift-register generator of Matsumoto and Nishimura. There is
1.896 - /// two different implementation in the Lemon library, one for
1.897 - /// 32-bit processors and one for 64-bit processors. The period of
1.898 - /// the generated sequence is \f$ 2^{19937} - 1 \f$, the generated
1.899 - /// sequence of 32-bit random numbers is equi-distributed in 623
1.900 - /// dimensions. The time performance of this generator is comparable
1.901 - /// to the commonly used generators.
1.902 - class Random {
1.903 - public:
1.904 -
1.905 - static const unsigned long min = 0ul;
1.906 - static const unsigned long max = ~0ul;
1.907 -
1.908 - /// \brief Constructor
1.909 - ///
1.910 - /// Constructor with time dependent seeding.
1.911 - Random() {}
1.912 -
1.913 - /// \brief Constructor
1.914 - ///
1.915 - /// Constructor
1.916 - Random(unsigned long seed) {}
1.917 -
1.918 - /// \brief Copy constructor
1.919 - ///
1.920 - /// Copy constructor. The generated sequence will be identical to
1.921 - /// the other sequence.
1.922 - Random(const Random& other) {}
1.923 -
1.924 - /// \brief Assign operator
1.925 - ///
1.926 - /// Assign operator. The generated sequence will be identical to
1.927 - /// the other sequence.
1.928 - Random& operator=(const Random& other) { return *this; }
1.929 -
1.930 - /// \brief Returns an unsigned random number
1.931 - ///
1.932 - /// It returns an unsigned integer random number from the range
1.933 - /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$ for 64-bit processors and from
1.934 - /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$ for 32-bit processors.
1.935 - unsigned long getUnsigned() { return 0ul; }
1.936 -
1.937 - /// \brief Returns a random number
1.938 - ///
1.939 - /// It returns an integer random number from the range
1.940 - /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$ for 64-bit processors and from
1.941 - /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$ for 32-bit processors.
1.942 - long getInt() { return 0l; }
1.943 -
1.944 - /// \brief Returns an unsigned integer random number
1.945 - ///
1.946 - /// It returns an unsigned integer random number from the range
1.947 - /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$ for 64-bit processors and
1.948 - /// from \f$ \{0, 1, \dots, 2^{31} - 1\} \f$ for 32-bit processors.
1.949 - long getNatural() { return 0l; }
1.950 -
1.951 - /// \brief Returns a random bool
1.952 - ///
1.953 - /// It returns a random bool.
1.954 - bool getBool() { return false; }
1.955 -
1.956 - /// \brief Returns a real random number
1.957 - ///
1.958 - /// It returns a real random number from the range
1.959 - /// \f$ [0, 1) \f$. For 32-bit processors the generated random
1.960 - /// number will just have 32 significant bits.
1.961 - double getReal() { return 0.0; }
1.962 -
1.963 - /// \brief Returns a real random number
1.964 - ///
1.965 - /// It returns a real random number from the range
1.966 - /// \f$ [0, 1) \f$. This function returns random numbers with 53
1.967 - /// significant bits for 32-bit processors. For 64-bit processors
1.968 - /// it is identical to the \c getReal().
1.969 - double getPrecReal() { return 0.0; }
1.970 -
1.971 - /// \brief Returns an unsigned random number from a given range
1.972 - ///
1.973 - /// It returns an unsigned integer random number from the range
1.974 - /// \f$ \{0, 1, \dots, n - 1\} \f$.
1.975 - unsigned long getUnsigned(unsigned long n) { return 0; }
1.976 -
1.977 - /// \brief Returns a random number from a given range
1.978 - ///
1.979 - /// It returns an unsigned integer random number from the range
1.980 - /// \f$ \{0, 1, \dots, n - 1\} \f$.
1.981 - long getInt(long n) { return 0; }
1.982 -
1.983 - };
1.984 -
1.985 -#endif
1.986
1.987 /// \brief Global random number generator instance
1.988 ///