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