alpar@10: /* -*- C++ -*- alpar@10: * alpar@10: * This file is a part of LEMON, a generic C++ optimization library alpar@10: * alpar@39: * Copyright (C) 2003-2008 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@10: * alpar@10: * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, alpar@10: * 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@10: * 3. The names of its contributors may not be used to endorse or promote alpar@10: * 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@10: #include alpar@10: #include alpar@10: #include alpar@10: alpar@10: #include alpar@10: alpar@68: #include alpar@10: #include alpar@68: 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@10: 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@10: 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: 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@10: 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@10: 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@10: num = length > end - begin ? length : end - begin; alpar@10: while (num--) { alpar@10: 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@10: num = length - 1; cnt = 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@10: alpar@10: state[length - 1] = Word(1) << (bits - 1); alpar@10: } alpar@10: 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: 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@10: current = state + length; alpar@10: alpar@10: register Word *curr = state + length - 1; alpar@10: register long num; alpar@10: 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: alpar@10: Word *current; alpar@10: Word state[length]; alpar@10: alpar@10: }; alpar@10: alpar@10: alpar@10: 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@10: 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@10: template ::digits, int shift = 0, alpar@10: bool last = rest <= std::numeric_limits::digits> alpar@10: struct IntConversion { alpar@10: static const int bits = std::numeric_limits::digits; alpar@10: alpar@10: static Result convert(RandomCore& rnd) { alpar@10: return static_cast(rnd() >> (bits - rest)) << shift; alpar@10: } alpar@10: alpar@10: }; alpar@10: alpar@10: 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@10: return (static_cast(rnd()) << shift) | alpar@10: IntConversion::convert(rnd); alpar@10: } alpar@10: }; alpar@10: alpar@10: alpar@10: template ::digits < alpar@10: 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@10: 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@10: return num; alpar@10: } alpar@10: }; alpar@10: alpar@10: template = 0)> 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(2.0); alpar@10: return res; alpar@10: } alpar@10: }; alpar@10: alpar@10: 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@10: return res; alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@10: struct ShiftMultiplier { alpar@10: static const Result multiplier() { alpar@10: return static_cast(1.0); alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@10: struct ShiftMultiplier { alpar@10: static const Result multiplier() { alpar@10: return static_cast(1.0/1048576.0); alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@10: struct ShiftMultiplier { alpar@10: static const Result multiplier() { alpar@10: return static_cast(1.0/424967296.0); alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@10: struct ShiftMultiplier { alpar@10: static const Result multiplier() { alpar@10: return static_cast(1.0/9007199254740992.0); alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@10: struct ShiftMultiplier { alpar@10: static const Result multiplier() { alpar@10: 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@10: struct RealConversion{ alpar@10: static const int bits = std::numeric_limits::digits; alpar@10: alpar@10: static Result convert(RandomCore& rnd) { alpar@10: return Shifting:: alpar@10: shift(static_cast(rnd() >> (bits - rest))); alpar@10: } alpar@10: }; alpar@10: alpar@10: template alpar@10: struct RealConversion { alpar@10: static const int bits = std::numeric_limits::digits; alpar@10: alpar@10: static Result convert(RandomCore& rnd) { alpar@10: 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@10: 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@10: } alpar@10: alpar@10: /// \ingroup misc alpar@10: /// alpar@10: /// \brief Mersenne Twister random number generator alpar@10: /// alpar@10: /// The Mersenne Twister is a twisted generalized feedback alpar@10: /// shift-register generator of Matsumoto and Nishimura. The period alpar@10: /// of this generator is \f$ 2^{19937} - 1 \f$ and it is alpar@10: /// equi-distributed in 623 dimensions for 32-bit numbers. The time alpar@10: /// performance of this generator is comparable to the commonly used alpar@10: /// generators. alpar@10: /// alpar@10: /// This implementation is specialized for both 32-bit and 64-bit alpar@10: /// architectures. The generators differ sligthly in the alpar@10: /// initialization and generation phase so they produce two alpar@10: /// completly different sequences. alpar@10: /// alpar@10: /// The generator gives back random numbers of serveral types. To alpar@10: /// get a random number from a range of a floating point type you alpar@10: /// can use one form of the \c operator() or the \c real() member alpar@10: /// function. If you want to get random number from the {0, 1, ..., alpar@10: /// n-1} integer range use the \c operator[] or the \c integer() alpar@10: /// method. And to get random number from the whole range of an alpar@10: /// integer type you can use the argumentless \c integer() or \c alpar@10: /// uinteger() functions. After all you can get random bool with alpar@10: /// equal chance of true and false or given probability of true alpar@10: /// result with the \c boolean() member functions. alpar@10: /// alpar@10: ///\code alpar@10: /// // The commented code is identical to the other alpar@10: /// double a = rnd(); // [0.0, 1.0) alpar@10: /// // double a = rnd.real(); // [0.0, 1.0) alpar@10: /// double b = rnd(100.0); // [0.0, 100.0) alpar@10: /// // double b = rnd.real(100.0); // [0.0, 100.0) alpar@10: /// double c = rnd(1.0, 2.0); // [1.0, 2.0) alpar@10: /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0) alpar@10: /// int d = rnd[100000]; // 0..99999 alpar@10: /// // int d = rnd.integer(100000); // 0..99999 alpar@10: /// int e = rnd[6] + 1; // 1..6 alpar@10: /// // int e = rnd.integer(1, 1 + 6); // 1..6 alpar@10: /// int b = rnd.uinteger(); // 0 .. 2^31 - 1 alpar@10: /// int c = rnd.integer(); // - 2^31 .. 2^31 - 1 alpar@10: /// bool g = rnd.boolean(); // P(g = true) = 0.5 alpar@10: /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 alpar@10: ///\endcode alpar@10: /// kpeter@49: /// LEMON provides a global instance of the random number alpar@10: /// generator which name is \ref lemon::rnd "rnd". Usually it is a alpar@10: /// good programming convenience to use this global generator to get alpar@10: /// random numbers. alpar@10: class Random { alpar@10: private: alpar@10: kpeter@16: // Architecture word alpar@10: typedef unsigned long Word; alpar@10: alpar@10: _random_bits::RandomCore core; alpar@10: _random_bits::BoolProducer bool_producer; alpar@10: alpar@10: alpar@10: public: alpar@10: kpeter@49: /// \brief Default constructor alpar@10: /// alpar@10: /// Constructor with constant seeding. alpar@10: Random() { core.initState(); } alpar@10: kpeter@49: /// \brief Constructor with seed alpar@10: /// alpar@10: /// Constructor with seed. The current number type will be converted alpar@10: /// to the architecture word type. alpar@10: template alpar@10: Random(Number seed) { alpar@10: _random_bits::Initializer::init(core, seed); alpar@10: } alpar@10: kpeter@49: /// \brief Constructor with array seeding alpar@10: /// alpar@10: /// Constructor with array seeding. The given range should contain alpar@10: /// any number type and the numbers will be converted to the alpar@10: /// architecture word type. alpar@10: template alpar@10: Random(Iterator begin, Iterator end) { alpar@10: typedef typename std::iterator_traits::value_type Number; alpar@10: _random_bits::Initializer::init(core, begin, end); alpar@10: } alpar@10: alpar@10: /// \brief Copy constructor alpar@10: /// alpar@10: /// Copy constructor. The generated sequence will be identical to alpar@10: /// the other sequence. It can be used to save the current state alpar@10: /// of the generator and later use it to generate the same alpar@10: /// sequence. alpar@10: Random(const Random& other) { alpar@10: core.copyState(other.core); alpar@10: } alpar@10: alpar@10: /// \brief Assign operator alpar@10: /// alpar@10: /// Assign operator. The generated sequence will be identical to alpar@10: /// the other sequence. It can be used to save the current state alpar@10: /// of the generator and later use it to generate the same alpar@10: /// sequence. alpar@10: Random& operator=(const Random& other) { alpar@10: if (&other != this) { alpar@10: core.copyState(other.core); alpar@10: } alpar@10: return *this; alpar@10: } alpar@10: alpar@10: /// \brief Returns a random real number from the range [0, 1) alpar@10: /// alpar@10: /// It returns a random real number from the range [0, 1). The kpeter@49: /// default Number type is \c double. alpar@10: template alpar@10: Number real() { alpar@10: return _random_bits::RealConversion::convert(core); alpar@10: } alpar@10: alpar@10: double real() { alpar@10: return real(); alpar@10: } alpar@10: alpar@10: /// \brief Returns a random real number the range [0, b) alpar@10: /// alpar@10: /// It returns a random real number from the range [0, b). alpar@10: template alpar@10: Number real(Number b) { alpar@10: return real() * b; alpar@10: } alpar@10: alpar@10: /// \brief Returns a random real number from the range [a, b) alpar@10: /// alpar@10: /// It returns a random real number from the range [a, b). alpar@10: template alpar@10: Number real(Number a, Number b) { alpar@10: return real() * (b - a) + a; alpar@10: } alpar@10: alpar@10: /// \brief Returns a random real number from the range [0, 1) alpar@10: /// alpar@10: /// It returns a random double from the range [0, 1). alpar@10: double operator()() { alpar@10: return real(); alpar@10: } alpar@10: alpar@10: /// \brief Returns a random real number from the range [0, b) alpar@10: /// alpar@10: /// It returns a random real number from the range [0, b). alpar@10: template alpar@10: Number operator()(Number b) { alpar@10: return real() * b; alpar@10: } alpar@10: alpar@10: /// \brief Returns a random real number from the range [a, b) alpar@10: /// alpar@10: /// It returns a random real number from the range [a, b). alpar@10: template alpar@10: Number operator()(Number a, Number b) { alpar@10: return real() * (b - a) + a; alpar@10: } alpar@10: alpar@10: /// \brief Returns a random integer from a range alpar@10: /// alpar@10: /// It returns a random integer from the range {0, 1, ..., b - 1}. alpar@10: template alpar@10: Number integer(Number b) { alpar@10: return _random_bits::Mapping::map(core, b); alpar@10: } alpar@10: alpar@10: /// \brief Returns a random integer from a range alpar@10: /// alpar@10: /// It returns a random integer from the range {a, a + 1, ..., b - 1}. alpar@10: template alpar@10: Number integer(Number a, Number b) { alpar@10: return _random_bits::Mapping::map(core, b - a) + a; alpar@10: } alpar@10: alpar@10: /// \brief Returns a random integer from a range alpar@10: /// alpar@10: /// It returns a random integer from the range {0, 1, ..., b - 1}. alpar@10: template alpar@10: Number operator[](Number b) { alpar@10: return _random_bits::Mapping::map(core, b); alpar@10: } alpar@10: alpar@10: /// \brief Returns a random non-negative integer alpar@10: /// alpar@10: /// It returns a random non-negative integer uniformly from the kpeter@49: /// whole range of the current \c Number type. The default result kpeter@49: /// type of this function is unsigned int. alpar@10: template alpar@10: Number uinteger() { alpar@10: return _random_bits::IntConversion::convert(core); alpar@10: } alpar@10: alpar@10: unsigned int uinteger() { alpar@10: return uinteger(); alpar@10: } alpar@10: alpar@10: /// \brief Returns a random integer alpar@10: /// alpar@10: /// It returns a random integer uniformly from the whole range of alpar@10: /// the current \c Number type. The default result type of this kpeter@49: /// function is \c int. alpar@10: template alpar@10: Number integer() { alpar@10: static const int nb = std::numeric_limits::digits + alpar@10: (std::numeric_limits::is_signed ? 1 : 0); alpar@10: return _random_bits::IntConversion::convert(core); alpar@10: } alpar@10: alpar@10: int integer() { alpar@10: return integer(); alpar@10: } alpar@10: alpar@10: /// \brief Returns a random bool alpar@10: /// alpar@10: /// It returns a random bool. The generator holds a buffer for alpar@10: /// random bits. Every time when it become empty the generator makes alpar@10: /// a new random word and fill the buffer up. alpar@10: bool boolean() { alpar@10: return bool_producer.convert(core); alpar@10: } alpar@10: kpeter@49: ///\name Non-uniform distributions alpar@10: /// alpar@10: alpar@10: ///@{ alpar@10: alpar@10: /// \brief Returns a random bool alpar@10: /// kpeter@23: /// It returns a random bool with given probability of true result. alpar@10: bool boolean(double p) { alpar@10: return operator()() < p; alpar@10: } alpar@10: alpar@10: /// Standard Gauss distribution alpar@10: alpar@10: /// Standard Gauss distribution. alpar@10: /// \note The Cartesian form of the Box-Muller alpar@10: /// transformation is used to generate a random normal distribution. alpar@10: /// \todo Consider using the "ziggurat" method instead. alpar@10: double gauss() alpar@10: { alpar@10: double V1,V2,S; alpar@10: do { alpar@10: V1=2*real()-1; alpar@10: V2=2*real()-1; alpar@10: S=V1*V1+V2*V2; alpar@10: } while(S>=1); alpar@10: return std::sqrt(-2*std::log(S)/S)*V1; alpar@10: } alpar@10: /// Gauss distribution with given mean and standard deviation alpar@10: kpeter@23: /// Gauss distribution with given mean and standard deviation. alpar@10: /// \sa gauss() alpar@10: double gauss(double mean,double std_dev) alpar@10: { alpar@10: return gauss()*std_dev+mean; alpar@10: } alpar@10: alpar@10: /// Exponential distribution with given mean alpar@10: alpar@10: /// This function generates an exponential distribution random number alpar@10: /// with mean 1/lambda. alpar@10: /// alpar@10: double exponential(double lambda=1.0) alpar@10: { alpar@11: return -std::log(1.0-real())/lambda; alpar@10: } alpar@10: alpar@10: /// Gamma distribution with given integer shape alpar@10: alpar@10: /// This function generates a gamma distribution random number. alpar@10: /// alpar@10: ///\param k shape parameter (k>0 integer) alpar@10: double gamma(int k) alpar@10: { alpar@10: double s = 0; alpar@10: for(int i=0;i()); alpar@10: return s; alpar@10: } alpar@10: alpar@10: /// Gamma distribution with given shape and scale parameter alpar@10: alpar@10: /// This function generates a gamma distribution random number. alpar@10: /// alpar@10: ///\param k shape parameter (k>0) alpar@10: ///\param theta scale parameter alpar@10: /// alpar@10: double gamma(double k,double theta=1.0) alpar@10: { alpar@10: double xi,nu; alpar@10: const double delta = k-std::floor(k); alpar@68: const double v0=E/(E-delta); alpar@10: do { alpar@10: double V0=1.0-real(); alpar@10: double V1=1.0-real(); alpar@10: double V2=1.0-real(); alpar@10: if(V2<=v0) alpar@10: { alpar@10: xi=std::pow(V1,1.0/delta); alpar@10: nu=V0*std::pow(xi,delta-1.0); alpar@10: } alpar@10: else alpar@10: { alpar@10: xi=1.0-std::log(V1); alpar@10: nu=V0*std::exp(-xi); alpar@10: } alpar@10: } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); alpar@10: return theta*(xi-gamma(int(std::floor(k)))); alpar@10: } alpar@10: alpar@11: /// Weibull distribution alpar@11: alpar@11: /// This function generates a Weibull distribution random number. alpar@11: /// alpar@11: ///\param k shape parameter (k>0) alpar@11: ///\param lambda scale parameter (lambda>0) alpar@11: /// alpar@11: double weibull(double k,double lambda) alpar@11: { alpar@11: return lambda*pow(-std::log(1.0-real()),1.0/k); alpar@11: } alpar@11: alpar@11: /// Pareto distribution alpar@11: alpar@11: /// This function generates a Pareto distribution random number. alpar@11: /// alpar@12: ///\param k shape parameter (k>0) alpar@11: ///\param x_min location parameter (x_min>0) alpar@11: /// alpar@12: double pareto(double k,double x_min) alpar@11: { alpar@12: return exponential(gamma(k,1.0/x_min)); alpar@11: } alpar@10: alpar@92: /// Poisson distribution alpar@92: alpar@92: /// This function generates a Poisson distribution random number with alpar@92: /// parameter \c lambda. alpar@92: /// alpar@92: /// The probability mass function of this distribusion is alpar@92: /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f] alpar@92: /// \note The algorithm is taken from the book of Donald E. Knuth titled alpar@92: /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the alpar@92: /// return value. alpar@92: alpar@92: int poisson(double lambda) alpar@92: { alpar@92: const double l = std::exp(-lambda); alpar@92: int k=0; alpar@92: double p = 1.0; alpar@92: do { alpar@92: k++; alpar@92: p*=real(); alpar@92: } while (p>=l); alpar@92: return k-1; alpar@92: } alpar@92: alpar@10: ///@} alpar@10: alpar@10: ///\name Two dimensional distributions alpar@10: /// alpar@10: alpar@10: ///@{ alpar@10: kpeter@23: /// Uniform distribution on the full unit circle kpeter@16: kpeter@16: /// Uniform distribution on the full unit circle. kpeter@16: /// alpar@10: dim2::Point disc() alpar@10: { alpar@10: double V1,V2; alpar@10: do { alpar@10: V1=2*real()-1; alpar@10: V2=2*real()-1; alpar@10: alpar@10: } while(V1*V1+V2*V2>=1); alpar@10: return dim2::Point(V1,V2); alpar@10: } alpar@10: /// A kind of two dimensional Gauss distribution alpar@10: alpar@10: /// This function provides a turning symmetric two-dimensional distribution. alpar@10: /// Both coordinates are of standard normal distribution, but they are not alpar@10: /// independent. alpar@10: /// alpar@10: /// \note The coordinates are the two random variables provided by alpar@10: /// the Box-Muller method. alpar@10: dim2::Point gauss2() alpar@10: { alpar@10: double V1,V2,S; alpar@10: do { alpar@10: V1=2*real()-1; alpar@10: V2=2*real()-1; alpar@10: S=V1*V1+V2*V2; alpar@10: } while(S>=1); alpar@10: double W=std::sqrt(-2*std::log(S)/S); alpar@10: return dim2::Point(W*V1,W*V2); alpar@10: } alpar@10: /// A kind of two dimensional exponential distribution alpar@10: alpar@10: /// This function provides a turning symmetric two-dimensional distribution. alpar@10: /// The x-coordinate is of conditionally exponential distribution alpar@10: /// with the condition that x is positive and y=0. If x is negative and alpar@10: /// y=0 then, -x is of exponential distribution. The same is true for the alpar@10: /// y-coordinate. alpar@10: dim2::Point exponential2() alpar@10: { alpar@10: double V1,V2,S; alpar@10: do { alpar@10: V1=2*real()-1; alpar@10: V2=2*real()-1; alpar@10: S=V1*V1+V2*V2; alpar@10: } while(S>=1); alpar@10: double W=-std::log(S)/S; alpar@10: return dim2::Point(W*V1,W*V2); alpar@10: } alpar@10: alpar@10: ///@} alpar@10: }; alpar@10: alpar@10: alpar@10: extern Random rnd; alpar@10: alpar@10: } alpar@10: alpar@10: #endif