RhodeCode

#include <algorithm>

#include <iterator>

#include <vector>

#include <ctime>

#include <cmath>

#include <lemon/dim2.h>

///\ingroup misc

///\file

///\brief Mersenne Twister random number generator

namespace lemon {

  namespace _random_bits {

    template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>

    struct RandomTraits {};

    template <typename _Word>

    struct RandomTraits<_Word, 32> {

      typedef _Word Word;

      static const int bits = 32;

      static const int length = 624;

      static const int shift = 397;

      static const Word mul = 0x6c078965u;

      static const Word arrayInit = 0x012BD6AAu;

      static const Word arrayMul1 = 0x0019660Du;

      static const Word arrayMul2 = 0x5D588B65u;

      static const Word mask = 0x9908B0DFu;

      static const Word loMask = (1u << 31) - 1;

      static const Word hiMask = ~loMask;

      static Word tempering(Word rnd) {

        rnd ^= (rnd >> 11);

        rnd ^= (rnd << 7) & 0x9D2C5680u;

        rnd ^= (rnd << 15) & 0xEFC60000u;

        rnd ^= (rnd >> 18);

        return rnd;

      }

    };

    template <typename _Word>

    struct RandomTraits<_Word, 64> {

      typedef _Word Word;

      static const int bits = 64;

      static const int length = 312;

      static const int shift = 156;

      static const Word mul = Word(0x5851F42Du) << 32 | Word(0x4C957F2Du);

      static const Word arrayInit = Word(0x00000000u) << 32 |Word(0x012BD6AAu);

      static const Word arrayMul1 = Word(0x369DEA0Fu) << 32 |Word(0x31A53F85u);

      static const Word arrayMul2 = Word(0x27BB2EE6u) << 32 |Word(0x87B0B0FDu);

      static const Word mask = Word(0xB5026F5Au) << 32 | Word(0xA96619E9u);

      static const Word loMask = (Word(1u) << 31) - 1;

      static const Word hiMask = ~loMask;

      static Word tempering(Word rnd) {

        rnd ^= (rnd >> 29) & (Word(0x55555555u) << 32 | Word(0x55555555u));

        rnd ^= (rnd << 17) & (Word(0x71D67FFFu) << 32 | Word(0xEDA60000u));

        rnd ^= (rnd << 37) & (Word(0xFFF7EEE0u) << 32 | Word(0x00000000u));

        rnd ^= (rnd >> 43);

        return rnd;

      }

    };

    template <typename _Word>

    class RandomCore {

    public:

      typedef _Word Word;

    private:

      static const int bits = RandomTraits<Word>::bits;

      static const int length = RandomTraits<Word>::length;

      static const int shift = RandomTraits<Word>::shift;

    public:

      void initState() {

        static const Word seedArray[4] = {

          0x12345u, 0x23456u, 0x34567u, 0x45678u

        };

        initState(seedArray, seedArray + 4);

      }

      void initState(Word seed) {

        static const Word mul = RandomTraits<Word>::mul;

        current = state; 

        Word *curr = state + length - 1;

        curr[0] = seed; --curr;

        for (int i = 1; i < length; ++i) {

          curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i);

          --curr;

        }

      }

      template <typename Iterator>

      void initState(Iterator begin, Iterator end) {

        static const Word init = RandomTraits<Word>::arrayInit;

        static const Word mul1 = RandomTraits<Word>::arrayMul1;

        static const Word mul2 = RandomTraits<Word>::arrayMul2;

        Word *curr = state + length - 1; --curr;

        Iterator it = begin; int cnt = 0;

        int num;

        initState(init);

        num = length > end - begin ? length : end - begin;

        while (num--) {

          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1)) 

            + *it + cnt;

          ++it; ++cnt;

          if (it == end) {

            it = begin; cnt = 0;

          }

          if (curr == state) {

            curr = state + length - 1; curr[0] = state[0];

          }

          --curr;

        }

        num = length - 1; cnt = length - (curr - state) - 1;

