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

      }