deba@2229: /* -*- C++ -*- deba@2229: * deba@2229: * This file is a part of LEMON, a generic C++ optimization library deba@2229: * alpar@2391: * Copyright (C) 2003-2007 deba@2229: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport deba@2229: * (Egervary Research Group on Combinatorial Optimization, EGRES). deba@2229: * deba@2229: * Permission to use, modify and distribute this software is granted deba@2229: * provided that this copyright notice appears in all copies. For deba@2229: * precise terms see the accompanying LICENSE file. deba@2229: * deba@2229: * This software is provided "AS IS" with no warranty of any kind, deba@2229: * express or implied, and with no claim as to its suitability for any deba@2229: * purpose. deba@2229: * deba@2229: */ deba@2229: deba@2372: /* deba@2372: * This file contains the reimplemented version of the Mersenne Twister deba@2372: * Generator of Matsumoto and Nishimura. deba@2372: * deba@2372: * See the appropriate copyright notice below. deba@2372: * deba@2372: * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, deba@2372: * All rights reserved. deba@2372: * deba@2372: * Redistribution and use in source and binary forms, with or without deba@2372: * modification, are permitted provided that the following conditions deba@2372: * are met: deba@2372: * deba@2372: * 1. Redistributions of source code must retain the above copyright deba@2372: * notice, this list of conditions and the following disclaimer. deba@2372: * deba@2372: * 2. Redistributions in binary form must reproduce the above copyright deba@2372: * notice, this list of conditions and the following disclaimer in the deba@2372: * documentation and/or other materials provided with the distribution. deba@2372: * deba@2372: * 3. The names of its contributors may not be used to endorse or promote deba@2372: * products derived from this software without specific prior written deba@2372: * permission. deba@2372: * deba@2372: * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS deba@2372: * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT deba@2372: * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS deba@2372: * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE deba@2372: * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, deba@2372: * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES deba@2372: * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR deba@2372: * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) deba@2372: * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, deba@2372: * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) deba@2372: * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED deba@2372: * OF THE POSSIBILITY OF SUCH DAMAGE. deba@2372: * deba@2372: * deba@2372: * Any feedback is very welcome. deba@2372: * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html deba@2372: * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) deba@2372: */ deba@2372: deba@2229: #ifndef LEMON_RANDOM_H deba@2229: #define LEMON_RANDOM_H deba@2229: deba@2229: #include deba@2242: #include deba@2242: #include deba@2229: deba@2229: #include deba@2229: #include deba@2229: alpar@2374: #include deba@2229: ///\ingroup misc deba@2229: ///\file deba@2229: ///\brief Mersenne Twister random number generator deba@2229: /// deba@2229: ///\author Balazs Dezso deba@2229: deba@2229: namespace lemon { deba@2229: deba@2242: namespace _random_bits { deba@2242: deba@2242: template ::digits> deba@2242: struct RandomTraits {}; deba@2242: deba@2242: template deba@2242: struct RandomTraits<_Word, 32> { deba@2242: deba@2242: typedef _Word Word; deba@2242: static const int bits = 32; deba@2242: deba@2242: static const int length = 624; deba@2242: static const int shift = 397; deba@2242: deba@2242: static const Word mul = 0x6c078965u; deba@2242: static const Word arrayInit = 0x012BD6AAu; deba@2242: static const Word arrayMul1 = 0x0019660Du; deba@2242: static const Word arrayMul2 = 0x5D588B65u; deba@2242: deba@2242: static const Word mask = 0x9908B0DFu; deba@2242: static const Word loMask = (1u << 31) - 1; deba@2242: static const Word hiMask = ~loMask; deba@2242: deba@2242: deba@2242: static Word tempering(Word rnd) { deba@2242: rnd ^= (rnd >> 11); deba@2242: rnd ^= (rnd << 7) & 0x9D2C5680u; deba@2242: rnd ^= (rnd << 15) & 0xEFC60000u; deba@2242: rnd ^= (rnd >> 18); deba@2242: return rnd; deba@2242: } deba@2242: deba@2242: }; deba@2242: deba@2242: template deba@2242: struct RandomTraits<_Word, 64> { deba@2242: deba@2242: typedef _Word Word; deba@2242: static const int bits = 64; deba@2242: deba@2242: static const int length = 312; deba@2242: static const int shift = 156; deba@2242: deba@2386: static const Word mul = Word(0x5851F42Du) << 32 | Word(0x4C957F2Du); deba@2386: static const Word arrayInit = Word(0x00000000u) << 32 |Word(0x012BD6AAu); deba@2386: static const Word arrayMul1 = Word(0x369DEA0Fu) << 32 |Word(0x31A53F85u); deba@2386: static const Word arrayMul2 = Word(0x27BB2EE6u) << 32 |Word(0x87B0B0FDu); deba@2242: deba@2386: static const Word mask = Word(0xB5026F5Au) << 32 | Word(0xA96619E9u); deba@2386: static const Word loMask = (Word(1u) << 31) - 1; deba@2242: static const Word hiMask = ~loMask; deba@2242: deba@2242: static Word tempering(Word rnd) { deba@2386: rnd ^= (rnd >> 29) & (Word(0x55555555u) << 32 | Word(0x55555555u)); deba@2386: rnd ^= (rnd << 17) & (Word(0x71D67FFFu) << 32 | Word(0xEDA60000u)); deba@2386: rnd ^= (rnd << 37) & (Word(0xFFF7EEE0u) << 32 | Word(0x00000000u)); deba@2242: rnd ^= (rnd >> 43); deba@2242: return rnd; deba@2242: } deba@2242: deba@2242: }; deba@2242: deba@2242: template deba@2242: class RandomCore { deba@2242: public: deba@2242: deba@2242: typedef _Word Word; deba@2242: deba@2242: private: deba@2242: deba@2242: static const int bits = RandomTraits::bits; deba@2242: deba@2242: static const int length = RandomTraits::length; deba@2242: static const int shift = RandomTraits::shift; deba@2242: deba@2242: public: deba@2242: deba@2242: void initState() { deba@2242: static const Word seedArray[4] = { deba@2242: 0x12345u, 0x23456u, 0x34567u, 0x45678u deba@2242: }; deba@2242: deba@2242: initState(seedArray, seedArray + 4); deba@2242: } deba@2242: deba@2242: void initState(Word seed) { deba@2242: deba@2242: static const Word mul = RandomTraits::mul; deba@2242: deba@2242: current = state; deba@2242: deba@2242: Word *curr = state + length - 1; deba@2242: curr[0] = seed; --curr; deba@2242: for (int i = 1; i < length; ++i) { deba@2242: curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i); deba@2242: --curr; deba@2242: } deba@2242: } deba@2242: deba@2242: template deba@2242: void initState(Iterator begin, Iterator end) { deba@2242: deba@2242: static const Word init = RandomTraits::arrayInit; deba@2242: static const Word mul1 = RandomTraits::arrayMul1; deba@2242: static const Word mul2 = RandomTraits::arrayMul2; deba@2242: deba@2242: deba@2242: Word *curr = state + length - 1; --curr; deba@2242: Iterator it = begin; int cnt = 0; deba@2242: int num; deba@2242: deba@2242: initState(init); deba@2242: deba@2242: num = length > end - begin ? length : end - begin; deba@2242: while (num--) { deba@2242: curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1)) deba@2242: + *it + cnt; deba@2242: ++it; ++cnt; deba@2242: if (it == end) { deba@2242: it = begin; cnt = 0; deba@2242: } deba@2242: if (curr == state) { deba@2242: curr = state + length - 1; curr[0] = state[0]; deba@2242: } deba@2242: --curr; deba@2242: } deba@2242: deba@2242: num = length - 1; cnt = length - (curr - state) - 1; deba@2242: while (num--) { deba@2242: curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2)) deba@2242: - cnt; deba@2242: --curr; ++cnt; deba@2242: if (curr == state) { deba@2242: curr = state + length - 1; curr[0] = state[0]; --curr; deba@2242: cnt = 1; deba@2242: } deba@2242: } deba@2242: deba@2386: state[length - 1] = Word(1) << (bits - 1); deba@2242: } deba@2242: deba@2242: void copyState(const RandomCore& other) { deba@2242: std::copy(other.state, other.state + length, state); deba@2242: current = state + (other.current - other.state); deba@2242: } deba@2242: deba@2242: Word operator()() { deba@2242: if (current == state) fillState(); deba@2242: --current; deba@2242: Word rnd = *current; deba@2242: return RandomTraits::tempering(rnd); deba@2242: } deba@2242: deba@2242: private: deba@2242: deba@2242: deba@2242: void fillState() { deba@2242: static const Word mask[2] = { 0x0ul, RandomTraits::mask }; deba@2242: static const Word loMask = RandomTraits::loMask; deba@2242: static const Word hiMask = RandomTraits::hiMask; deba@2242: deba@2242: current = state + length; deba@2242: deba@2242: register Word *curr = state + length - 1; deba@2242: register long num; deba@2242: deba@2242: num = length - shift; deba@2242: while (num--) { deba@2242: curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ deba@2242: curr[- shift] ^ mask[curr[-1] & 1ul]; deba@2242: --curr; deba@2242: } deba@2242: num = shift - 1; deba@2242: while (num--) { deba@2242: curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ deba@2242: curr[length - shift] ^ mask[curr[-1] & 1ul]; deba@2242: --curr; deba@2242: } deba@2242: curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^ deba@2242: curr[length - shift] ^ mask[curr[length - 1] & 1ul]; deba@2242: deba@2242: } deba@2242: deba@2242: deba@2242: Word *current; deba@2242: Word state[length]; deba@2242: deba@2242: }; deba@2242: deba@2242: deba@2242: template ::digits + 1) / 2> deba@2242: struct Masker { deba@2242: static Result mask(const Result& result) { deba@2242: return Masker:: deba@2386: mask(static_cast(result | (result >> shift))); deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct Masker { deba@2242: static Result mask(const Result& result) { deba@2386: return static_cast(result | (result >> 1)); deba@2242: } deba@2242: }; deba@2242: deba@2242: template ::digits, int shift = 0, deba@2242: bool last = rest <= std::numeric_limits::digits> deba@2242: struct IntConversion { deba@2242: static const int bits = std::numeric_limits::digits; deba@2242: deba@2242: static Result convert(RandomCore& rnd) { deba@2386: return static_cast(rnd() >> (bits - rest)) << shift; deba@2242: } deba@2242: deba@2242: }; deba@2242: deba@2242: template deba@2242: struct IntConversion { deba@2242: static const int bits = std::numeric_limits::digits; deba@2242: deba@2242: static Result convert(RandomCore& rnd) { deba@2386: return (static_cast(rnd()) << shift) | deba@2242: IntConversion::convert(rnd); deba@2242: } deba@2242: }; deba@2242: deba@2242: deba@2242: template ::digits < alpar@2387: std::numeric_limits::digits) > deba@2242: struct Mapping { deba@2242: static Result map(RandomCore& rnd, const Result& bound) { deba@2386: Word max = Word(bound - 1); deba@2242: Result mask = Masker::mask(bound - 1); deba@2242: Result num; deba@2242: do { deba@2242: num = IntConversion::convert(rnd) & mask; deba@2242: } while (num > max); deba@2242: return num; deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct Mapping { deba@2242: static Result map(RandomCore& rnd, const Result& bound) { deba@2386: Word max = Word(bound - 1); deba@2242: Word mask = Masker::digits + 1) / 2> deba@2242: ::mask(max); deba@2242: Word num; deba@2242: do { deba@2242: num = rnd() & mask; deba@2242: } while (num > max); deba@2242: return num; deba@2242: } deba@2242: }; deba@2242: deba@2242: template = 0)> deba@2242: struct ShiftMultiplier { deba@2242: static const Result multiplier() { deba@2242: Result res = ShiftMultiplier::multiplier(); deba@2242: res *= res; deba@2386: if ((exp & 1) == 1) res *= static_cast(2.0); deba@2242: return res; deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct ShiftMultiplier { deba@2242: static const Result multiplier() { deba@2242: Result res = ShiftMultiplier::multiplier(); deba@2242: res *= res; deba@2386: if ((exp & 1) == 1) res *= static_cast(0.5); deba@2242: return res; deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct ShiftMultiplier { deba@2242: static const Result multiplier() { deba@2386: return static_cast(1.0); deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct ShiftMultiplier { deba@2242: static const Result multiplier() { deba@2386: return static_cast(1.0/1048576.0); deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct ShiftMultiplier { deba@2242: static const Result multiplier() { deba@2386: return static_cast(1.0/424967296.0); deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct ShiftMultiplier { deba@2242: static const Result multiplier() { deba@2386: return static_cast(1.0/9007199254740992.0); deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct ShiftMultiplier { deba@2242: static const Result multiplier() { deba@2386: return static_cast(1.0/18446744073709551616.0); deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct Shifting { deba@2242: static Result shift(const Result& result) { deba@2242: return result * ShiftMultiplier::multiplier(); deba@2242: } deba@2242: }; deba@2242: deba@2242: template ::digits, int shift = 0, deba@2242: bool last = rest <= std::numeric_limits::digits> deba@2242: struct RealConversion{ deba@2242: static const int bits = std::numeric_limits::digits; deba@2242: deba@2242: static Result convert(RandomCore& rnd) { deba@2242: return Shifting:: deba@2386: shift(static_cast(rnd() >> (bits - rest))); deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct RealConversion { deba@2242: static const int bits = std::numeric_limits::digits; deba@2242: deba@2242: static Result convert(RandomCore& rnd) { deba@2386: return Shifting:: deba@2386: shift(static_cast(rnd())) + deba@2386: RealConversion:: deba@2386: convert(rnd); deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct Initializer { deba@2242: deba@2242: template deba@2242: static void init(RandomCore& rnd, Iterator begin, Iterator end) { deba@2242: std::vector ws; deba@2242: for (Iterator it = begin; it != end; ++it) { deba@2386: ws.push_back(Word(*it)); deba@2242: } deba@2242: rnd.initState(ws.begin(), ws.end()); deba@2242: } deba@2242: deba@2242: static void init(RandomCore& rnd, Result seed) { deba@2242: rnd.initState(seed); deba@2242: } deba@2242: }; deba@2242: deba@2242: template deba@2242: struct BoolConversion { deba@2242: static bool convert(RandomCore& rnd) { deba@2242: return (rnd() & 1) == 1; deba@2242: } deba@2242: }; deba@2242: deba@2372: template deba@2372: struct BoolProducer { deba@2372: Word buffer; deba@2372: int num; deba@2372: deba@2372: BoolProducer() : num(0) {} deba@2372: deba@2372: bool convert(RandomCore& rnd) { deba@2372: if (num == 0) { deba@2372: buffer = rnd(); deba@2372: num = RandomTraits::bits; deba@2372: } deba@2372: bool r = (buffer & 1); deba@2372: buffer >>= 1; deba@2372: --num; deba@2372: return r; deba@2372: } deba@2372: }; deba@2372: deba@2242: } deba@2229: deba@2229: /// \ingroup misc deba@2229: /// deba@2229: /// \brief Mersenne Twister random number generator deba@2229: /// deba@2229: /// The Mersenne Twister is a twisted generalized feedback deba@2242: /// shift-register generator of Matsumoto and Nishimura. The period deba@2242: /// of this generator is \f$ 2^{19937} - 1 \f$ and it is deba@2242: /// equi-distributed in 623 dimensions for 32-bit numbers. The time deba@2242: /// performance of this generator is comparable to the commonly used deba@2242: /// generators. deba@2242: /// deba@2242: /// This implementation is specialized for both 32-bit and 64-bit deba@2242: /// architectures. The generators differ sligthly in the deba@2242: /// initialization and generation phase so they produce two deba@2242: /// completly different sequences. deba@2242: /// deba@2242: /// The generator gives back random numbers of serveral types. To deba@2242: /// get a random number from a range of a floating point type you deba@2245: /// can use one form of the \c operator() or the \c real() member deba@2245: /// function. If you want to get random number from the {0, 1, ..., deba@2245: /// n-1} integer range use the \c operator[] or the \c integer() deba@2245: /// method. And to get random number from the whole range of an deba@2245: /// integer type you can use the argumentless \c integer() or \c deba@2245: /// uinteger() functions. After all you can get random bool with deba@2245: /// equal chance of true and false or given probability of true deba@2245: /// result with the \c boolean() member functions. deba@2242: /// deba@2242: ///\code deba@2245: /// // The commented code is identical to the other deba@2245: /// double a = rnd(); // [0.0, 1.0) deba@2245: /// // double a = rnd.real(); // [0.0, 1.0) deba@2245: /// double b = rnd(100.0); // [0.0, 100.0) deba@2245: /// // double b = rnd.real(100.0); // [0.0, 100.0) deba@2245: /// double c = rnd(1.0, 2.0); // [1.0, 2.0) deba@2245: /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0) deba@2245: /// int d = rnd[100000]; // 0..99999 deba@2245: /// // int d = rnd.integer(100000); // 0..99999 deba@2245: /// int e = rnd[6] + 1; // 1..6 deba@2245: /// // int e = rnd.integer(1, 1 + 6); // 1..6 deba@2242: /// int b = rnd.uinteger(); // 0 .. 2^31 - 1 deba@2242: /// int c = rnd.integer(); // - 2^31 .. 2^31 - 1 deba@2242: /// bool g = rnd.boolean(); // P(g = true) = 0.5 deba@2242: /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 deba@2242: ///\endcode deba@2242: /// deba@2245: /// The lemon provides a global instance of the random number deba@2245: /// generator which name is \ref lemon::rnd "rnd". Usually it is a deba@2245: /// good programming convenience to use this global generator to get deba@2245: /// random numbers. deba@2229: /// deba@2229: /// \author Balazs Dezso deba@2229: class Random { deba@2242: private: deba@2229: deba@2245: // architecture word deba@2242: typedef unsigned long Word; deba@2242: deba@2242: _random_bits::RandomCore core; deba@2372: _random_bits::BoolProducer bool_producer; deba@2372: deba@2229: deba@2229: public: deba@2229: deba@2229: /// \brief Constructor deba@2229: /// deba@2242: /// Constructor with constant seeding. deba@2242: Random() { core.initState(); } deba@2229: deba@2229: /// \brief Constructor deba@2229: /// deba@2242: /// Constructor with seed. The current number type will be converted deba@2242: /// to the architecture word type. deba@2242: template deba@2242: Random(Number seed) { deba@2242: _random_bits::Initializer::init(core, seed); deba@2242: } deba@2242: deba@2242: /// \brief Constructor deba@2242: /// deba@2242: /// Constructor with array seeding. The given range should contain deba@2242: /// any number type and the numbers will be converted to the deba@2242: /// architecture word type. deba@2242: template deba@2242: Random(Iterator begin, Iterator end) { deba@2242: typedef typename std::iterator_traits::value_type Number; deba@2242: _random_bits::Initializer::init(core, begin, end); deba@2242: } deba@2229: deba@2229: /// \brief Copy constructor deba@2229: /// deba@2229: /// Copy constructor. The generated sequence will be identical to deba@2245: /// the other sequence. It can be used to save the current state deba@2245: /// of the generator and later use it to generate the same deba@2245: /// sequence. deba@2242: Random(const Random& other) { deba@2242: core.copyState(other.core); deba@2229: } deba@2229: deba@2229: /// \brief Assign operator deba@2229: /// deba@2229: /// Assign operator. The generated sequence will be identical to deba@2245: /// the other sequence. It can be used to save the current state deba@2245: /// of the generator and later use it to generate the same deba@2245: /// sequence. deba@2229: Random& operator=(const Random& other) { deba@2229: if (&other != this) { deba@2242: core.copyState(other.core); deba@2229: } deba@2229: return *this; deba@2229: } deba@2229: alpar@2257: /// \brief Returns a random real number from the range [0, 1) deba@2229: /// deba@2245: /// It returns a random real number from the range [0, 1). The deba@2245: /// default Number type is double. deba@2242: template deba@2245: Number real() { deba@2242: return _random_bits::RealConversion::convert(core); deba@2229: } deba@2229: deba@2245: double real() { deba@2245: return real(); deba@2245: } deba@2245: alpar@2257: /// \brief Returns a random real number the range [0, b) deba@2245: /// deba@2245: /// It returns a random real number from the range [0, b). deba@2245: template deba@2245: Number real(Number b) { deba@2245: return real() * b; deba@2245: } deba@2245: alpar@2257: /// \brief Returns a random real number from the range [a, b) deba@2245: /// deba@2245: /// It returns a random real number from the range [a, b). deba@2245: template deba@2245: Number real(Number a, Number b) { deba@2245: return real() * (b - a) + a; deba@2245: } deba@2245: alpar@2257: /// \brief Returns a random real number from the range [0, 1) deba@2245: /// deba@2245: /// It returns a random double from the range [0, 1). deba@2242: double operator()() { deba@2245: return real(); deba@2242: } deba@2242: alpar@2257: /// \brief Returns a random real number from the range [0, b) deba@2229: /// deba@2242: /// It returns a random real number from the range [0, b). deba@2242: template deba@2242: Number operator()(Number b) { deba@2245: return real() * b; deba@2242: } deba@2242: alpar@2257: /// \brief Returns a random real number from the range [a, b) deba@2242: /// deba@2242: /// It returns a random real number from the range [a, b). deba@2242: template deba@2242: Number operator()(Number a, Number b) { deba@2245: return real() * (b - a) + a; deba@2242: } deba@2242: deba@2242: /// \brief Returns a random integer from a range deba@2242: /// deba@2245: /// It returns a random integer from the range {0, 1, ..., b - 1}. deba@2242: template deba@2245: Number integer(Number b) { deba@2245: return _random_bits::Mapping::map(core, b); deba@2245: } deba@2245: deba@2245: /// \brief Returns a random integer from a range deba@2245: /// deba@2245: /// It returns a random integer from the range {a, a + 1, ..., b - 1}. deba@2245: template deba@2245: Number integer(Number a, Number b) { deba@2245: return _random_bits::Mapping::map(core, b - a) + a; deba@2245: } deba@2245: deba@2245: /// \brief Returns a random integer from a range deba@2245: /// deba@2245: /// It returns a random integer from the range {0, 1, ..., b - 1}. deba@2245: template deba@2245: Number operator[](Number b) { deba@2245: return _random_bits::Mapping::map(core, b); deba@2242: } deba@2242: deba@2242: /// \brief Returns a random non-negative integer deba@2242: /// deba@2242: /// It returns a random non-negative integer uniformly from the deba@2242: /// whole range of the current \c Number type. The default result deba@2242: /// type of this function is unsigned int. deba@2242: template deba@2242: Number uinteger() { deba@2242: return _random_bits::IntConversion::convert(core); deba@2242: } deba@2242: deba@2242: unsigned int uinteger() { deba@2242: return uinteger(); deba@2242: } deba@2242: deba@2242: /// \brief Returns a random integer deba@2242: /// deba@2242: /// It returns a random integer uniformly from the whole range of deba@2242: /// the current \c Number type. The default result type of this deba@2242: /// function is int. deba@2242: template deba@2242: Number integer() { deba@2242: static const int nb = std::numeric_limits::digits + deba@2242: (std::numeric_limits::is_signed ? 1 : 0); deba@2242: return _random_bits::IntConversion::convert(core); deba@2242: } deba@2242: deba@2242: int integer() { deba@2242: return integer(); deba@2229: } deba@2229: deba@2242: /// \brief Returns a random bool deba@2229: /// deba@2372: /// It returns a random bool. The generator holds a buffer for deba@2372: /// random bits. Every time when it become empty the generator makes deba@2372: /// a new random word and fill the buffer up. deba@2242: bool boolean() { deba@2372: return bool_producer.convert(core); deba@2229: } deba@2229: alpar@2356: ///\name Nonuniform distributions alpar@2356: /// alpar@2356: alpar@2356: ///@{ alpar@2356: deba@2229: /// \brief Returns a random bool deba@2229: /// deba@2242: /// It returns a random bool with given probability of true result deba@2242: bool boolean(double p) { deba@2242: return operator()() < p; deba@2229: } alpar@2355: alpar@2355: /// Standard Gauss distribution alpar@2355: alpar@2355: /// Standard Gauss distribution. alpar@2356: /// \note The Cartesian form of the Box-Muller alpar@2356: /// transformation is used to generate a random normal distribution. alpar@2356: /// \todo Consider using the "ziggurat" method instead. alpar@2355: double gauss() alpar@2355: { alpar@2355: double V1,V2,S; alpar@2355: do { alpar@2355: V1=2*real()-1; alpar@2355: V2=2*real()-1; alpar@2355: S=V1*V1+V2*V2; alpar@2355: } while(S>=1); alpar@2355: return std::sqrt(-2*std::log(S)/S)*V1; alpar@2355: } alpar@2356: /// Gauss distribution with given standard deviation and mean 0 alpar@2356: alpar@2356: /// \sa gauss() alpar@2356: /// alpar@2356: double gauss(double std_dev) alpar@2355: { alpar@2356: return gauss()*std_dev; alpar@2355: } alpar@2356: /// Gauss distribution with given mean and standard deviation alpar@2356: alpar@2356: /// \sa gauss() alpar@2356: /// alpar@2356: double gauss(double mean,double std_dev) alpar@2355: { alpar@2356: return gauss()*std_dev+mean; alpar@2355: } alpar@2355: alpar@2356: /// Exponential distribution with given mean alpar@2356: alpar@2356: /// This function generates an exponential distribution random number alpar@2356: /// with mean 1/lambda. alpar@2356: /// alpar@2356: double exponential(double lambda=1.0) alpar@2355: { alpar@2374: return -std::log(real())/lambda; alpar@2355: } alpar@2356: alpar@2483: double gamma(int k) alpar@2483: { alpar@2483: double s = 0; alpar@2483: for(int i=0;i()); alpar@2483: return s; alpar@2483: } alpar@2483: alpar@2483: /// Gamma distribution with given shape and scale parameter alpar@2483: alpar@2483: /// This function generates a gamma distribution random number. alpar@2483: /// alpar@2483: ///\param k shape parameter (k>0) alpar@2483: ///\param theta scale parameter alpar@2483: /// alpar@2483: double gamma(double k,double theta=1.0) alpar@2483: { alpar@2483: double xi,nu; alpar@2483: const double delta = k-std::floor(k); alpar@2483: const double v0=M_E/(M_E-delta); alpar@2483: do { alpar@2483: double V0=1.0-real(); alpar@2483: double V1=1.0-real(); alpar@2483: double V2=1.0-real(); alpar@2483: if(V2<=v0) alpar@2483: { alpar@2483: xi=std::pow(V1,1.0/delta); alpar@2483: nu=V0*std::pow(xi,delta-1.0); alpar@2483: } alpar@2483: else alpar@2483: { alpar@2483: xi=1.0-std::log(V1); alpar@2483: nu=V0*std::exp(-xi); alpar@2483: } alpar@2483: } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); alpar@2483: return theta*(xi-gamma(int(std::floor(k)))); alpar@2483: } alpar@2483: alpar@2483: alpar@2356: ///@} deba@2229: alpar@2374: ///\name Two dimensional distributions alpar@2374: /// alpar@2374: alpar@2374: ///@{ alpar@2374: alpar@2374: /// Uniform distribution on the full unit circle. alpar@2380: dim2::Point disc() alpar@2374: { alpar@2374: double V1,V2; alpar@2374: do { alpar@2374: V1=2*real()-1; alpar@2374: V2=2*real()-1; alpar@2374: alpar@2374: } while(V1*V1+V2*V2>=1); alpar@2374: return dim2::Point(V1,V2); alpar@2374: } alpar@2374: /// A kind of two dimensional Gauss distribution alpar@2374: alpar@2374: /// This function provides a turning symmetric two-dimensional distribution. alpar@2374: /// Both coordinates are of standard normal distribution, but they are not alpar@2374: /// independent. alpar@2374: /// alpar@2374: /// \note The coordinates are the two random variables provided by alpar@2374: /// the Box-Muller method. alpar@2374: dim2::Point gauss2() alpar@2374: { alpar@2374: double V1,V2,S; alpar@2374: do { alpar@2374: V1=2*real()-1; alpar@2374: V2=2*real()-1; alpar@2374: S=V1*V1+V2*V2; alpar@2374: } while(S>=1); alpar@2374: double W=std::sqrt(-2*std::log(S)/S); alpar@2374: return dim2::Point(W*V1,W*V2); alpar@2374: } alpar@2374: /// A kind of two dimensional exponential distribution alpar@2374: alpar@2374: /// This function provides a turning symmetric two-dimensional distribution. alpar@2374: /// The x-coordinate is of conditionally exponential distribution alpar@2374: /// with the condition that x is positive and y=0. If x is negative and alpar@2374: /// y=0 then, -x is of exponential distribution. The same is true for the alpar@2374: /// y-coordinate. alpar@2374: dim2::Point exponential2() alpar@2374: { alpar@2374: double V1,V2,S; alpar@2374: do { alpar@2374: V1=2*real()-1; alpar@2374: V2=2*real()-1; alpar@2374: S=V1*V1+V2*V2; alpar@2374: } while(S>=1); alpar@2374: double W=-std::log(S)/S; alpar@2374: return dim2::Point(W*V1,W*V2); alpar@2374: } alpar@2374: alpar@2374: ///@} deba@2229: }; deba@2229: deba@2229: deba@2229: extern Random rnd; deba@2229: deba@2229: } deba@2229: deba@2229: #endif