        while (num--) {

          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2))

            - cnt;

          --curr; ++cnt;

          if (curr == state) {

            curr = state + length - 1; curr[0] = state[0]; --curr;

            cnt = 1;

          }

        }

        state[length - 1] = Word(1) << (bits - 1);

      }

      void copyState(const RandomCore& other) {

        std::copy(other.state, other.state + length, state);

        current = state + (other.current - other.state);

      }

      Word operator()() {

        if (current == state) fillState();

        --current;

        Word rnd = *current;

        return RandomTraits<Word>::tempering(rnd);

      }

    private:

      void fillState() {

        static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };

        static const Word loMask = RandomTraits<Word>::loMask;

        static const Word hiMask = RandomTraits<Word>::hiMask;

        current = state + length; 

        register Word *curr = state + length - 1;

        register long num;

        num = length - shift;

        while (num--) {

          curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^

            curr[- shift] ^ mask[curr[-1] & 1ul];

          --curr;

        }

        num = shift - 1;

        while (num--) {

          curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^

            curr[length - shift] ^ mask[curr[-1] & 1ul];

          --curr;

        }

          curr[length - shift] ^ mask[curr[length - 1] & 1ul];

      }

      Word *current;

      Word state[length];

    };

    template <typename Result, 

              int shift = (std::numeric_limits<Result>::digits + 1) / 2>

    struct Masker {

      static Result mask(const Result& result) {

        return Masker<Result, (shift + 1) / 2>::

          mask(static_cast<Result>(result | (result >> shift)));

      }

    };

    template <typename Result>

    struct Masker<Result, 1> {

      static Result mask(const Result& result) {

        return static_cast<Result>(result | (result >> 1));

      }

    };

    template <typename Result, typename Word, 

              int rest = std::numeric_limits<Result>::digits, int shift = 0, 

              bool last = rest <= std::numeric_limits<Word>::digits>

    struct IntConversion {

      static const int bits = std::numeric_limits<Word>::digits;

      static Result convert(RandomCore<Word>& rnd) {

        return static_cast<Result>(rnd() >> (bits - rest)) << shift;

      }

    }; 

    template <typename Result, typename Word, int rest, int shift> 

    struct IntConversion<Result, Word, rest, shift, false> {

      static const int bits = std::numeric_limits<Word>::digits;

      static Result convert(RandomCore<Word>& rnd) {

        return (static_cast<Result>(rnd()) << shift) | 

          IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);

      }

    };

    template <typename Result, typename Word,

              bool one_word = (std::numeric_limits<Word>::digits < 

			       std::numeric_limits<Result>::digits) >

    struct Mapping {

      static Result map(RandomCore<Word>& rnd, const Result& bound) {

        Word max = Word(bound - 1);

        Result mask = Masker<Result>::mask(bound - 1);

        Result num;

        do {

          num = IntConversion<Result, Word>::convert(rnd) & mask; 

        } while (num > max);

        return num;

      }

    };

    template <typename Result, typename Word>

    struct Mapping<Result, Word, false> {

      static Result map(RandomCore<Word>& rnd, const Result& bound) {

        Word max = Word(bound - 1);

        Word mask = Masker<Word, (std::numeric_limits<Result>::digits + 1) / 2>

          ::mask(max);

        Word num;

        do {

          num = rnd() & mask;

        } while (num > max);

        return num;

      }

    };

    template <typename Result, int exp, bool pos = (exp >= 0)>

    struct ShiftMultiplier {

      static const Result multiplier() {

        Result res = ShiftMultiplier<Result, exp / 2>::multiplier();

        res *= res;

        if ((exp & 1) == 1) res *= static_cast<Result>(2.0);

        return res; 

      }

    };

    template <typename Result, int exp>

    struct ShiftMultiplier<Result, exp, false> {

      static const Result multiplier() {

        Result res = ShiftMultiplier<Result, exp / 2>::multiplier();

        res *= res;

        if ((exp & 1) == 1) res *= static_cast<Result>(0.5);

        return res; 

      }

    };

    template <typename Result>

    struct ShiftMultiplier<Result, 0, true> {

      static const Result multiplier() {

        return static_cast<Result>(1.0); 

      }

    };

    template <typename Result>

    struct ShiftMultiplier<Result, -20, true> {

      static const Result multiplier() {

        return static_cast<Result>(1.0/1048576.0); 

      }

    };

    template <typename Result>

    struct ShiftMultiplier<Result, -32, true> {

      static const Result multiplier() {

        return static_cast<Result>(1.0/424967296.0); 

      }

    };

    template <typename Result>

    struct ShiftMultiplier<Result, -53, true> {

      static const Result multiplier() {

        return static_cast<Result>(1.0/9007199254740992.0); 

      }

    };

    template <typename Result>

    struct ShiftMultiplier<Result, -64, true> {

      static const Result multiplier() {

        return static_cast<Result>(1.0/18446744073709551616.0); 

      }

    };

    template <typename Result, int exp>

    struct Shifting {

      static Result shift(const Result& result) {

        return result * ShiftMultiplier<Result, exp>::multiplier();

      }

    };

    template <typename Result, typename Word,

              int rest = std::numeric_limits<Result>::digits, int shift = 0, 

              bool last = rest <= std::numeric_limits<Word>::digits>

    struct RealConversion{ 

      static const int bits = std::numeric_limits<Word>::digits;

      static Result convert(RandomCore<Word>& rnd) {

        return Shifting<Result, - shift - rest>::

          shift(static_cast<Result>(rnd() >> (bits - rest)));

      }

    };

    template <typename Result, typename Word, int rest, int shift>

    struct RealConversion<Result, Word, rest, shift, false> { 

      static const int bits = std::numeric_limits<Word>::digits;

      static Result convert(RandomCore<Word>& rnd) {

        return Shifting<Result, - shift - bits>::

          shift(static_cast<Result>(rnd())) +

          RealConversion<Result, Word, rest-bits, shift + bits>::

          convert(rnd);

      }

    };

    template <typename Result, typename Word>

    struct Initializer {

      template <typename Iterator>

      static void init(RandomCore<Word>& rnd, Iterator begin, Iterator end) {

        std::vector<Word> ws;

        for (Iterator it = begin; it != end; ++it) {

          ws.push_back(Word(*it));

        }

        rnd.initState(ws.begin(), ws.end());

      }

      static void init(RandomCore<Word>& rnd, Result seed) {

        rnd.initState(seed);

      }

    };

    template <typename Word>

    struct BoolConversion {

      static bool convert(RandomCore<Word>& rnd) {

        return (rnd() & 1) == 1;

      }

    };

    template <typename Word>

    struct BoolProducer {

      Word buffer;