alpar@209: /* -*- mode: C++; indent-tabs-mode: nil; -*- alpar@10: * alpar@209: * This file is a part of LEMON, a generic C++ optimization library. alpar@10: * alpar@440: * Copyright (C) 2003-2009 alpar@10: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport alpar@10: * (Egervary Research Group on Combinatorial Optimization, EGRES). alpar@10: * alpar@10: * Permission to use, modify and distribute this software is granted alpar@10: * provided that this copyright notice appears in all copies. For alpar@10: * precise terms see the accompanying LICENSE file. alpar@10: * alpar@10: * This software is provided "AS IS" with no warranty of any kind, alpar@10: * express or implied, and with no claim as to its suitability for any alpar@10: * purpose. alpar@10: * alpar@10: */ alpar@10: alpar@10: /* alpar@10: * This file contains the reimplemented version of the Mersenne Twister alpar@10: * Generator of Matsumoto and Nishimura. alpar@10: * alpar@10: * See the appropriate copyright notice below. alpar@209: * alpar@10: * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, alpar@209: * All rights reserved. alpar@10: * alpar@10: * Redistribution and use in source and binary forms, with or without alpar@10: * modification, are permitted provided that the following conditions alpar@10: * are met: alpar@10: * alpar@10: * 1. Redistributions of source code must retain the above copyright alpar@10: * notice, this list of conditions and the following disclaimer. alpar@10: * alpar@10: * 2. Redistributions in binary form must reproduce the above copyright alpar@10: * notice, this list of conditions and the following disclaimer in the alpar@10: * documentation and/or other materials provided with the distribution. alpar@10: * alpar@209: * 3. The names of its contributors may not be used to endorse or promote alpar@209: * products derived from this software without specific prior written alpar@10: * permission. alpar@10: * alpar@10: * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS alpar@10: * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT alpar@10: * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS alpar@10: * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE alpar@10: * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, alpar@10: * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES alpar@10: * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR alpar@10: * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) alpar@10: * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, alpar@10: * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) alpar@10: * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED alpar@10: * OF THE POSSIBILITY OF SUCH DAMAGE. alpar@10: * alpar@10: * alpar@10: * Any feedback is very welcome. alpar@10: * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html alpar@10: * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) alpar@10: */ alpar@10: alpar@10: #ifndef LEMON_RANDOM_H alpar@10: #define LEMON_RANDOM_H alpar@10: alpar@1134: #include alpar@1134: alpar@10: #include alpar@10: #include alpar@10: #include deba@110: #include deba@177: #include alpar@10: alpar@68: #include alpar@10: #include alpar@68: alpar@1134: #ifndef LEMON_WIN32 deba@177: #include deba@177: #include deba@177: #include deba@177: #include deba@177: #else alpar@491: #include deba@177: #endif deba@177: alpar@10: ///\ingroup misc alpar@10: ///\file alpar@10: ///\brief Mersenne Twister random number generator alpar@10: alpar@10: namespace lemon { alpar@10: alpar@10: namespace _random_bits { alpar@209: alpar@10: template ::digits> alpar@10: struct RandomTraits {}; alpar@10: alpar@10: template alpar@10: struct RandomTraits<_Word, 32> { alpar@10: alpar@10: typedef _Word Word; alpar@10: static const int bits = 32; alpar@10: alpar@10: static const int length = 624; alpar@10: static const int shift = 397; alpar@209: alpar@10: static const Word mul = 0x6c078965u; alpar@10: static const Word arrayInit = 0x012BD6AAu; alpar@10: static const Word arrayMul1 = 0x0019660Du; alpar@10: static const Word arrayMul2 = 0x5D588B65u; alpar@10: alpar@10: static const Word mask = 0x9908B0DFu; alpar@10: static const Word loMask = (1u << 31) - 1; alpar@10: static const Word hiMask = ~loMask; alpar@10: alpar@10: static Word tempering(Word rnd) { alpar@10: rnd ^= (rnd >> 11); alpar@10: rnd ^= (rnd << 7) & 0x9D2C5680u; alpar@10: rnd ^= (rnd << 15) & 0xEFC60000u; alpar@10: rnd ^= (rnd >> 18); alpar@10: return rnd; alpar@10: } alpar@10: alpar@10: }; alpar@10: alpar@10: template alpar@10: struct RandomTraits<_Word, 64> { alpar@10: alpar@10: typedef _Word Word; alpar@10: static const int bits = 64; alpar@10: alpar@10: static const int length = 312; alpar@10: static const int shift = 156; alpar@10: alpar@10: static const Word mul = Word(0x5851F42Du) << 32 | Word(0x4C957F2Du); alpar@10: static const Word arrayInit = Word(0x00000000u) << 32 |Word(0x012BD6AAu); alpar@10: static const Word arrayMul1 = Word(0x369DEA0Fu) << 32 |Word(0x31A53F85u); alpar@10: static const Word arrayMul2 = Word(0x27BB2EE6u) << 32 |Word(0x87B0B0FDu); alpar@10: alpar@10: static const Word mask = Word(0xB5026F5Au) << 32 | Word(0xA96619E9u); alpar@10: static const Word loMask = (Word(1u) << 31) - 1; alpar@10: static const Word hiMask = ~loMask; alpar@10: alpar@10: static Word tempering(Word rnd) { alpar@10: rnd ^= (rnd >> 29) & (Word(0x55555555u) << 32 | Word(0x55555555u)); alpar@10: rnd ^= (rnd << 17) & (Word(0x71D67FFFu) << 32 | Word(0xEDA60000u)); alpar@10: rnd ^= (rnd << 37) & (Word(0xFFF7EEE0u) << 32 | Word(0x00000000u)); alpar@10: rnd ^= (rnd >> 43); alpar@10: return rnd; alpar@10: } alpar@10: alpar@10: }; alpar@10: alpar@10: template alpar@10: class RandomCore { alpar@10: public: alpar@10: alpar@10: typedef _Word Word; alpar@10: alpar@10: private: alpar@10: alpar@10: static const int bits = RandomTraits::bits; alpar@10: alpar@10: static const int length = RandomTraits::length; alpar@10: static const int shift = RandomTraits::shift; alpar@10: alpar@10: public: alpar@10: alpar@10: void initState() { alpar@10: static const Word seedArray[4] = { alpar@10: 0x12345u, 0x23456u, 0x34567u, 0x45678u alpar@10: }; alpar@209: alpar@10: initState(seedArray, seedArray + 4); alpar@10: } alpar@10: alpar@10: void initState(Word seed) { alpar@10: alpar@10: static const Word mul = RandomTraits::mul; alpar@10: alpar@209: current = state; alpar@10: alpar@10: Word *curr = state + length - 1; alpar@10: curr[0] = seed; --curr; alpar@10: for (int i = 1; i < length; ++i) { alpar@10: curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i); alpar@10: --curr; alpar@10: } alpar@10: } alpar@10: alpar@10: template alpar@10: void initState(Iterator begin, Iterator end) { alpar@10: alpar@10: static const Word init = RandomTraits::arrayInit; alpar@10: static const Word mul1 = RandomTraits::arrayMul1; alpar@10: static const Word mul2 = RandomTraits::arrayMul2; alpar@10: alpar@10: alpar@10: Word *curr = state + length - 1; --curr; alpar@10: Iterator it = begin; int cnt = 0; alpar@10: int num; alpar@10: alpar@10: initState(init); alpar@10: alpar@1124: num = static_cast(length > end - begin ? length : end - begin); alpar@10: while (num--) { alpar@209: curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1)) alpar@10: + *it + cnt; alpar@10: ++it; ++cnt; alpar@10: if (it == end) { alpar@10: it = begin; cnt = 0; alpar@10: } alpar@10: if (curr == state) { alpar@10: curr = state + length - 1; curr[0] = state[0]; alpar@10: } alpar@10: --curr; alpar@10: } alpar@10: alpar@1124: num = length - 1; cnt = static_cast(length - (curr - state) - 1); alpar@10: while (num--) { alpar@10: curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2)) alpar@10: - cnt; alpar@10: --curr; ++cnt; alpar@10: if (curr == state) { alpar@10: curr = state + length - 1; curr[0] = state[0]; --curr; alpar@10: cnt = 1; alpar@10: } alpar@10: } alpar@209: alpar@10: state[length - 1] = Word(1) << (bits - 1); alpar@10: } alpar@209: alpar@10: void copyState(const RandomCore& other) { alpar@10: std::copy(other.state, other.state + length, state); alpar@10: current = state + (other.current - other.state); alpar@10: } alpar@10: alpar@10: Word operator()() { alpar@10: if (current == state) fillState(); alpar@10: --current; alpar@10: Word rnd = *current; alpar@10: return RandomTraits::tempering(rnd); alpar@10: } alpar@10: alpar@10: private: alpar@10: alpar@10: void fillState() { alpar@10: static const Word mask[2] = { 0x0ul, RandomTraits::mask }; alpar@10: static const Word loMask = RandomTraits::loMask; alpar@10: static const Word hiMask = RandomTraits::hiMask; alpar@10: alpar@209: current = state + length; alpar@10: alpar@1132: Word *curr = state + length - 1; alpar@1132: long num; alpar@209: alpar@10: num = length - shift; alpar@10: while (num--) { alpar@10: curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ alpar@10: curr[- shift] ^ mask[curr[-1] & 1ul]; alpar@10: --curr; alpar@10: } alpar@10: num = shift - 1; alpar@10: while (num--) { alpar@10: curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ alpar@10: curr[length - shift] ^ mask[curr[-1] & 1ul]; alpar@10: --curr; alpar@10: } deba@62: state[0] = (((state[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^ alpar@10: curr[length - shift] ^ mask[curr[length - 1] & 1ul]; alpar@10: alpar@10: } alpar@10: alpar@10: Word *current; alpar@10: Word state[length]; alpar@209: alpar@10: }; alpar@10: alpar@10: alpar@209: template ::digits + 1) / 2> alpar@10: struct Masker { alpar@10: static Result mask(const Result& result) { alpar@10: return Masker:: alpar@10: mask(static_cast(result | (result >> shift))); alpar@10: } alpar@10: }; alpar@209: alpar@10: template alpar@10: struct Masker { alpar@10: static Result mask(const Result& result) { alpar@10: return static_cast(result | (result >> 1)); alpar@10: } alpar@10: }; alpar@10: alpar@209: template ::digits, int shift = 0, alpar@1163: bool last = (rest <= std::numeric_limits::digits)> alpar@10: struct IntConversion { alpar@10: static const int bits = std::numeric_limits::digits; alpar@209: alpar@10: static Result convert(RandomCore& rnd) { alpar@10: return static_cast(rnd() >> (bits - rest)) << shift; alpar@10: } alpar@10: alpar@209: }; alpar@209: alpar@209: template alpar@10: struct IntConversion { alpar@10: static const int bits = std::numeric_limits::digits; alpar@10: alpar@10: static Result convert(RandomCore& rnd) { alpar@209: return (static_cast(rnd()) << shift) | alpar@10: IntConversion::convert(rnd); alpar@10: } alpar@10: }; alpar@10: alpar@10: alpar@10: template ::digits < alpar@209: std::numeric_limits::digits) > alpar@10: struct Mapping { alpar@10: static Result map(RandomCore& rnd, const Result& bound) { alpar@10: Word max = Word(bound - 1); alpar@10: Result mask = Masker::mask(bound - 1); alpar@10: Result num; alpar@10: do { alpar@209: num = IntConversion::convert(rnd) & mask; alpar@10: } while (num > max); alpar@10: return num; alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@10: struct Mapping { alpar@10: static Result map(RandomCore& rnd, const Result& bound) { alpar@10: Word max = Word(bound - 1); alpar@10: Word mask = Masker::digits + 1) / 2> alpar@10: ::mask(max); alpar@10: Word num; alpar@10: do { alpar@10: num = rnd() & mask; alpar@10: } while (num > max); alpar@1178: return static_cast(num); alpar@10: } alpar@10: }; alpar@10: kpeter@498: template alpar@10: struct ShiftMultiplier { alpar@10: static const Result multiplier() { alpar@10: Result res = ShiftMultiplier::multiplier(); alpar@10: res *= res; alpar@10: if ((exp & 1) == 1) res *= static_cast(0.5); alpar@209: return res; alpar@10: } alpar@10: }; alpar@10: alpar@10: template kpeter@498: struct ShiftMultiplier { alpar@10: static const Result multiplier() { alpar@209: return static_cast(1.0); alpar@10: } alpar@10: }; alpar@10: alpar@10: template kpeter@498: struct ShiftMultiplier { alpar@10: static const Result multiplier() { alpar@209: return static_cast(1.0/1048576.0); alpar@10: } alpar@10: }; alpar@209: alpar@10: template kpeter@498: struct ShiftMultiplier { alpar@10: static const Result multiplier() { kpeter@498: return static_cast(1.0/4294967296.0); alpar@10: } alpar@10: }; alpar@10: alpar@10: template kpeter@498: struct ShiftMultiplier { alpar@10: static const Result multiplier() { alpar@209: return static_cast(1.0/9007199254740992.0); alpar@10: } alpar@10: }; alpar@10: alpar@10: template kpeter@498: struct ShiftMultiplier { alpar@10: static const Result multiplier() { alpar@209: return static_cast(1.0/18446744073709551616.0); alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@10: struct Shifting { alpar@10: static Result shift(const Result& result) { alpar@10: return result * ShiftMultiplier::multiplier(); alpar@10: } alpar@10: }; alpar@10: alpar@10: template ::digits, int shift = 0, alpar@10: bool last = rest <= std::numeric_limits::digits> alpar@209: struct RealConversion{ alpar@10: static const int bits = std::numeric_limits::digits; alpar@10: alpar@10: static Result convert(RandomCore& rnd) { kpeter@498: return Shifting:: alpar@10: shift(static_cast(rnd() >> (bits - rest))); alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@209: struct RealConversion { alpar@10: static const int bits = std::numeric_limits::digits; alpar@10: alpar@10: static Result convert(RandomCore& rnd) { kpeter@498: return Shifting:: alpar@10: shift(static_cast(rnd())) + alpar@10: RealConversion:: alpar@10: convert(rnd); alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@10: struct Initializer { alpar@10: alpar@10: template alpar@10: static void init(RandomCore& rnd, Iterator begin, Iterator end) { alpar@10: std::vector ws; alpar@10: for (Iterator it = begin; it != end; ++it) { alpar@10: ws.push_back(Word(*it)); alpar@10: } alpar@10: rnd.initState(ws.begin(), ws.end()); alpar@10: } alpar@10: alpar@10: static void init(RandomCore& rnd, Result seed) { alpar@10: rnd.initState(seed); alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@10: struct BoolConversion { alpar@10: static bool convert(RandomCore& rnd) { alpar@10: return (rnd() & 1) == 1; alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@10: struct BoolProducer { alpar@10: Word buffer; alpar@10: int num; alpar@209: alpar@10: BoolProducer() : num(0) {} alpar@10: alpar@10: bool convert(RandomCore& rnd) { alpar@10: if (num == 0) { alpar@10: buffer = rnd(); alpar@10: num = RandomTraits::bits; alpar@10: } alpar@10: bool r = (buffer & 1); alpar@10: buffer >>= 1; alpar@10: --num; alpar@10: return r; alpar@10: } alpar@10: }; alpar@10: alpar@1163: /// \ingroup misc alpar@1163: /// alpar@1163: /// \brief Mersenne Twister random number generator alpar@1163: /// alpar@1163: /// The Mersenne Twister is a twisted generalized feedback alpar@1163: /// shift-register generator of Matsumoto and Nishimura. The period kpeter@1164: /// of this generator is \f$ 2^{19937} - 1\f$ and it is alpar@1163: /// equi-distributed in 623 dimensions for 32-bit numbers. The time alpar@1163: /// performance of this generator is comparable to the commonly used alpar@1163: /// generators. alpar@1163: /// kpeter@1164: /// This is a template implementation of both 32-bit and alpar@1163: /// 64-bit architecture optimized versions. The generators differ alpar@1163: /// sligthly in the initialization and generation phase so they alpar@1163: /// produce two completly different sequences. alpar@1163: /// alpar@1163: /// \alert Do not use this class directly, but instead one of \ref alpar@1163: /// Random, \ref Random32 or \ref Random64. alpar@1163: /// alpar@1163: /// The generator gives back random numbers of serveral types. To kpeter@1164: /// get a random number from a range of a floating point type, you alpar@1163: /// can use one form of the \c operator() or the \c real() member alpar@1163: /// function. If you want to get random number from the {0, 1, ..., kpeter@1164: /// n-1} integer range, use the \c operator[] or the \c integer() alpar@1163: /// method. And to get random number from the whole range of an kpeter@1164: /// integer type, you can use the argumentless \c integer() or kpeter@1164: /// \c uinteger() functions. Finally, you can get random bool with kpeter@1164: /// equal chance of true and false or with given probability of true kpeter@1164: /// result using the \c boolean() member functions. kpeter@1164: /// kpeter@1164: /// Various non-uniform distributions are also supported: normal (Gauss), kpeter@1164: /// exponential, gamma, Poisson, etc.; and a few two-dimensional kpeter@1164: /// distributions, too. alpar@1163: /// alpar@1163: ///\code alpar@1163: /// // The commented code is identical to the other alpar@1163: /// double a = rnd(); // [0.0, 1.0) alpar@1163: /// // double a = rnd.real(); // [0.0, 1.0) alpar@1163: /// double b = rnd(100.0); // [0.0, 100.0) alpar@1163: /// // double b = rnd.real(100.0); // [0.0, 100.0) alpar@1163: /// double c = rnd(1.0, 2.0); // [1.0, 2.0) alpar@1163: /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0) alpar@1163: /// int d = rnd[100000]; // 0..99999 alpar@1163: /// // int d = rnd.integer(100000); // 0..99999 alpar@1163: /// int e = rnd[6] + 1; // 1..6 alpar@1163: /// // int e = rnd.integer(1, 1 + 6); // 1..6 alpar@1163: /// int b = rnd.uinteger(); // 0 .. 2^31 - 1 alpar@1163: /// int c = rnd.integer(); // - 2^31 .. 2^31 - 1 alpar@1163: /// bool g = rnd.boolean(); // P(g = true) = 0.5 alpar@1163: /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 alpar@1163: ///\endcode alpar@1163: /// kpeter@1164: /// LEMON provides a global instance of the random number generator: kpeter@1164: /// \ref lemon::rnd "rnd". In most cases, it is a good practice kpeter@1164: /// to use this global generator to get random numbers. alpar@1163: /// alpar@1163: /// \sa \ref Random, \ref Random32 or \ref Random64. alpar@1163: template alpar@1163: class Random { alpar@1163: private: alpar@1163: alpar@1163: _random_bits::RandomCore core; alpar@1163: _random_bits::BoolProducer bool_producer; alpar@1163: alpar@1163: alpar@1163: public: alpar@1163: alpar@1163: ///\name Initialization alpar@1163: /// alpar@1163: /// @{ alpar@1163: alpar@1163: /// \brief Default constructor alpar@1163: /// alpar@1163: /// Constructor with constant seeding. alpar@1163: Random() { core.initState(); } alpar@1163: alpar@1163: /// \brief Constructor with seed alpar@1163: /// alpar@1163: /// Constructor with seed. The current number type will be converted alpar@1163: /// to the architecture word type. alpar@1163: template alpar@1163: Random(Number seed) { alpar@1163: _random_bits::Initializer::init(core, seed); alpar@1163: } alpar@1163: alpar@1163: /// \brief Constructor with array seeding alpar@1163: /// alpar@1163: /// Constructor with array seeding. The given range should contain alpar@1163: /// any number type and the numbers will be converted to the alpar@1163: /// architecture word type. alpar@1163: template alpar@1163: Random(Iterator begin, Iterator end) { alpar@1163: typedef typename std::iterator_traits::value_type Number; alpar@1163: _random_bits::Initializer::init(core, begin, end); alpar@1163: } alpar@1163: alpar@1163: /// \brief Copy constructor alpar@1163: /// alpar@1163: /// Copy constructor. The generated sequence will be identical to alpar@1163: /// the other sequence. It can be used to save the current state alpar@1163: /// of the generator and later use it to generate the same alpar@1163: /// sequence. alpar@1163: Random(const Random& other) { alpar@1163: core.copyState(other.core); alpar@1163: } alpar@1163: alpar@1163: /// \brief Assign operator alpar@1163: /// alpar@1163: /// Assign operator. The generated sequence will be identical to alpar@1163: /// the other sequence. It can be used to save the current state alpar@1163: /// of the generator and later use it to generate the same alpar@1163: /// sequence. alpar@1163: Random& operator=(const Random& other) { alpar@1163: if (&other != this) { alpar@1163: core.copyState(other.core); alpar@1163: } alpar@1163: return *this; alpar@1163: } alpar@1163: alpar@1163: /// \brief Seeding random sequence alpar@1163: /// alpar@1163: /// Seeding the random sequence. The current number type will be alpar@1163: /// converted to the architecture word type. alpar@1163: template alpar@1163: void seed(Number seed) { alpar@1163: _random_bits::Initializer::init(core, seed); alpar@1163: } alpar@1163: alpar@1163: /// \brief Seeding random sequence alpar@1163: /// alpar@1163: /// Seeding the random sequence. The given range should contain alpar@1163: /// any number type and the numbers will be converted to the alpar@1163: /// architecture word type. alpar@1163: template alpar@1163: void seed(Iterator begin, Iterator end) { alpar@1163: typedef typename std::iterator_traits::value_type Number; alpar@1163: _random_bits::Initializer::init(core, begin, end); alpar@1163: } alpar@1163: alpar@1163: /// \brief Seeding from file or from process id and time alpar@1163: /// alpar@1163: /// By default, this function calls the \c seedFromFile() member alpar@1163: /// function with the /dev/urandom file. If it does not success, alpar@1163: /// it uses the \c seedFromTime(). alpar@1163: /// \return Currently always \c true. alpar@1163: bool seed() { alpar@1163: #ifndef LEMON_WIN32 alpar@1163: if (seedFromFile("/dev/urandom", 0)) return true; alpar@1163: #endif alpar@1163: if (seedFromTime()) return true; alpar@1163: return false; alpar@1163: } alpar@1163: alpar@1163: /// \brief Seeding from file alpar@1163: /// alpar@1163: /// Seeding the random sequence from file. The linux kernel has two alpar@1163: /// devices, /dev/random and /dev/urandom which alpar@1163: /// could give good seed values for pseudo random generators (The alpar@1163: /// difference between two devices is that the random may alpar@1163: /// block the reading operation while the kernel can give good alpar@1163: /// source of randomness, while the urandom does not alpar@1163: /// block the input, but it could give back bytes with worse alpar@1163: /// entropy). alpar@1163: /// \param file The source file alpar@1163: /// \param offset The offset, from the file read. alpar@1163: /// \return \c true when the seeding successes. alpar@1163: #ifndef LEMON_WIN32 alpar@1163: bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0) alpar@1163: #else alpar@1163: bool seedFromFile(const std::string& file = "", int offset = 0) alpar@1163: #endif alpar@1163: { alpar@1163: std::ifstream rs(file.c_str()); alpar@1163: const int size = 4; alpar@1163: Word buf[size]; alpar@1163: if (offset != 0 && !rs.seekg(offset)) return false; alpar@1163: if (!rs.read(reinterpret_cast(buf), sizeof(buf))) return false; alpar@1163: seed(buf, buf + size); alpar@1163: return true; alpar@1163: } alpar@1163: kpeter@1164: /// \brief Seeding from process id and time alpar@1163: /// kpeter@1164: /// Seeding from process id and time. This function uses the alpar@1163: /// current process id and the current time for initialize the alpar@1163: /// random sequence. alpar@1163: /// \return Currently always \c true. alpar@1163: bool seedFromTime() { alpar@1163: #ifndef LEMON_WIN32 alpar@1163: timeval tv; alpar@1163: gettimeofday(&tv, 0); alpar@1163: seed(getpid() + tv.tv_sec + tv.tv_usec); alpar@1163: #else alpar@1163: seed(bits::getWinRndSeed()); alpar@1163: #endif alpar@1163: return true; alpar@1163: } alpar@1163: alpar@1163: /// @} alpar@1163: alpar@1163: ///\name Uniform Distributions alpar@1163: /// alpar@1163: /// @{ alpar@1163: alpar@1163: /// \brief Returns a random real number from the range [0, 1) alpar@1163: /// alpar@1163: /// It returns a random real number from the range [0, 1). The alpar@1163: /// default Number type is \c double. alpar@1163: template alpar@1163: Number real() { alpar@1163: return _random_bits::RealConversion::convert(core); alpar@1163: } alpar@1163: alpar@1163: double real() { alpar@1163: return real(); alpar@1163: } alpar@1163: alpar@1163: /// \brief Returns a random real number from the range [0, 1) alpar@1163: /// alpar@1163: /// It returns a random double from the range [0, 1). alpar@1163: double operator()() { alpar@1163: return real(); alpar@1163: } alpar@1163: alpar@1163: /// \brief Returns a random real number from the range [0, b) alpar@1163: /// alpar@1163: /// It returns a random real number from the range [0, b). alpar@1163: double operator()(double b) { alpar@1163: return real() * b; alpar@1163: } alpar@1163: alpar@1163: /// \brief Returns a random real number from the range [a, b) alpar@1163: /// alpar@1163: /// It returns a random real number from the range [a, b). alpar@1163: double operator()(double a, double b) { alpar@1163: return real() * (b - a) + a; alpar@1163: } alpar@1163: alpar@1163: /// \brief Returns a random integer from a range alpar@1163: /// alpar@1163: /// It returns a random integer from the range {0, 1, ..., b - 1}. alpar@1163: template alpar@1163: Number integer(Number b) { alpar@1163: return _random_bits::Mapping::map(core, b); alpar@1163: } alpar@1163: alpar@1163: /// \brief Returns a random integer from a range alpar@1163: /// alpar@1163: /// It returns a random integer from the range {a, a + 1, ..., b - 1}. alpar@1163: template alpar@1163: Number integer(Number a, Number b) { alpar@1163: return _random_bits::Mapping::map(core, b - a) + a; alpar@1163: } alpar@1163: alpar@1163: /// \brief Returns a random integer from a range alpar@1163: /// alpar@1163: /// It returns a random integer from the range {0, 1, ..., b - 1}. alpar@1163: template alpar@1163: Number operator[](Number b) { alpar@1163: return _random_bits::Mapping::map(core, b); alpar@1163: } alpar@1163: alpar@1163: /// \brief Returns a random non-negative integer alpar@1163: /// alpar@1163: /// It returns a random non-negative integer uniformly from the alpar@1163: /// whole range of the current \c Number type. The default result alpar@1163: /// type of this function is unsigned int. alpar@1163: template alpar@1163: Number uinteger() { alpar@1163: return _random_bits::IntConversion::convert(core); alpar@1163: } alpar@1163: alpar@1163: unsigned int uinteger() { alpar@1163: return uinteger(); alpar@1163: } alpar@1163: alpar@1163: /// \brief Returns a random integer alpar@1163: /// alpar@1163: /// It returns a random integer uniformly from the whole range of alpar@1163: /// the current \c Number type. The default result type of this alpar@1163: /// function is \c int. alpar@1163: template alpar@1163: Number integer() { alpar@1163: static const int nb = std::numeric_limits::digits + alpar@1163: (std::numeric_limits::is_signed ? 1 : 0); alpar@1163: return _random_bits::IntConversion::convert(core); alpar@1163: } alpar@1163: alpar@1163: int integer() { alpar@1163: return integer(); alpar@1163: } alpar@1163: alpar@1163: /// \brief Returns a random bool alpar@1163: /// alpar@1163: /// It returns a random bool. The generator holds a buffer for alpar@1163: /// random bits. Every time when it become empty the generator makes alpar@1163: /// a new random word and fill the buffer up. alpar@1163: bool boolean() { alpar@1163: return bool_producer.convert(core); alpar@1163: } alpar@1163: alpar@1163: /// @} alpar@1163: alpar@1163: ///\name Non-uniform Distributions alpar@1163: /// alpar@1163: ///@{ alpar@1163: alpar@1163: /// \brief Returns a random bool with given probability of true result. alpar@1163: /// alpar@1163: /// It returns a random bool with given probability of true result. alpar@1163: bool boolean(double p) { alpar@1163: return operator()() < p; alpar@1163: } alpar@1163: alpar@1163: /// Standard normal (Gauss) distribution alpar@1163: alpar@1163: /// Standard normal (Gauss) distribution. alpar@1163: /// \note The Cartesian form of the Box-Muller alpar@1163: /// transformation is used to generate a random normal distribution. alpar@1163: double gauss() alpar@1163: { alpar@1163: double V1,V2,S; alpar@1163: do { alpar@1163: V1=2*real()-1; alpar@1163: V2=2*real()-1; alpar@1163: S=V1*V1+V2*V2; alpar@1163: } while(S>=1); alpar@1163: return std::sqrt(-2*std::log(S)/S)*V1; alpar@1163: } alpar@1163: /// Normal (Gauss) distribution with given mean and standard deviation alpar@1163: alpar@1163: /// Normal (Gauss) distribution with given mean and standard deviation. alpar@1163: /// \sa gauss() alpar@1163: double gauss(double mean,double std_dev) alpar@1163: { alpar@1163: return gauss()*std_dev+mean; alpar@1163: } alpar@1163: alpar@1163: /// Lognormal distribution alpar@1163: alpar@1163: /// Lognormal distribution. The parameters are the mean and the standard alpar@1163: /// deviation of exp(X). alpar@1163: /// alpar@1163: double lognormal(double n_mean,double n_std_dev) alpar@1163: { alpar@1163: return std::exp(gauss(n_mean,n_std_dev)); alpar@1163: } alpar@1163: /// Lognormal distribution alpar@1163: alpar@1163: /// Lognormal distribution. The parameter is an std::pair of alpar@1163: /// the mean and the standard deviation of exp(X). alpar@1163: /// alpar@1163: double lognormal(const std::pair ¶ms) alpar@1163: { alpar@1163: return std::exp(gauss(params.first,params.second)); alpar@1163: } alpar@1163: /// Compute the lognormal parameters from mean and standard deviation alpar@1163: alpar@1163: /// This function computes the lognormal parameters from mean and alpar@1163: /// standard deviation. The return value can direcly be passed to alpar@1163: /// lognormal(). alpar@1163: std::pair lognormalParamsFromMD(double mean, alpar@1163: double std_dev) alpar@1163: { alpar@1163: double fr=std_dev/mean; alpar@1163: fr*=fr; alpar@1163: double lg=std::log(1+fr); alpar@1163: return std::pair(std::log(mean)-lg/2.0,std::sqrt(lg)); alpar@1163: } alpar@1163: /// Lognormal distribution with given mean and standard deviation alpar@1163: alpar@1163: /// Lognormal distribution with given mean and standard deviation. alpar@1163: /// alpar@1163: double lognormalMD(double mean,double std_dev) alpar@1163: { alpar@1163: return lognormal(lognormalParamsFromMD(mean,std_dev)); alpar@1163: } alpar@1163: alpar@1163: /// Exponential distribution with given mean alpar@1163: alpar@1163: /// This function generates an exponential distribution random number alpar@1163: /// with mean 1/lambda. alpar@1163: /// alpar@1163: double exponential(double lambda=1.0) alpar@1163: { alpar@1163: return -std::log(1.0-real())/lambda; alpar@1163: } alpar@1163: alpar@1163: /// Gamma distribution with given integer shape alpar@1163: alpar@1163: /// This function generates a gamma distribution random number. alpar@1163: /// alpar@1163: ///\param k shape parameter (k>0 integer) alpar@1163: double gamma(int k) alpar@1163: { alpar@1163: double s = 0; alpar@1163: for(int i=0;i()); alpar@1163: return s; alpar@1163: } alpar@1163: alpar@1163: /// Gamma distribution with given shape and scale parameter alpar@1163: alpar@1163: /// This function generates a gamma distribution random number. alpar@1163: /// alpar@1163: ///\param k shape parameter (k>0) alpar@1163: ///\param theta scale parameter alpar@1163: /// alpar@1163: double gamma(double k,double theta=1.0) alpar@1163: { alpar@1163: double xi,nu; alpar@1163: const double delta = k-std::floor(k); alpar@1163: const double v0=E/(E-delta); alpar@1163: do { alpar@1163: double V0=1.0-real(); alpar@1163: double V1=1.0-real(); alpar@1163: double V2=1.0-real(); alpar@1163: if(V2<=v0) alpar@1163: { alpar@1163: xi=std::pow(V1,1.0/delta); alpar@1163: nu=V0*std::pow(xi,delta-1.0); alpar@1163: } alpar@1163: else alpar@1163: { alpar@1163: xi=1.0-std::log(V1); alpar@1163: nu=V0*std::exp(-xi); alpar@1163: } alpar@1163: } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); alpar@1163: return theta*(xi+gamma(int(std::floor(k)))); alpar@1163: } alpar@1163: alpar@1163: /// Weibull distribution alpar@1163: alpar@1163: /// This function generates a Weibull distribution random number. alpar@1163: /// alpar@1163: ///\param k shape parameter (k>0) alpar@1163: ///\param lambda scale parameter (lambda>0) alpar@1163: /// alpar@1163: double weibull(double k,double lambda) alpar@1163: { alpar@1163: return lambda*pow(-std::log(1.0-real()),1.0/k); alpar@1163: } alpar@1163: alpar@1163: /// Pareto distribution alpar@1163: alpar@1163: /// This function generates a Pareto distribution random number. alpar@1163: /// alpar@1163: ///\param k shape parameter (k>0) alpar@1163: ///\param x_min location parameter (x_min>0) alpar@1163: /// alpar@1163: double pareto(double k,double x_min) alpar@1163: { alpar@1163: return exponential(gamma(k,1.0/x_min))+x_min; alpar@1163: } alpar@1163: alpar@1163: /// Poisson distribution alpar@1163: alpar@1163: /// This function generates a Poisson distribution random number with alpar@1163: /// parameter \c lambda. alpar@1163: /// alpar@1163: /// The probability mass function of this distribusion is alpar@1163: /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f] alpar@1163: /// \note The algorithm is taken from the book of Donald E. Knuth titled alpar@1163: /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the alpar@1163: /// return value. alpar@1163: alpar@1163: int poisson(double lambda) alpar@1163: { alpar@1163: const double l = std::exp(-lambda); alpar@1163: int k=0; alpar@1163: double p = 1.0; alpar@1163: do { alpar@1163: k++; alpar@1163: p*=real(); alpar@1163: } while (p>=l); alpar@1163: return k-1; alpar@1163: } alpar@1163: alpar@1163: ///@} alpar@1163: kpeter@1164: ///\name Two-Dimensional Distributions alpar@1163: /// alpar@1163: ///@{ alpar@1163: alpar@1163: /// Uniform distribution on the full unit circle alpar@1163: alpar@1163: /// Uniform distribution on the full unit circle. alpar@1163: /// alpar@1163: dim2::Point disc() alpar@1163: { alpar@1163: double V1,V2; alpar@1163: do { alpar@1163: V1=2*real()-1; alpar@1163: V2=2*real()-1; alpar@1163: alpar@1163: } while(V1*V1+V2*V2>=1); alpar@1163: return dim2::Point(V1,V2); alpar@1163: } kpeter@1164: /// A kind of two-dimensional normal (Gauss) distribution alpar@1163: alpar@1163: /// This function provides a turning symmetric two-dimensional distribution. alpar@1163: /// Both coordinates are of standard normal distribution, but they are not alpar@1163: /// independent. alpar@1163: /// alpar@1163: /// \note The coordinates are the two random variables provided by alpar@1163: /// the Box-Muller method. alpar@1163: dim2::Point gauss2() alpar@1163: { alpar@1163: double V1,V2,S; alpar@1163: do { alpar@1163: V1=2*real()-1; alpar@1163: V2=2*real()-1; alpar@1163: S=V1*V1+V2*V2; alpar@1163: } while(S>=1); alpar@1163: double W=std::sqrt(-2*std::log(S)/S); alpar@1163: return dim2::Point(W*V1,W*V2); alpar@1163: } kpeter@1164: /// A kind of two-dimensional exponential distribution alpar@1163: alpar@1163: /// This function provides a turning symmetric two-dimensional distribution. alpar@1163: /// The x-coordinate is of conditionally exponential distribution alpar@1163: /// with the condition that x is positive and y=0. If x is negative and alpar@1163: /// y=0 then, -x is of exponential distribution. The same is true for the alpar@1163: /// y-coordinate. alpar@1163: dim2::Point exponential2() alpar@1163: { alpar@1163: double V1,V2,S; alpar@1163: do { alpar@1163: V1=2*real()-1; alpar@1163: V2=2*real()-1; alpar@1163: S=V1*V1+V2*V2; alpar@1163: } while(S>=1); alpar@1163: double W=-std::log(S)/S; alpar@1163: return dim2::Point(W*V1,W*V2); alpar@1163: } alpar@1163: alpar@1163: ///@} alpar@1163: }; alpar@1163: alpar@1163: alpar@1163: }; alpar@10: alpar@10: /// \ingroup misc alpar@10: /// alpar@10: /// \brief Mersenne Twister random number generator alpar@10: /// kpeter@1164: /// This class implements either the 32-bit or the 64-bit version of alpar@1163: /// the Mersenne Twister random number generator algorithm kpeter@1164: /// depending on the system architecture. alpar@1163: /// kpeter@1164: /// For the API description, see its base class kpeter@1164: /// \ref _random_bits::Random. alpar@10: /// alpar@1163: /// \sa \ref _random_bits::Random alpar@1163: typedef _random_bits::Random Random; kpeter@1164: alpar@1163: /// \ingroup misc alpar@10: /// kpeter@1164: /// \brief Mersenne Twister random number generator (32-bit version) alpar@10: /// kpeter@1164: /// This class implements the 32-bit version of the Mersenne Twister alpar@1163: /// random number generator algorithm. It is recommended to be used alpar@1163: /// when someone wants to make sure that the \e same pseudo random alpar@1163: /// sequence will be generated on every platfrom. alpar@10: /// kpeter@1164: /// For the API description, see its base class kpeter@1164: /// \ref _random_bits::Random. alpar@1163: /// alpar@1163: /// \sa \ref _random_bits::Random kpeter@1164: typedef _random_bits::Random Random32; alpar@10: alpar@1163: /// \ingroup misc alpar@1163: /// kpeter@1164: /// \brief Mersenne Twister random number generator (64-bit version) alpar@1163: /// kpeter@1164: /// This class implements the 64-bit version of the Mersenne Twister kpeter@1164: /// random number generator algorithm. (Even though it runs kpeter@1164: /// on 32-bit architectures, too.) It is recommended to be used when alpar@1163: /// someone wants to make sure that the \e same pseudo random sequence alpar@1163: /// will be generated on every platfrom. alpar@1163: /// kpeter@1164: /// For the API description, see its base class kpeter@1164: /// \ref _random_bits::Random. alpar@1163: /// alpar@1163: /// \sa \ref _random_bits::Random alpar@1163: typedef _random_bits::Random Random64; alpar@10: alpar@10: extern Random rnd; alpar@1163: alpar@10: } alpar@10: alpar@10: #endif