deba@2229: /* -*- C++ -*-
deba@2229:  *
deba@2229:  * This file is a part of LEMON, a generic C++ optimization library
deba@2229:  *
deba@2229:  * Copyright (C) 2003-2006
deba@2229:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@2229:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@2229:  *
deba@2229:  * Permission to use, modify and distribute this software is granted
deba@2229:  * provided that this copyright notice appears in all copies. For
deba@2229:  * precise terms see the accompanying LICENSE file.
deba@2229:  *
deba@2229:  * This software is provided "AS IS" with no warranty of any kind,
deba@2229:  * express or implied, and with no claim as to its suitability for any
deba@2229:  * purpose.
deba@2229:  *
deba@2229:  */
deba@2229: 
deba@2229: #ifndef LEMON_RANDOM_H
deba@2229: #define LEMON_RANDOM_H
deba@2229: 
deba@2229: #include <algorithm>
deba@2229: 
deba@2229: #include <ctime>
deba@2229: #include <cmath>
deba@2229: #include <cmath>
deba@2229: 
deba@2229: ///\ingroup misc
deba@2229: ///\file
deba@2229: ///\brief Mersenne Twister random number generator
deba@2229: ///
deba@2229: ///\author Balazs Dezso
deba@2229: 
deba@2229: namespace lemon {
deba@2229: 
deba@2229: #if WORD_BIT == 32
deba@2229: 
deba@2229:   /// \ingroup misc
deba@2229:   ///
deba@2229:   /// \brief Mersenne Twister random number generator
deba@2229:   ///
deba@2229:   /// The Mersenne Twister is a twisted generalized feedback
deba@2229:   /// shift-register generator of Matsumoto and Nishimura. The period of
deba@2229:   /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in
deba@2229:   /// 623 dimensions. The time performance of this generator is comparable
deba@2229:   /// to the commonly used generators.
deba@2229:   ///
deba@2229:   /// \author Balazs Dezso
deba@2229:   class Random {
deba@2229: 
deba@2229:     static const int length = 624;
deba@2229:     static const int shift = 397;
deba@2229: 
deba@2229:   public:
deba@2229: 
deba@2229:     static const unsigned long min = 0ul;
deba@2229:     static const unsigned long max = ~0ul;
deba@2229:   
deba@2229:     /// \brief Constructor
deba@2229:     ///
deba@2229:     /// Constructor with time dependent seeding.
deba@2229:     Random() { initState(std::time(0)); }
deba@2229: 
deba@2229:     /// \brief Constructor
deba@2229:     ///
deba@2229:     /// Constructor
deba@2229:     Random(unsigned long seed) { initState(seed); }
deba@2229: 
deba@2229:     /// \brief Copy constructor
deba@2229:     ///
deba@2229:     /// Copy constructor. The generated sequence will be identical to
deba@2229:     /// the other sequence.
deba@2229:     Random(const Random& other) { 
deba@2229:       std::copy(other.state, other.state + length, state);
deba@2229:       current = state + (other.current - other.state);
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Assign operator
deba@2229:     ///
deba@2229:     /// Assign operator. The generated sequence will be identical to
deba@2229:     /// the other sequence.
deba@2229:     Random& operator=(const Random& other) {
deba@2229:       if (&other != this) {
deba@2229:         std::copy(other.state, other.state + length, state);
deba@2229:         current = state + (other.current - other.state);
deba@2229:       }
deba@2229:       return *this;
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns an unsigned random number
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$.
deba@2229:     unsigned long getUnsigned() {
deba@2229:       if (current == state) fillState();
deba@2229:       --current;
deba@2229:       unsigned long rnd = *current;
deba@2229:       rnd ^= (rnd >> 11);
deba@2229:       rnd ^= (rnd << 7) & 0x9D2C5680ul;
deba@2229:       rnd ^= (rnd << 15) & 0xEFC60000ul;
deba@2229:       rnd ^= (rnd >> 18);
deba@2229:       return rnd;
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns a random number
deba@2229:     ///
deba@2229:     /// It returns an integer random number from the range 
deba@2229:     /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$.
deba@2229:     long getInt() {
deba@2229:       return (long)getUnsigned();
deba@2229:     }
deba@2229:     
deba@2229:     /// \brief Returns an unsigned integer random number
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, 2^{31} - 1\} \f$.
deba@2229:     long getNatural() {
deba@2229:       return (long)(getUnsigned() >> 1);
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns a random bool
deba@2229:     ///
deba@2229:     /// It returns a random bool.
deba@2229:     bool getBool() {
deba@2229:       return (bool)(getUnsigned() & 1);
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns a real random number
deba@2229:     ///
deba@2229:     /// It returns a real random number from the range 
deba@2229:     /// \f$ [0, 1) \f$. The double will have 32 significant bits.
deba@2229:     double getReal() {
deba@2229:       return std::ldexp((double)getUnsigned(), -32);
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns a real random number
deba@2229:     ///
deba@2229:     /// It returns a real random number from the range 
deba@2229:     /// \f$ [0, 1) \f$. The double will have 53 significant bits,
deba@2229:     /// but it is slower than the \c getReal().
deba@2229:     double getPrecReal() {
deba@2229:       return std::ldexp((double)(getUnsigned() >> 5), -27) + 
deba@2229:         std::ldexp((double)(getUnsigned() >> 6), -53);
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns an unsigned random number from a given range
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, n - 1\} \f$.
deba@2229:     unsigned long getUnsigned(unsigned long n) {
deba@2229:       unsigned long mask = n - 1, rnd;
deba@2229:       mask |= mask >> 1;
deba@2229:       mask |= mask >> 2;
deba@2229:       mask |= mask >> 4;
deba@2229:       mask |= mask >> 8;
deba@2229:       mask |= mask >> 16;
deba@2229:       do rnd = getUnsigned() & mask; while (rnd >= n);
deba@2229:       return rnd;
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns a random number from a given range
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, n - 1\} \f$.
deba@2229:     long getInt(long n) {
deba@2229:       long mask = n - 1, rnd;
deba@2229:       mask |= mask >> 1;
deba@2229:       mask |= mask >> 2;
deba@2229:       mask |= mask >> 4;
deba@2229:       mask |= mask >> 8;
deba@2229:       mask |= mask >> 16;
deba@2229:       do rnd = getUnsigned() & mask; while (rnd >= n);
deba@2229:       return rnd;
deba@2229:     }
deba@2229: 
deba@2229:   private:
deba@2229: 
deba@2229:     void initState(unsigned long seed) {
deba@2229:       static const unsigned long mul = 0x6c078965ul;
deba@2229: 
deba@2229:       current = state; 
deba@2229: 
deba@2229:       unsigned long *curr = state + length - 1;
deba@2229:       curr[0] = seed; --curr;
deba@2229:       for (int i = 1; i < length; ++i) {
deba@2229:         curr[0] = (mul * ( curr[1] ^ (curr[1] >> 30) ) + i);
deba@2229:         --curr;
deba@2229:       }
deba@2229:     }
deba@2229:   
deba@2229:     void fillState() {
deba@2229:       static const unsigned long mask[2] = { 0x0ul, 0x9908B0DFul };
deba@2229:       static const unsigned long loMask = (1ul << 31) - 1;
deba@2229:       static const unsigned long hiMask = ~loMask;
deba@2229: 
deba@2229:       current = state + length; 
deba@2229: 
deba@2229:       register unsigned long *curr = state + length - 1;
deba@2229:       register long num;
deba@2229:       
deba@2229:       num = length - shift;
deba@2229:       while (num--) {
deba@2229:         curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
deba@2229:           curr[- shift] ^ mask[curr[-1] & 1ul];
deba@2229:         --curr;
deba@2229:       }
deba@2229:       num = shift - 1;
deba@2229:       while (num--) {
deba@2229:         curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
deba@2229:           curr[length - shift] ^ mask[curr[-1] & 1ul];
deba@2229:         --curr;
deba@2229:       }
deba@2229:       curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
deba@2229:         curr[length - shift] ^ mask[curr[length - 1] & 1ul];
deba@2229: 
deba@2229:     }
deba@2229:   
deba@2229:     unsigned long *current;
deba@2229:     unsigned long state[length];
deba@2229:     
deba@2229:   };
deba@2229: 
deba@2229: #elif WORD_BIT == 64
deba@2229: 
deba@2230:   /// \ingroup misc
deba@2230:   ///
deba@2229:   /// \brief Mersenne Twister random number generator
deba@2229:   ///
deba@2229:   /// The Mersenne Twister is a twisted generalized feedback
deba@2229:   /// shift-register generator of Matsumoto and Nishimura. The period of
deba@2229:   /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in
deba@2229:   /// 311 dimensions. The time performance of this generator is comparable
deba@2229:   /// to the commonly used generators.
deba@2229:   class Random {
deba@2229: 
deba@2229:     static const int length = 312;
deba@2229:     static const int shift = 156;
deba@2229: 
deba@2229:   public:
deba@2229: 
deba@2229:     static const unsigned long min = 0ul;
deba@2229:     static const unsigned long max = ~0ul;
deba@2229:   
deba@2229:     /// \brief Constructor
deba@2229:     ///
deba@2229:     /// Constructor with time dependent seeding.
deba@2229:     Random() { initState(std::time(0)); }
deba@2229: 
deba@2229:     /// \brief Constructor
deba@2229:     ///
deba@2229:     /// Constructor
deba@2229:     Random(unsigned long seed) { initState(seed); }
deba@2229: 
deba@2229:     /// \brief Copy constructor
deba@2229:     ///
deba@2229:     /// Copy constructor. The generated sequence will be identical to
deba@2229:     /// the other sequence.
deba@2229:     Random(const Random& other) { 
deba@2229:       std::copy(other.state, other.state + length, state);
deba@2229:       current = state + (other.current - other.state);
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Assign operator
deba@2229:     ///
deba@2229:     /// Assign operator. The generated sequence will be identical to
deba@2229:     /// the other sequence.
deba@2229:     Random& operator=(const Random& other) {
deba@2229:       if (&other != this) {
deba@2229:         std::copy(other.state, other.state + length, state);
deba@2229:         current = state + (other.current - other.state);
deba@2229:       }
deba@2229:       return *this;
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns an unsigned random number
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$.
deba@2229:     unsigned long getUnsigned() {
deba@2229:       if (current == state) fillState();
deba@2229:       --current;
deba@2229:       unsigned long rnd = *current;
deba@2229:       rnd ^= (rnd >> 29) & 0x5555555555555555ul;
deba@2229:       rnd ^= (rnd << 17) & 0x71D67FFFEDA60000ul;
deba@2229:       rnd ^= (rnd << 37) & 0xFFF7EEE000000000ul;
deba@2229:       rnd ^= (rnd >> 43);
deba@2229:       return rnd;
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns a random number
deba@2229:     ///
deba@2229:     /// It returns an integer random number from the range 
deba@2229:     /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$.
deba@2229:     long getInt() {
deba@2229:       return (long)getUnsigned();
deba@2229:     }
deba@2229:     
deba@2229:     /// \brief Returns an unsigned integer random number
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$.
deba@2229:     long getNatural() {
deba@2229:       return (long)(getUnsigned() >> 1);
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns a random bool
deba@2229:     ///
deba@2229:     /// It returns a random bool.
deba@2229:     bool getBool() {
deba@2229:       return (bool)(getUnsigned() & 1);
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns a real random number
deba@2229:     ///
deba@2229:     /// It returns a real random number from the range 
deba@2229:     /// \f$ [0, 1) \f$.
deba@2229:     double getReal() {
deba@2229:       return std::ldexp((double)(getUnsigned() >> 11), -53);
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns a real random number
deba@2229:     ///
deba@2229:     /// It returns a real random number from the range 
deba@2229:     /// \f$ [0, 1) \f$. This function is identical to the \c getReal().
deba@2229:     double getPrecReal() {
deba@2229:       return getReal();
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns an unsigned random number from a given range
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, n - 1\} \f$.
deba@2229:     unsigned long getUnsigned(unsigned long n) {
deba@2229:       unsigned long mask = n - 1, rnd;
deba@2229:       mask |= mask >> 1;
deba@2229:       mask |= mask >> 2;
deba@2229:       mask |= mask >> 4;
deba@2229:       mask |= mask >> 8;
deba@2229:       mask |= mask >> 16;
deba@2229:       mask |= mask >> 32;
deba@2229:       do rnd = getUnsigned() & mask; while (rnd >= n);
deba@2229:       return rnd;
deba@2229:     }
deba@2229: 
deba@2229:     /// \brief Returns a random number from a given range
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, n - 1\} \f$.
deba@2229:     long getInt(long n) {
deba@2229:       long mask = n - 1, rnd;
deba@2229:       mask |= mask >> 1;
deba@2229:       mask |= mask >> 2;
deba@2229:       mask |= mask >> 4;
deba@2229:       mask |= mask >> 8;
deba@2229:       mask |= mask >> 16;
deba@2229:       mask |= mask >> 32;
deba@2229:       do rnd = getUnsigned() & mask; while (rnd >= n);
deba@2229:       return rnd;
deba@2229:     }
deba@2229: 
deba@2229:   private:
deba@2229: 
deba@2229:     void initState(unsigned long seed) {
deba@2229: 
deba@2229:       static const unsigned long mul = 0x5851F42D4C957F2Dul;
deba@2229: 
deba@2229:       current = state; 
deba@2229: 
deba@2229:       unsigned long *curr = state + length - 1;
deba@2229:       curr[0] = seed; --curr;
deba@2229:       for (int i = 1; i < length; ++i) {
deba@2229:         curr[0] = (mul * ( curr[1] ^ (curr[1] >> 62) ) + i);
deba@2229:         --curr;
deba@2229:       }
deba@2229:     }
deba@2229:   
deba@2229:     void fillState() {
deba@2229:       static const unsigned long mask[2] = { 0x0ul, 0xB5026F5AA96619E9ul };
deba@2229:       static const unsigned long loMask = (1ul << 31) - 1;
deba@2229:       static const unsigned long hiMask = ~loMask;
deba@2229: 
deba@2229:       current = state + length; 
deba@2229: 
deba@2229:       register unsigned long *curr = state + length - 1;
deba@2229:       register int num;
deba@2229:       
deba@2229:       num = length - shift;
deba@2229:       while (num--) {
deba@2229:         curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
deba@2229:           curr[- shift] ^ mask[curr[-1] & 1ul];
deba@2229:         --curr;
deba@2229:       }
deba@2229:       num = shift - 1;
deba@2229:       while (num--) {
deba@2229:         curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
deba@2229:           curr[length - shift] ^ mask[curr[-1] & 1ul];
deba@2229:         --curr;
deba@2229:       }
deba@2229:       curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
deba@2229:         curr[length - shift] ^ mask[curr[length - 1] & 1ul];
deba@2229: 
deba@2229:     }
deba@2229:   
deba@2229:     unsigned long *current;
deba@2229:     unsigned long state[length];
deba@2229:     
deba@2229:   };
deba@2229: 
deba@2229: #else
deba@2229: 
deba@2230:   /// \ingroup misc
deba@2230:   ///
deba@2229:   /// \brief Mersenne Twister random number generator
deba@2229:   ///
deba@2229:   /// The Mersenne Twister is a twisted generalized feedback
deba@2229:   /// shift-register generator of Matsumoto and Nishimura. There is
deba@2229:   /// two different implementation in the Lemon library, one for
deba@2229:   /// 32-bit processors and one for 64-bit processors. The period of
deba@2229:   /// the generated sequence is \f$ 2^{19937} - 1 \f$, the generated
deba@2229:   /// sequence of 32-bit random numbers is equi-distributed in 623
deba@2229:   /// dimensions. The time performance of this generator is comparable
deba@2229:   /// to the commonly used generators.
deba@2229:   class Random {
deba@2229:   public:
deba@2229: 
deba@2229:     static const unsigned long min = 0ul;
deba@2229:     static const unsigned long max = ~0ul;
deba@2229:   
deba@2229:     /// \brief Constructor
deba@2229:     ///
deba@2229:     /// Constructor with time dependent seeding.
deba@2229:     Random() {}
deba@2229: 
deba@2229:     /// \brief Constructor
deba@2229:     ///
deba@2229:     /// Constructor
deba@2229:     Random(unsigned long seed) {}
deba@2229: 
deba@2229:     /// \brief Copy constructor
deba@2229:     ///
deba@2229:     /// Copy constructor. The generated sequence will be identical to
deba@2229:     /// the other sequence.
deba@2229:     Random(const Random& other) {}
deba@2229: 
deba@2229:     /// \brief Assign operator
deba@2229:     ///
deba@2229:     /// Assign operator. The generated sequence will be identical to
deba@2229:     /// the other sequence.
deba@2229:     Random& operator=(const Random& other) { return *this; }
deba@2229: 
deba@2229:     /// \brief Returns an unsigned random number
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$ for 64-bit processors and from
deba@2229:     /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$ for 32-bit processors.
deba@2229:     unsigned long getUnsigned() { return 0ul; }
deba@2229: 
deba@2229:     /// \brief Returns a random number
deba@2229:     ///
deba@2229:     /// It returns an integer random number from the range 
deba@2229:     /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$ for 64-bit processors and from
deba@2229:     /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$ for 32-bit processors.
deba@2229:     long getInt() { return 0l; }
deba@2229:     
deba@2229:     /// \brief Returns an unsigned integer random number
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$ for 64-bit processors and
deba@2229:     /// from \f$ \{0, 1, \dots, 2^{31} - 1\} \f$ for 32-bit processors.
deba@2229:     long getNatural() { return 0l; }
deba@2229: 
deba@2229:     /// \brief Returns a random bool
deba@2229:     ///
deba@2229:     /// It returns a random bool.
deba@2229:     bool getBool() { return false; }
deba@2229: 
deba@2229:     /// \brief Returns a real random number
deba@2229:     ///
deba@2229:     /// It returns a real random number from the range 
deba@2229:     /// \f$ [0, 1) \f$. For 32-bit processors the generated random
deba@2229:     /// number will just have 32 significant bits.
deba@2229:     double getReal() { return 0.0; }
deba@2229: 
deba@2229:     /// \brief Returns a real random number
deba@2229:     ///
deba@2229:     /// It returns a real random number from the range 
deba@2229:     /// \f$ [0, 1) \f$. This function returns random numbers with 53
deba@2229:     /// significant bits for 32-bit processors. For 64-bit processors
deba@2229:     /// it is identical to the \c getReal().
deba@2229:     double getPrecReal() { return 0.0; }
deba@2229: 
deba@2229:     /// \brief Returns an unsigned random number from a given range
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, n - 1\} \f$.
deba@2229:     unsigned long getUnsigned(unsigned long n) { return 0; }
deba@2229: 
deba@2229:     /// \brief Returns a random number from a given range
deba@2229:     ///
deba@2229:     /// It returns an unsigned integer random number from the range 
deba@2229:     /// \f$ \{0, 1, \dots, n - 1\} \f$.
deba@2229:     long getInt(long n) { return 0; }
deba@2229: 
deba@2229:   };
deba@2229: 
deba@2229: #endif
deba@2229: 
deba@2229:   /// \brief Global random number generator instance
deba@2229:   ///
deba@2229:   /// A global mersenne twister random number generator instance
deba@2229:   extern Random rnd;
deba@2229: 
deba@2229: }
deba@2229: 
deba@2229: #endif