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