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

  /// \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

  /// \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