# HG changeset patch # User Alpar Juttner # Date 1521816210 -3600 # Node ID 57167d92e96c6c52005817dbd523ae36805a60eb # Parent 0fdf84c79bc192625fef343812123141c8c64c2d# Parent 04f57dad1b07e03dc3da650fad31c9b08622bc06 Merge #602 diff -r 0fdf84c79bc1 -r 57167d92e96c lemon/random.h --- a/lemon/random.h Fri Mar 23 15:39:54 2018 +0100 +++ b/lemon/random.h Fri Mar 23 15:43:30 2018 +0100 @@ -111,7 +111,6 @@ static const Word loMask = (1u << 31) - 1; static const Word hiMask = ~loMask; - static Word tempering(Word rnd) { rnd ^= (rnd >> 11); rnd ^= (rnd << 7) & 0x9D2C5680u; @@ -243,7 +242,6 @@ private: - void fillState() { static const Word mask[2] = { 0x0ul, RandomTraits::mask }; static const Word loMask = RandomTraits::loMask; @@ -271,7 +269,6 @@ } - Word *current; Word state[length]; @@ -296,7 +293,7 @@ template ::digits, int shift = 0, - bool last = rest <= std::numeric_limits::digits> + bool last = (rest <= std::numeric_limits::digits)> struct IntConversion { static const int bits = std::numeric_limits::digits; @@ -465,543 +462,594 @@ } }; - } + /// \ingroup misc + /// + /// \brief Mersenne Twister random number generator + /// + /// The Mersenne Twister is a twisted generalized feedback + /// shift-register generator of Matsumoto and Nishimura. The period + /// of this generator is \f$ 2^{19937} - 1\f$ and it is + /// equi-distributed in 623 dimensions for 32-bit numbers. The time + /// performance of this generator is comparable to the commonly used + /// generators. + /// + /// This is a template implementation of both 32-bit and + /// 64-bit architecture optimized versions. The generators differ + /// sligthly in the initialization and generation phase so they + /// produce two completly different sequences. + /// + /// \alert Do not use this class directly, but instead one of \ref + /// Random, \ref Random32 or \ref Random64. + /// + /// The generator gives back random numbers of serveral types. To + /// get a random number from a range of a floating point type, you + /// can use one form of the \c operator() or the \c real() member + /// function. If you want to get random number from the {0, 1, ..., + /// n-1} integer range, use the \c operator[] or the \c integer() + /// method. And to get random number from the whole range of an + /// integer type, you can use the argumentless \c integer() or + /// \c uinteger() functions. Finally, you can get random bool with + /// equal chance of true and false or with given probability of true + /// result using the \c boolean() member functions. + /// + /// Various non-uniform distributions are also supported: normal (Gauss), + /// exponential, gamma, Poisson, etc.; and a few two-dimensional + /// distributions, too. + /// + ///\code + /// // The commented code is identical to the other + /// double a = rnd(); // [0.0, 1.0) + /// // double a = rnd.real(); // [0.0, 1.0) + /// double b = rnd(100.0); // [0.0, 100.0) + /// // double b = rnd.real(100.0); // [0.0, 100.0) + /// double c = rnd(1.0, 2.0); // [1.0, 2.0) + /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0) + /// int d = rnd[100000]; // 0..99999 + /// // int d = rnd.integer(100000); // 0..99999 + /// int e = rnd[6] + 1; // 1..6 + /// // int e = rnd.integer(1, 1 + 6); // 1..6 + /// int b = rnd.uinteger(); // 0 .. 2^31 - 1 + /// int c = rnd.integer(); // - 2^31 .. 2^31 - 1 + /// bool g = rnd.boolean(); // P(g = true) = 0.5 + /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 + ///\endcode + /// + /// LEMON provides a global instance of the random number generator: + /// \ref lemon::rnd "rnd". In most cases, it is a good practice + /// to use this global generator to get random numbers. + /// + /// \sa \ref Random, \ref Random32 or \ref Random64. + template + class Random { + private: + + _random_bits::RandomCore core; + _random_bits::BoolProducer bool_producer; + + + public: + + ///\name Initialization + /// + /// @{ + + /// \brief Default constructor + /// + /// Constructor with constant seeding. + Random() { core.initState(); } + + /// \brief Constructor with seed + /// + /// Constructor with seed. The current number type will be converted + /// to the architecture word type. + template + Random(Number seed) { + _random_bits::Initializer::init(core, seed); + } + + /// \brief Constructor with array seeding + /// + /// Constructor with array seeding. The given range should contain + /// any number type and the numbers will be converted to the + /// architecture word type. + template + Random(Iterator begin, Iterator end) { + typedef typename std::iterator_traits::value_type Number; + _random_bits::Initializer::init(core, begin, end); + } + + /// \brief Copy constructor + /// + /// Copy constructor. The generated sequence will be identical to + /// the other sequence. It can be used to save the current state + /// of the generator and later use it to generate the same + /// sequence. + Random(const Random& other) { + core.copyState(other.core); + } + + /// \brief Assign operator + /// + /// Assign operator. The generated sequence will be identical to + /// the other sequence. It can be used to save the current state + /// of the generator and later use it to generate the same + /// sequence. + Random& operator=(const Random& other) { + if (&other != this) { + core.copyState(other.core); + } + return *this; + } + + /// \brief Seeding random sequence + /// + /// Seeding the random sequence. The current number type will be + /// converted to the architecture word type. + template + void seed(Number seed) { + _random_bits::Initializer::init(core, seed); + } + + /// \brief Seeding random sequence + /// + /// Seeding the random sequence. The given range should contain + /// any number type and the numbers will be converted to the + /// architecture word type. + template + void seed(Iterator begin, Iterator end) { + typedef typename std::iterator_traits::value_type Number; + _random_bits::Initializer::init(core, begin, end); + } + + /// \brief Seeding from file or from process id and time + /// + /// By default, this function calls the \c seedFromFile() member + /// function with the /dev/urandom file. If it does not success, + /// it uses the \c seedFromTime(). + /// \return Currently always \c true. + bool seed() { +#ifndef LEMON_WIN32 + if (seedFromFile("/dev/urandom", 0)) return true; +#endif + if (seedFromTime()) return true; + return false; + } + + /// \brief Seeding from file + /// + /// Seeding the random sequence from file. The linux kernel has two + /// devices, /dev/random and /dev/urandom which + /// could give good seed values for pseudo random generators (The + /// difference between two devices is that the random may + /// block the reading operation while the kernel can give good + /// source of randomness, while the urandom does not + /// block the input, but it could give back bytes with worse + /// entropy). + /// \param file The source file + /// \param offset The offset, from the file read. + /// \return \c true when the seeding successes. +#ifndef LEMON_WIN32 + bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0) +#else + bool seedFromFile(const std::string& file = "", int offset = 0) +#endif + { + std::ifstream rs(file.c_str()); + const int size = 4; + Word buf[size]; + if (offset != 0 && !rs.seekg(offset)) return false; + if (!rs.read(reinterpret_cast(buf), sizeof(buf))) return false; + seed(buf, buf + size); + return true; + } + + /// \brief Seeding from process id and time + /// + /// Seeding from process id and time. This function uses the + /// current process id and the current time for initialize the + /// random sequence. + /// \return Currently always \c true. + bool seedFromTime() { +#ifndef LEMON_WIN32 + timeval tv; + gettimeofday(&tv, 0); + seed(getpid() + tv.tv_sec + tv.tv_usec); +#else + seed(bits::getWinRndSeed()); +#endif + return true; + } + + /// @} + + ///\name Uniform Distributions + /// + /// @{ + + /// \brief Returns a random real number from the range [0, 1) + /// + /// It returns a random real number from the range [0, 1). The + /// default Number type is \c double. + template + Number real() { + return _random_bits::RealConversion::convert(core); + } + + double real() { + return real(); + } + + /// \brief Returns a random real number from the range [0, 1) + /// + /// It returns a random double from the range [0, 1). + double operator()() { + return real(); + } + + /// \brief Returns a random real number from the range [0, b) + /// + /// It returns a random real number from the range [0, b). + double operator()(double b) { + return real() * b; + } + + /// \brief Returns a random real number from the range [a, b) + /// + /// It returns a random real number from the range [a, b). + double operator()(double a, double b) { + return real() * (b - a) + a; + } + + /// \brief Returns a random integer from a range + /// + /// It returns a random integer from the range {0, 1, ..., b - 1}. + template + Number integer(Number b) { + return _random_bits::Mapping::map(core, b); + } + + /// \brief Returns a random integer from a range + /// + /// It returns a random integer from the range {a, a + 1, ..., b - 1}. + template + Number integer(Number a, Number b) { + return _random_bits::Mapping::map(core, b - a) + a; + } + + /// \brief Returns a random integer from a range + /// + /// It returns a random integer from the range {0, 1, ..., b - 1}. + template + Number operator[](Number b) { + return _random_bits::Mapping::map(core, b); + } + + /// \brief Returns a random non-negative integer + /// + /// It returns a random non-negative integer uniformly from the + /// whole range of the current \c Number type. The default result + /// type of this function is unsigned int. + template + Number uinteger() { + return _random_bits::IntConversion::convert(core); + } + + unsigned int uinteger() { + return uinteger(); + } + + /// \brief Returns a random integer + /// + /// It returns a random integer uniformly from the whole range of + /// the current \c Number type. The default result type of this + /// function is \c int. + template + Number integer() { + static const int nb = std::numeric_limits::digits + + (std::numeric_limits::is_signed ? 1 : 0); + return _random_bits::IntConversion::convert(core); + } + + int integer() { + return integer(); + } + + /// \brief Returns a random bool + /// + /// It returns a random bool. The generator holds a buffer for + /// random bits. Every time when it become empty the generator makes + /// a new random word and fill the buffer up. + bool boolean() { + return bool_producer.convert(core); + } + + /// @} + + ///\name Non-uniform Distributions + /// + ///@{ + + /// \brief Returns a random bool with given probability of true result. + /// + /// It returns a random bool with given probability of true result. + bool boolean(double p) { + return operator()() < p; + } + + /// Standard normal (Gauss) distribution + + /// Standard normal (Gauss) distribution. + /// \note The Cartesian form of the Box-Muller + /// transformation is used to generate a random normal distribution. + double gauss() + { + double V1,V2,S; + do { + V1=2*real()-1; + V2=2*real()-1; + S=V1*V1+V2*V2; + } while(S>=1); + return std::sqrt(-2*std::log(S)/S)*V1; + } + /// Normal (Gauss) distribution with given mean and standard deviation + + /// Normal (Gauss) distribution with given mean and standard deviation. + /// \sa gauss() + double gauss(double mean,double std_dev) + { + return gauss()*std_dev+mean; + } + + /// Lognormal distribution + + /// Lognormal distribution. The parameters are the mean and the standard + /// deviation of exp(X). + /// + double lognormal(double n_mean,double n_std_dev) + { + return std::exp(gauss(n_mean,n_std_dev)); + } + /// Lognormal distribution + + /// Lognormal distribution. The parameter is an std::pair of + /// the mean and the standard deviation of exp(X). + /// + double lognormal(const std::pair ¶ms) + { + return std::exp(gauss(params.first,params.second)); + } + /// Compute the lognormal parameters from mean and standard deviation + + /// This function computes the lognormal parameters from mean and + /// standard deviation. The return value can direcly be passed to + /// lognormal(). + std::pair lognormalParamsFromMD(double mean, + double std_dev) + { + double fr=std_dev/mean; + fr*=fr; + double lg=std::log(1+fr); + return std::pair(std::log(mean)-lg/2.0,std::sqrt(lg)); + } + /// Lognormal distribution with given mean and standard deviation + + /// Lognormal distribution with given mean and standard deviation. + /// + double lognormalMD(double mean,double std_dev) + { + return lognormal(lognormalParamsFromMD(mean,std_dev)); + } + + /// Exponential distribution with given mean + + /// This function generates an exponential distribution random number + /// with mean 1/lambda. + /// + double exponential(double lambda=1.0) + { + return -std::log(1.0-real())/lambda; + } + + /// Gamma distribution with given integer shape + + /// This function generates a gamma distribution random number. + /// + ///\param k shape parameter (k>0 integer) + double gamma(int k) + { + double s = 0; + for(int i=0;i()); + return s; + } + + /// Gamma distribution with given shape and scale parameter + + /// This function generates a gamma distribution random number. + /// + ///\param k shape parameter (k>0) + ///\param theta scale parameter + /// + double gamma(double k,double theta=1.0) + { + double xi,nu; + const double delta = k-std::floor(k); + const double v0=E/(E-delta); + do { + double V0=1.0-real(); + double V1=1.0-real(); + double V2=1.0-real(); + if(V2<=v0) + { + xi=std::pow(V1,1.0/delta); + nu=V0*std::pow(xi,delta-1.0); + } + else + { + xi=1.0-std::log(V1); + nu=V0*std::exp(-xi); + } + } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); + return theta*(xi+gamma(int(std::floor(k)))); + } + + /// Weibull distribution + + /// This function generates a Weibull distribution random number. + /// + ///\param k shape parameter (k>0) + ///\param lambda scale parameter (lambda>0) + /// + double weibull(double k,double lambda) + { + return lambda*pow(-std::log(1.0-real()),1.0/k); + } + + /// Pareto distribution + + /// This function generates a Pareto distribution random number. + /// + ///\param k shape parameter (k>0) + ///\param x_min location parameter (x_min>0) + /// + double pareto(double k,double x_min) + { + return exponential(gamma(k,1.0/x_min))+x_min; + } + + /// Poisson distribution + + /// This function generates a Poisson distribution random number with + /// parameter \c lambda. + /// + /// The probability mass function of this distribusion is + /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f] + /// \note The algorithm is taken from the book of Donald E. Knuth titled + /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the + /// return value. + + int poisson(double lambda) + { + const double l = std::exp(-lambda); + int k=0; + double p = 1.0; + do { + k++; + p*=real(); + } while (p>=l); + return k-1; + } + + ///@} + + ///\name Two-Dimensional Distributions + /// + ///@{ + + /// Uniform distribution on the full unit circle + + /// Uniform distribution on the full unit circle. + /// + dim2::Point disc() + { + double V1,V2; + do { + V1=2*real()-1; + V2=2*real()-1; + + } while(V1*V1+V2*V2>=1); + return dim2::Point(V1,V2); + } + /// A kind of two-dimensional normal (Gauss) distribution + + /// This function provides a turning symmetric two-dimensional distribution. + /// Both coordinates are of standard normal distribution, but they are not + /// independent. + /// + /// \note The coordinates are the two random variables provided by + /// the Box-Muller method. + dim2::Point gauss2() + { + double V1,V2,S; + do { + V1=2*real()-1; + V2=2*real()-1; + S=V1*V1+V2*V2; + } while(S>=1); + double W=std::sqrt(-2*std::log(S)/S); + return dim2::Point(W*V1,W*V2); + } + /// A kind of two-dimensional exponential distribution + + /// This function provides a turning symmetric two-dimensional distribution. + /// The x-coordinate is of conditionally exponential distribution + /// with the condition that x is positive and y=0. If x is negative and + /// y=0 then, -x is of exponential distribution. The same is true for the + /// y-coordinate. + dim2::Point exponential2() + { + double V1,V2,S; + do { + V1=2*real()-1; + V2=2*real()-1; + S=V1*V1+V2*V2; + } while(S>=1); + double W=-std::log(S)/S; + return dim2::Point(W*V1,W*V2); + } + + ///@} + }; + + + }; /// \ingroup misc /// /// \brief Mersenne Twister random number generator /// - /// The Mersenne Twister is a twisted generalized feedback - /// shift-register generator of Matsumoto and Nishimura. The period - /// of this generator is \f$ 2^{19937} - 1 \f$ and it is - /// equi-distributed in 623 dimensions for 32-bit numbers. The time - /// performance of this generator is comparable to the commonly used - /// generators. + /// This class implements either the 32-bit or the 64-bit version of + /// the Mersenne Twister random number generator algorithm + /// depending on the system architecture. + /// + /// For the API description, see its base class + /// \ref _random_bits::Random. /// - /// This implementation is specialized for both 32-bit and 64-bit - /// architectures. The generators differ sligthly in the - /// initialization and generation phase so they produce two - /// completly different sequences. + /// \sa \ref _random_bits::Random + typedef _random_bits::Random Random; + + /// \ingroup misc /// - /// The generator gives back random numbers of serveral types. To - /// get a random number from a range of a floating point type you - /// can use one form of the \c operator() or the \c real() member - /// function. If you want to get random number from the {0, 1, ..., - /// n-1} integer range use the \c operator[] or the \c integer() - /// method. And to get random number from the whole range of an - /// integer type you can use the argumentless \c integer() or \c - /// uinteger() functions. After all you can get random bool with - /// equal chance of true and false or given probability of true - /// result with the \c boolean() member functions. + /// \brief Mersenne Twister random number generator (32-bit version) /// - ///\code - /// // The commented code is identical to the other - /// double a = rnd(); // [0.0, 1.0) - /// // double a = rnd.real(); // [0.0, 1.0) - /// double b = rnd(100.0); // [0.0, 100.0) - /// // double b = rnd.real(100.0); // [0.0, 100.0) - /// double c = rnd(1.0, 2.0); // [1.0, 2.0) - /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0) - /// int d = rnd[100000]; // 0..99999 - /// // int d = rnd.integer(100000); // 0..99999 - /// int e = rnd[6] + 1; // 1..6 - /// // int e = rnd.integer(1, 1 + 6); // 1..6 - /// int b = rnd.uinteger(); // 0 .. 2^31 - 1 - /// int c = rnd.integer(); // - 2^31 .. 2^31 - 1 - /// bool g = rnd.boolean(); // P(g = true) = 0.5 - /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 - ///\endcode + /// This class implements the 32-bit version of the Mersenne Twister + /// random number generator algorithm. It is recommended to be used + /// when someone wants to make sure that the \e same pseudo random + /// sequence will be generated on every platfrom. /// - /// LEMON provides a global instance of the random number - /// generator which name is \ref lemon::rnd "rnd". Usually it is a - /// good programming convenience to use this global generator to get - /// random numbers. - class Random { - private: + /// For the API description, see its base class + /// \ref _random_bits::Random. + /// + /// \sa \ref _random_bits::Random + typedef _random_bits::Random Random32; - // Architecture word - typedef unsigned long Word; - - _random_bits::RandomCore core; - _random_bits::BoolProducer bool_producer; - - - public: - - ///\name Initialization - /// - /// @{ - - /// \brief Default constructor - /// - /// Constructor with constant seeding. - Random() { core.initState(); } - - /// \brief Constructor with seed - /// - /// Constructor with seed. The current number type will be converted - /// to the architecture word type. - template - Random(Number seed) { - _random_bits::Initializer::init(core, seed); - } - - /// \brief Constructor with array seeding - /// - /// Constructor with array seeding. The given range should contain - /// any number type and the numbers will be converted to the - /// architecture word type. - template - Random(Iterator begin, Iterator end) { - typedef typename std::iterator_traits::value_type Number; - _random_bits::Initializer::init(core, begin, end); - } - - /// \brief Copy constructor - /// - /// Copy constructor. The generated sequence will be identical to - /// the other sequence. It can be used to save the current state - /// of the generator and later use it to generate the same - /// sequence. - Random(const Random& other) { - core.copyState(other.core); - } - - /// \brief Assign operator - /// - /// Assign operator. The generated sequence will be identical to - /// the other sequence. It can be used to save the current state - /// of the generator and later use it to generate the same - /// sequence. - Random& operator=(const Random& other) { - if (&other != this) { - core.copyState(other.core); - } - return *this; - } - - /// \brief Seeding random sequence - /// - /// Seeding the random sequence. The current number type will be - /// converted to the architecture word type. - template - void seed(Number seed) { - _random_bits::Initializer::init(core, seed); - } - - /// \brief Seeding random sequence - /// - /// Seeding the random sequence. The given range should contain - /// any number type and the numbers will be converted to the - /// architecture word type. - template - void seed(Iterator begin, Iterator end) { - typedef typename std::iterator_traits::value_type Number; - _random_bits::Initializer::init(core, begin, end); - } - - /// \brief Seeding from file or from process id and time - /// - /// By default, this function calls the \c seedFromFile() member - /// function with the /dev/urandom file. If it does not success, - /// it uses the \c seedFromTime(). - /// \return Currently always \c true. - bool seed() { -#ifndef LEMON_WIN32 - if (seedFromFile("/dev/urandom", 0)) return true; -#endif - if (seedFromTime()) return true; - return false; - } - - /// \brief Seeding from file - /// - /// Seeding the random sequence from file. The linux kernel has two - /// devices, /dev/random and /dev/urandom which - /// could give good seed values for pseudo random generators (The - /// difference between two devices is that the random may - /// block the reading operation while the kernel can give good - /// source of randomness, while the urandom does not - /// block the input, but it could give back bytes with worse - /// entropy). - /// \param file The source file - /// \param offset The offset, from the file read. - /// \return \c true when the seeding successes. -#ifndef LEMON_WIN32 - bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0) -#else - bool seedFromFile(const std::string& file = "", int offset = 0) -#endif - { - std::ifstream rs(file.c_str()); - const int size = 4; - Word buf[size]; - if (offset != 0 && !rs.seekg(offset)) return false; - if (!rs.read(reinterpret_cast(buf), sizeof(buf))) return false; - seed(buf, buf + size); - return true; - } - - /// \brief Seding from process id and time - /// - /// Seding from process id and time. This function uses the - /// current process id and the current time for initialize the - /// random sequence. - /// \return Currently always \c true. - bool seedFromTime() { -#ifndef LEMON_WIN32 - timeval tv; - gettimeofday(&tv, 0); - seed(getpid() + tv.tv_sec + tv.tv_usec); -#else - seed(bits::getWinRndSeed()); -#endif - return true; - } - - /// @} - - ///\name Uniform Distributions - /// - /// @{ - - /// \brief Returns a random real number from the range [0, 1) - /// - /// It returns a random real number from the range [0, 1). The - /// default Number type is \c double. - template - Number real() { - return _random_bits::RealConversion::convert(core); - } - - double real() { - return real(); - } - - /// \brief Returns a random real number from the range [0, 1) - /// - /// It returns a random double from the range [0, 1). - double operator()() { - return real(); - } - - /// \brief Returns a random real number from the range [0, b) - /// - /// It returns a random real number from the range [0, b). - double operator()(double b) { - return real() * b; - } - - /// \brief Returns a random real number from the range [a, b) - /// - /// It returns a random real number from the range [a, b). - double operator()(double a, double b) { - return real() * (b - a) + a; - } - - /// \brief Returns a random integer from a range - /// - /// It returns a random integer from the range {0, 1, ..., b - 1}. - template - Number integer(Number b) { - return _random_bits::Mapping::map(core, b); - } - - /// \brief Returns a random integer from a range - /// - /// It returns a random integer from the range {a, a + 1, ..., b - 1}. - template - Number integer(Number a, Number b) { - return _random_bits::Mapping::map(core, b - a) + a; - } - - /// \brief Returns a random integer from a range - /// - /// It returns a random integer from the range {0, 1, ..., b - 1}. - template - Number operator[](Number b) { - return _random_bits::Mapping::map(core, b); - } - - /// \brief Returns a random non-negative integer - /// - /// It returns a random non-negative integer uniformly from the - /// whole range of the current \c Number type. The default result - /// type of this function is unsigned int. - template - Number uinteger() { - return _random_bits::IntConversion::convert(core); - } - - unsigned int uinteger() { - return uinteger(); - } - - /// \brief Returns a random integer - /// - /// It returns a random integer uniformly from the whole range of - /// the current \c Number type. The default result type of this - /// function is \c int. - template - Number integer() { - static const int nb = std::numeric_limits::digits + - (std::numeric_limits::is_signed ? 1 : 0); - return _random_bits::IntConversion::convert(core); - } - - int integer() { - return integer(); - } - - /// \brief Returns a random bool - /// - /// It returns a random bool. The generator holds a buffer for - /// random bits. Every time when it become empty the generator makes - /// a new random word and fill the buffer up. - bool boolean() { - return bool_producer.convert(core); - } - - /// @} - - ///\name Non-uniform Distributions - /// - ///@{ - - /// \brief Returns a random bool with given probability of true result. - /// - /// It returns a random bool with given probability of true result. - bool boolean(double p) { - return operator()() < p; - } - - /// Standard normal (Gauss) distribution - - /// Standard normal (Gauss) distribution. - /// \note The Cartesian form of the Box-Muller - /// transformation is used to generate a random normal distribution. - double gauss() - { - double V1,V2,S; - do { - V1=2*real()-1; - V2=2*real()-1; - S=V1*V1+V2*V2; - } while(S>=1); - return std::sqrt(-2*std::log(S)/S)*V1; - } - /// Normal (Gauss) distribution with given mean and standard deviation - - /// Normal (Gauss) distribution with given mean and standard deviation. - /// \sa gauss() - double gauss(double mean,double std_dev) - { - return gauss()*std_dev+mean; - } - - /// Lognormal distribution - - /// Lognormal distribution. The parameters are the mean and the standard - /// deviation of exp(X). - /// - double lognormal(double n_mean,double n_std_dev) - { - return std::exp(gauss(n_mean,n_std_dev)); - } - /// Lognormal distribution - - /// Lognormal distribution. The parameter is an std::pair of - /// the mean and the standard deviation of exp(X). - /// - double lognormal(const std::pair ¶ms) - { - return std::exp(gauss(params.first,params.second)); - } - /// Compute the lognormal parameters from mean and standard deviation - - /// This function computes the lognormal parameters from mean and - /// standard deviation. The return value can direcly be passed to - /// lognormal(). - std::pair lognormalParamsFromMD(double mean, - double std_dev) - { - double fr=std_dev/mean; - fr*=fr; - double lg=std::log(1+fr); - return std::pair(std::log(mean)-lg/2.0,std::sqrt(lg)); - } - /// Lognormal distribution with given mean and standard deviation - - /// Lognormal distribution with given mean and standard deviation. - /// - double lognormalMD(double mean,double std_dev) - { - return lognormal(lognormalParamsFromMD(mean,std_dev)); - } - - /// Exponential distribution with given mean - - /// This function generates an exponential distribution random number - /// with mean 1/lambda. - /// - double exponential(double lambda=1.0) - { - return -std::log(1.0-real())/lambda; - } - - /// Gamma distribution with given integer shape - - /// This function generates a gamma distribution random number. - /// - ///\param k shape parameter (k>0 integer) - double gamma(int k) - { - double s = 0; - for(int i=0;i()); - return s; - } - - /// Gamma distribution with given shape and scale parameter - - /// This function generates a gamma distribution random number. - /// - ///\param k shape parameter (k>0) - ///\param theta scale parameter - /// - double gamma(double k,double theta=1.0) - { - double xi,nu; - const double delta = k-std::floor(k); - const double v0=E/(E-delta); - do { - double V0=1.0-real(); - double V1=1.0-real(); - double V2=1.0-real(); - if(V2<=v0) - { - xi=std::pow(V1,1.0/delta); - nu=V0*std::pow(xi,delta-1.0); - } - else - { - xi=1.0-std::log(V1); - nu=V0*std::exp(-xi); - } - } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); - return theta*(xi+gamma(int(std::floor(k)))); - } - - /// Weibull distribution - - /// This function generates a Weibull distribution random number. - /// - ///\param k shape parameter (k>0) - ///\param lambda scale parameter (lambda>0) - /// - double weibull(double k,double lambda) - { - return lambda*pow(-std::log(1.0-real()),1.0/k); - } - - /// Pareto distribution - - /// This function generates a Pareto distribution random number. - /// - ///\param k shape parameter (k>0) - ///\param x_min location parameter (x_min>0) - /// - double pareto(double k,double x_min) - { - return exponential(gamma(k,1.0/x_min))+x_min; - } - - /// Poisson distribution - - /// This function generates a Poisson distribution random number with - /// parameter \c lambda. - /// - /// The probability mass function of this distribusion is - /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f] - /// \note The algorithm is taken from the book of Donald E. Knuth titled - /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the - /// return value. - - int poisson(double lambda) - { - const double l = std::exp(-lambda); - int k=0; - double p = 1.0; - do { - k++; - p*=real(); - } while (p>=l); - return k-1; - } - - ///@} - - ///\name Two Dimensional Distributions - /// - ///@{ - - /// Uniform distribution on the full unit circle - - /// Uniform distribution on the full unit circle. - /// - dim2::Point disc() - { - double V1,V2; - do { - V1=2*real()-1; - V2=2*real()-1; - - } while(V1*V1+V2*V2>=1); - return dim2::Point(V1,V2); - } - /// A kind of two dimensional normal (Gauss) distribution - - /// This function provides a turning symmetric two-dimensional distribution. - /// Both coordinates are of standard normal distribution, but they are not - /// independent. - /// - /// \note The coordinates are the two random variables provided by - /// the Box-Muller method. - dim2::Point gauss2() - { - double V1,V2,S; - do { - V1=2*real()-1; - V2=2*real()-1; - S=V1*V1+V2*V2; - } while(S>=1); - double W=std::sqrt(-2*std::log(S)/S); - return dim2::Point(W*V1,W*V2); - } - /// A kind of two dimensional exponential distribution - - /// This function provides a turning symmetric two-dimensional distribution. - /// The x-coordinate is of conditionally exponential distribution - /// with the condition that x is positive and y=0. If x is negative and - /// y=0 then, -x is of exponential distribution. The same is true for the - /// y-coordinate. - dim2::Point exponential2() - { - double V1,V2,S; - do { - V1=2*real()-1; - V2=2*real()-1; - S=V1*V1+V2*V2; - } while(S>=1); - double W=-std::log(S)/S; - return dim2::Point(W*V1,W*V2); - } - - ///@} - }; - + /// \ingroup misc + /// + /// \brief Mersenne Twister random number generator (64-bit version) + /// + /// This class implements the 64-bit version of the Mersenne Twister + /// random number generator algorithm. (Even though it runs + /// on 32-bit architectures, too.) It is recommended to be used when + /// someone wants to make sure that the \e same pseudo random sequence + /// will be generated on every platfrom. + /// + /// For the API description, see its base class + /// \ref _random_bits::Random. + /// + /// \sa \ref _random_bits::Random + typedef _random_bits::Random Random64; extern Random rnd; - + } #endif diff -r 0fdf84c79bc1 -r 57167d92e96c test/random_test.cc --- a/test/random_test.cc Fri Mar 23 15:39:54 2018 +0100 +++ b/test/random_test.cc Fri Mar 23 15:43:30 2018 +0100 @@ -21,6 +21,33 @@ int seed_array[] = {1, 2}; +int rnd_seq32[] = { +2732, 43567, 42613, 52416, 45891, 21243, 30403, 32103, +62501, 33003, 12172, 5192, 32511, 50057, 43723, 7813, +23720, 35343, 6637, 30280, 44566, 31019, 18898, 33867, +5994, 1688, 11513, 59011, 48056, 25544, 39168, 25365, +17530, 8366, 27063, 49861, 55169, 63848, 11863, 49608 +}; +int rnd_seq64[] = { +56382, 63883, 59577, 64750, 9644, 59886, 57647, 18152, +28520, 64078, 17818, 49294, 26424, 26697, 53684, 19209, +35404, 12121, 12837, 11827, 32156, 58333, 62553, 7907, +64427, 39399, 21971, 48789, 46981, 15716, 53335, 65256, +12999, 15308, 10906, 42162, 47587, 43006, 53921, 18716 +}; + +void seq_test() { + for(int i=0;i<5;i++) { + lemon::Random32 r32(i); + lemon::Random64 r64(i); + for(int j=0;j<8;j++) { + check(r32[65536]==rnd_seq32[i*8+j], "Wrong random sequence"); + check(r64[65536]==rnd_seq64[i*8+j], "Wrong random sequence"); + } + } +} + + int main() { double a=lemon::rnd(); @@ -36,5 +63,6 @@ lemon::rnd.seed(seed_array, seed_array + (sizeof(seed_array) / sizeof(seed_array[0]))); + seq_test(); return 0; }