Merge #602
authorAlpar Juttner <alpar@cs.elte.hu>
Fri, 23 Mar 2018 15:43:30 +0100 (2018-03-23)
changeset 117357167d92e96c
parent 1172 0fdf84c79bc1
parent 1164 04f57dad1b07
child 1185 c8d0179a32a2
Merge #602
     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> &params)
   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> &params)
   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  }