lemon/random.h
author Balazs Dezso <deba@inf.elte.hu>
Sun, 17 Aug 2008 15:02:03 +0200
changeset 1199 15282595e6f4
parent 1164 04f57dad1b07
parent 1178 61fdd06833a6
permissions -rw-r--r--
Using Arc instead of ArcIt in Dfs (#32)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2009
     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 /*
    20  * This file contains the reimplemented version of the Mersenne Twister
    21  * Generator of Matsumoto and Nishimura.
    22  *
    23  * See the appropriate copyright notice below.
    24  *
    25  * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
    26  * All rights reserved.
    27  *
    28  * Redistribution and use in source and binary forms, with or without
    29  * modification, are permitted provided that the following conditions
    30  * are met:
    31  *
    32  * 1. Redistributions of source code must retain the above copyright
    33  *    notice, this list of conditions and the following disclaimer.
    34  *
    35  * 2. Redistributions in binary form must reproduce the above copyright
    36  *    notice, this list of conditions and the following disclaimer in the
    37  *    documentation and/or other materials provided with the distribution.
    38  *
    39  * 3. The names of its contributors may not be used to endorse or promote
    40  *    products derived from this software without specific prior written
    41  *    permission.
    42  *
    43  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
    44  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
    45  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
    46  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
    47  * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
    48  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
    49  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
    50  * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
    51  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
    52  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
    53  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
    54  * OF THE POSSIBILITY OF SUCH DAMAGE.
    55  *
    56  *
    57  * Any feedback is very welcome.
    58  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
    59  * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
    60  */
    61 
    62 #ifndef LEMON_RANDOM_H
    63 #define LEMON_RANDOM_H
    64 
    65 #include <lemon/config.h>
    66 
    67 #include <algorithm>
    68 #include <iterator>
    69 #include <vector>
    70 #include <limits>
    71 #include <fstream>
    72 
    73 #include <lemon/math.h>
    74 #include <lemon/dim2.h>
    75 
    76 #ifndef LEMON_WIN32
    77 #include <sys/time.h>
    78 #include <ctime>
    79 #include <sys/types.h>
    80 #include <unistd.h>
    81 #else
    82 #include <lemon/bits/windows.h>
    83 #endif
    84 
    85 ///\ingroup misc
    86 ///\file
    87 ///\brief Mersenne Twister random number generator
    88 
    89 namespace lemon {
    90 
    91   namespace _random_bits {
    92 
    93     template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
    94     struct RandomTraits {};
    95 
    96     template <typename _Word>
    97     struct RandomTraits<_Word, 32> {
    98 
    99       typedef _Word Word;
   100       static const int bits = 32;
   101 
   102       static const int length = 624;
   103       static const int shift = 397;
   104 
   105       static const Word mul = 0x6c078965u;
   106       static const Word arrayInit = 0x012BD6AAu;
   107       static const Word arrayMul1 = 0x0019660Du;
   108       static const Word arrayMul2 = 0x5D588B65u;
   109 
   110       static const Word mask = 0x9908B0DFu;
   111       static const Word loMask = (1u << 31) - 1;
   112       static const Word hiMask = ~loMask;
   113 
   114       static Word tempering(Word rnd) {
   115         rnd ^= (rnd >> 11);
   116         rnd ^= (rnd << 7) & 0x9D2C5680u;
   117         rnd ^= (rnd << 15) & 0xEFC60000u;
   118         rnd ^= (rnd >> 18);
   119         return rnd;
   120       }
   121 
   122     };
   123 
   124     template <typename _Word>
   125     struct RandomTraits<_Word, 64> {
   126 
   127       typedef _Word Word;
   128       static const int bits = 64;
   129 
   130       static const int length = 312;
   131       static const int shift = 156;
   132 
   133       static const Word mul = Word(0x5851F42Du) << 32 | Word(0x4C957F2Du);
   134       static const Word arrayInit = Word(0x00000000u) << 32 |Word(0x012BD6AAu);
   135       static const Word arrayMul1 = Word(0x369DEA0Fu) << 32 |Word(0x31A53F85u);
   136       static const Word arrayMul2 = Word(0x27BB2EE6u) << 32 |Word(0x87B0B0FDu);
   137 
   138       static const Word mask = Word(0xB5026F5Au) << 32 | Word(0xA96619E9u);
   139       static const Word loMask = (Word(1u) << 31) - 1;
   140       static const Word hiMask = ~loMask;
   141 
   142       static Word tempering(Word rnd) {
   143         rnd ^= (rnd >> 29) & (Word(0x55555555u) << 32 | Word(0x55555555u));
   144         rnd ^= (rnd << 17) & (Word(0x71D67FFFu) << 32 | Word(0xEDA60000u));
   145         rnd ^= (rnd << 37) & (Word(0xFFF7EEE0u) << 32 | Word(0x00000000u));
   146         rnd ^= (rnd >> 43);
   147         return rnd;
   148       }
   149 
   150     };
   151 
   152     template <typename _Word>
   153     class RandomCore {
   154     public:
   155 
   156       typedef _Word Word;
   157 
   158     private:
   159 
   160       static const int bits = RandomTraits<Word>::bits;
   161 
   162       static const int length = RandomTraits<Word>::length;
   163       static const int shift = RandomTraits<Word>::shift;
   164 
   165     public:
   166 
   167       void initState() {
   168         static const Word seedArray[4] = {
   169           0x12345u, 0x23456u, 0x34567u, 0x45678u
   170         };
   171 
   172         initState(seedArray, seedArray + 4);
   173       }
   174 
   175       void initState(Word seed) {
   176 
   177         static const Word mul = RandomTraits<Word>::mul;
   178 
   179         current = state;
   180 
   181         Word *curr = state + length - 1;
   182         curr[0] = seed; --curr;
   183         for (int i = 1; i < length; ++i) {
   184           curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i);
   185           --curr;
   186         }
   187       }
   188 
   189       template <typename Iterator>
   190       void initState(Iterator begin, Iterator end) {
   191 
   192         static const Word init = RandomTraits<Word>::arrayInit;
   193         static const Word mul1 = RandomTraits<Word>::arrayMul1;
   194         static const Word mul2 = RandomTraits<Word>::arrayMul2;
   195 
   196 
   197         Word *curr = state + length - 1; --curr;
   198         Iterator it = begin; int cnt = 0;
   199         int num;
   200 
   201         initState(init);
   202 
   203         num = static_cast<int>(length > end - begin ? length : end - begin);
   204         while (num--) {
   205           curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1))
   206             + *it + cnt;
   207           ++it; ++cnt;
   208           if (it == end) {
   209             it = begin; cnt = 0;
   210           }
   211           if (curr == state) {
   212             curr = state + length - 1; curr[0] = state[0];
   213           }
   214           --curr;
   215         }
   216 
   217         num = length - 1; cnt = static_cast<int>(length - (curr - state) - 1);
   218         while (num--) {
   219           curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2))
   220             - cnt;
   221           --curr; ++cnt;
   222           if (curr == state) {
   223             curr = state + length - 1; curr[0] = state[0]; --curr;
   224             cnt = 1;
   225           }
   226         }
   227 
   228         state[length - 1] = Word(1) << (bits - 1);
   229       }
   230 
   231       void copyState(const RandomCore& other) {
   232         std::copy(other.state, other.state + length, state);
   233         current = state + (other.current - other.state);
   234       }
   235 
   236       Word operator()() {
   237         if (current == state) fillState();
   238         --current;
   239         Word rnd = *current;
   240         return RandomTraits<Word>::tempering(rnd);
   241       }
   242 
   243     private:
   244 
   245       void fillState() {
   246         static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
   247         static const Word loMask = RandomTraits<Word>::loMask;
   248         static const Word hiMask = RandomTraits<Word>::hiMask;
   249 
   250         current = state + length;
   251 
   252         Word *curr = state + length - 1;
   253         long num;
   254 
   255         num = length - shift;
   256         while (num--) {
   257           curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   258             curr[- shift] ^ mask[curr[-1] & 1ul];
   259           --curr;
   260         }
   261         num = shift - 1;
   262         while (num--) {
   263           curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   264             curr[length - shift] ^ mask[curr[-1] & 1ul];
   265           --curr;
   266         }
   267         state[0] = (((state[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
   268           curr[length - shift] ^ mask[curr[length - 1] & 1ul];
   269 
   270       }
   271 
   272       Word *current;
   273       Word state[length];
   274 
   275     };
   276 
   277 
   278     template <typename Result,
   279               int shift = (std::numeric_limits<Result>::digits + 1) / 2>
   280     struct Masker {
   281       static Result mask(const Result& result) {
   282         return Masker<Result, (shift + 1) / 2>::
   283           mask(static_cast<Result>(result | (result >> shift)));
   284       }
   285     };
   286 
   287     template <typename Result>
   288     struct Masker<Result, 1> {
   289       static Result mask(const Result& result) {
   290         return static_cast<Result>(result | (result >> 1));
   291       }
   292     };
   293 
   294     template <typename Result, typename Word,
   295               int rest = std::numeric_limits<Result>::digits, int shift = 0,
   296               bool last = (rest <= std::numeric_limits<Word>::digits)>
   297     struct IntConversion {
   298       static const int bits = std::numeric_limits<Word>::digits;
   299 
   300       static Result convert(RandomCore<Word>& rnd) {
   301         return static_cast<Result>(rnd() >> (bits - rest)) << shift;
   302       }
   303 
   304     };
   305 
   306     template <typename Result, typename Word, int rest, int shift>
   307     struct IntConversion<Result, Word, rest, shift, false> {
   308       static const int bits = std::numeric_limits<Word>::digits;
   309 
   310       static Result convert(RandomCore<Word>& rnd) {
   311         return (static_cast<Result>(rnd()) << shift) |
   312           IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);
   313       }
   314     };
   315 
   316 
   317     template <typename Result, typename Word,
   318               bool one_word = (std::numeric_limits<Word>::digits <
   319                                std::numeric_limits<Result>::digits) >
   320     struct Mapping {
   321       static Result map(RandomCore<Word>& rnd, const Result& bound) {
   322         Word max = Word(bound - 1);
   323         Result mask = Masker<Result>::mask(bound - 1);
   324         Result num;
   325         do {
   326           num = IntConversion<Result, Word>::convert(rnd) & mask;
   327         } while (num > max);
   328         return num;
   329       }
   330     };
   331 
   332     template <typename Result, typename Word>
   333     struct Mapping<Result, Word, false> {
   334       static Result map(RandomCore<Word>& rnd, const Result& bound) {
   335         Word max = Word(bound - 1);
   336         Word mask = Masker<Word, (std::numeric_limits<Result>::digits + 1) / 2>
   337           ::mask(max);
   338         Word num;
   339         do {
   340           num = rnd() & mask;
   341         } while (num > max);
   342         return static_cast<Result>(num);
   343       }
   344     };
   345 
   346     template <typename Result, int exp>
   347     struct ShiftMultiplier {
   348       static const Result multiplier() {
   349         Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
   350         res *= res;
   351         if ((exp & 1) == 1) res *= static_cast<Result>(0.5);
   352         return res;
   353       }
   354     };
   355 
   356     template <typename Result>
   357     struct ShiftMultiplier<Result, 0> {
   358       static const Result multiplier() {
   359         return static_cast<Result>(1.0);
   360       }
   361     };
   362 
   363     template <typename Result>
   364     struct ShiftMultiplier<Result, 20> {
   365       static const Result multiplier() {
   366         return static_cast<Result>(1.0/1048576.0);
   367       }
   368     };
   369 
   370     template <typename Result>
   371     struct ShiftMultiplier<Result, 32> {
   372       static const Result multiplier() {
   373         return static_cast<Result>(1.0/4294967296.0);
   374       }
   375     };
   376 
   377     template <typename Result>
   378     struct ShiftMultiplier<Result, 53> {
   379       static const Result multiplier() {
   380         return static_cast<Result>(1.0/9007199254740992.0);
   381       }
   382     };
   383 
   384     template <typename Result>
   385     struct ShiftMultiplier<Result, 64> {
   386       static const Result multiplier() {
   387         return static_cast<Result>(1.0/18446744073709551616.0);
   388       }
   389     };
   390 
   391     template <typename Result, int exp>
   392     struct Shifting {
   393       static Result shift(const Result& result) {
   394         return result * ShiftMultiplier<Result, exp>::multiplier();
   395       }
   396     };
   397 
   398     template <typename Result, typename Word,
   399               int rest = std::numeric_limits<Result>::digits, int shift = 0,
   400               bool last = rest <= std::numeric_limits<Word>::digits>
   401     struct RealConversion{
   402       static const int bits = std::numeric_limits<Word>::digits;
   403 
   404       static Result convert(RandomCore<Word>& rnd) {
   405         return Shifting<Result, shift + rest>::
   406           shift(static_cast<Result>(rnd() >> (bits - rest)));
   407       }
   408     };
   409 
   410     template <typename Result, typename Word, int rest, int shift>
   411     struct RealConversion<Result, Word, rest, shift, false> {
   412       static const int bits = std::numeric_limits<Word>::digits;
   413 
   414       static Result convert(RandomCore<Word>& rnd) {
   415         return Shifting<Result, shift + bits>::
   416           shift(static_cast<Result>(rnd())) +
   417           RealConversion<Result, Word, rest-bits, shift + bits>::
   418           convert(rnd);
   419       }
   420     };
   421 
   422     template <typename Result, typename Word>
   423     struct Initializer {
   424 
   425       template <typename Iterator>
   426       static void init(RandomCore<Word>& rnd, Iterator begin, Iterator end) {
   427         std::vector<Word> ws;
   428         for (Iterator it = begin; it != end; ++it) {
   429           ws.push_back(Word(*it));
   430         }
   431         rnd.initState(ws.begin(), ws.end());
   432       }
   433 
   434       static void init(RandomCore<Word>& rnd, Result seed) {
   435         rnd.initState(seed);
   436       }
   437     };
   438 
   439     template <typename Word>
   440     struct BoolConversion {
   441       static bool convert(RandomCore<Word>& rnd) {
   442         return (rnd() & 1) == 1;
   443       }
   444     };
   445 
   446     template <typename Word>
   447     struct BoolProducer {
   448       Word buffer;
   449       int num;
   450 
   451       BoolProducer() : num(0) {}
   452 
   453       bool convert(RandomCore<Word>& rnd) {
   454         if (num == 0) {
   455           buffer = rnd();
   456           num = RandomTraits<Word>::bits;
   457         }
   458         bool r = (buffer & 1);
   459         buffer >>= 1;
   460         --num;
   461         return r;
   462       }
   463     };
   464 
   465     /// \ingroup misc
   466     ///
   467     /// \brief Mersenne Twister random number generator
   468     ///
   469     /// The Mersenne Twister is a twisted generalized feedback
   470     /// shift-register generator of Matsumoto and Nishimura. The period
   471     /// of this generator is \f$ 2^{19937} - 1\f$ and it is
   472     /// equi-distributed in 623 dimensions for 32-bit numbers. The time
   473     /// performance of this generator is comparable to the commonly used
   474     /// generators.
   475     ///
   476     /// This is a template implementation of both 32-bit and
   477     /// 64-bit architecture optimized versions. The generators differ
   478     /// sligthly in the initialization and generation phase so they
   479     /// produce two completly different sequences.
   480     ///
   481     /// \alert Do not use this class directly, but instead one of \ref
   482     /// Random, \ref Random32 or \ref Random64.
   483     ///
   484     /// The generator gives back random numbers of serveral types. To
   485     /// get a random number from a range of a floating point type, you
   486     /// can use one form of the \c operator() or the \c real() member
   487     /// function. If you want to get random number from the {0, 1, ...,
   488     /// n-1} integer range, use the \c operator[] or the \c integer()
   489     /// method. And to get random number from the whole range of an
   490     /// integer type, you can use the argumentless \c integer() or
   491     /// \c uinteger() functions. Finally, you can get random bool with
   492     /// equal chance of true and false or with given probability of true
   493     /// result using the \c boolean() member functions.
   494     ///
   495     /// Various non-uniform distributions are also supported: normal (Gauss),
   496     /// exponential, gamma, Poisson, etc.; and a few two-dimensional
   497     /// distributions, too.
   498     ///
   499     ///\code
   500     /// // The commented code is identical to the other
   501     /// double a = rnd();                     // [0.0, 1.0)
   502     /// // double a = rnd.real();             // [0.0, 1.0)
   503     /// double b = rnd(100.0);                // [0.0, 100.0)
   504     /// // double b = rnd.real(100.0);        // [0.0, 100.0)
   505     /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)
   506     /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)
   507     /// int d = rnd[100000];                  // 0..99999
   508     /// // int d = rnd.integer(100000);       // 0..99999
   509     /// int e = rnd[6] + 1;                   // 1..6
   510     /// // int e = rnd.integer(1, 1 + 6);     // 1..6
   511     /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
   512     /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
   513     /// bool g = rnd.boolean();               // P(g = true) = 0.5
   514     /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
   515     ///\endcode
   516     ///
   517     /// LEMON provides a global instance of the random number generator:
   518     /// \ref lemon::rnd "rnd". In most cases, it is a good practice
   519     /// to use this global generator to get random numbers.
   520     ///
   521     /// \sa \ref Random, \ref Random32 or \ref Random64.
   522     template<class Word>
   523     class Random {
   524     private:
   525 
   526       _random_bits::RandomCore<Word> core;
   527       _random_bits::BoolProducer<Word> bool_producer;
   528 
   529 
   530     public:
   531 
   532       ///\name Initialization
   533       ///
   534       /// @{
   535 
   536       /// \brief Default constructor
   537       ///
   538       /// Constructor with constant seeding.
   539       Random() { core.initState(); }
   540 
   541       /// \brief Constructor with seed
   542       ///
   543       /// Constructor with seed. The current number type will be converted
   544       /// to the architecture word type.
   545       template <typename Number>
   546       Random(Number seed) {
   547         _random_bits::Initializer<Number, Word>::init(core, seed);
   548       }
   549 
   550       /// \brief Constructor with array seeding
   551       ///
   552       /// Constructor with array seeding. The given range should contain
   553       /// any number type and the numbers will be converted to the
   554       /// architecture word type.
   555       template <typename Iterator>
   556       Random(Iterator begin, Iterator end) {
   557         typedef typename std::iterator_traits<Iterator>::value_type Number;
   558         _random_bits::Initializer<Number, Word>::init(core, begin, end);
   559       }
   560 
   561       /// \brief Copy constructor
   562       ///
   563       /// Copy constructor. The generated sequence will be identical to
   564       /// the other sequence. It can be used to save the current state
   565       /// of the generator and later use it to generate the same
   566       /// sequence.
   567       Random(const Random& other) {
   568         core.copyState(other.core);
   569       }
   570 
   571       /// \brief Assign operator
   572       ///
   573       /// Assign operator. The generated sequence will be identical to
   574       /// the other sequence. It can be used to save the current state
   575       /// of the generator and later use it to generate the same
   576       /// sequence.
   577       Random& operator=(const Random& other) {
   578         if (&other != this) {
   579           core.copyState(other.core);
   580         }
   581         return *this;
   582       }
   583 
   584       /// \brief Seeding random sequence
   585       ///
   586       /// Seeding the random sequence. The current number type will be
   587       /// converted to the architecture word type.
   588       template <typename Number>
   589       void seed(Number seed) {
   590         _random_bits::Initializer<Number, Word>::init(core, seed);
   591       }
   592 
   593       /// \brief Seeding random sequence
   594       ///
   595       /// Seeding the random sequence. The given range should contain
   596       /// any number type and the numbers will be converted to the
   597       /// architecture word type.
   598       template <typename Iterator>
   599       void seed(Iterator begin, Iterator end) {
   600         typedef typename std::iterator_traits<Iterator>::value_type Number;
   601         _random_bits::Initializer<Number, Word>::init(core, begin, end);
   602       }
   603 
   604       /// \brief Seeding from file or from process id and time
   605       ///
   606       /// By default, this function calls the \c seedFromFile() member
   607       /// function with the <tt>/dev/urandom</tt> file. If it does not success,
   608       /// it uses the \c seedFromTime().
   609       /// \return Currently always \c true.
   610       bool seed() {
   611 #ifndef LEMON_WIN32
   612         if (seedFromFile("/dev/urandom", 0)) return true;
   613 #endif
   614         if (seedFromTime()) return true;
   615         return false;
   616       }
   617 
   618       /// \brief Seeding from file
   619       ///
   620       /// Seeding the random sequence from file. The linux kernel has two
   621       /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which
   622       /// could give good seed values for pseudo random generators (The
   623       /// difference between two devices is that the <tt>random</tt> may
   624       /// block the reading operation while the kernel can give good
   625       /// source of randomness, while the <tt>urandom</tt> does not
   626       /// block the input, but it could give back bytes with worse
   627       /// entropy).
   628       /// \param file The source file
   629       /// \param offset The offset, from the file read.
   630       /// \return \c true when the seeding successes.
   631 #ifndef LEMON_WIN32
   632       bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
   633 #else
   634         bool seedFromFile(const std::string& file = "", int offset = 0)
   635 #endif
   636       {
   637         std::ifstream rs(file.c_str());
   638         const int size = 4;
   639         Word buf[size];
   640         if (offset != 0 && !rs.seekg(offset)) return false;
   641         if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false;
   642         seed(buf, buf + size);
   643         return true;
   644       }
   645 
   646       /// \brief Seeding from process id and time
   647       ///
   648       /// Seeding from process id and time. This function uses the
   649       /// current process id and the current time for initialize the
   650       /// random sequence.
   651       /// \return Currently always \c true.
   652       bool seedFromTime() {
   653 #ifndef LEMON_WIN32
   654         timeval tv;
   655         gettimeofday(&tv, 0);
   656         seed(getpid() + tv.tv_sec + tv.tv_usec);
   657 #else
   658         seed(bits::getWinRndSeed());
   659 #endif
   660         return true;
   661       }
   662 
   663       /// @}
   664 
   665       ///\name Uniform Distributions
   666       ///
   667       /// @{
   668 
   669       /// \brief Returns a random real number from the range [0, 1)
   670       ///
   671       /// It returns a random real number from the range [0, 1). The
   672       /// default Number type is \c double.
   673       template <typename Number>
   674       Number real() {
   675         return _random_bits::RealConversion<Number, Word>::convert(core);
   676       }
   677 
   678       double real() {
   679         return real<double>();
   680       }
   681 
   682       /// \brief Returns a random real number from the range [0, 1)
   683       ///
   684       /// It returns a random double from the range [0, 1).
   685       double operator()() {
   686         return real<double>();
   687       }
   688 
   689       /// \brief Returns a random real number from the range [0, b)
   690       ///
   691       /// It returns a random real number from the range [0, b).
   692       double operator()(double b) {
   693         return real<double>() * b;
   694       }
   695 
   696       /// \brief Returns a random real number from the range [a, b)
   697       ///
   698       /// It returns a random real number from the range [a, b).
   699       double operator()(double a, double b) {
   700         return real<double>() * (b - a) + a;
   701       }
   702 
   703       /// \brief Returns a random integer from a range
   704       ///
   705       /// It returns a random integer from the range {0, 1, ..., b - 1}.
   706       template <typename Number>
   707       Number integer(Number b) {
   708         return _random_bits::Mapping<Number, Word>::map(core, b);
   709       }
   710 
   711       /// \brief Returns a random integer from a range
   712       ///
   713       /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
   714       template <typename Number>
   715       Number integer(Number a, Number b) {
   716         return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
   717       }
   718 
   719       /// \brief Returns a random integer from a range
   720       ///
   721       /// It returns a random integer from the range {0, 1, ..., b - 1}.
   722       template <typename Number>
   723       Number operator[](Number b) {
   724         return _random_bits::Mapping<Number, Word>::map(core, b);
   725       }
   726 
   727       /// \brief Returns a random non-negative integer
   728       ///
   729       /// It returns a random non-negative integer uniformly from the
   730       /// whole range of the current \c Number type. The default result
   731       /// type of this function is <tt>unsigned int</tt>.
   732       template <typename Number>
   733       Number uinteger() {
   734         return _random_bits::IntConversion<Number, Word>::convert(core);
   735       }
   736 
   737       unsigned int uinteger() {
   738         return uinteger<unsigned int>();
   739       }
   740 
   741       /// \brief Returns a random integer
   742       ///
   743       /// It returns a random integer uniformly from the whole range of
   744       /// the current \c Number type. The default result type of this
   745       /// function is \c int.
   746       template <typename Number>
   747       Number integer() {
   748         static const int nb = std::numeric_limits<Number>::digits +
   749           (std::numeric_limits<Number>::is_signed ? 1 : 0);
   750         return _random_bits::IntConversion<Number, Word, nb>::convert(core);
   751       }
   752 
   753       int integer() {
   754         return integer<int>();
   755       }
   756 
   757       /// \brief Returns a random bool
   758       ///
   759       /// It returns a random bool. The generator holds a buffer for
   760       /// random bits. Every time when it become empty the generator makes
   761       /// a new random word and fill the buffer up.
   762       bool boolean() {
   763         return bool_producer.convert(core);
   764       }
   765 
   766       /// @}
   767 
   768       ///\name Non-uniform Distributions
   769       ///
   770       ///@{
   771 
   772       /// \brief Returns a random bool with given probability of true result.
   773       ///
   774       /// It returns a random bool with given probability of true result.
   775       bool boolean(double p) {
   776         return operator()() < p;
   777       }
   778 
   779       /// Standard normal (Gauss) distribution
   780 
   781       /// Standard normal (Gauss) distribution.
   782       /// \note The Cartesian form of the Box-Muller
   783       /// transformation is used to generate a random normal distribution.
   784       double gauss()
   785       {
   786         double V1,V2,S;
   787         do {
   788           V1=2*real<double>()-1;
   789           V2=2*real<double>()-1;
   790           S=V1*V1+V2*V2;
   791         } while(S>=1);
   792         return std::sqrt(-2*std::log(S)/S)*V1;
   793       }
   794       /// Normal (Gauss) distribution with given mean and standard deviation
   795 
   796       /// Normal (Gauss) distribution with given mean and standard deviation.
   797       /// \sa gauss()
   798       double gauss(double mean,double std_dev)
   799       {
   800         return gauss()*std_dev+mean;
   801       }
   802 
   803       /// Lognormal distribution
   804 
   805       /// Lognormal distribution. The parameters are the mean and the standard
   806       /// deviation of <tt>exp(X)</tt>.
   807       ///
   808       double lognormal(double n_mean,double n_std_dev)
   809       {
   810         return std::exp(gauss(n_mean,n_std_dev));
   811       }
   812       /// Lognormal distribution
   813 
   814       /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of
   815       /// the mean and the standard deviation of <tt>exp(X)</tt>.
   816       ///
   817       double lognormal(const std::pair<double,double> &params)
   818       {
   819         return std::exp(gauss(params.first,params.second));
   820       }
   821       /// Compute the lognormal parameters from mean and standard deviation
   822 
   823       /// This function computes the lognormal parameters from mean and
   824       /// standard deviation. The return value can direcly be passed to
   825       /// lognormal().
   826       std::pair<double,double> lognormalParamsFromMD(double mean,
   827                                                      double std_dev)
   828       {
   829         double fr=std_dev/mean;
   830         fr*=fr;
   831         double lg=std::log(1+fr);
   832         return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg));
   833       }
   834       /// Lognormal distribution with given mean and standard deviation
   835 
   836       /// Lognormal distribution with given mean and standard deviation.
   837       ///
   838       double lognormalMD(double mean,double std_dev)
   839       {
   840         return lognormal(lognormalParamsFromMD(mean,std_dev));
   841       }
   842 
   843       /// Exponential distribution with given mean
   844 
   845       /// This function generates an exponential distribution random number
   846       /// with mean <tt>1/lambda</tt>.
   847       ///
   848       double exponential(double lambda=1.0)
   849       {
   850         return -std::log(1.0-real<double>())/lambda;
   851       }
   852 
   853       /// Gamma distribution with given integer shape
   854 
   855       /// This function generates a gamma distribution random number.
   856       ///
   857       ///\param k shape parameter (<tt>k>0</tt> integer)
   858       double gamma(int k)
   859       {
   860         double s = 0;
   861         for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
   862         return s;
   863       }
   864 
   865       /// Gamma distribution with given shape and scale parameter
   866 
   867       /// This function generates a gamma distribution random number.
   868       ///
   869       ///\param k shape parameter (<tt>k>0</tt>)
   870       ///\param theta scale parameter
   871       ///
   872       double gamma(double k,double theta=1.0)
   873       {
   874         double xi,nu;
   875         const double delta = k-std::floor(k);
   876         const double v0=E/(E-delta);
   877         do {
   878           double V0=1.0-real<double>();
   879           double V1=1.0-real<double>();
   880           double V2=1.0-real<double>();
   881           if(V2<=v0)
   882             {
   883               xi=std::pow(V1,1.0/delta);
   884               nu=V0*std::pow(xi,delta-1.0);
   885             }
   886           else
   887             {
   888               xi=1.0-std::log(V1);
   889               nu=V0*std::exp(-xi);
   890             }
   891         } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
   892         return theta*(xi+gamma(int(std::floor(k))));
   893       }
   894 
   895       /// Weibull distribution
   896 
   897       /// This function generates a Weibull distribution random number.
   898       ///
   899       ///\param k shape parameter (<tt>k>0</tt>)
   900       ///\param lambda scale parameter (<tt>lambda>0</tt>)
   901       ///
   902       double weibull(double k,double lambda)
   903       {
   904         return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
   905       }
   906 
   907       /// Pareto distribution
   908 
   909       /// This function generates a Pareto distribution random number.
   910       ///
   911       ///\param k shape parameter (<tt>k>0</tt>)
   912       ///\param x_min location parameter (<tt>x_min>0</tt>)
   913       ///
   914       double pareto(double k,double x_min)
   915       {
   916         return exponential(gamma(k,1.0/x_min))+x_min;
   917       }
   918 
   919       /// Poisson distribution
   920 
   921       /// This function generates a Poisson distribution random number with
   922       /// parameter \c lambda.
   923       ///
   924       /// The probability mass function of this distribusion is
   925       /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
   926       /// \note The algorithm is taken from the book of Donald E. Knuth titled
   927       /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
   928       /// return value.
   929 
   930       int poisson(double lambda)
   931       {
   932         const double l = std::exp(-lambda);
   933         int k=0;
   934         double p = 1.0;
   935         do {
   936           k++;
   937           p*=real<double>();
   938         } while (p>=l);
   939         return k-1;
   940       }
   941 
   942       ///@}
   943 
   944       ///\name Two-Dimensional Distributions
   945       ///
   946       ///@{
   947 
   948       /// Uniform distribution on the full unit circle
   949 
   950       /// Uniform distribution on the full unit circle.
   951       ///
   952       dim2::Point<double> disc()
   953       {
   954         double V1,V2;
   955         do {
   956           V1=2*real<double>()-1;
   957           V2=2*real<double>()-1;
   958 
   959         } while(V1*V1+V2*V2>=1);
   960         return dim2::Point<double>(V1,V2);
   961       }
   962       /// A kind of two-dimensional normal (Gauss) distribution
   963 
   964       /// This function provides a turning symmetric two-dimensional distribution.
   965       /// Both coordinates are of standard normal distribution, but they are not
   966       /// independent.
   967       ///
   968       /// \note The coordinates are the two random variables provided by
   969       /// the Box-Muller method.
   970       dim2::Point<double> gauss2()
   971       {
   972         double V1,V2,S;
   973         do {
   974           V1=2*real<double>()-1;
   975           V2=2*real<double>()-1;
   976           S=V1*V1+V2*V2;
   977         } while(S>=1);
   978         double W=std::sqrt(-2*std::log(S)/S);
   979         return dim2::Point<double>(W*V1,W*V2);
   980       }
   981       /// A kind of two-dimensional exponential distribution
   982 
   983       /// This function provides a turning symmetric two-dimensional distribution.
   984       /// The x-coordinate is of conditionally exponential distribution
   985       /// with the condition that x is positive and y=0. If x is negative and
   986       /// y=0 then, -x is of exponential distribution. The same is true for the
   987       /// y-coordinate.
   988       dim2::Point<double> exponential2()
   989       {
   990         double V1,V2,S;
   991         do {
   992           V1=2*real<double>()-1;
   993           V2=2*real<double>()-1;
   994           S=V1*V1+V2*V2;
   995         } while(S>=1);
   996         double W=-std::log(S)/S;
   997         return dim2::Point<double>(W*V1,W*V2);
   998       }
   999 
  1000       ///@}
  1001     };
  1002 
  1003 
  1004   };
  1005 
  1006   /// \ingroup misc
  1007   ///
  1008   /// \brief Mersenne Twister random number generator
  1009   ///
  1010   /// This class implements either the 32-bit or the 64-bit version of
  1011   /// the Mersenne Twister random number generator algorithm
  1012   /// depending on the system architecture.
  1013   /// 
  1014   /// For the API description, see its base class
  1015   /// \ref _random_bits::Random.
  1016   ///
  1017   /// \sa \ref _random_bits::Random
  1018   typedef _random_bits::Random<unsigned long> Random;
  1019 
  1020   /// \ingroup misc
  1021   ///
  1022   /// \brief Mersenne Twister random number generator (32-bit version)
  1023   ///
  1024   /// This class implements the 32-bit version of the Mersenne Twister
  1025   /// random number generator algorithm. It is recommended to be used
  1026   /// when someone wants to make sure that the \e same pseudo random
  1027   /// sequence will be generated on every platfrom.
  1028   ///
  1029   /// For the API description, see its base class
  1030   /// \ref _random_bits::Random.
  1031   ///
  1032   /// \sa \ref _random_bits::Random
  1033   typedef _random_bits::Random<unsigned int> Random32;
  1034 
  1035   /// \ingroup misc
  1036   ///
  1037   /// \brief Mersenne Twister random number generator (64-bit version)
  1038   ///
  1039   /// This class implements the 64-bit version of the Mersenne Twister
  1040   /// random number generator algorithm. (Even though it runs
  1041   /// on 32-bit architectures, too.) It is recommended to be used when
  1042   /// someone wants to make sure that the \e same pseudo random sequence
  1043   /// will be generated on every platfrom.
  1044   ///
  1045   /// For the API description, see its base class
  1046   /// \ref _random_bits::Random.
  1047   ///
  1048   /// \sa \ref _random_bits::Random
  1049   typedef _random_bits::Random<unsigned long long> Random64;
  1050 
  1051   extern Random rnd;
  1052   
  1053 }
  1054 
  1055 #endif