lemon/random.h
author deba
Thu, 15 Feb 2007 13:06:23 +0000
changeset 2362 eb37b9774ef6
parent 2355 ac0d843b8873
child 2372 7fcc0179fb21
permissions -rw-r--r--
Small changes
     1 /* -*- C++ -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library
     4  *
     5  * Copyright (C) 2003-2006
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_RANDOM_H
    20 #define LEMON_RANDOM_H
    21 
    22 #include <algorithm>
    23 #include <iterator>
    24 #include <vector>
    25 
    26 #include <ctime>
    27 #include <cmath>
    28 
    29 ///\ingroup misc
    30 ///\file
    31 ///\brief Mersenne Twister random number generator
    32 ///
    33 ///\author Balazs Dezso
    34 
    35 namespace lemon {
    36 
    37   namespace _random_bits {
    38     
    39     template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
    40     struct RandomTraits {};
    41 
    42     template <typename _Word>
    43     struct RandomTraits<_Word, 32> {
    44 
    45       typedef _Word Word;
    46       static const int bits = 32;
    47 
    48       static const int length = 624;
    49       static const int shift = 397;
    50       
    51       static const Word mul = 0x6c078965u;
    52       static const Word arrayInit = 0x012BD6AAu;
    53       static const Word arrayMul1 = 0x0019660Du;
    54       static const Word arrayMul2 = 0x5D588B65u;
    55 
    56       static const Word mask = 0x9908B0DFu;
    57       static const Word loMask = (1u << 31) - 1;
    58       static const Word hiMask = ~loMask;
    59 
    60 
    61       static Word tempering(Word rnd) {
    62         rnd ^= (rnd >> 11);
    63         rnd ^= (rnd << 7) & 0x9D2C5680u;
    64         rnd ^= (rnd << 15) & 0xEFC60000u;
    65         rnd ^= (rnd >> 18);
    66         return rnd;
    67       }
    68 
    69     };
    70 
    71     template <typename _Word>
    72     struct RandomTraits<_Word, 64> {
    73 
    74       typedef _Word Word;
    75       static const int bits = 64;
    76 
    77       static const int length = 312;
    78       static const int shift = 156;
    79 
    80       static const Word mul = (Word)0x5851F42Du << 32 | (Word)0x4C957F2Du;
    81       static const Word arrayInit = (Word)0x00000000u << 32 |(Word)0x012BD6AAu;
    82       static const Word arrayMul1 = (Word)0x369DEA0Fu << 32 |(Word)0x31A53F85u;
    83       static const Word arrayMul2 = (Word)0x27BB2EE6u << 32 |(Word)0x87B0B0FDu;
    84 
    85       static const Word mask = (Word)0xB5026F5Au << 32 | (Word)0xA96619E9u;
    86       static const Word loMask = ((Word)1u << 31) - 1;
    87       static const Word hiMask = ~loMask;
    88 
    89       static Word tempering(Word rnd) {
    90         rnd ^= (rnd >> 29) & ((Word)0x55555555u << 32 | (Word)0x55555555u);
    91         rnd ^= (rnd << 17) & ((Word)0x71D67FFFu << 32 | (Word)0xEDA60000u);
    92         rnd ^= (rnd << 37) & ((Word)0xFFF7EEE0u << 32 | (Word)0x00000000u);
    93         rnd ^= (rnd >> 43);
    94         return rnd;
    95       }
    96 
    97     };
    98 
    99     template <typename _Word>
   100     class RandomCore {
   101     public:
   102 
   103       typedef _Word Word;
   104 
   105     private:
   106 
   107       static const int bits = RandomTraits<Word>::bits;
   108 
   109       static const int length = RandomTraits<Word>::length;
   110       static const int shift = RandomTraits<Word>::shift;
   111 
   112     public:
   113 
   114       void initState() {
   115         static const Word seedArray[4] = {
   116           0x12345u, 0x23456u, 0x34567u, 0x45678u
   117         };
   118     
   119         initState(seedArray, seedArray + 4);
   120       }
   121 
   122       void initState(Word seed) {
   123 
   124         static const Word mul = RandomTraits<Word>::mul;
   125 
   126         current = state; 
   127 
   128         Word *curr = state + length - 1;
   129         curr[0] = seed; --curr;
   130         for (int i = 1; i < length; ++i) {
   131           curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i);
   132           --curr;
   133         }
   134       }
   135 
   136       template <typename Iterator>
   137       void initState(Iterator begin, Iterator end) {
   138 
   139         static const Word init = RandomTraits<Word>::arrayInit;
   140         static const Word mul1 = RandomTraits<Word>::arrayMul1;
   141         static const Word mul2 = RandomTraits<Word>::arrayMul2;
   142 
   143 
   144         Word *curr = state + length - 1; --curr;
   145         Iterator it = begin; int cnt = 0;
   146         int num;
   147 
   148         initState(init);
   149 
   150         num = length > end - begin ? length : end - begin;
   151         while (num--) {
   152           curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1)) 
   153             + *it + cnt;
   154           ++it; ++cnt;
   155           if (it == end) {
   156             it = begin; cnt = 0;
   157           }
   158           if (curr == state) {
   159             curr = state + length - 1; curr[0] = state[0];
   160           }
   161           --curr;
   162         }
   163 
   164         num = length - 1; cnt = length - (curr - state) - 1;
   165         while (num--) {
   166           curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2))
   167             - cnt;
   168           --curr; ++cnt;
   169           if (curr == state) {
   170             curr = state + length - 1; curr[0] = state[0]; --curr;
   171             cnt = 1;
   172           }
   173         }
   174         
   175         state[length - 1] = (Word)1 << (bits - 1);
   176       }
   177       
   178       void copyState(const RandomCore& other) {
   179         std::copy(other.state, other.state + length, state);
   180         current = state + (other.current - other.state);
   181       }
   182 
   183       Word operator()() {
   184         if (current == state) fillState();
   185         --current;
   186         Word rnd = *current;
   187         return RandomTraits<Word>::tempering(rnd);
   188       }
   189 
   190     private:
   191 
   192   
   193       void fillState() {
   194         static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
   195         static const Word loMask = RandomTraits<Word>::loMask;
   196         static const Word hiMask = RandomTraits<Word>::hiMask;
   197 
   198         current = state + length; 
   199 
   200         register Word *curr = state + length - 1;
   201         register long num;
   202       
   203         num = length - shift;
   204         while (num--) {
   205           curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   206             curr[- shift] ^ mask[curr[-1] & 1ul];
   207           --curr;
   208         }
   209         num = shift - 1;
   210         while (num--) {
   211           curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   212             curr[length - shift] ^ mask[curr[-1] & 1ul];
   213           --curr;
   214         }
   215         curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
   216           curr[length - shift] ^ mask[curr[length - 1] & 1ul];
   217 
   218       }
   219 
   220   
   221       Word *current;
   222       Word state[length];
   223       
   224     };
   225 
   226 
   227     template <typename Result, 
   228               int shift = (std::numeric_limits<Result>::digits + 1) / 2>
   229     struct Masker {
   230       static Result mask(const Result& result) {
   231         return Masker<Result, (shift + 1) / 2>::
   232           mask((Result)(result | (result >> shift)));
   233       }
   234     };
   235     
   236     template <typename Result>
   237     struct Masker<Result, 1> {
   238       static Result mask(const Result& result) {
   239         return (Result)(result | (result >> 1));
   240       }
   241     };
   242 
   243     template <typename Result, typename Word, 
   244               int rest = std::numeric_limits<Result>::digits, int shift = 0, 
   245               bool last = rest <= std::numeric_limits<Word>::digits>
   246     struct IntConversion {
   247       static const int bits = std::numeric_limits<Word>::digits;
   248     
   249       static Result convert(RandomCore<Word>& rnd) {
   250         return (Result)(rnd() >> (bits - rest)) << shift;
   251       }
   252       
   253     }; 
   254 
   255     template <typename Result, typename Word, int rest, int shift> 
   256     struct IntConversion<Result, Word, rest, shift, false> {
   257       static const int bits = std::numeric_limits<Word>::digits;
   258 
   259       static Result convert(RandomCore<Word>& rnd) {
   260         return ((Result)rnd() << shift) | 
   261           IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);
   262       }
   263     };
   264 
   265 
   266     template <typename Result, typename Word,
   267               bool one_word = std::numeric_limits<Word>::digits < 
   268                               std::numeric_limits<Result>::digits>
   269     struct Mapping {
   270       static Result map(RandomCore<Word>& rnd, const Result& bound) {
   271         Word max = (Word)(bound - 1);
   272         Result mask = Masker<Result>::mask(bound - 1);
   273         Result num;
   274         do {
   275           num = IntConversion<Result, Word>::convert(rnd) & mask; 
   276         } while (num > max);
   277         return num;
   278       }
   279     };
   280 
   281     template <typename Result, typename Word>
   282     struct Mapping<Result, Word, false> {
   283       static Result map(RandomCore<Word>& rnd, const Result& bound) {
   284         Word max = (Word)(bound - 1);
   285         Word mask = Masker<Word, (std::numeric_limits<Result>::digits + 1) / 2>
   286           ::mask(max);
   287         Word num;
   288         do {
   289           num = rnd() & mask;
   290         } while (num > max);
   291         return num;
   292       }
   293     };
   294 
   295     template <typename Result, int exp, bool pos = (exp >= 0)>
   296     struct ShiftMultiplier {
   297       static const Result multiplier() {
   298         Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
   299         res *= res;
   300         if ((exp & 1) == 1) res *= (Result)2.0;
   301         return res; 
   302       }
   303     };
   304 
   305     template <typename Result, int exp>
   306     struct ShiftMultiplier<Result, exp, false> {
   307       static const Result multiplier() {
   308         Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
   309         res *= res;
   310         if ((exp & 1) == 1) res *= (Result)0.5;
   311         return res; 
   312       }
   313     };
   314 
   315     template <typename Result>
   316     struct ShiftMultiplier<Result, 0, true> {
   317       static const Result multiplier() {
   318         return (Result)1.0; 
   319       }
   320     };
   321 
   322     template <typename Result>
   323     struct ShiftMultiplier<Result, -20, true> {
   324       static const Result multiplier() {
   325         return (Result)(1.0/1048576.0); 
   326       }
   327     };
   328     
   329     template <typename Result>
   330     struct ShiftMultiplier<Result, -32, true> {
   331       static const Result multiplier() {
   332         return (Result)(1.0/424967296.0); 
   333       }
   334     };
   335 
   336     template <typename Result>
   337     struct ShiftMultiplier<Result, -53, true> {
   338       static const Result multiplier() {
   339         return (Result)(1.0/9007199254740992.0); 
   340       }
   341     };
   342 
   343     template <typename Result>
   344     struct ShiftMultiplier<Result, -64, true> {
   345       static const Result multiplier() {
   346         return (Result)(1.0/18446744073709551616.0); 
   347       }
   348     };
   349 
   350     template <typename Result, int exp>
   351     struct Shifting {
   352       static Result shift(const Result& result) {
   353         return result * ShiftMultiplier<Result, exp>::multiplier();
   354       }
   355     };
   356 
   357     template <typename Result, typename Word,
   358               int rest = std::numeric_limits<Result>::digits, int shift = 0, 
   359               bool last = rest <= std::numeric_limits<Word>::digits>
   360     struct RealConversion{ 
   361       static const int bits = std::numeric_limits<Word>::digits;
   362 
   363       static Result convert(RandomCore<Word>& rnd) {
   364         return Shifting<Result, - shift - rest>::
   365           shift((Result)(rnd() >> (bits - rest)));
   366       }
   367     };
   368 
   369     template <typename Result, typename Word, int rest, int shift>
   370     struct RealConversion<Result, Word, rest, shift, false> { 
   371       static const int bits = std::numeric_limits<Word>::digits;
   372 
   373       static Result convert(RandomCore<Word>& rnd) {
   374         return Shifting<Result, - shift - bits>::shift((Result)rnd()) +
   375           RealConversion<Result, Word, rest-bits, shift + bits>::convert(rnd);
   376       }
   377     };
   378 
   379     template <typename Result, typename Word>
   380     struct Initializer {
   381 
   382       template <typename Iterator>
   383       static void init(RandomCore<Word>& rnd, Iterator begin, Iterator end) {
   384         std::vector<Word> ws;
   385         for (Iterator it = begin; it != end; ++it) {
   386           ws.push_back((Word)*it);
   387         }
   388         rnd.initState(ws.begin(), ws.end());
   389       }
   390 
   391       static void init(RandomCore<Word>& rnd, Result seed) {
   392         rnd.initState(seed);
   393       }
   394     };
   395 
   396     template <typename Word>
   397     struct BoolConversion {
   398       static bool convert(RandomCore<Word>& rnd) {
   399         return (rnd() & 1) == 1;
   400       }
   401     };
   402 
   403   }
   404 
   405   /// \ingroup misc
   406   ///
   407   /// \brief Mersenne Twister random number generator
   408   ///
   409   /// The Mersenne Twister is a twisted generalized feedback
   410   /// shift-register generator of Matsumoto and Nishimura. The period
   411   /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
   412   /// equi-distributed in 623 dimensions for 32-bit numbers. The time
   413   /// performance of this generator is comparable to the commonly used
   414   /// generators.
   415   ///
   416   /// This implementation is specialized for both 32-bit and 64-bit
   417   /// architectures. The generators differ sligthly in the
   418   /// initialization and generation phase so they produce two
   419   /// completly different sequences.
   420   ///
   421   /// The generator gives back random numbers of serveral types. To
   422   /// get a random number from a range of a floating point type you
   423   /// can use one form of the \c operator() or the \c real() member
   424   /// function. If you want to get random number from the {0, 1, ...,
   425   /// n-1} integer range use the \c operator[] or the \c integer()
   426   /// method. And to get random number from the whole range of an
   427   /// integer type you can use the argumentless \c integer() or \c
   428   /// uinteger() functions. After all you can get random bool with
   429   /// equal chance of true and false or given probability of true
   430   /// result with the \c boolean() member functions.
   431   ///
   432   ///\code
   433   /// // The commented code is identical to the other
   434   /// double a = rnd();                     // [0.0, 1.0)
   435   /// // double a = rnd.real();             // [0.0, 1.0)
   436   /// double b = rnd(100.0);                // [0.0, 100.0)
   437   /// // double b = rnd.real(100.0);        // [0.0, 100.0)
   438   /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)
   439   /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)
   440   /// int d = rnd[100000];                  // 0..99999
   441   /// // int d = rnd.integer(100000);       // 0..99999
   442   /// int e = rnd[6] + 1;                   // 1..6
   443   /// // int e = rnd.integer(1, 1 + 6);     // 1..6
   444   /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
   445   /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
   446   /// bool g = rnd.boolean();               // P(g = true) = 0.5
   447   /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
   448   ///\endcode
   449   ///
   450   /// The lemon provides a global instance of the random number
   451   /// generator which name is \ref lemon::rnd "rnd". Usually it is a
   452   /// good programming convenience to use this global generator to get
   453   /// random numbers.
   454   ///
   455   /// \author Balazs Dezso
   456   class Random {
   457   private:
   458 
   459     // architecture word
   460     typedef unsigned long Word;
   461     
   462     _random_bits::RandomCore<Word> core;
   463 
   464   public:
   465 
   466     /// \brief Constructor
   467     ///
   468     /// Constructor with constant seeding.
   469     Random() { core.initState(); }
   470 
   471     /// \brief Constructor
   472     ///
   473     /// Constructor with seed. The current number type will be converted
   474     /// to the architecture word type.
   475     template <typename Number>
   476     Random(Number seed) { 
   477       _random_bits::Initializer<Number, Word>::init(core, seed);
   478     }
   479 
   480     /// \brief Constructor
   481     ///
   482     /// Constructor with array seeding. The given range should contain
   483     /// any number type and the numbers will be converted to the
   484     /// architecture word type.
   485     template <typename Iterator>
   486     Random(Iterator begin, Iterator end) { 
   487       typedef typename std::iterator_traits<Iterator>::value_type Number;
   488       _random_bits::Initializer<Number, Word>::init(core, begin, end);
   489     }
   490 
   491     /// \brief Copy constructor
   492     ///
   493     /// Copy constructor. The generated sequence will be identical to
   494     /// the other sequence. It can be used to save the current state
   495     /// of the generator and later use it to generate the same
   496     /// sequence.
   497     Random(const Random& other) {
   498       core.copyState(other.core);
   499     }
   500 
   501     /// \brief Assign operator
   502     ///
   503     /// Assign operator. The generated sequence will be identical to
   504     /// the other sequence. It can be used to save the current state
   505     /// of the generator and later use it to generate the same
   506     /// sequence.
   507     Random& operator=(const Random& other) {
   508       if (&other != this) {
   509         core.copyState(other.core);
   510       }
   511       return *this;
   512     }
   513 
   514     /// \brief Returns a random real number from the range [0, 1)
   515     ///
   516     /// It returns a random real number from the range [0, 1). The
   517     /// default Number type is double.
   518     template <typename Number>
   519     Number real() {
   520       return _random_bits::RealConversion<Number, Word>::convert(core);
   521     }
   522 
   523     double real() {
   524       return real<double>();
   525     }
   526 
   527     /// \brief Returns a random real number the range [0, b)
   528     ///
   529     /// It returns a random real number from the range [0, b).
   530     template <typename Number>
   531     Number real(Number b) { 
   532       return real<Number>() * b; 
   533     }
   534 
   535     /// \brief Returns a random real number from the range [a, b)
   536     ///
   537     /// It returns a random real number from the range [a, b).
   538     template <typename Number>
   539     Number real(Number a, Number b) { 
   540       return real<Number>() * (b - a) + a; 
   541     }
   542 
   543     /// \brief Returns a random real number from the range [0, 1)
   544     ///
   545     /// It returns a random double from the range [0, 1).
   546     double operator()() {
   547       return real<double>();
   548     }
   549 
   550     /// \brief Returns a random real number from the range [0, b)
   551     ///
   552     /// It returns a random real number from the range [0, b).
   553     template <typename Number>
   554     Number operator()(Number b) { 
   555       return real<Number>() * b; 
   556     }
   557 
   558     /// \brief Returns a random real number from the range [a, b)
   559     ///
   560     /// It returns a random real number from the range [a, b).
   561     template <typename Number>
   562     Number operator()(Number a, Number b) { 
   563       return real<Number>() * (b - a) + a; 
   564     }
   565 
   566     /// \brief Returns a random integer from a range
   567     ///
   568     /// It returns a random integer from the range {0, 1, ..., b - 1}.
   569     template <typename Number>
   570     Number integer(Number b) {
   571       return _random_bits::Mapping<Number, Word>::map(core, b);
   572     }
   573 
   574     /// \brief Returns a random integer from a range
   575     ///
   576     /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
   577     template <typename Number>
   578     Number integer(Number a, Number b) {
   579       return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
   580     }
   581 
   582     /// \brief Returns a random integer from a range
   583     ///
   584     /// It returns a random integer from the range {0, 1, ..., b - 1}.
   585     template <typename Number>
   586     Number operator[](Number b) {
   587       return _random_bits::Mapping<Number, Word>::map(core, b);
   588     }
   589 
   590     /// \brief Returns a random non-negative integer
   591     ///
   592     /// It returns a random non-negative integer uniformly from the
   593     /// whole range of the current \c Number type.  The default result
   594     /// type of this function is unsigned int.
   595     template <typename Number>
   596     Number uinteger() {
   597       return _random_bits::IntConversion<Number, Word>::convert(core);
   598     }
   599 
   600     unsigned int uinteger() {
   601       return uinteger<unsigned int>();
   602     }
   603 
   604     /// \brief Returns a random integer
   605     ///
   606     /// It returns a random integer uniformly from the whole range of
   607     /// the current \c Number type. The default result type of this
   608     /// function is int.
   609     template <typename Number>
   610     Number integer() {
   611       static const int nb = std::numeric_limits<Number>::digits + 
   612         (std::numeric_limits<Number>::is_signed ? 1 : 0);
   613       return _random_bits::IntConversion<Number, Word, nb>::convert(core);
   614     }
   615 
   616     int integer() {
   617       return integer<int>();
   618     }
   619     
   620     /// \brief Returns a random bool
   621     ///
   622     /// It returns a random bool
   623     bool boolean() {
   624       return _random_bits::BoolConversion<Word>::convert(core);
   625     }
   626 
   627     ///\name Nonuniform distributions
   628     ///
   629     
   630     ///@{
   631     
   632     /// \brief Returns a random bool
   633     ///
   634     /// It returns a random bool with given probability of true result
   635     bool boolean(double p) {
   636       return operator()() < p;
   637     }
   638 
   639     /// Standard Gauss distribution
   640 
   641     /// Standard Gauss distribution.
   642     /// \note The Cartesian form of the Box-Muller
   643     /// transformation is used to generate a random normal distribution.
   644     /// \todo Consider using the "ziggurat" method instead.
   645     double gauss() 
   646     {
   647       double V1,V2,S;
   648       do {
   649 	V1=2*real<double>()-1;
   650 	V2=2*real<double>()-1;
   651 	S=V1*V1+V2*V2;
   652       } while(S>=1);
   653       return std::sqrt(-2*std::log(S)/S)*V1;
   654     }
   655     /// Gauss distribution with given standard deviation and mean 0
   656 
   657     /// \sa gauss()
   658     ///
   659     double gauss(double std_dev) 
   660     {
   661       return gauss()*std_dev;
   662     }
   663     /// Gauss distribution with given mean and standard deviation
   664 
   665     /// \sa gauss()
   666     ///
   667     double gauss(double mean,double std_dev)
   668     {
   669       return gauss()*std_dev+mean;
   670     }
   671 
   672     /// Exponential distribution with given mean
   673 
   674     /// This function generates an exponential distribution random number
   675     /// with mean <tt>1/lambda</tt>.
   676     ///
   677     double exponential(double lambda=1.0)
   678     {
   679       return -log(real<double>())/lambda;
   680     }
   681 
   682     ///@}
   683     
   684   };
   685 
   686 
   687   extern Random rnd;
   688 
   689 }
   690 
   691 #endif