/* -*- C++ -*-

 *

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

 *

 * Copyright (C) 2003-2006

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

 *

 */

#ifndef LEMON_RANDOM_H

#define LEMON_RANDOM_H

#include <algorithm>

#include <ctime>

#include <cmath>

#include <cmath>

///\ingroup misc

///\file

///\brief Mersenne Twister random number generator

///

///\author Balazs Dezso

namespace lemon {

#if WORD_BIT == 32

  /// \ingroup misc

  ///

  /// \brief Mersenne Twister random number generator

  ///

  /// The Mersenne Twister is a twisted generalized feedback

  /// shift-register generator of Matsumoto and Nishimura. The period of

  /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in

  /// 623 dimensions. The time performance of this generator is comparable

  /// to the commonly used generators.

  ///

  /// \author Balazs Dezso

  class Random {

    static const int length = 624;

    static const int shift = 397;

  public:

    static const unsigned long min = 0ul;

    static const unsigned long max = ~0ul;

    /// \brief Constructor

    ///

    /// Constructor with time dependent seeding.

    Random() { initState(std::time(0)); }

    /// \brief Constructor

    ///

    /// Constructor

    Random(unsigned long seed) { initState(seed); }

    /// \brief Copy constructor

    ///

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

    /// the other sequence.

    Random(const Random& other) { 

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

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

    }

    /// \brief Assign operator

    ///

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

    /// the other sequence.

    Random& operator=(const Random& other) {

      if (&other != this) {

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

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

      }

      return *this;

    }

    /// \brief Returns an unsigned random number

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$.

    unsigned long getUnsigned() {

      if (current == state) fillState();

      --current;

      unsigned long rnd = *current;

      rnd ^= (rnd >> 11);

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

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

      rnd ^= (rnd >> 18);

      return rnd;

    }

    /// \brief Returns a random number

    ///

    /// It returns an integer random number from the range 

    /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$.

    long getInt() {

      return (long)getUnsigned();

    }

    /// \brief Returns an unsigned integer random number

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, 2^{31} - 1\} \f$.

    long getNatural() {

      return (long)(getUnsigned() >> 1);

    }

    /// \brief Returns a random bool

    ///

    /// It returns a random bool.

    bool getBool() {

      return (bool)(getUnsigned() & 1);

    }

    /// \brief Returns a real random number

    ///

    /// It returns a real random number from the range 

    /// \f$ [0, 1) \f$. The double will have 32 significant bits.

    double getReal() {

      return std::ldexp((double)getUnsigned(), -32);

    }

    /// \brief Returns a real random number

    ///

    /// It returns a real random number from the range 

    /// \f$ [0, 1) \f$. The double will have 53 significant bits,

    /// but it is slower than the \c getReal().

    double getPrecReal() {

      return std::ldexp((double)(getUnsigned() >> 5), -27) + 

        std::ldexp((double)(getUnsigned() >> 6), -53);

    }

    /// \brief Returns an unsigned random number from a given range

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, n - 1\} \f$.

    unsigned long getUnsigned(unsigned long n) {

      unsigned long mask = n - 1, rnd;

      mask |= mask >> 1;

      mask |= mask >> 2;

      mask |= mask >> 4;

      mask |= mask >> 8;

      mask |= mask >> 16;

      do rnd = getUnsigned() & mask; while (rnd >= n);

      return rnd;

    }

    /// \brief Returns a random number from a given range

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, n - 1\} \f$.

    long getInt(long n) {

      long mask = n - 1, rnd;

      mask |= mask >> 1;

      mask |= mask >> 2;

      mask |= mask >> 4;

      mask |= mask >> 8;

      mask |= mask >> 16;

      do rnd = getUnsigned() & mask; while (rnd >= n);

      return rnd;

    }

  private:

    void initState(unsigned long seed) {

      static const unsigned long mul = 0x6c078965ul;

      current = state; 

      unsigned long *curr = state + length - 1;

      curr[0] = seed; --curr;

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

        curr[0] = (mul * ( curr[1] ^ (curr[1] >> 30) ) + i);

        --curr;

      }

    }

    void fillState() {

      static const unsigned long mask[2] = { 0x0ul, 0x9908B0DFul };

      static const unsigned long loMask = (1ul << 31) - 1;

      static const unsigned long hiMask = ~loMask;

      current = state + length; 

      register unsigned long *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[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^

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

    }

    unsigned long *current;

    unsigned long state[length];

  };

