1.1 --- a/lemon/random.h Sun Jul 13 16:46:56 2008 +0100
1.2 +++ b/lemon/random.h Sun Jul 13 19:51:02 2008 +0100
1.3 @@ -1,6 +1,6 @@
1.4 -/* -*- C++ -*-
1.5 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.6 *
1.7 - * This file is a part of LEMON, a generic C++ optimization library
1.8 + * This file is a part of LEMON, a generic C++ optimization library.
1.9 *
1.10 * Copyright (C) 2003-2008
1.11 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.12 @@ -21,9 +21,9 @@
1.13 * Generator of Matsumoto and Nishimura.
1.14 *
1.15 * See the appropriate copyright notice below.
1.16 - *
1.17 + *
1.18 * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
1.19 - * All rights reserved.
1.20 + * All rights reserved.
1.21 *
1.22 * Redistribution and use in source and binary forms, with or without
1.23 * modification, are permitted provided that the following conditions
1.24 @@ -36,8 +36,8 @@
1.25 * notice, this list of conditions and the following disclaimer in the
1.26 * documentation and/or other materials provided with the distribution.
1.27 *
1.28 - * 3. The names of its contributors may not be used to endorse or promote
1.29 - * products derived from this software without specific prior written
1.30 + * 3. The names of its contributors may not be used to endorse or promote
1.31 + * products derived from this software without specific prior written
1.32 * permission.
1.33 *
1.34 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
1.35 @@ -87,7 +87,7 @@
1.36 namespace lemon {
1.37
1.38 namespace _random_bits {
1.39 -
1.40 +
1.41 template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
1.42 struct RandomTraits {};
1.43
1.44 @@ -99,7 +99,7 @@
1.45
1.46 static const int length = 624;
1.47 static const int shift = 397;
1.48 -
1.49 +
1.50 static const Word mul = 0x6c078965u;
1.51 static const Word arrayInit = 0x012BD6AAu;
1.52 static const Word arrayMul1 = 0x0019660Du;
1.53 @@ -167,7 +167,7 @@
1.54 static const Word seedArray[4] = {
1.55 0x12345u, 0x23456u, 0x34567u, 0x45678u
1.56 };
1.57 -
1.58 +
1.59 initState(seedArray, seedArray + 4);
1.60 }
1.61
1.62 @@ -175,7 +175,7 @@
1.63
1.64 static const Word mul = RandomTraits<Word>::mul;
1.65
1.66 - current = state;
1.67 + current = state;
1.68
1.69 Word *curr = state + length - 1;
1.70 curr[0] = seed; --curr;
1.71 @@ -201,7 +201,7 @@
1.72
1.73 num = length > end - begin ? length : end - begin;
1.74 while (num--) {
1.75 - curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1))
1.76 + curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1))
1.77 + *it + cnt;
1.78 ++it; ++cnt;
1.79 if (it == end) {
1.80 @@ -223,10 +223,10 @@
1.81 cnt = 1;
1.82 }
1.83 }
1.84 -
1.85 +
1.86 state[length - 1] = Word(1) << (bits - 1);
1.87 }
1.88 -
1.89 +
1.90 void copyState(const RandomCore& other) {
1.91 std::copy(other.state, other.state + length, state);
1.92 current = state + (other.current - other.state);
1.93 @@ -241,17 +241,17 @@
1.94
1.95 private:
1.96
1.97 -
1.98 +
1.99 void fillState() {
1.100 static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
1.101 static const Word loMask = RandomTraits<Word>::loMask;
1.102 static const Word hiMask = RandomTraits<Word>::hiMask;
1.103
1.104 - current = state + length;
1.105 + current = state + length;
1.106
1.107 register Word *curr = state + length - 1;
1.108 register long num;
1.109 -
1.110 +
1.111 num = length - shift;
1.112 while (num--) {
1.113 curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
1.114 @@ -269,14 +269,14 @@
1.115
1.116 }
1.117
1.118 -
1.119 +
1.120 Word *current;
1.121 Word state[length];
1.122 -
1.123 +
1.124 };
1.125
1.126
1.127 - template <typename Result,
1.128 + template <typename Result,
1.129 int shift = (std::numeric_limits<Result>::digits + 1) / 2>
1.130 struct Masker {
1.131 static Result mask(const Result& result) {
1.132 @@ -284,7 +284,7 @@
1.133 mask(static_cast<Result>(result | (result >> shift)));
1.134 }
1.135 };
1.136 -
1.137 +
1.138 template <typename Result>
1.139 struct Masker<Result, 1> {
1.140 static Result mask(const Result& result) {
1.141 @@ -292,39 +292,39 @@
1.142 }
1.143 };
1.144
1.145 - template <typename Result, typename Word,
1.146 - int rest = std::numeric_limits<Result>::digits, int shift = 0,
1.147 + template <typename Result, typename Word,
1.148 + int rest = std::numeric_limits<Result>::digits, int shift = 0,
1.149 bool last = rest <= std::numeric_limits<Word>::digits>
1.150 struct IntConversion {
1.151 static const int bits = std::numeric_limits<Word>::digits;
1.152 -
1.153 +
1.154 static Result convert(RandomCore<Word>& rnd) {
1.155 return static_cast<Result>(rnd() >> (bits - rest)) << shift;
1.156 }
1.157 -
1.158 - };
1.159
1.160 - template <typename Result, typename Word, int rest, int shift>
1.161 + };
1.162 +
1.163 + template <typename Result, typename Word, int rest, int shift>
1.164 struct IntConversion<Result, Word, rest, shift, false> {
1.165 static const int bits = std::numeric_limits<Word>::digits;
1.166
1.167 static Result convert(RandomCore<Word>& rnd) {
1.168 - return (static_cast<Result>(rnd()) << shift) |
1.169 + return (static_cast<Result>(rnd()) << shift) |
1.170 IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);
1.171 }
1.172 };
1.173
1.174
1.175 template <typename Result, typename Word,
1.176 - bool one_word = (std::numeric_limits<Word>::digits <
1.177 - std::numeric_limits<Result>::digits) >
1.178 + bool one_word = (std::numeric_limits<Word>::digits <
1.179 + std::numeric_limits<Result>::digits) >
1.180 struct Mapping {
1.181 static Result map(RandomCore<Word>& rnd, const Result& bound) {
1.182 Word max = Word(bound - 1);
1.183 Result mask = Masker<Result>::mask(bound - 1);
1.184 Result num;
1.185 do {
1.186 - num = IntConversion<Result, Word>::convert(rnd) & mask;
1.187 + num = IntConversion<Result, Word>::convert(rnd) & mask;
1.188 } while (num > max);
1.189 return num;
1.190 }
1.191 @@ -350,7 +350,7 @@
1.192 Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
1.193 res *= res;
1.194 if ((exp & 1) == 1) res *= static_cast<Result>(2.0);
1.195 - return res;
1.196 + return res;
1.197 }
1.198 };
1.199
1.200 @@ -360,42 +360,42 @@
1.201 Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
1.202 res *= res;
1.203 if ((exp & 1) == 1) res *= static_cast<Result>(0.5);
1.204 - return res;
1.205 + return res;
1.206 }
1.207 };
1.208
1.209 template <typename Result>
1.210 struct ShiftMultiplier<Result, 0, true> {
1.211 static const Result multiplier() {
1.212 - return static_cast<Result>(1.0);
1.213 + return static_cast<Result>(1.0);
1.214 }
1.215 };
1.216
1.217 template <typename Result>
1.218 struct ShiftMultiplier<Result, -20, true> {
1.219 static const Result multiplier() {
1.220 - return static_cast<Result>(1.0/1048576.0);
1.221 + return static_cast<Result>(1.0/1048576.0);
1.222 }
1.223 };
1.224 -
1.225 +
1.226 template <typename Result>
1.227 struct ShiftMultiplier<Result, -32, true> {
1.228 static const Result multiplier() {
1.229 - return static_cast<Result>(1.0/424967296.0);
1.230 + return static_cast<Result>(1.0/424967296.0);
1.231 }
1.232 };
1.233
1.234 template <typename Result>
1.235 struct ShiftMultiplier<Result, -53, true> {
1.236 static const Result multiplier() {
1.237 - return static_cast<Result>(1.0/9007199254740992.0);
1.238 + return static_cast<Result>(1.0/9007199254740992.0);
1.239 }
1.240 };
1.241
1.242 template <typename Result>
1.243 struct ShiftMultiplier<Result, -64, true> {
1.244 static const Result multiplier() {
1.245 - return static_cast<Result>(1.0/18446744073709551616.0);
1.246 + return static_cast<Result>(1.0/18446744073709551616.0);
1.247 }
1.248 };
1.249
1.250 @@ -407,9 +407,9 @@
1.251 };
1.252
1.253 template <typename Result, typename Word,
1.254 - int rest = std::numeric_limits<Result>::digits, int shift = 0,
1.255 + int rest = std::numeric_limits<Result>::digits, int shift = 0,
1.256 bool last = rest <= std::numeric_limits<Word>::digits>
1.257 - struct RealConversion{
1.258 + struct RealConversion{
1.259 static const int bits = std::numeric_limits<Word>::digits;
1.260
1.261 static Result convert(RandomCore<Word>& rnd) {
1.262 @@ -419,7 +419,7 @@
1.263 };
1.264
1.265 template <typename Result, typename Word, int rest, int shift>
1.266 - struct RealConversion<Result, Word, rest, shift, false> {
1.267 + struct RealConversion<Result, Word, rest, shift, false> {
1.268 static const int bits = std::numeric_limits<Word>::digits;
1.269
1.270 static Result convert(RandomCore<Word>& rnd) {
1.271 @@ -458,7 +458,7 @@
1.272 struct BoolProducer {
1.273 Word buffer;
1.274 int num;
1.275 -
1.276 +
1.277 BoolProducer() : num(0) {}
1.278
1.279 bool convert(RandomCore<Word>& rnd) {
1.280 @@ -529,10 +529,10 @@
1.281
1.282 // Architecture word
1.283 typedef unsigned long Word;
1.284 -
1.285 +
1.286 _random_bits::RandomCore<Word> core;
1.287 _random_bits::BoolProducer<Word> bool_producer;
1.288 -
1.289 +
1.290
1.291 public:
1.292
1.293 @@ -554,7 +554,7 @@
1.294 /// Constructor with seed. The current number type will be converted
1.295 /// to the architecture word type.
1.296 template <typename Number>
1.297 - Random(Number seed) {
1.298 + Random(Number seed) {
1.299 _random_bits::Initializer<Number, Word>::init(core, seed);
1.300 }
1.301
1.302 @@ -564,7 +564,7 @@
1.303 /// any number type and the numbers will be converted to the
1.304 /// architecture word type.
1.305 template <typename Iterator>
1.306 - Random(Iterator begin, Iterator end) {
1.307 + Random(Iterator begin, Iterator end) {
1.308 typedef typename std::iterator_traits<Iterator>::value_type Number;
1.309 _random_bits::Initializer<Number, Word>::init(core, begin, end);
1.310 }
1.311 @@ -597,7 +597,7 @@
1.312 /// Seeding the random sequence. The current number type will be
1.313 /// converted to the architecture word type.
1.314 template <typename Number>
1.315 - void seed(Number seed) {
1.316 + void seed(Number seed) {
1.317 _random_bits::Initializer<Number, Word>::init(core, seed);
1.318 }
1.319
1.320 @@ -607,7 +607,7 @@
1.321 /// any number type and the numbers will be converted to the
1.322 /// architecture word type.
1.323 template <typename Iterator>
1.324 - void seed(Iterator begin, Iterator end) {
1.325 + void seed(Iterator begin, Iterator end) {
1.326 typedef typename std::iterator_traits<Iterator>::value_type Number;
1.327 _random_bits::Initializer<Number, Word>::init(core, begin, end);
1.328 }
1.329 @@ -625,7 +625,7 @@
1.330 if (seedFromTime()) return true;
1.331 return false;
1.332 }
1.333 -
1.334 +
1.335 /// \brief Seeding from file
1.336 ///
1.337 /// Seeding the random sequence from file. The linux kernel has two
1.338 @@ -640,9 +640,9 @@
1.339 /// \param offset The offset, from the file read.
1.340 /// \return True when the seeding successes.
1.341 #ifndef WIN32
1.342 - bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
1.343 + bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
1.344 #else
1.345 - bool seedFromFile(const std::string& file = "", int offset = 0)
1.346 + bool seedFromFile(const std::string& file = "", int offset = 0)
1.347 #endif
1.348 {
1.349 std::ifstream rs(file.c_str());
1.350 @@ -660,7 +660,7 @@
1.351 /// current process id and the current time for initialize the
1.352 /// random sequence.
1.353 /// \return Currently always true.
1.354 - bool seedFromTime() {
1.355 + bool seedFromTime() {
1.356 #ifndef WIN32
1.357 timeval tv;
1.358 gettimeofday(&tv, 0);
1.359 @@ -696,16 +696,16 @@
1.360 ///
1.361 /// It returns a random real number from the range [0, b).
1.362 template <typename Number>
1.363 - Number real(Number b) {
1.364 - return real<Number>() * b;
1.365 + Number real(Number b) {
1.366 + return real<Number>() * b;
1.367 }
1.368
1.369 /// \brief Returns a random real number from the range [a, b)
1.370 ///
1.371 /// It returns a random real number from the range [a, b).
1.372 template <typename Number>
1.373 - Number real(Number a, Number b) {
1.374 - return real<Number>() * (b - a) + a;
1.375 + Number real(Number a, Number b) {
1.376 + return real<Number>() * (b - a) + a;
1.377 }
1.378
1.379 /// @}
1.380 @@ -725,16 +725,16 @@
1.381 ///
1.382 /// It returns a random real number from the range [0, b).
1.383 template <typename Number>
1.384 - Number operator()(Number b) {
1.385 - return real<Number>() * b;
1.386 + Number operator()(Number b) {
1.387 + return real<Number>() * b;
1.388 }
1.389
1.390 /// \brief Returns a random real number from the range [a, b)
1.391 ///
1.392 /// It returns a random real number from the range [a, b).
1.393 template <typename Number>
1.394 - Number operator()(Number a, Number b) {
1.395 - return real<Number>() * (b - a) + a;
1.396 + Number operator()(Number a, Number b) {
1.397 + return real<Number>() * (b - a) + a;
1.398 }
1.399
1.400 /// \brief Returns a random integer from a range
1.401 @@ -784,7 +784,7 @@
1.402 /// function is \c int.
1.403 template <typename Number>
1.404 Number integer() {
1.405 - static const int nb = std::numeric_limits<Number>::digits +
1.406 + static const int nb = std::numeric_limits<Number>::digits +
1.407 (std::numeric_limits<Number>::is_signed ? 1 : 0);
1.408 return _random_bits::IntConversion<Number, Word, nb>::convert(core);
1.409 }
1.410 @@ -792,7 +792,7 @@
1.411 int integer() {
1.412 return integer<int>();
1.413 }
1.414 -
1.415 +
1.416 /// \brief Returns a random bool
1.417 ///
1.418 /// It returns a random bool. The generator holds a buffer for
1.419 @@ -806,9 +806,9 @@
1.420
1.421 ///\name Non-uniform distributions
1.422 ///
1.423 -
1.424 +
1.425 ///@{
1.426 -
1.427 +
1.428 /// \brief Returns a random bool
1.429 ///
1.430 /// It returns a random bool with given probability of true result.
1.431 @@ -822,13 +822,13 @@
1.432 /// \note The Cartesian form of the Box-Muller
1.433 /// transformation is used to generate a random normal distribution.
1.434 /// \todo Consider using the "ziggurat" method instead.
1.435 - double gauss()
1.436 + double gauss()
1.437 {
1.438 double V1,V2,S;
1.439 do {
1.440 - V1=2*real<double>()-1;
1.441 - V2=2*real<double>()-1;
1.442 - S=V1*V1+V2*V2;
1.443 + V1=2*real<double>()-1;
1.444 + V2=2*real<double>()-1;
1.445 + S=V1*V1+V2*V2;
1.446 } while(S>=1);
1.447 return std::sqrt(-2*std::log(S)/S)*V1;
1.448 }
1.449 @@ -854,19 +854,19 @@
1.450 /// Gamma distribution with given integer shape
1.451
1.452 /// This function generates a gamma distribution random number.
1.453 - ///
1.454 + ///
1.455 ///\param k shape parameter (<tt>k>0</tt> integer)
1.456 - double gamma(int k)
1.457 + double gamma(int k)
1.458 {
1.459 double s = 0;
1.460 for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
1.461 return s;
1.462 }
1.463 -
1.464 +
1.465 /// Gamma distribution with given shape and scale parameter
1.466
1.467 /// This function generates a gamma distribution random number.
1.468 - ///
1.469 + ///
1.470 ///\param k shape parameter (<tt>k>0</tt>)
1.471 ///\param theta scale parameter
1.472 ///
1.473 @@ -876,88 +876,88 @@
1.474 const double delta = k-std::floor(k);
1.475 const double v0=E/(E-delta);
1.476 do {
1.477 - double V0=1.0-real<double>();
1.478 - double V1=1.0-real<double>();
1.479 - double V2=1.0-real<double>();
1.480 - if(V2<=v0)
1.481 - {
1.482 - xi=std::pow(V1,1.0/delta);
1.483 - nu=V0*std::pow(xi,delta-1.0);
1.484 - }
1.485 - else
1.486 - {
1.487 - xi=1.0-std::log(V1);
1.488 - nu=V0*std::exp(-xi);
1.489 - }
1.490 + double V0=1.0-real<double>();
1.491 + double V1=1.0-real<double>();
1.492 + double V2=1.0-real<double>();
1.493 + if(V2<=v0)
1.494 + {
1.495 + xi=std::pow(V1,1.0/delta);
1.496 + nu=V0*std::pow(xi,delta-1.0);
1.497 + }
1.498 + else
1.499 + {
1.500 + xi=1.0-std::log(V1);
1.501 + nu=V0*std::exp(-xi);
1.502 + }
1.503 } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
1.504 return theta*(xi+gamma(int(std::floor(k))));
1.505 }
1.506 -
1.507 +
1.508 /// Weibull distribution
1.509
1.510 /// This function generates a Weibull distribution random number.
1.511 - ///
1.512 + ///
1.513 ///\param k shape parameter (<tt>k>0</tt>)
1.514 ///\param lambda scale parameter (<tt>lambda>0</tt>)
1.515 ///
1.516 double weibull(double k,double lambda)
1.517 {
1.518 return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
1.519 - }
1.520 -
1.521 + }
1.522 +
1.523 /// Pareto distribution
1.524
1.525 /// This function generates a Pareto distribution random number.
1.526 - ///
1.527 + ///
1.528 ///\param k shape parameter (<tt>k>0</tt>)
1.529 ///\param x_min location parameter (<tt>x_min>0</tt>)
1.530 ///
1.531 double pareto(double k,double x_min)
1.532 {
1.533 return exponential(gamma(k,1.0/x_min))+x_min;
1.534 - }
1.535 -
1.536 + }
1.537 +
1.538 /// Poisson distribution
1.539
1.540 /// This function generates a Poisson distribution random number with
1.541 /// parameter \c lambda.
1.542 - ///
1.543 + ///
1.544 /// The probability mass function of this distribusion is
1.545 /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
1.546 /// \note The algorithm is taken from the book of Donald E. Knuth titled
1.547 /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
1.548 /// return value.
1.549 -
1.550 +
1.551 int poisson(double lambda)
1.552 {
1.553 const double l = std::exp(-lambda);
1.554 int k=0;
1.555 double p = 1.0;
1.556 do {
1.557 - k++;
1.558 - p*=real<double>();
1.559 + k++;
1.560 + p*=real<double>();
1.561 } while (p>=l);
1.562 return k-1;
1.563 - }
1.564 -
1.565 + }
1.566 +
1.567 ///@}
1.568 -
1.569 +
1.570 ///\name Two dimensional distributions
1.571 ///
1.572
1.573 ///@{
1.574 -
1.575 +
1.576 /// Uniform distribution on the full unit circle
1.577
1.578 /// Uniform distribution on the full unit circle.
1.579 ///
1.580 - dim2::Point<double> disc()
1.581 + dim2::Point<double> disc()
1.582 {
1.583 double V1,V2;
1.584 do {
1.585 - V1=2*real<double>()-1;
1.586 - V2=2*real<double>()-1;
1.587 -
1.588 + V1=2*real<double>()-1;
1.589 + V2=2*real<double>()-1;
1.590 +
1.591 } while(V1*V1+V2*V2>=1);
1.592 return dim2::Point<double>(V1,V2);
1.593 }
1.594 @@ -973,9 +973,9 @@
1.595 {
1.596 double V1,V2,S;
1.597 do {
1.598 - V1=2*real<double>()-1;
1.599 - V2=2*real<double>()-1;
1.600 - S=V1*V1+V2*V2;
1.601 + V1=2*real<double>()-1;
1.602 + V2=2*real<double>()-1;
1.603 + S=V1*V1+V2*V2;
1.604 } while(S>=1);
1.605 double W=std::sqrt(-2*std::log(S)/S);
1.606 return dim2::Point<double>(W*V1,W*V2);
1.607 @@ -984,22 +984,22 @@
1.608
1.609 /// This function provides a turning symmetric two-dimensional distribution.
1.610 /// The x-coordinate is of conditionally exponential distribution
1.611 - /// with the condition that x is positive and y=0. If x is negative and
1.612 + /// with the condition that x is positive and y=0. If x is negative and
1.613 /// y=0 then, -x is of exponential distribution. The same is true for the
1.614 /// y-coordinate.
1.615 - dim2::Point<double> exponential2()
1.616 + dim2::Point<double> exponential2()
1.617 {
1.618 double V1,V2,S;
1.619 do {
1.620 - V1=2*real<double>()-1;
1.621 - V2=2*real<double>()-1;
1.622 - S=V1*V1+V2*V2;
1.623 + V1=2*real<double>()-1;
1.624 + V2=2*real<double>()-1;
1.625 + S=V1*V1+V2*V2;
1.626 } while(S>=1);
1.627 double W=-std::log(S)/S;
1.628 return dim2::Point<double>(W*V1,W*V2);
1.629 }
1.630
1.631 - ///@}
1.632 + ///@}
1.633 };
1.634
1.635