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