#elif WORD_BIT == 64

  /// \ingroup misc

  ///

  /// \brief Mersenne Twister random number generator

  ///

  /// The Mersenne Twister is a twisted generalized feedback

  /// shift-register generator of Matsumoto and Nishimura. The period of

  /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in

  /// 311 dimensions. The time performance of this generator is comparable

  /// to the commonly used generators.

  class Random {

    static const int length = 312;

    static const int shift = 156;

  public:

    static const unsigned long min = 0ul;

    static const unsigned long max = ~0ul;

    /// \brief Constructor

    ///

    /// Constructor with time dependent seeding.

    Random() { initState(std::time(0)); }

    /// \brief Constructor

    ///

    /// Constructor

    Random(unsigned long seed) { initState(seed); }

    /// \brief Copy constructor

    ///

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

    /// the other sequence.

    Random(const Random& other) { 

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

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

    }

    /// \brief Assign operator

    ///

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

    /// the other sequence.

    Random& operator=(const Random& other) {

      if (&other != this) {

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

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

      }

      return *this;

    }

    /// \brief Returns an unsigned random number

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$.

    unsigned long getUnsigned() {

      if (current == state) fillState();

      --current;

      unsigned long rnd = *current;

      rnd ^= (rnd >> 29) & 0x5555555555555555ul;

      rnd ^= (rnd << 17) & 0x71D67FFFEDA60000ul;

      rnd ^= (rnd << 37) & 0xFFF7EEE000000000ul;

      rnd ^= (rnd >> 43);

      return rnd;

    }

    /// \brief Returns a random number

    ///

    /// It returns an integer random number from the range 

    /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$.

    long getInt() {

      return (long)getUnsigned();

    }

    /// \brief Returns an unsigned integer random number

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$.

    long getNatural() {

      return (long)(getUnsigned() >> 1);

    }

    /// \brief Returns a random bool

    ///

    /// It returns a random bool.

    bool getBool() {

      return (bool)(getUnsigned() & 1);

    }

    /// \brief Returns a real random number

    ///

    /// It returns a real random number from the range 

    /// \f$ [0, 1) \f$.

    double getReal() {

      return std::ldexp((double)(getUnsigned() >> 11), -53);

    }

    /// \brief Returns a real random number

    ///

    /// It returns a real random number from the range 

    /// \f$ [0, 1) \f$. This function is identical to the \c getReal().

    double getPrecReal() {

      return getReal();

    }

    /// \brief Returns an unsigned random number from a given range

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, n - 1\} \f$.

    unsigned long getUnsigned(unsigned long n) {

      unsigned long mask = n - 1, rnd;

      mask |= mask >> 1;

      mask |= mask >> 2;

      mask |= mask >> 4;

      mask |= mask >> 8;

      mask |= mask >> 16;

      mask |= mask >> 32;

      do rnd = getUnsigned() & mask; while (rnd >= n);

      return rnd;

    }

    /// \brief Returns a random number from a given range

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, n - 1\} \f$.

    long getInt(long n) {

      long mask = n - 1, rnd;

      mask |= mask >> 1;

      mask |= mask >> 2;

      mask |= mask >> 4;

      mask |= mask >> 8;

      mask |= mask >> 16;

      mask |= mask >> 32;

      do rnd = getUnsigned() & mask; while (rnd >= n);

      return rnd;

    }

  private:

    void initState(unsigned long seed) {

      static const unsigned long mul = 0x5851F42D4C957F2Dul;

      current = state; 

      unsigned long *curr = state + length - 1;

      curr[0] = seed; --curr;

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

        curr[0] = (mul * ( curr[1] ^ (curr[1] >> 62) ) + i);

        --curr;

      }

    }

    void fillState() {

      static const unsigned long mask[2] = { 0x0ul, 0xB5026F5AA96619E9ul };

      static const unsigned long loMask = (1ul << 31) - 1;

      static const unsigned long hiMask = ~loMask;

      current = state + length; 

      register unsigned long *curr = state + length - 1;

      register int 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[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^

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

    }

    unsigned long *current;

    unsigned long state[length];

  };

#else

  /// \ingroup misc

  ///

  /// \brief Mersenne Twister random number generator

  ///

  /// The Mersenne Twister is a twisted generalized feedback

  /// shift-register generator of Matsumoto and Nishimura. There is

  /// two different implementation in the Lemon library, one for

  /// 32-bit processors and one for 64-bit processors. The period of

  /// the generated sequence is \f$ 2^{19937} - 1 \f$, the generated

  /// sequence of 32-bit random numbers is equi-distributed in 623

  /// dimensions. The time performance of this generator is comparable

  /// to the commonly used generators.

  class Random {

  public:

    static const unsigned long min = 0ul;

    static const unsigned long max = ~0ul;

    /// \brief Constructor

    ///

    /// Constructor with time dependent seeding.

    Random() {}

    /// \brief Constructor

    ///

    /// Constructor

    Random(unsigned long seed) {}

    /// \brief Copy constructor

    ///

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

    /// the other sequence.

    Random(const Random& other) {}

    /// \brief Assign operator

    ///

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

    /// the other sequence.

    Random& operator=(const Random& other) { return *this; }

    /// \brief Returns an unsigned random number

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$ for 64-bit processors and from

    /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$ for 32-bit processors.

    unsigned long getUnsigned() { return 0ul; }

    /// \brief Returns a random number

    ///

    /// It returns an integer random number from the range 

    /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$ for 64-bit processors and from

    /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$ for 32-bit processors.

    long getInt() { return 0l; }

    /// \brief Returns an unsigned integer random number

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$ for 64-bit processors and

    /// from \f$ \{0, 1, \dots, 2^{31} - 1\} \f$ for 32-bit processors.

    long getNatural() { return 0l; }

    /// \brief Returns a random bool

    ///

    /// It returns a random bool.

    bool getBool() { return false; }

    /// \brief Returns a real random number

    ///

    /// It returns a real random number from the range 

    /// \f$ [0, 1) \f$. For 32-bit processors the generated random

    /// number will just have 32 significant bits.

    double getReal() { return 0.0; }

    /// \brief Returns a real random number

    ///

    /// It returns a real random number from the range 

    /// \f$ [0, 1) \f$. This function returns random numbers with 53

    /// significant bits for 32-bit processors. For 64-bit processors

    /// it is identical to the \c getReal().

    double getPrecReal() { return 0.0; }

    /// \brief Returns an unsigned random number from a given range

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, n - 1\} \f$.

    unsigned long getUnsigned(unsigned long n) { return 0; }

    /// \brief Returns a random number from a given range

    ///

    /// It returns an unsigned integer random number from the range 

    /// \f$ \{0, 1, \dots, n - 1\} \f$.

    long getInt(long n) { return 0; }

  };

#endif

  /// \brief Global random number generator instance

  ///

  /// A global mersenne twister random number generator instance

  extern Random rnd;

}

#endif