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