RhodeCode

/* -*- C++ -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library

 *

 * Copyright (C) 2003-2008

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

/*

 * This file contains the reimplemented version of the Mersenne Twister

 * Generator of Matsumoto and Nishimura.

 *

 * See the appropriate copyright notice below.

 * 

 * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,

 * All rights reserved.                          

 *

 * Redistribution and use in source and binary forms, with or without

 * modification, are permitted provided that the following conditions

 * are met:

 *

 * 1. Redistributions of source code must retain the above copyright

 *    notice, this list of conditions and the following disclaimer.

 *

 * 2. Redistributions in binary form must reproduce the above copyright

 *    notice, this list of conditions and the following disclaimer in the

 *    documentation and/or other materials provided with the distribution.

 *

 * 3. The names of its contributors may not be used to endorse or promote 

 *    products derived from this software without specific prior written 

 *    permission.

 *

 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS

 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS

 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE

 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,

 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES

 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR

 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)

 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,

 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)

 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED

 * OF THE POSSIBILITY OF SUCH DAMAGE.

 *

 *

 * Any feedback is very welcome.

 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html

 * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)

 */

#ifndef LEMON_RANDOM_H

#define LEMON_RANDOM_H

#include <algorithm>

#include <iterator>

#include <vector>

#include <lemon/math.h>

#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;

        }

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

          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;

      int num;

      BoolProducer() : num(0) {}