RhodeCode

 * 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 <limits>

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

  /// initialization and generation phase so they produce two

  /// completly different sequences.

  ///

  /// The generator gives back random numbers of serveral types. To

  /// get a random number from a range of a floating point type you

  /// can use one form of the \c operator() or the \c real() member

  /// function. If you want to get random number from the {0, 1, ...,

  /// n-1} integer range use the \c operator[] or the \c integer()

  /// method. And to get random number from the whole range of an

  /// integer type you can use the argumentless \c integer() or \c

  /// uinteger() functions. After all you can get random bool with

  /// equal chance of true and false or given probability of true

  /// result with the \c boolean() member functions.

  ///

  ///\code

  /// // The commented code is identical to the other

  /// double a = rnd();                     // [0.0, 1.0)

  /// // double a = rnd.real();             // [0.0, 1.0)

  /// double b = rnd(100.0);                // [0.0, 100.0)

  /// // double b = rnd.real(100.0);        // [0.0, 100.0)

  /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)

  /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)

  /// int d = rnd[100000];                  // 0..99999

  /// // int d = rnd.integer(100000);       // 0..99999

  /// int e = rnd[6] + 1;                   // 1..6

  /// // int e = rnd.integer(1, 1 + 6);     // 1..6

  /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1

  /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1

  /// bool g = rnd.boolean();               // P(g = true) = 0.5

  /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8

  ///\endcode

  ///

  /// LEMON provides a global instance of the random number

  /// generator which name is \ref lemon::rnd "rnd". Usually it is a

  /// good programming convenience to use this global generator to get

  /// random numbers.

  class Random {

  private:

    // Architecture word

    typedef unsigned long Word;

    _random_bits::RandomCore<Word> core;

    _random_bits::BoolProducer<Word> bool_producer;

  public:

    /// \brief Default constructor

    ///

    /// Constructor with constant seeding.

    Random() { core.initState(); }

    /// \brief Constructor with seed

    ///

    /// Constructor with seed. The current number type will be converted

    /// to the architecture word type.

    template <typename Number>

    Random(Number seed) { 

      _random_bits::Initializer<Number, Word>::init(core, seed);

    }

    /// \brief Constructor with array seeding

    ///

    /// Constructor with array seeding. The given range should contain

    /// any number type and the numbers will be converted to the

    /// architecture word type.

    template <typename Iterator>

    Random(Iterator begin, Iterator end) { 

      typedef typename std::iterator_traits<Iterator>::value_type Number;

      _random_bits::Initializer<Number, Word>::init(core, begin, end);

    }

    /// \brief Copy constructor

    ///

    /// Copy constructor. The generated sequence will be identical to

    /// the other sequence. It can be used to save the current state

    /// of the generator and later use it to generate the same

    /// sequence.

    Random(const Random& other) {

      core.copyState(other.core);

    }

    /// \brief Assign operator

    ///

    /// Assign operator. The generated sequence will be identical to

    /// the other sequence. It can be used to save the current state

    /// of the generator and later use it to generate the same

    /// sequence.

    Random& operator=(const Random& other) {

      if (&other != this) {

        core.copyState(other.core);

      }

      return *this;

    }

    /// \brief Seeding random sequence

    ///

    /// Seeding the random sequence. The current number type will be

    /// converted to the architecture word type.

    template <typename Number>

    void seed(Number seed) { 

      _random_bits::Initializer<Number, Word>::init(core, seed);

    }

    /// \brief Seeding random sequence

    ///

    /// Seeding the random sequence. The given range should contain

    /// any number type and the numbers will be converted to the

    /// architecture word type.

    template <typename Iterator>

    void seed(Iterator begin, Iterator end) { 

      typedef typename std::iterator_traits<Iterator>::value_type Number;

      _random_bits::Initializer<Number, Word>::init(core, begin, end);

    }

    /// \brief Returns a random real number from the range [0, 1)

    ///

    /// It returns a random real number from the range [0, 1). The

    /// default Number type is \c double.

    template <typename Number>

    Number real() {

      return _random_bits::RealConversion<Number, Word>::convert(core);

    }

    double real() {

      return real<double>();

    }

    /// \brief Returns a random real number the range [0, b)

    ///

    /// It returns a random real number from the range [0, b).

    template <typename Number>

    Number real(Number b) { 

      return real<Number>() * b; 

    }

    /// \brief Returns a random real number from the range [a, b)

    ///

    /// It returns a random real number from the range [a, b).

    template <typename Number>

    Number real(Number a, Number b) { 

      return real<Number>() * (b - a) + a; 

    }

    /// \brief Returns a random real number from the range [0, 1)

    ///

    /// It returns a random double from the range [0, 1).

    double operator()() {

      return real<double>();

    }

    /// \brief Returns a random real number from the range [0, b)

    ///

    /// It returns a random real number from the range [0, b).

    template <typename Number>

    Number operator()(Number b) { 

      return real<Number>() * b; 

    }

    /// \brief Returns a random real number from the range [a, b)

    ///

    /// It returns a random real number from the range [a, b).

    template <typename Number>

    /// \brief Returns a random integer from a range

    ///

    /// It returns a random integer from the range {0, 1, ..., b - 1}.

    template <typename Number>

    Number operator[](Number b) {

      return _random_bits::Mapping<Number, Word>::map(core, b);

    }

    /// \brief Returns a random non-negative integer

    ///

    /// It returns a random non-negative integer uniformly from the

    /// whole range of the current \c Number type. The default result

    /// type of this function is <tt>unsigned int</tt>.

    template <typename Number>

    Number uinteger() {

      return _random_bits::IntConversion<Number, Word>::convert(core);

    }

    unsigned int uinteger() {

      return uinteger<unsigned int>();

    }

    /// \brief Returns a random integer

    ///

    /// It returns a random integer uniformly from the whole range of

    /// the current \c Number type. The default result type of this

    /// function is \c int.

    template <typename Number>

    Number integer() {

      static const int nb = std::numeric_limits<Number>::digits + 

        (std::numeric_limits<Number>::is_signed ? 1 : 0);

      return _random_bits::IntConversion<Number, Word, nb>::convert(core);

    }

    int integer() {

      return integer<int>();

    }

    /// \brief Returns a random bool

    ///

    /// It returns a random bool. The generator holds a buffer for

    /// random bits. Every time when it become empty the generator makes

    /// a new random word and fill the buffer up.

    bool boolean() {

      return bool_producer.convert(core);

    }

    ///\name Non-uniform distributions

    ///

    ///@{

    /// \brief Returns a random bool

    ///

    /// It returns a random bool with given probability of true result.

    bool boolean(double p) {

      return operator()() < p;

    }

    /// Standard Gauss distribution

    /// Standard Gauss distribution.

    /// \note The Cartesian form of the Box-Muller

    /// transformation is used to generate a random normal distribution.

    /// \todo Consider using the "ziggurat" method instead.

    double gauss() 

    {

      double V1,V2,S;

      do {

	V1=2*real<double>()-1;

	V2=2*real<double>()-1;

	S=V1*V1+V2*V2;

      } while(S>=1);

      return std::sqrt(-2*std::log(S)/S)*V1;

    }

    /// Gauss distribution with given mean and standard deviation

    /// Gauss distribution with given mean and standard deviation.

    /// \sa gauss()

    double gauss(double mean,double std_dev)

    {

      return gauss()*std_dev+mean;

    }

    /// Exponential distribution with given mean

    /// This function generates an exponential distribution random number

    /// with mean <tt>1/lambda</tt>.

    ///

    double exponential(double lambda=1.0)

    {

      return -std::log(1.0-real<double>())/lambda;

    }

    /// Gamma distribution with given integer shape