1.1 --- a/lemon/random.h Thu Oct 08 10:13:24 2015 +0200
1.2 +++ b/lemon/random.h Thu Oct 08 13:48:09 2015 +0200
1.3 @@ -296,7 +296,7 @@
1.4
1.5 template <typename Result, typename Word,
1.6 int rest = std::numeric_limits<Result>::digits, int shift = 0,
1.7 - bool last = rest <= std::numeric_limits<Word>::digits>
1.8 + bool last = (rest <= std::numeric_limits<Word>::digits)>
1.9 struct IntConversion {
1.10 static const int bits = std::numeric_limits<Word>::digits;
1.11
1.12 @@ -465,543 +465,593 @@
1.13 }
1.14 };
1.15
1.16 - }
1.17 + /// \ingroup misc
1.18 + ///
1.19 + /// \brief Mersenne Twister random number generator
1.20 + ///
1.21 + /// The Mersenne Twister is a twisted generalized feedback
1.22 + /// shift-register generator of Matsumoto and Nishimura. The period
1.23 + /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
1.24 + /// equi-distributed in 623 dimensions for 32-bit numbers. The time
1.25 + /// performance of this generator is comparable to the commonly used
1.26 + /// generators.
1.27 + ///
1.28 + /// This is a template version implementation both 32-bit and
1.29 + /// 64-bit architecture optimized versions. The generators differ
1.30 + /// sligthly in the initialization and generation phase so they
1.31 + /// produce two completly different sequences.
1.32 + ///
1.33 + /// \alert Do not use this class directly, but instead one of \ref
1.34 + /// Random, \ref Random32 or \ref Random64.
1.35 + ///
1.36 + /// The generator gives back random numbers of serveral types. To
1.37 + /// get a random number from a range of a floating point type you
1.38 + /// can use one form of the \c operator() or the \c real() member
1.39 + /// function. If you want to get random number from the {0, 1, ...,
1.40 + /// n-1} integer range use the \c operator[] or the \c integer()
1.41 + /// method. And to get random number from the whole range of an
1.42 + /// integer type you can use the argumentless \c integer() or \c
1.43 + /// uinteger() functions. After all you can get random bool with
1.44 + /// equal chance of true and false or given probability of true
1.45 + /// result with the \c boolean() member functions.
1.46 + ///
1.47 + ///\code
1.48 + /// // The commented code is identical to the other
1.49 + /// double a = rnd(); // [0.0, 1.0)
1.50 + /// // double a = rnd.real(); // [0.0, 1.0)
1.51 + /// double b = rnd(100.0); // [0.0, 100.0)
1.52 + /// // double b = rnd.real(100.0); // [0.0, 100.0)
1.53 + /// double c = rnd(1.0, 2.0); // [1.0, 2.0)
1.54 + /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0)
1.55 + /// int d = rnd[100000]; // 0..99999
1.56 + /// // int d = rnd.integer(100000); // 0..99999
1.57 + /// int e = rnd[6] + 1; // 1..6
1.58 + /// // int e = rnd.integer(1, 1 + 6); // 1..6
1.59 + /// int b = rnd.uinteger<int>(); // 0 .. 2^31 - 1
1.60 + /// int c = rnd.integer<int>(); // - 2^31 .. 2^31 - 1
1.61 + /// bool g = rnd.boolean(); // P(g = true) = 0.5
1.62 + /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8
1.63 + ///\endcode
1.64 + ///
1.65 + /// LEMON provides a global instance of the random number
1.66 + /// generator which name is \ref lemon::rnd "rnd". Usually it is a
1.67 + /// good programming convenience to use this global generator to get
1.68 + /// random numbers.
1.69 + ///
1.70 + /// \sa \ref Random, \ref Random32 or \ref Random64.
1.71 + ///
1.72 + template<class Word>
1.73 + class Random {
1.74 + private:
1.75 +
1.76 + _random_bits::RandomCore<Word> core;
1.77 + _random_bits::BoolProducer<Word> bool_producer;
1.78 +
1.79 +
1.80 + public:
1.81 +
1.82 + ///\name Initialization
1.83 + ///
1.84 + /// @{
1.85 +
1.86 + /// \brief Default constructor
1.87 + ///
1.88 + /// Constructor with constant seeding.
1.89 + Random() { core.initState(); }
1.90 +
1.91 + /// \brief Constructor with seed
1.92 + ///
1.93 + /// Constructor with seed. The current number type will be converted
1.94 + /// to the architecture word type.
1.95 + template <typename Number>
1.96 + Random(Number seed) {
1.97 + _random_bits::Initializer<Number, Word>::init(core, seed);
1.98 + }
1.99 +
1.100 + /// \brief Constructor with array seeding
1.101 + ///
1.102 + /// Constructor with array seeding. The given range should contain
1.103 + /// any number type and the numbers will be converted to the
1.104 + /// architecture word type.
1.105 + template <typename Iterator>
1.106 + Random(Iterator begin, Iterator end) {
1.107 + typedef typename std::iterator_traits<Iterator>::value_type Number;
1.108 + _random_bits::Initializer<Number, Word>::init(core, begin, end);
1.109 + }
1.110 +
1.111 + /// \brief Copy constructor
1.112 + ///
1.113 + /// Copy constructor. The generated sequence will be identical to
1.114 + /// the other sequence. It can be used to save the current state
1.115 + /// of the generator and later use it to generate the same
1.116 + /// sequence.
1.117 + Random(const Random& other) {
1.118 + core.copyState(other.core);
1.119 + }
1.120 +
1.121 + /// \brief Assign operator
1.122 + ///
1.123 + /// Assign operator. The generated sequence will be identical to
1.124 + /// the other sequence. It can be used to save the current state
1.125 + /// of the generator and later use it to generate the same
1.126 + /// sequence.
1.127 + Random& operator=(const Random& other) {
1.128 + if (&other != this) {
1.129 + core.copyState(other.core);
1.130 + }
1.131 + return *this;
1.132 + }
1.133 +
1.134 + /// \brief Seeding random sequence
1.135 + ///
1.136 + /// Seeding the random sequence. The current number type will be
1.137 + /// converted to the architecture word type.
1.138 + template <typename Number>
1.139 + void seed(Number seed) {
1.140 + _random_bits::Initializer<Number, Word>::init(core, seed);
1.141 + }
1.142 +
1.143 + /// \brief Seeding random sequence
1.144 + ///
1.145 + /// Seeding the random sequence. The given range should contain
1.146 + /// any number type and the numbers will be converted to the
1.147 + /// architecture word type.
1.148 + template <typename Iterator>
1.149 + void seed(Iterator begin, Iterator end) {
1.150 + typedef typename std::iterator_traits<Iterator>::value_type Number;
1.151 + _random_bits::Initializer<Number, Word>::init(core, begin, end);
1.152 + }
1.153 +
1.154 + /// \brief Seeding from file or from process id and time
1.155 + ///
1.156 + /// By default, this function calls the \c seedFromFile() member
1.157 + /// function with the <tt>/dev/urandom</tt> file. If it does not success,
1.158 + /// it uses the \c seedFromTime().
1.159 + /// \return Currently always \c true.
1.160 + bool seed() {
1.161 +#ifndef LEMON_WIN32
1.162 + if (seedFromFile("/dev/urandom", 0)) return true;
1.163 +#endif
1.164 + if (seedFromTime()) return true;
1.165 + return false;
1.166 + }
1.167 +
1.168 + /// \brief Seeding from file
1.169 + ///
1.170 + /// Seeding the random sequence from file. The linux kernel has two
1.171 + /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which
1.172 + /// could give good seed values for pseudo random generators (The
1.173 + /// difference between two devices is that the <tt>random</tt> may
1.174 + /// block the reading operation while the kernel can give good
1.175 + /// source of randomness, while the <tt>urandom</tt> does not
1.176 + /// block the input, but it could give back bytes with worse
1.177 + /// entropy).
1.178 + /// \param file The source file
1.179 + /// \param offset The offset, from the file read.
1.180 + /// \return \c true when the seeding successes.
1.181 +#ifndef LEMON_WIN32
1.182 + bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
1.183 +#else
1.184 + bool seedFromFile(const std::string& file = "", int offset = 0)
1.185 +#endif
1.186 + {
1.187 + std::ifstream rs(file.c_str());
1.188 + const int size = 4;
1.189 + Word buf[size];
1.190 + if (offset != 0 && !rs.seekg(offset)) return false;
1.191 + if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false;
1.192 + seed(buf, buf + size);
1.193 + return true;
1.194 + }
1.195 +
1.196 + /// \brief Seding from process id and time
1.197 + ///
1.198 + /// Seding from process id and time. This function uses the
1.199 + /// current process id and the current time for initialize the
1.200 + /// random sequence.
1.201 + /// \return Currently always \c true.
1.202 + bool seedFromTime() {
1.203 +#ifndef LEMON_WIN32
1.204 + timeval tv;
1.205 + gettimeofday(&tv, 0);
1.206 + seed(getpid() + tv.tv_sec + tv.tv_usec);
1.207 +#else
1.208 + seed(bits::getWinRndSeed());
1.209 +#endif
1.210 + return true;
1.211 + }
1.212 +
1.213 + /// @}
1.214 +
1.215 + ///\name Uniform Distributions
1.216 + ///
1.217 + /// @{
1.218 +
1.219 + /// \brief Returns a random real number from the range [0, 1)
1.220 + ///
1.221 + /// It returns a random real number from the range [0, 1). The
1.222 + /// default Number type is \c double.
1.223 + template <typename Number>
1.224 + Number real() {
1.225 + return _random_bits::RealConversion<Number, Word>::convert(core);
1.226 + }
1.227 +
1.228 + double real() {
1.229 + return real<double>();
1.230 + }
1.231 +
1.232 + /// \brief Returns a random real number from the range [0, 1)
1.233 + ///
1.234 + /// It returns a random double from the range [0, 1).
1.235 + double operator()() {
1.236 + return real<double>();
1.237 + }
1.238 +
1.239 + /// \brief Returns a random real number from the range [0, b)
1.240 + ///
1.241 + /// It returns a random real number from the range [0, b).
1.242 + double operator()(double b) {
1.243 + return real<double>() * b;
1.244 + }
1.245 +
1.246 + /// \brief Returns a random real number from the range [a, b)
1.247 + ///
1.248 + /// It returns a random real number from the range [a, b).
1.249 + double operator()(double a, double b) {
1.250 + return real<double>() * (b - a) + a;
1.251 + }
1.252 +
1.253 + /// \brief Returns a random integer from a range
1.254 + ///
1.255 + /// It returns a random integer from the range {0, 1, ..., b - 1}.
1.256 + template <typename Number>
1.257 + Number integer(Number b) {
1.258 + return _random_bits::Mapping<Number, Word>::map(core, b);
1.259 + }
1.260 +
1.261 + /// \brief Returns a random integer from a range
1.262 + ///
1.263 + /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
1.264 + template <typename Number>
1.265 + Number integer(Number a, Number b) {
1.266 + return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
1.267 + }
1.268 +
1.269 + /// \brief Returns a random integer from a range
1.270 + ///
1.271 + /// It returns a random integer from the range {0, 1, ..., b - 1}.
1.272 + template <typename Number>
1.273 + Number operator[](Number b) {
1.274 + return _random_bits::Mapping<Number, Word>::map(core, b);
1.275 + }
1.276 +
1.277 + /// \brief Returns a random non-negative integer
1.278 + ///
1.279 + /// It returns a random non-negative integer uniformly from the
1.280 + /// whole range of the current \c Number type. The default result
1.281 + /// type of this function is <tt>unsigned int</tt>.
1.282 + template <typename Number>
1.283 + Number uinteger() {
1.284 + return _random_bits::IntConversion<Number, Word>::convert(core);
1.285 + }
1.286 +
1.287 + unsigned int uinteger() {
1.288 + return uinteger<unsigned int>();
1.289 + }
1.290 +
1.291 + /// \brief Returns a random integer
1.292 + ///
1.293 + /// It returns a random integer uniformly from the whole range of
1.294 + /// the current \c Number type. The default result type of this
1.295 + /// function is \c int.
1.296 + template <typename Number>
1.297 + Number integer() {
1.298 + static const int nb = std::numeric_limits<Number>::digits +
1.299 + (std::numeric_limits<Number>::is_signed ? 1 : 0);
1.300 + return _random_bits::IntConversion<Number, Word, nb>::convert(core);
1.301 + }
1.302 +
1.303 + int integer() {
1.304 + return integer<int>();
1.305 + }
1.306 +
1.307 + /// \brief Returns a random bool
1.308 + ///
1.309 + /// It returns a random bool. The generator holds a buffer for
1.310 + /// random bits. Every time when it become empty the generator makes
1.311 + /// a new random word and fill the buffer up.
1.312 + bool boolean() {
1.313 + return bool_producer.convert(core);
1.314 + }
1.315 +
1.316 + /// @}
1.317 +
1.318 + ///\name Non-uniform Distributions
1.319 + ///
1.320 + ///@{
1.321 +
1.322 + /// \brief Returns a random bool with given probability of true result.
1.323 + ///
1.324 + /// It returns a random bool with given probability of true result.
1.325 + bool boolean(double p) {
1.326 + return operator()() < p;
1.327 + }
1.328 +
1.329 + /// Standard normal (Gauss) distribution
1.330 +
1.331 + /// Standard normal (Gauss) distribution.
1.332 + /// \note The Cartesian form of the Box-Muller
1.333 + /// transformation is used to generate a random normal distribution.
1.334 + double gauss()
1.335 + {
1.336 + double V1,V2,S;
1.337 + do {
1.338 + V1=2*real<double>()-1;
1.339 + V2=2*real<double>()-1;
1.340 + S=V1*V1+V2*V2;
1.341 + } while(S>=1);
1.342 + return std::sqrt(-2*std::log(S)/S)*V1;
1.343 + }
1.344 + /// Normal (Gauss) distribution with given mean and standard deviation
1.345 +
1.346 + /// Normal (Gauss) distribution with given mean and standard deviation.
1.347 + /// \sa gauss()
1.348 + double gauss(double mean,double std_dev)
1.349 + {
1.350 + return gauss()*std_dev+mean;
1.351 + }
1.352 +
1.353 + /// Lognormal distribution
1.354 +
1.355 + /// Lognormal distribution. The parameters are the mean and the standard
1.356 + /// deviation of <tt>exp(X)</tt>.
1.357 + ///
1.358 + double lognormal(double n_mean,double n_std_dev)
1.359 + {
1.360 + return std::exp(gauss(n_mean,n_std_dev));
1.361 + }
1.362 + /// Lognormal distribution
1.363 +
1.364 + /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of
1.365 + /// the mean and the standard deviation of <tt>exp(X)</tt>.
1.366 + ///
1.367 + double lognormal(const std::pair<double,double> ¶ms)
1.368 + {
1.369 + return std::exp(gauss(params.first,params.second));
1.370 + }
1.371 + /// Compute the lognormal parameters from mean and standard deviation
1.372 +
1.373 + /// This function computes the lognormal parameters from mean and
1.374 + /// standard deviation. The return value can direcly be passed to
1.375 + /// lognormal().
1.376 + std::pair<double,double> lognormalParamsFromMD(double mean,
1.377 + double std_dev)
1.378 + {
1.379 + double fr=std_dev/mean;
1.380 + fr*=fr;
1.381 + double lg=std::log(1+fr);
1.382 + return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg));
1.383 + }
1.384 + /// Lognormal distribution with given mean and standard deviation
1.385 +
1.386 + /// Lognormal distribution with given mean and standard deviation.
1.387 + ///
1.388 + double lognormalMD(double mean,double std_dev)
1.389 + {
1.390 + return lognormal(lognormalParamsFromMD(mean,std_dev));
1.391 + }
1.392 +
1.393 + /// Exponential distribution with given mean
1.394 +
1.395 + /// This function generates an exponential distribution random number
1.396 + /// with mean <tt>1/lambda</tt>.
1.397 + ///
1.398 + double exponential(double lambda=1.0)
1.399 + {
1.400 + return -std::log(1.0-real<double>())/lambda;
1.401 + }
1.402 +
1.403 + /// Gamma distribution with given integer shape
1.404 +
1.405 + /// This function generates a gamma distribution random number.
1.406 + ///
1.407 + ///\param k shape parameter (<tt>k>0</tt> integer)
1.408 + double gamma(int k)
1.409 + {
1.410 + double s = 0;
1.411 + for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
1.412 + return s;
1.413 + }
1.414 +
1.415 + /// Gamma distribution with given shape and scale parameter
1.416 +
1.417 + /// This function generates a gamma distribution random number.
1.418 + ///
1.419 + ///\param k shape parameter (<tt>k>0</tt>)
1.420 + ///\param theta scale parameter
1.421 + ///
1.422 + double gamma(double k,double theta=1.0)
1.423 + {
1.424 + double xi,nu;
1.425 + const double delta = k-std::floor(k);
1.426 + const double v0=E/(E-delta);
1.427 + do {
1.428 + double V0=1.0-real<double>();
1.429 + double V1=1.0-real<double>();
1.430 + double V2=1.0-real<double>();
1.431 + if(V2<=v0)
1.432 + {
1.433 + xi=std::pow(V1,1.0/delta);
1.434 + nu=V0*std::pow(xi,delta-1.0);
1.435 + }
1.436 + else
1.437 + {
1.438 + xi=1.0-std::log(V1);
1.439 + nu=V0*std::exp(-xi);
1.440 + }
1.441 + } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
1.442 + return theta*(xi+gamma(int(std::floor(k))));
1.443 + }
1.444 +
1.445 + /// Weibull distribution
1.446 +
1.447 + /// This function generates a Weibull distribution random number.
1.448 + ///
1.449 + ///\param k shape parameter (<tt>k>0</tt>)
1.450 + ///\param lambda scale parameter (<tt>lambda>0</tt>)
1.451 + ///
1.452 + double weibull(double k,double lambda)
1.453 + {
1.454 + return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
1.455 + }
1.456 +
1.457 + /// Pareto distribution
1.458 +
1.459 + /// This function generates a Pareto distribution random number.
1.460 + ///
1.461 + ///\param k shape parameter (<tt>k>0</tt>)
1.462 + ///\param x_min location parameter (<tt>x_min>0</tt>)
1.463 + ///
1.464 + double pareto(double k,double x_min)
1.465 + {
1.466 + return exponential(gamma(k,1.0/x_min))+x_min;
1.467 + }
1.468 +
1.469 + /// Poisson distribution
1.470 +
1.471 + /// This function generates a Poisson distribution random number with
1.472 + /// parameter \c lambda.
1.473 + ///
1.474 + /// The probability mass function of this distribusion is
1.475 + /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
1.476 + /// \note The algorithm is taken from the book of Donald E. Knuth titled
1.477 + /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
1.478 + /// return value.
1.479 +
1.480 + int poisson(double lambda)
1.481 + {
1.482 + const double l = std::exp(-lambda);
1.483 + int k=0;
1.484 + double p = 1.0;
1.485 + do {
1.486 + k++;
1.487 + p*=real<double>();
1.488 + } while (p>=l);
1.489 + return k-1;
1.490 + }
1.491 +
1.492 + ///@}
1.493 +
1.494 + ///\name Two Dimensional Distributions
1.495 + ///
1.496 + ///@{
1.497 +
1.498 + /// Uniform distribution on the full unit circle
1.499 +
1.500 + /// Uniform distribution on the full unit circle.
1.501 + ///
1.502 + dim2::Point<double> disc()
1.503 + {
1.504 + double V1,V2;
1.505 + do {
1.506 + V1=2*real<double>()-1;
1.507 + V2=2*real<double>()-1;
1.508 +
1.509 + } while(V1*V1+V2*V2>=1);
1.510 + return dim2::Point<double>(V1,V2);
1.511 + }
1.512 + /// A kind of two dimensional normal (Gauss) distribution
1.513 +
1.514 + /// This function provides a turning symmetric two-dimensional distribution.
1.515 + /// Both coordinates are of standard normal distribution, but they are not
1.516 + /// independent.
1.517 + ///
1.518 + /// \note The coordinates are the two random variables provided by
1.519 + /// the Box-Muller method.
1.520 + dim2::Point<double> gauss2()
1.521 + {
1.522 + double V1,V2,S;
1.523 + do {
1.524 + V1=2*real<double>()-1;
1.525 + V2=2*real<double>()-1;
1.526 + S=V1*V1+V2*V2;
1.527 + } while(S>=1);
1.528 + double W=std::sqrt(-2*std::log(S)/S);
1.529 + return dim2::Point<double>(W*V1,W*V2);
1.530 + }
1.531 + /// A kind of two dimensional exponential distribution
1.532 +
1.533 + /// This function provides a turning symmetric two-dimensional distribution.
1.534 + /// The x-coordinate is of conditionally exponential distribution
1.535 + /// with the condition that x is positive and y=0. If x is negative and
1.536 + /// y=0 then, -x is of exponential distribution. The same is true for the
1.537 + /// y-coordinate.
1.538 + dim2::Point<double> exponential2()
1.539 + {
1.540 + double V1,V2,S;
1.541 + do {
1.542 + V1=2*real<double>()-1;
1.543 + V2=2*real<double>()-1;
1.544 + S=V1*V1+V2*V2;
1.545 + } while(S>=1);
1.546 + double W=-std::log(S)/S;
1.547 + return dim2::Point<double>(W*V1,W*V2);
1.548 + }
1.549 +
1.550 + ///@}
1.551 + };
1.552 +
1.553 +
1.554 + };
1.555
1.556 /// \ingroup misc
1.557 ///
1.558 /// \brief Mersenne Twister random number generator
1.559 ///
1.560 - /// The Mersenne Twister is a twisted generalized feedback
1.561 - /// shift-register generator of Matsumoto and Nishimura. The period
1.562 - /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
1.563 - /// equi-distributed in 623 dimensions for 32-bit numbers. The time
1.564 - /// performance of this generator is comparable to the commonly used
1.565 - /// generators.
1.566 + /// This class implements either the 32 bit or the 64 bit version of
1.567 + /// the Mersenne Twister random number generator algorithm
1.568 + /// depending the the system architecture.
1.569 + ///
1.570 + /// For the API description, see its base class \ref
1.571 + /// _random_bits::Random
1.572 ///
1.573 - /// This implementation is specialized for both 32-bit and 64-bit
1.574 - /// architectures. The generators differ sligthly in the
1.575 - /// initialization and generation phase so they produce two
1.576 - /// completly different sequences.
1.577 + /// \sa \ref _random_bits::Random
1.578 + typedef _random_bits::Random<unsigned long> Random;
1.579 + /// \ingroup misc
1.580 ///
1.581 - /// The generator gives back random numbers of serveral types. To
1.582 - /// get a random number from a range of a floating point type you
1.583 - /// can use one form of the \c operator() or the \c real() member
1.584 - /// function. If you want to get random number from the {0, 1, ...,
1.585 - /// n-1} integer range use the \c operator[] or the \c integer()
1.586 - /// method. And to get random number from the whole range of an
1.587 - /// integer type you can use the argumentless \c integer() or \c
1.588 - /// uinteger() functions. After all you can get random bool with
1.589 - /// equal chance of true and false or given probability of true
1.590 - /// result with the \c boolean() member functions.
1.591 + /// \brief Mersenne Twister random number generator (32 bit version)
1.592 ///
1.593 - ///\code
1.594 - /// // The commented code is identical to the other
1.595 - /// double a = rnd(); // [0.0, 1.0)
1.596 - /// // double a = rnd.real(); // [0.0, 1.0)
1.597 - /// double b = rnd(100.0); // [0.0, 100.0)
1.598 - /// // double b = rnd.real(100.0); // [0.0, 100.0)
1.599 - /// double c = rnd(1.0, 2.0); // [1.0, 2.0)
1.600 - /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0)
1.601 - /// int d = rnd[100000]; // 0..99999
1.602 - /// // int d = rnd.integer(100000); // 0..99999
1.603 - /// int e = rnd[6] + 1; // 1..6
1.604 - /// // int e = rnd.integer(1, 1 + 6); // 1..6
1.605 - /// int b = rnd.uinteger<int>(); // 0 .. 2^31 - 1
1.606 - /// int c = rnd.integer<int>(); // - 2^31 .. 2^31 - 1
1.607 - /// bool g = rnd.boolean(); // P(g = true) = 0.5
1.608 - /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8
1.609 - ///\endcode
1.610 + /// This class implements the 32 bit version of the Mersenne Twister
1.611 + /// random number generator algorithm. It is recommended to be used
1.612 + /// when someone wants to make sure that the \e same pseudo random
1.613 + /// sequence will be generated on every platfrom.
1.614 ///
1.615 - /// LEMON provides a global instance of the random number
1.616 - /// generator which name is \ref lemon::rnd "rnd". Usually it is a
1.617 - /// good programming convenience to use this global generator to get
1.618 - /// random numbers.
1.619 - class Random {
1.620 - private:
1.621 + /// For the API description, see its base class \ref
1.622 + /// _random_bits::Random
1.623 + ///
1.624 + /// \sa \ref _random_bits::Random
1.625
1.626 - // Architecture word
1.627 - typedef unsigned long Word;
1.628 -
1.629 - _random_bits::RandomCore<Word> core;
1.630 - _random_bits::BoolProducer<Word> bool_producer;
1.631 -
1.632 -
1.633 - public:
1.634 -
1.635 - ///\name Initialization
1.636 - ///
1.637 - /// @{
1.638 -
1.639 - /// \brief Default constructor
1.640 - ///
1.641 - /// Constructor with constant seeding.
1.642 - Random() { core.initState(); }
1.643 -
1.644 - /// \brief Constructor with seed
1.645 - ///
1.646 - /// Constructor with seed. The current number type will be converted
1.647 - /// to the architecture word type.
1.648 - template <typename Number>
1.649 - Random(Number seed) {
1.650 - _random_bits::Initializer<Number, Word>::init(core, seed);
1.651 - }
1.652 -
1.653 - /// \brief Constructor with array seeding
1.654 - ///
1.655 - /// Constructor with array seeding. The given range should contain
1.656 - /// any number type and the numbers will be converted to the
1.657 - /// architecture word type.
1.658 - template <typename Iterator>
1.659 - Random(Iterator begin, Iterator end) {
1.660 - typedef typename std::iterator_traits<Iterator>::value_type Number;
1.661 - _random_bits::Initializer<Number, Word>::init(core, begin, end);
1.662 - }
1.663 -
1.664 - /// \brief Copy constructor
1.665 - ///
1.666 - /// Copy constructor. The generated sequence will be identical to
1.667 - /// the other sequence. It can be used to save the current state
1.668 - /// of the generator and later use it to generate the same
1.669 - /// sequence.
1.670 - Random(const Random& other) {
1.671 - core.copyState(other.core);
1.672 - }
1.673 -
1.674 - /// \brief Assign operator
1.675 - ///
1.676 - /// Assign operator. The generated sequence will be identical to
1.677 - /// the other sequence. It can be used to save the current state
1.678 - /// of the generator and later use it to generate the same
1.679 - /// sequence.
1.680 - Random& operator=(const Random& other) {
1.681 - if (&other != this) {
1.682 - core.copyState(other.core);
1.683 - }
1.684 - return *this;
1.685 - }
1.686 -
1.687 - /// \brief Seeding random sequence
1.688 - ///
1.689 - /// Seeding the random sequence. The current number type will be
1.690 - /// converted to the architecture word type.
1.691 - template <typename Number>
1.692 - void seed(Number seed) {
1.693 - _random_bits::Initializer<Number, Word>::init(core, seed);
1.694 - }
1.695 -
1.696 - /// \brief Seeding random sequence
1.697 - ///
1.698 - /// Seeding the random sequence. The given range should contain
1.699 - /// any number type and the numbers will be converted to the
1.700 - /// architecture word type.
1.701 - template <typename Iterator>
1.702 - void seed(Iterator begin, Iterator end) {
1.703 - typedef typename std::iterator_traits<Iterator>::value_type Number;
1.704 - _random_bits::Initializer<Number, Word>::init(core, begin, end);
1.705 - }
1.706 -
1.707 - /// \brief Seeding from file or from process id and time
1.708 - ///
1.709 - /// By default, this function calls the \c seedFromFile() member
1.710 - /// function with the <tt>/dev/urandom</tt> file. If it does not success,
1.711 - /// it uses the \c seedFromTime().
1.712 - /// \return Currently always \c true.
1.713 - bool seed() {
1.714 -#ifndef LEMON_WIN32
1.715 - if (seedFromFile("/dev/urandom", 0)) return true;
1.716 -#endif
1.717 - if (seedFromTime()) return true;
1.718 - return false;
1.719 - }
1.720 -
1.721 - /// \brief Seeding from file
1.722 - ///
1.723 - /// Seeding the random sequence from file. The linux kernel has two
1.724 - /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which
1.725 - /// could give good seed values for pseudo random generators (The
1.726 - /// difference between two devices is that the <tt>random</tt> may
1.727 - /// block the reading operation while the kernel can give good
1.728 - /// source of randomness, while the <tt>urandom</tt> does not
1.729 - /// block the input, but it could give back bytes with worse
1.730 - /// entropy).
1.731 - /// \param file The source file
1.732 - /// \param offset The offset, from the file read.
1.733 - /// \return \c true when the seeding successes.
1.734 -#ifndef LEMON_WIN32
1.735 - bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
1.736 -#else
1.737 - bool seedFromFile(const std::string& file = "", int offset = 0)
1.738 -#endif
1.739 - {
1.740 - std::ifstream rs(file.c_str());
1.741 - const int size = 4;
1.742 - Word buf[size];
1.743 - if (offset != 0 && !rs.seekg(offset)) return false;
1.744 - if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false;
1.745 - seed(buf, buf + size);
1.746 - return true;
1.747 - }
1.748 -
1.749 - /// \brief Seding from process id and time
1.750 - ///
1.751 - /// Seding from process id and time. This function uses the
1.752 - /// current process id and the current time for initialize the
1.753 - /// random sequence.
1.754 - /// \return Currently always \c true.
1.755 - bool seedFromTime() {
1.756 -#ifndef LEMON_WIN32
1.757 - timeval tv;
1.758 - gettimeofday(&tv, 0);
1.759 - seed(getpid() + tv.tv_sec + tv.tv_usec);
1.760 -#else
1.761 - seed(bits::getWinRndSeed());
1.762 -#endif
1.763 - return true;
1.764 - }
1.765 -
1.766 - /// @}
1.767 -
1.768 - ///\name Uniform Distributions
1.769 - ///
1.770 - /// @{
1.771 -
1.772 - /// \brief Returns a random real number from the range [0, 1)
1.773 - ///
1.774 - /// It returns a random real number from the range [0, 1). The
1.775 - /// default Number type is \c double.
1.776 - template <typename Number>
1.777 - Number real() {
1.778 - return _random_bits::RealConversion<Number, Word>::convert(core);
1.779 - }
1.780 -
1.781 - double real() {
1.782 - return real<double>();
1.783 - }
1.784 -
1.785 - /// \brief Returns a random real number from the range [0, 1)
1.786 - ///
1.787 - /// It returns a random double from the range [0, 1).
1.788 - double operator()() {
1.789 - return real<double>();
1.790 - }
1.791 -
1.792 - /// \brief Returns a random real number from the range [0, b)
1.793 - ///
1.794 - /// It returns a random real number from the range [0, b).
1.795 - double operator()(double b) {
1.796 - return real<double>() * b;
1.797 - }
1.798 -
1.799 - /// \brief Returns a random real number from the range [a, b)
1.800 - ///
1.801 - /// It returns a random real number from the range [a, b).
1.802 - double operator()(double a, double b) {
1.803 - return real<double>() * (b - a) + a;
1.804 - }
1.805 -
1.806 - /// \brief Returns a random integer from a range
1.807 - ///
1.808 - /// It returns a random integer from the range {0, 1, ..., b - 1}.
1.809 - template <typename Number>
1.810 - Number integer(Number b) {
1.811 - return _random_bits::Mapping<Number, Word>::map(core, b);
1.812 - }
1.813 -
1.814 - /// \brief Returns a random integer from a range
1.815 - ///
1.816 - /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
1.817 - template <typename Number>
1.818 - Number integer(Number a, Number b) {
1.819 - return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
1.820 - }
1.821 -
1.822 - /// \brief Returns a random integer from a range
1.823 - ///
1.824 - /// It returns a random integer from the range {0, 1, ..., b - 1}.
1.825 - template <typename Number>
1.826 - Number operator[](Number b) {
1.827 - return _random_bits::Mapping<Number, Word>::map(core, b);
1.828 - }
1.829 -
1.830 - /// \brief Returns a random non-negative integer
1.831 - ///
1.832 - /// It returns a random non-negative integer uniformly from the
1.833 - /// whole range of the current \c Number type. The default result
1.834 - /// type of this function is <tt>unsigned int</tt>.
1.835 - template <typename Number>
1.836 - Number uinteger() {
1.837 - return _random_bits::IntConversion<Number, Word>::convert(core);
1.838 - }
1.839 -
1.840 - unsigned int uinteger() {
1.841 - return uinteger<unsigned int>();
1.842 - }
1.843 -
1.844 - /// \brief Returns a random integer
1.845 - ///
1.846 - /// It returns a random integer uniformly from the whole range of
1.847 - /// the current \c Number type. The default result type of this
1.848 - /// function is \c int.
1.849 - template <typename Number>
1.850 - Number integer() {
1.851 - static const int nb = std::numeric_limits<Number>::digits +
1.852 - (std::numeric_limits<Number>::is_signed ? 1 : 0);
1.853 - return _random_bits::IntConversion<Number, Word, nb>::convert(core);
1.854 - }
1.855 -
1.856 - int integer() {
1.857 - return integer<int>();
1.858 - }
1.859 -
1.860 - /// \brief Returns a random bool
1.861 - ///
1.862 - /// It returns a random bool. The generator holds a buffer for
1.863 - /// random bits. Every time when it become empty the generator makes
1.864 - /// a new random word and fill the buffer up.
1.865 - bool boolean() {
1.866 - return bool_producer.convert(core);
1.867 - }
1.868 -
1.869 - /// @}
1.870 -
1.871 - ///\name Non-uniform Distributions
1.872 - ///
1.873 - ///@{
1.874 -
1.875 - /// \brief Returns a random bool with given probability of true result.
1.876 - ///
1.877 - /// It returns a random bool with given probability of true result.
1.878 - bool boolean(double p) {
1.879 - return operator()() < p;
1.880 - }
1.881 -
1.882 - /// Standard normal (Gauss) distribution
1.883 -
1.884 - /// Standard normal (Gauss) distribution.
1.885 - /// \note The Cartesian form of the Box-Muller
1.886 - /// transformation is used to generate a random normal distribution.
1.887 - double gauss()
1.888 - {
1.889 - double V1,V2,S;
1.890 - do {
1.891 - V1=2*real<double>()-1;
1.892 - V2=2*real<double>()-1;
1.893 - S=V1*V1+V2*V2;
1.894 - } while(S>=1);
1.895 - return std::sqrt(-2*std::log(S)/S)*V1;
1.896 - }
1.897 - /// Normal (Gauss) distribution with given mean and standard deviation
1.898 -
1.899 - /// Normal (Gauss) distribution with given mean and standard deviation.
1.900 - /// \sa gauss()
1.901 - double gauss(double mean,double std_dev)
1.902 - {
1.903 - return gauss()*std_dev+mean;
1.904 - }
1.905 -
1.906 - /// Lognormal distribution
1.907 -
1.908 - /// Lognormal distribution. The parameters are the mean and the standard
1.909 - /// deviation of <tt>exp(X)</tt>.
1.910 - ///
1.911 - double lognormal(double n_mean,double n_std_dev)
1.912 - {
1.913 - return std::exp(gauss(n_mean,n_std_dev));
1.914 - }
1.915 - /// Lognormal distribution
1.916 -
1.917 - /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of
1.918 - /// the mean and the standard deviation of <tt>exp(X)</tt>.
1.919 - ///
1.920 - double lognormal(const std::pair<double,double> ¶ms)
1.921 - {
1.922 - return std::exp(gauss(params.first,params.second));
1.923 - }
1.924 - /// Compute the lognormal parameters from mean and standard deviation
1.925 -
1.926 - /// This function computes the lognormal parameters from mean and
1.927 - /// standard deviation. The return value can direcly be passed to
1.928 - /// lognormal().
1.929 - std::pair<double,double> lognormalParamsFromMD(double mean,
1.930 - double std_dev)
1.931 - {
1.932 - double fr=std_dev/mean;
1.933 - fr*=fr;
1.934 - double lg=std::log(1+fr);
1.935 - return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg));
1.936 - }
1.937 - /// Lognormal distribution with given mean and standard deviation
1.938 -
1.939 - /// Lognormal distribution with given mean and standard deviation.
1.940 - ///
1.941 - double lognormalMD(double mean,double std_dev)
1.942 - {
1.943 - return lognormal(lognormalParamsFromMD(mean,std_dev));
1.944 - }
1.945 -
1.946 - /// Exponential distribution with given mean
1.947 -
1.948 - /// This function generates an exponential distribution random number
1.949 - /// with mean <tt>1/lambda</tt>.
1.950 - ///
1.951 - double exponential(double lambda=1.0)
1.952 - {
1.953 - return -std::log(1.0-real<double>())/lambda;
1.954 - }
1.955 -
1.956 - /// Gamma distribution with given integer shape
1.957 -
1.958 - /// This function generates a gamma distribution random number.
1.959 - ///
1.960 - ///\param k shape parameter (<tt>k>0</tt> integer)
1.961 - double gamma(int k)
1.962 - {
1.963 - double s = 0;
1.964 - for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
1.965 - return s;
1.966 - }
1.967 -
1.968 - /// Gamma distribution with given shape and scale parameter
1.969 -
1.970 - /// This function generates a gamma distribution random number.
1.971 - ///
1.972 - ///\param k shape parameter (<tt>k>0</tt>)
1.973 - ///\param theta scale parameter
1.974 - ///
1.975 - double gamma(double k,double theta=1.0)
1.976 - {
1.977 - double xi,nu;
1.978 - const double delta = k-std::floor(k);
1.979 - const double v0=E/(E-delta);
1.980 - do {
1.981 - double V0=1.0-real<double>();
1.982 - double V1=1.0-real<double>();
1.983 - double V2=1.0-real<double>();
1.984 - if(V2<=v0)
1.985 - {
1.986 - xi=std::pow(V1,1.0/delta);
1.987 - nu=V0*std::pow(xi,delta-1.0);
1.988 - }
1.989 - else
1.990 - {
1.991 - xi=1.0-std::log(V1);
1.992 - nu=V0*std::exp(-xi);
1.993 - }
1.994 - } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
1.995 - return theta*(xi+gamma(int(std::floor(k))));
1.996 - }
1.997 -
1.998 - /// Weibull distribution
1.999 -
1.1000 - /// This function generates a Weibull distribution random number.
1.1001 - ///
1.1002 - ///\param k shape parameter (<tt>k>0</tt>)
1.1003 - ///\param lambda scale parameter (<tt>lambda>0</tt>)
1.1004 - ///
1.1005 - double weibull(double k,double lambda)
1.1006 - {
1.1007 - return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
1.1008 - }
1.1009 -
1.1010 - /// Pareto distribution
1.1011 -
1.1012 - /// This function generates a Pareto distribution random number.
1.1013 - ///
1.1014 - ///\param k shape parameter (<tt>k>0</tt>)
1.1015 - ///\param x_min location parameter (<tt>x_min>0</tt>)
1.1016 - ///
1.1017 - double pareto(double k,double x_min)
1.1018 - {
1.1019 - return exponential(gamma(k,1.0/x_min))+x_min;
1.1020 - }
1.1021 -
1.1022 - /// Poisson distribution
1.1023 -
1.1024 - /// This function generates a Poisson distribution random number with
1.1025 - /// parameter \c lambda.
1.1026 - ///
1.1027 - /// The probability mass function of this distribusion is
1.1028 - /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
1.1029 - /// \note The algorithm is taken from the book of Donald E. Knuth titled
1.1030 - /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
1.1031 - /// return value.
1.1032 -
1.1033 - int poisson(double lambda)
1.1034 - {
1.1035 - const double l = std::exp(-lambda);
1.1036 - int k=0;
1.1037 - double p = 1.0;
1.1038 - do {
1.1039 - k++;
1.1040 - p*=real<double>();
1.1041 - } while (p>=l);
1.1042 - return k-1;
1.1043 - }
1.1044 -
1.1045 - ///@}
1.1046 -
1.1047 - ///\name Two Dimensional Distributions
1.1048 - ///
1.1049 - ///@{
1.1050 -
1.1051 - /// Uniform distribution on the full unit circle
1.1052 -
1.1053 - /// Uniform distribution on the full unit circle.
1.1054 - ///
1.1055 - dim2::Point<double> disc()
1.1056 - {
1.1057 - double V1,V2;
1.1058 - do {
1.1059 - V1=2*real<double>()-1;
1.1060 - V2=2*real<double>()-1;
1.1061 -
1.1062 - } while(V1*V1+V2*V2>=1);
1.1063 - return dim2::Point<double>(V1,V2);
1.1064 - }
1.1065 - /// A kind of two dimensional normal (Gauss) distribution
1.1066 -
1.1067 - /// This function provides a turning symmetric two-dimensional distribution.
1.1068 - /// Both coordinates are of standard normal distribution, but they are not
1.1069 - /// independent.
1.1070 - ///
1.1071 - /// \note The coordinates are the two random variables provided by
1.1072 - /// the Box-Muller method.
1.1073 - dim2::Point<double> gauss2()
1.1074 - {
1.1075 - double V1,V2,S;
1.1076 - do {
1.1077 - V1=2*real<double>()-1;
1.1078 - V2=2*real<double>()-1;
1.1079 - S=V1*V1+V2*V2;
1.1080 - } while(S>=1);
1.1081 - double W=std::sqrt(-2*std::log(S)/S);
1.1082 - return dim2::Point<double>(W*V1,W*V2);
1.1083 - }
1.1084 - /// A kind of two dimensional exponential distribution
1.1085 -
1.1086 - /// This function provides a turning symmetric two-dimensional distribution.
1.1087 - /// The x-coordinate is of conditionally exponential distribution
1.1088 - /// with the condition that x is positive and y=0. If x is negative and
1.1089 - /// y=0 then, -x is of exponential distribution. The same is true for the
1.1090 - /// y-coordinate.
1.1091 - dim2::Point<double> exponential2()
1.1092 - {
1.1093 - double V1,V2,S;
1.1094 - do {
1.1095 - V1=2*real<double>()-1;
1.1096 - V2=2*real<double>()-1;
1.1097 - S=V1*V1+V2*V2;
1.1098 - } while(S>=1);
1.1099 - double W=-std::log(S)/S;
1.1100 - return dim2::Point<double>(W*V1,W*V2);
1.1101 - }
1.1102 -
1.1103 - ///@}
1.1104 - };
1.1105 + typedef _random_bits::Random<unsigned int> Random32;
1.1106 + /// \ingroup misc
1.1107 + ///
1.1108 + /// \brief Mersenne Twister random number generator (64 bit version)
1.1109 + ///
1.1110 + /// This class implements the 64 bit version of the Mersenne Twister
1.1111 + /// random number generator algorithm. (Even though it will run
1.1112 + /// on 32 bit architectures, too.) It is recommended to ber used when
1.1113 + /// someone wants to make sure that the \e same pseudo random sequence
1.1114 + /// will be generated on every platfrom.
1.1115 + ///
1.1116 + /// For the API description, see its base class \ref
1.1117 + /// _random_bits::Random
1.1118 + ///
1.1119 + /// \sa \ref _random_bits::Random
1.1120 + typedef _random_bits::Random<unsigned long long> Random64;
1.1121
1.1122
1.1123 extern Random rnd;
1.1124
1.1125 +
1.1126 }
1.1127
1.1128 #endif
2.1 --- a/test/random_test.cc Thu Oct 08 10:13:24 2015 +0200
2.2 +++ b/test/random_test.cc Thu Oct 08 13:48:09 2015 +0200
2.3 @@ -21,6 +21,33 @@
2.4
2.5 int seed_array[] = {1, 2};
2.6
2.7 +int rnd_seq32[] = {
2.8 +2732, 43567, 42613, 52416, 45891, 21243, 30403, 32103,
2.9 +62501, 33003, 12172, 5192, 32511, 50057, 43723, 7813,
2.10 +23720, 35343, 6637, 30280, 44566, 31019, 18898, 33867,
2.11 +5994, 1688, 11513, 59011, 48056, 25544, 39168, 25365,
2.12 +17530, 8366, 27063, 49861, 55169, 63848, 11863, 49608
2.13 +};
2.14 +int rnd_seq64[] = {
2.15 +56382, 63883, 59577, 64750, 9644, 59886, 57647, 18152,
2.16 +28520, 64078, 17818, 49294, 26424, 26697, 53684, 19209,
2.17 +35404, 12121, 12837, 11827, 32156, 58333, 62553, 7907,
2.18 +64427, 39399, 21971, 48789, 46981, 15716, 53335, 65256,
2.19 +12999, 15308, 10906, 42162, 47587, 43006, 53921, 18716
2.20 +};
2.21 +
2.22 +void seq_test() {
2.23 + for(int i=0;i<5;i++) {
2.24 + lemon::Random32 r32(i);
2.25 + lemon::Random64 r64(i);
2.26 + for(int j=0;j<8;j++) {
2.27 + check(r32[65536]==rnd_seq32[i*8+j], "Wrong random sequence");
2.28 + check(r64[65536]==rnd_seq64[i*8+j], "Wrong random sequence");
2.29 + }
2.30 + }
2.31 +}
2.32 +
2.33 +
2.34 int main()
2.35 {
2.36 double a=lemon::rnd();
2.37 @@ -36,5 +63,6 @@
2.38 lemon::rnd.seed(seed_array, seed_array +
2.39 (sizeof(seed_array) / sizeof(seed_array[0])));
2.40
2.41 + seq_test();
2.42 return 0;
2.43 }