COIN-OR::LEMON - Graph Library

Changeset 2242:16523135943d in lemon-0.x for lemon


Ignore:
Timestamp:
10/14/06 17:26:05 (13 years ago)
Author:
Balazs Dezso
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@2991
Message:

New random interface
Switching to the new interface

Location:
lemon
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • lemon/hypercube_graph.h

    r2229 r2242  
    261261    /// dim2::Point<double> base[DIM];
    262262    /// for (int k = 0; k < DIM; ++k) {
    263     ///   base[k].x = random.getReal();
    264     ///   base[k].y = random.getReal();
     263    ///   base[k].x = rnd();
     264    ///   base[k].y = rnd();
    265265    /// }
    266266    /// HyperCubeGraph::HyperMap<dim2::Point<double> >
  • lemon/random.h

    r2230 r2242  
    2121
    2222#include <algorithm>
     23#include <iterator>
     24#include <vector>
    2325
    2426#include <ctime>
     
    3436namespace lemon {
    3537
    36 #if WORD_BIT == 32
     38  namespace _random_bits {
     39   
     40    template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
     41    struct RandomTraits {};
     42
     43    template <typename _Word>
     44    struct RandomTraits<_Word, 32> {
     45
     46      typedef _Word Word;
     47      static const int bits = 32;
     48
     49      static const int length = 624;
     50      static const int shift = 397;
     51     
     52      static const Word mul = 0x6c078965u;
     53      static const Word arrayInit = 0x012BD6AAu;
     54      static const Word arrayMul1 = 0x0019660Du;
     55      static const Word arrayMul2 = 0x5D588B65u;
     56
     57      static const Word mask = 0x9908B0DFu;
     58      static const Word loMask = (1u << 31) - 1;
     59      static const Word hiMask = ~loMask;
     60
     61
     62      static Word tempering(Word rnd) {
     63        rnd ^= (rnd >> 11);
     64        rnd ^= (rnd << 7) & 0x9D2C5680u;
     65        rnd ^= (rnd << 15) & 0xEFC60000u;
     66        rnd ^= (rnd >> 18);
     67        return rnd;
     68      }
     69
     70    };
     71
     72    template <typename _Word>
     73    struct RandomTraits<_Word, 64> {
     74
     75      typedef _Word Word;
     76      static const int bits = 64;
     77
     78      static const int length = 312;
     79      static const int shift = 156;
     80
     81      static const Word mul = (Word)0x5851F42Du << 32 | (Word)0x4C957F2Du;
     82      static const Word arrayInit = (Word)0x00000000u << 32 |(Word)0x012BD6AAu;
     83      static const Word arrayMul1 = (Word)0x369DEA0Fu << 32 |(Word)0x31A53F85u;
     84      static const Word arrayMul2 = (Word)0x27BB2EE6u << 32 |(Word)0x87B0B0FDu;
     85
     86      static const Word mask = (Word)0xB5026F5Au << 32 | (Word)0xA96619E9u;
     87      static const Word loMask = ((Word)1u << 31) - 1;
     88      static const Word hiMask = ~loMask;
     89
     90      static Word tempering(Word rnd) {
     91        rnd ^= (rnd >> 29) & ((Word)0x55555555u << 32 | (Word)0x55555555u);
     92        rnd ^= (rnd << 17) & ((Word)0x71D67FFFu << 32 | (Word)0xEDA60000u);
     93        rnd ^= (rnd << 37) & ((Word)0xFFF7EEE0u << 32 | (Word)0x00000000u);
     94        rnd ^= (rnd >> 43);
     95        return rnd;
     96      }
     97
     98    };
     99
     100    template <typename _Word>
     101    class RandomCore {
     102    public:
     103
     104      typedef _Word Word;
     105
     106    private:
     107
     108      static const int bits = RandomTraits<Word>::bits;
     109
     110      static const int length = RandomTraits<Word>::length;
     111      static const int shift = RandomTraits<Word>::shift;
     112
     113    public:
     114
     115      void initState() {
     116        static const Word seedArray[4] = {
     117          0x12345u, 0x23456u, 0x34567u, 0x45678u
     118        };
     119   
     120        initState(seedArray, seedArray + 4);
     121      }
     122
     123      void initState(Word seed) {
     124
     125        static const Word mul = RandomTraits<Word>::mul;
     126
     127        current = state;
     128
     129        Word *curr = state + length - 1;
     130        curr[0] = seed; --curr;
     131        for (int i = 1; i < length; ++i) {
     132          curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i);
     133          --curr;
     134        }
     135      }
     136
     137      template <typename Iterator>
     138      void initState(Iterator begin, Iterator end) {
     139
     140        static const Word init = RandomTraits<Word>::arrayInit;
     141        static const Word mul1 = RandomTraits<Word>::arrayMul1;
     142        static const Word mul2 = RandomTraits<Word>::arrayMul2;
     143
     144
     145        Word *curr = state + length - 1; --curr;
     146        Iterator it = begin; int cnt = 0;
     147        int num;
     148
     149        initState(init);
     150
     151        num = length > end - begin ? length : end - begin;
     152        while (num--) {
     153          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1))
     154            + *it + cnt;
     155          ++it; ++cnt;
     156          if (it == end) {
     157            it = begin; cnt = 0;
     158          }
     159          if (curr == state) {
     160            curr = state + length - 1; curr[0] = state[0];
     161          }
     162          --curr;
     163        }
     164
     165        num = length - 1; cnt = length - (curr - state) - 1;
     166        while (num--) {
     167          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2))
     168            - cnt;
     169          --curr; ++cnt;
     170          if (curr == state) {
     171            curr = state + length - 1; curr[0] = state[0]; --curr;
     172            cnt = 1;
     173          }
     174        }
     175       
     176        state[length - 1] = (Word)1 << (bits - 1);
     177      }
     178     
     179      void copyState(const RandomCore& other) {
     180        std::copy(other.state, other.state + length, state);
     181        current = state + (other.current - other.state);
     182      }
     183
     184      Word operator()() {
     185        if (current == state) fillState();
     186        --current;
     187        Word rnd = *current;
     188        return RandomTraits<Word>::tempering(rnd);
     189      }
     190
     191    private:
     192
     193 
     194      void fillState() {
     195        static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
     196        static const Word loMask = RandomTraits<Word>::loMask;
     197        static const Word hiMask = RandomTraits<Word>::hiMask;
     198
     199        current = state + length;
     200
     201        register Word *curr = state + length - 1;
     202        register long num;
     203     
     204        num = length - shift;
     205        while (num--) {
     206          curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
     207            curr[- shift] ^ mask[curr[-1] & 1ul];
     208          --curr;
     209        }
     210        num = shift - 1;
     211        while (num--) {
     212          curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
     213            curr[length - shift] ^ mask[curr[-1] & 1ul];
     214          --curr;
     215        }
     216        curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
     217          curr[length - shift] ^ mask[curr[length - 1] & 1ul];
     218
     219      }
     220
     221 
     222      Word *current;
     223      Word state[length];
     224     
     225    };
     226
     227
     228    template <typename Result,
     229              int shift = (std::numeric_limits<Result>::digits + 1) / 2>
     230    struct Masker {
     231      static Result mask(const Result& result) {
     232        return Masker<Result, (shift + 1) / 2>::
     233          mask((Result)(result | (result >> shift)));
     234      }
     235    };
     236   
     237    template <typename Result>
     238    struct Masker<Result, 1> {
     239      static Result mask(const Result& result) {
     240        return (Result)(result | (result >> 1));
     241      }
     242    };
     243
     244    template <typename Result, typename Word,
     245              int rest = std::numeric_limits<Result>::digits, int shift = 0,
     246              bool last = rest <= std::numeric_limits<Word>::digits>
     247    struct IntConversion {
     248      static const int bits = std::numeric_limits<Word>::digits;
     249   
     250      static Result convert(RandomCore<Word>& rnd) {
     251        return (Result)(rnd() >> (bits - rest)) << shift;
     252      }
     253     
     254    };
     255
     256    template <typename Result, typename Word, int rest, int shift>
     257    struct IntConversion<Result, Word, rest, shift, false> {
     258      static const int bits = std::numeric_limits<Word>::digits;
     259
     260      static Result convert(RandomCore<Word>& rnd) {
     261        return ((Result)rnd() << shift) |
     262          IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);
     263      }
     264    };
     265
     266
     267    template <typename Result, typename Word,
     268              bool one_word = std::numeric_limits<Word>::digits <
     269                              std::numeric_limits<Result>::digits>
     270    struct Mapping {
     271      static Result map(RandomCore<Word>& rnd, const Result& bound) {
     272        Word max = (Word)(bound - 1);
     273        Result mask = Masker<Result>::mask(bound - 1);
     274        Result num;
     275        do {
     276          num = IntConversion<Result, Word>::convert(rnd) & mask;
     277        } while (num > max);
     278        return num;
     279      }
     280    };
     281
     282    template <typename Result, typename Word>
     283    struct Mapping<Result, Word, false> {
     284      static Result map(RandomCore<Word>& rnd, const Result& bound) {
     285        Word max = (Word)(bound - 1);
     286        Word mask = Masker<Word, (std::numeric_limits<Result>::digits + 1) / 2>
     287          ::mask(max);
     288        Word num;
     289        do {
     290          num = rnd() & mask;
     291        } while (num > max);
     292        return num;
     293      }
     294    };
     295
     296    template <typename Result, int exp, bool pos = (exp >= 0)>
     297    struct ShiftMultiplier {
     298      static const Result multiplier() {
     299        Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
     300        res *= res;
     301        if ((exp & 1) == 1) res *= (Result)2.0;
     302        return res;
     303      }
     304    };
     305
     306    template <typename Result, int exp>
     307    struct ShiftMultiplier<Result, exp, false> {
     308      static const Result multiplier() {
     309        Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
     310        res *= res;
     311        if ((exp & 1) == 1) res *= (Result)0.5;
     312        return res;
     313      }
     314    };
     315
     316    template <typename Result>
     317    struct ShiftMultiplier<Result, 0, true> {
     318      static const Result multiplier() {
     319        return (Result)1.0;
     320      }
     321    };
     322
     323    template <typename Result>
     324    struct ShiftMultiplier<Result, -20, true> {
     325      static const Result multiplier() {
     326        return (Result)(1.0/1048576.0);
     327      }
     328    };
     329   
     330    template <typename Result>
     331    struct ShiftMultiplier<Result, -32, true> {
     332      static const Result multiplier() {
     333        return (Result)(1.0/424967296.0);
     334      }
     335    };
     336
     337    template <typename Result>
     338    struct ShiftMultiplier<Result, -53, true> {
     339      static const Result multiplier() {
     340        return (Result)(1.0/9007199254740992.0);
     341      }
     342    };
     343
     344    template <typename Result>
     345    struct ShiftMultiplier<Result, -64, true> {
     346      static const Result multiplier() {
     347        return (Result)(1.0/18446744073709551616.0);
     348      }
     349    };
     350
     351    template <typename Result, int exp>
     352    struct Shifting {
     353      static Result shift(const Result& result) {
     354        return result * ShiftMultiplier<Result, exp>::multiplier();
     355      }
     356    };
     357
     358    template <typename Result, typename Word,
     359              int rest = std::numeric_limits<Result>::digits, int shift = 0,
     360              bool last = rest <= std::numeric_limits<Word>::digits>
     361    struct RealConversion{
     362      static const int bits = std::numeric_limits<Word>::digits;
     363
     364      static Result convert(RandomCore<Word>& rnd) {
     365        return Shifting<Result, - shift - rest>::
     366          shift((Result)(rnd() >> (bits - rest)));
     367      }
     368    };
     369
     370    template <typename Result, typename Word, int rest, int shift>
     371    struct RealConversion<Result, Word, rest, shift, false> {
     372      static const int bits = std::numeric_limits<Word>::digits;
     373
     374      static Result convert(RandomCore<Word>& rnd) {
     375        return Shifting<Result, - shift - bits>::shift((Result)rnd()) +
     376          RealConversion<Result, Word, rest-bits, shift + bits>::convert(rnd);
     377      }
     378    };
     379
     380    template <typename Result, typename Word>
     381    struct Initializer {
     382
     383      template <typename Iterator>
     384      static void init(RandomCore<Word>& rnd, Iterator begin, Iterator end) {
     385        std::vector<Word> ws;
     386        for (Iterator it = begin; it != end; ++it) {
     387          ws.push_back((Word)*it);
     388        }
     389        rnd.initState(ws.begin(), ws.end());
     390      }
     391
     392      template <typename Iterator>
     393      static void init(RandomCore<Word>& rnd, Result seed) {
     394        rnd.initState(seed);
     395      }
     396    };
     397
     398    template <typename Word>
     399    struct BoolConversion {
     400      static bool convert(RandomCore<Word>& rnd) {
     401        return (rnd() & 1) == 1;
     402      }
     403    };
     404
     405  }
    37406
    38407  /// \ingroup misc
     
    41410  ///
    42411  /// The Mersenne Twister is a twisted generalized feedback
    43   /// shift-register generator of Matsumoto and Nishimura. The period of
    44   /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in
    45   /// 623 dimensions. The time performance of this generator is comparable
    46   /// to the commonly used generators.
     412  /// shift-register generator of Matsumoto and Nishimura. The period
     413  /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
     414  /// equi-distributed in 623 dimensions for 32-bit numbers. The time
     415  /// performance of this generator is comparable to the commonly used
     416  /// generators.
     417  ///
     418  /// This implementation is specialized for both 32-bit and 64-bit
     419  /// architectures. The generators differ sligthly in the
     420  /// initialization and generation phase so they produce two
     421  /// completly different sequences.
     422  ///
     423  /// The generator gives back random numbers of serveral types. To
     424  /// get a random number from a range of a floating point type you
     425  /// can use one form of the \c operator(). If you want to get random
     426  /// number from a the {0, 1, ..., n-1} integer range use the \c
     427  /// operator[]. And to get random number from the whole range of an
     428  /// integer type you can use the \c integer() or \c uinteger()
     429  /// functions. After all you can get random bool with equal chance
     430  /// or given probability with the \c boolean() member function.
     431  ///
     432  ///\code
     433  /// int a = rnd[100000];                  // 0..99999
     434  /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
     435  /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
     436  /// double d = rnd();                     // [0.0, 1.0)
     437  /// double e = rnd(100.0);                // [0.0, 100.0)
     438  /// double f = rnd(1.0, 2.0);             // [1.0, 2.0)
     439  /// bool g = rnd.boolean();               // P(g = true) = 0.5
     440  /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
     441  ///\endcode
     442  ///
     443  /// The lemon provides a global instance of the random number generator
     444  /// which name is \c rnd. Usually it is a good programming convenience
     445  /// to use this global generator to get random numbers.
    47446  ///
    48447  /// \author Balazs Dezso
    49448  class Random {
    50 
    51     static const int length = 624;
    52     static const int shift = 397;
     449  private:
     450
     451#if WORD_BIT == 32
     452    typedef unsigned int Word;
     453#elif WORD_BIT == 64
     454    typedef unsigned long Word;
     455#endif
     456   
     457    _random_bits::RandomCore<Word> core;
    53458
    54459  public:
    55460
    56     static const unsigned long min = 0ul;
    57     static const unsigned long max = ~0ul;
    58  
    59461    /// \brief Constructor
    60462    ///
    61     /// Constructor with time dependent seeding.
    62     Random() { initState(std::time(0)); }
     463    /// Constructor with constant seeding.
     464    Random() { core.initState(); }
    63465
    64466    /// \brief Constructor
    65467    ///
    66     /// Constructor
    67     Random(unsigned long seed) { initState(seed); }
     468    /// Constructor with seed. The current number type will be converted
     469    /// to the architecture word type.
     470    template <typename Number>
     471    Random(Number seed) {
     472      _random_bits::Initializer<Number, Word>::init(core, seed);
     473    }
     474
     475    /// \brief Constructor
     476    ///
     477    /// Constructor with array seeding. The given range should contain
     478    /// any number type and the numbers will be converted to the
     479    /// architecture word type.
     480    template <typename Iterator>
     481    Random(Iterator begin, Iterator end) {
     482      typedef typename std::iterator_traits<Iterator>::value_type Number;
     483      _random_bits::Initializer<Number, Word>::init(core, begin, end);
     484    }
    68485
    69486    /// \brief Copy constructor
     
    71488    /// Copy constructor. The generated sequence will be identical to
    72489    /// the other sequence.
    73     Random(const Random& other) {
    74       std::copy(other.state, other.state + length, state);
    75       current = state + (other.current - other.state);
     490    Random(const Random& other) {
     491      core.copyState(other.core);
    76492    }
    77493
     
    82498    Random& operator=(const Random& other) {
    83499      if (&other != this) {
    84         std::copy(other.state, other.state + length, state);
    85         current = state + (other.current - other.state);
     500        core.copyState(other.core);
    86501      }
    87502      return *this;
    88503    }
    89504
    90     /// \brief Returns an unsigned random number
    91     ///
    92     /// It returns an unsigned integer random number from the range
    93     /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$.
    94     unsigned long getUnsigned() {
    95       if (current == state) fillState();
    96       --current;
    97       unsigned long rnd = *current;
    98       rnd ^= (rnd >> 11);
    99       rnd ^= (rnd << 7) & 0x9D2C5680ul;
    100       rnd ^= (rnd << 15) & 0xEFC60000ul;
    101       rnd ^= (rnd >> 18);
    102       return rnd;
    103     }
    104 
    105     /// \brief Returns a random number
    106     ///
    107     /// It returns an integer random number from the range
    108     /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$.
    109     long getInt() {
    110       return (long)getUnsigned();
    111     }
    112    
    113     /// \brief Returns an unsigned integer random number
    114     ///
    115     /// It returns an unsigned integer random number from the range
    116     /// \f$ \{0, 1, \dots, 2^{31} - 1\} \f$.
    117     long getNatural() {
    118       return (long)(getUnsigned() >> 1);
    119     }
    120 
     505    /// \brief Returns a random real number
     506    ///
     507    /// It returns a random real number from the range [0, 1). The default
     508    /// result type of this function is double.
     509    template <typename Number>
     510    Number operator()() {
     511      return _random_bits::RealConversion<Number, Word>::convert(core);
     512    }
     513
     514    double operator()() {
     515      return operator()<double>();
     516    }
     517
     518    /// \brief Returns a random real number
     519    ///
     520    /// It returns a random real number from the range [0, b).
     521    template <typename Number>
     522    Number operator()(Number b) {
     523      return operator()<Number>() * b;
     524    }
     525
     526    /// \brief Returns a random real number
     527    ///
     528    /// It returns a random real number from the range [a, b).
     529    template <typename Number>
     530    Number operator()(Number a, Number b) {
     531      return operator()<Number>() * (b - a) + a;
     532    }
     533
     534    /// \brief Returns a random integer from a range
     535    ///
     536    /// It returns a random integer from the range {0, 1, ..., bound - 1}.
     537    template <typename Number>
     538    Number operator[](const Number& bound) {
     539      return _random_bits::Mapping<Number, Word>::map(core, bound);
     540    }
     541
     542    /// \brief Returns a random non-negative integer
     543    ///
     544    /// It returns a random non-negative integer uniformly from the
     545    /// whole range of the current \c Number type.  The default result
     546    /// type of this function is unsigned int.
     547    template <typename Number>
     548    Number uinteger() {
     549      return _random_bits::IntConversion<Number, Word>::convert(core);
     550    }
     551
     552    unsigned int uinteger() {
     553      return uinteger<unsigned int>();
     554    }
     555
     556    /// \brief Returns a random integer
     557    ///
     558    /// It returns a random integer uniformly from the whole range of
     559    /// the current \c Number type. The default result type of this
     560    /// function is int.
     561    template <typename Number>
     562    Number integer() {
     563      static const int nb = std::numeric_limits<Number>::digits +
     564        (std::numeric_limits<Number>::is_signed ? 1 : 0);
     565      return _random_bits::IntConversion<Number, Word, nb>::convert(core);
     566    }
     567
     568    int integer() {
     569      return integer<int>();
     570    }
     571   
    121572    /// \brief Returns a random bool
    122573    ///
    123     /// It returns a random bool.
    124     bool getBool() {
    125       return (bool)(getUnsigned() & 1);
    126     }
    127 
    128     /// \brief Returns a real random number
    129     ///
    130     /// It returns a real random number from the range
    131     /// \f$ [0, 1) \f$. The double will have 32 significant bits.
    132     double getReal() {
    133       return std::ldexp((double)getUnsigned(), -32);
    134     }
    135 
    136     /// \brief Returns a real random number
    137     ///
    138     /// It returns a real random number from the range
    139     /// \f$ [0, 1) \f$. The double will have 53 significant bits,
    140     /// but it is slower than the \c getReal().
    141     double getPrecReal() {
    142       return std::ldexp((double)(getUnsigned() >> 5), -27) +
    143         std::ldexp((double)(getUnsigned() >> 6), -53);
    144     }
    145 
    146     /// \brief Returns an unsigned random number from a given range
    147     ///
    148     /// It returns an unsigned integer random number from the range
    149     /// \f$ \{0, 1, \dots, n - 1\} \f$.
    150     unsigned long getUnsigned(unsigned long n) {
    151       unsigned long mask = n - 1, rnd;
    152       mask |= mask >> 1;
    153       mask |= mask >> 2;
    154       mask |= mask >> 4;
    155       mask |= mask >> 8;
    156       mask |= mask >> 16;
    157       do rnd = getUnsigned() & mask; while (rnd >= n);
    158       return rnd;
    159     }
    160 
    161     /// \brief Returns a random number from a given range
    162     ///
    163     /// It returns an unsigned integer random number from the range
    164     /// \f$ \{0, 1, \dots, n - 1\} \f$.
    165     long getInt(long n) {
    166       long mask = n - 1, rnd;
    167       mask |= mask >> 1;
    168       mask |= mask >> 2;
    169       mask |= mask >> 4;
    170       mask |= mask >> 8;
    171       mask |= mask >> 16;
    172       do rnd = getUnsigned() & mask; while (rnd >= n);
    173       return rnd;
    174     }
    175 
    176   private:
    177 
    178     void initState(unsigned long seed) {
    179       static const unsigned long mul = 0x6c078965ul;
    180 
    181       current = state;
    182 
    183       unsigned long *curr = state + length - 1;
    184       curr[0] = seed; --curr;
    185       for (int i = 1; i < length; ++i) {
    186         curr[0] = (mul * ( curr[1] ^ (curr[1] >> 30) ) + i);
    187         --curr;
    188       }
    189     }
    190  
    191     void fillState() {
    192       static const unsigned long mask[2] = { 0x0ul, 0x9908B0DFul };
    193       static const unsigned long loMask = (1ul << 31) - 1;
    194       static const unsigned long hiMask = ~loMask;
    195 
    196       current = state + length;
    197 
    198       register unsigned long *curr = state + length - 1;
    199       register long num;
    200      
    201       num = length - shift;
    202       while (num--) {
    203         curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
    204           curr[- shift] ^ mask[curr[-1] & 1ul];
    205         --curr;
    206       }
    207       num = shift - 1;
    208       while (num--) {
    209         curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
    210           curr[length - shift] ^ mask[curr[-1] & 1ul];
    211         --curr;
    212       }
    213       curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
    214         curr[length - shift] ^ mask[curr[length - 1] & 1ul];
    215 
    216     }
    217  
    218     unsigned long *current;
    219     unsigned long state[length];
     574    /// It returns a random bool
     575    bool boolean() {
     576      return _random_bits::BoolConversion<Word>::convert(core);
     577    }
     578
     579    /// \brief Returns a random bool
     580    ///
     581    /// It returns a random bool with given probability of true result
     582    bool boolean(double p) {
     583      return operator()() < p;
     584    }
    220585   
    221586  };
    222587
    223 #elif WORD_BIT == 64
    224 
    225   /// \ingroup misc
    226   ///
    227   /// \brief Mersenne Twister random number generator
    228   ///
    229   /// The Mersenne Twister is a twisted generalized feedback
    230   /// shift-register generator of Matsumoto and Nishimura. The period of
    231   /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in
    232   /// 311 dimensions. The time performance of this generator is comparable
    233   /// to the commonly used generators.
    234   class Random {
    235 
    236     static const int length = 312;
    237     static const int shift = 156;
    238 
    239   public:
    240 
    241     static const unsigned long min = 0ul;
    242     static const unsigned long max = ~0ul;
    243  
    244     /// \brief Constructor
    245     ///
    246     /// Constructor with time dependent seeding.
    247     Random() { initState(std::time(0)); }
    248 
    249     /// \brief Constructor
    250     ///
    251     /// Constructor
    252     Random(unsigned long seed) { initState(seed); }
    253 
    254     /// \brief Copy constructor
    255     ///
    256     /// Copy constructor. The generated sequence will be identical to
    257     /// the other sequence.
    258     Random(const Random& other) {
    259       std::copy(other.state, other.state + length, state);
    260       current = state + (other.current - other.state);
    261     }
    262 
    263     /// \brief Assign operator
    264     ///
    265     /// Assign operator. The generated sequence will be identical to
    266     /// the other sequence.
    267     Random& operator=(const Random& other) {
    268       if (&other != this) {
    269         std::copy(other.state, other.state + length, state);
    270         current = state + (other.current - other.state);
    271       }
    272       return *this;
    273     }
    274 
    275     /// \brief Returns an unsigned random number
    276     ///
    277     /// It returns an unsigned integer random number from the range
    278     /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$.
    279     unsigned long getUnsigned() {
    280       if (current == state) fillState();
    281       --current;
    282       unsigned long rnd = *current;
    283       rnd ^= (rnd >> 29) & 0x5555555555555555ul;
    284       rnd ^= (rnd << 17) & 0x71D67FFFEDA60000ul;
    285       rnd ^= (rnd << 37) & 0xFFF7EEE000000000ul;
    286       rnd ^= (rnd >> 43);
    287       return rnd;
    288     }
    289 
    290     /// \brief Returns a random number
    291     ///
    292     /// It returns an integer random number from the range
    293     /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$.
    294     long getInt() {
    295       return (long)getUnsigned();
    296     }
    297    
    298     /// \brief Returns an unsigned integer random number
    299     ///
    300     /// It returns an unsigned integer random number from the range
    301     /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$.
    302     long getNatural() {
    303       return (long)(getUnsigned() >> 1);
    304     }
    305 
    306     /// \brief Returns a random bool
    307     ///
    308     /// It returns a random bool.
    309     bool getBool() {
    310       return (bool)(getUnsigned() & 1);
    311     }
    312 
    313     /// \brief Returns a real random number
    314     ///
    315     /// It returns a real random number from the range
    316     /// \f$ [0, 1) \f$.
    317     double getReal() {
    318       return std::ldexp((double)(getUnsigned() >> 11), -53);
    319     }
    320 
    321     /// \brief Returns a real random number
    322     ///
    323     /// It returns a real random number from the range
    324     /// \f$ [0, 1) \f$. This function is identical to the \c getReal().
    325     double getPrecReal() {
    326       return getReal();
    327     }
    328 
    329     /// \brief Returns an unsigned random number from a given range
    330     ///
    331     /// It returns an unsigned integer random number from the range
    332     /// \f$ \{0, 1, \dots, n - 1\} \f$.
    333     unsigned long getUnsigned(unsigned long n) {
    334       unsigned long mask = n - 1, rnd;
    335       mask |= mask >> 1;
    336       mask |= mask >> 2;
    337       mask |= mask >> 4;
    338       mask |= mask >> 8;
    339       mask |= mask >> 16;
    340       mask |= mask >> 32;
    341       do rnd = getUnsigned() & mask; while (rnd >= n);
    342       return rnd;
    343     }
    344 
    345     /// \brief Returns a random number from a given range
    346     ///
    347     /// It returns an unsigned integer random number from the range
    348     /// \f$ \{0, 1, \dots, n - 1\} \f$.
    349     long getInt(long n) {
    350       long mask = n - 1, rnd;
    351       mask |= mask >> 1;
    352       mask |= mask >> 2;
    353       mask |= mask >> 4;
    354       mask |= mask >> 8;
    355       mask |= mask >> 16;
    356       mask |= mask >> 32;
    357       do rnd = getUnsigned() & mask; while (rnd >= n);
    358       return rnd;
    359     }
    360 
    361   private:
    362 
    363     void initState(unsigned long seed) {
    364 
    365       static const unsigned long mul = 0x5851F42D4C957F2Dul;
    366 
    367       current = state;
    368 
    369       unsigned long *curr = state + length - 1;
    370       curr[0] = seed; --curr;
    371       for (int i = 1; i < length; ++i) {
    372         curr[0] = (mul * ( curr[1] ^ (curr[1] >> 62) ) + i);
    373         --curr;
    374       }
    375     }
    376  
    377     void fillState() {
    378       static const unsigned long mask[2] = { 0x0ul, 0xB5026F5AA96619E9ul };
    379       static const unsigned long loMask = (1ul << 31) - 1;
    380       static const unsigned long hiMask = ~loMask;
    381 
    382       current = state + length;
    383 
    384       register unsigned long *curr = state + length - 1;
    385       register int num;
    386      
    387       num = length - shift;
    388       while (num--) {
    389         curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
    390           curr[- shift] ^ mask[curr[-1] & 1ul];
    391         --curr;
    392       }
    393       num = shift - 1;
    394       while (num--) {
    395         curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
    396           curr[length - shift] ^ mask[curr[-1] & 1ul];
    397         --curr;
    398       }
    399       curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
    400         curr[length - shift] ^ mask[curr[length - 1] & 1ul];
    401 
    402     }
    403  
    404     unsigned long *current;
    405     unsigned long state[length];
    406    
    407   };
    408 
    409 #else
    410 
    411   /// \ingroup misc
    412   ///
    413   /// \brief Mersenne Twister random number generator
    414   ///
    415   /// The Mersenne Twister is a twisted generalized feedback
    416   /// shift-register generator of Matsumoto and Nishimura. There is
    417   /// two different implementation in the Lemon library, one for
    418   /// 32-bit processors and one for 64-bit processors. The period of
    419   /// the generated sequence is \f$ 2^{19937} - 1 \f$, the generated
    420   /// sequence of 32-bit random numbers is equi-distributed in 623
    421   /// dimensions. The time performance of this generator is comparable
    422   /// to the commonly used generators.
    423   class Random {
    424   public:
    425 
    426     static const unsigned long min = 0ul;
    427     static const unsigned long max = ~0ul;
    428  
    429     /// \brief Constructor
    430     ///
    431     /// Constructor with time dependent seeding.
    432     Random() {}
    433 
    434     /// \brief Constructor
    435     ///
    436     /// Constructor
    437     Random(unsigned long seed) {}
    438 
    439     /// \brief Copy constructor
    440     ///
    441     /// Copy constructor. The generated sequence will be identical to
    442     /// the other sequence.
    443     Random(const Random& other) {}
    444 
    445     /// \brief Assign operator
    446     ///
    447     /// Assign operator. The generated sequence will be identical to
    448     /// the other sequence.
    449     Random& operator=(const Random& other) { return *this; }
    450 
    451     /// \brief Returns an unsigned random number
    452     ///
    453     /// It returns an unsigned integer random number from the range
    454     /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$ for 64-bit processors and from
    455     /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$ for 32-bit processors.
    456     unsigned long getUnsigned() { return 0ul; }
    457 
    458     /// \brief Returns a random number
    459     ///
    460     /// It returns an integer random number from the range
    461     /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$ for 64-bit processors and from
    462     /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$ for 32-bit processors.
    463     long getInt() { return 0l; }
    464    
    465     /// \brief Returns an unsigned integer random number
    466     ///
    467     /// It returns an unsigned integer random number from the range
    468     /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$ for 64-bit processors and
    469     /// from \f$ \{0, 1, \dots, 2^{31} - 1\} \f$ for 32-bit processors.
    470     long getNatural() { return 0l; }
    471 
    472     /// \brief Returns a random bool
    473     ///
    474     /// It returns a random bool.
    475     bool getBool() { return false; }
    476 
    477     /// \brief Returns a real random number
    478     ///
    479     /// It returns a real random number from the range
    480     /// \f$ [0, 1) \f$. For 32-bit processors the generated random
    481     /// number will just have 32 significant bits.
    482     double getReal() { return 0.0; }
    483 
    484     /// \brief Returns a real random number
    485     ///
    486     /// It returns a real random number from the range
    487     /// \f$ [0, 1) \f$. This function returns random numbers with 53
    488     /// significant bits for 32-bit processors. For 64-bit processors
    489     /// it is identical to the \c getReal().
    490     double getPrecReal() { return 0.0; }
    491 
    492     /// \brief Returns an unsigned random number from a given range
    493     ///
    494     /// It returns an unsigned integer random number from the range
    495     /// \f$ \{0, 1, \dots, n - 1\} \f$.
    496     unsigned long getUnsigned(unsigned long n) { return 0; }
    497 
    498     /// \brief Returns a random number from a given range
    499     ///
    500     /// It returns an unsigned integer random number from the range
    501     /// \f$ \{0, 1, \dots, n - 1\} \f$.
    502     long getInt(long n) { return 0; }
    503 
    504   };
    505 
    506 #endif
    507588
    508589  /// \brief Global random number generator instance
  • lemon/simann.h

    r2229 r2242  
    258258    bool accept() {
    259259      double cost_diff = simann->getCurrCost() - simann->getPrevCost();
    260       return (rnd.getReal() <= exp(-(cost_diff / temp)));
     260      return (rnd() <= exp(-(cost_diff / temp)));
    261261    }
    262262    /// \brief Destructor.
     
    366366      else {
    367367        double cost_diff = simann->getCurrCost() - simann->getPrevCost();
    368         return (rnd.getReal() <= exp(-(cost_diff / temp)));
     368        return (rnd() <= exp(-(cost_diff / temp)));
    369369      }
    370370    }
Note: See TracChangeset for help on using the changeset viewer.