RhodeCode

  /// result with the \c boolean() member functions.

  ///

  ///\code

  /// // The commented code is identical to the other

  /// double a = rnd();                     // [0.0, 1.0)

  /// // double a = rnd.real();             // [0.0, 1.0)

  /// double b = rnd(100.0);                // [0.0, 100.0)

  /// // double b = rnd.real(100.0);        // [0.0, 100.0)

  /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)

  /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)

  /// int d = rnd[100000];                  // 0..99999

  /// // int d = rnd.integer(100000);       // 0..99999

  /// int e = rnd[6] + 1;                   // 1..6

  /// // int e = rnd.integer(1, 1 + 6);     // 1..6

  /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1

  /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1

  /// bool g = rnd.boolean();               // P(g = true) = 0.5

  /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8

  ///\endcode

  ///

  /// LEMON provides a global instance of the random number

  /// generator which name is \ref lemon::rnd "rnd". Usually it is a

  /// good programming convenience to use this global generator to get

  /// random numbers.

  class Random {

  private:

    // Architecture word

    typedef unsigned long Word;

    _random_bits::RandomCore<Word> core;

    _random_bits::BoolProducer<Word> bool_producer;

  public:

    ///\name Initialization

    ///

    /// @{

    ///\name Initialization

    ///

    /// @{

    /// \brief Default constructor

    ///

    /// Constructor with constant seeding.

    Random() { core.initState(); }

    /// \brief Constructor with seed

    ///

    /// Constructor with seed. The current number type will be converted

    /// to the architecture word type.

    template <typename Number>

    Random(Number seed) {

      _random_bits::Initializer<Number, Word>::init(core, seed);

    }

    /// \brief Constructor with array seeding

    ///

    /// Constructor with array seeding. The given range should contain

    /// any number type and the numbers will be converted to the

    /// architecture word type.

    template <typename Iterator>

    Random(Iterator begin, Iterator end) {

      typedef typename std::iterator_traits<Iterator>::value_type Number;

      _random_bits::Initializer<Number, Word>::init(core, begin, end);

    }

    /// \brief Copy constructor

    ///

    /// Copy constructor. The generated sequence will be identical to

    /// the other sequence. It can be used to save the current state

    /// of the generator and later use it to generate the same

    /// sequence.

    Random(const Random& other) {

      core.copyState(other.core);

    }

    /// \brief Assign operator

    ///

    /// Assign operator. The generated sequence will be identical to

    /// the other sequence. It can be used to save the current state

    /// of the generator and later use it to generate the same

    /// sequence.

    Random& operator=(const Random& other) {

      if (&other != this) {

        core.copyState(other.core);

      }

      return *this;

    }

    /// \brief Seeding random sequence

    ///

    /// Seeding the random sequence. The current number type will be

    /// converted to the architecture word type.

    template <typename Number>

    void seed(Number seed) {

      _random_bits::Initializer<Number, Word>::init(core, seed);

    }

    /// \brief Seeding random sequence

    ///

    /// Seeding the random sequence. The given range should contain

    /// any number type and the numbers will be converted to the

    /// architecture word type.

    template <typename Iterator>

    void seed(Iterator begin, Iterator end) {

      typedef typename std::iterator_traits<Iterator>::value_type Number;

      _random_bits::Initializer<Number, Word>::init(core, begin, end);

    }

    /// \brief Seeding from file or from process id and time

    ///

    /// By default, this function calls the \c seedFromFile() member

    /// function with the <tt>/dev/urandom</tt> file. If it does not success,

    /// it uses the \c seedFromTime().

    /// \return Currently always true.

    bool seed() {

#ifndef WIN32

      if (seedFromFile("/dev/urandom", 0)) return true;

#endif

      if (seedFromTime()) return true;

      return false;

    }

    /// \brief Seeding from file

    ///

    /// Seeding the random sequence from file. The linux kernel has two

    /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which

    /// could give good seed values for pseudo random generators (The

    /// difference between two devices is that the <tt>random</tt> may

    /// block the reading operation while the kernel can give good

    /// source of randomness, while the <tt>urandom</tt> does not

    /// block the input, but it could give back bytes with worse

    /// entropy).

    /// \param file The source file

    /// \param offset The offset, from the file read.

    /// \return True when the seeding successes.

#ifndef WIN32

    bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)

#else

    bool seedFromFile(const std::string& file = "", int offset = 0)

#endif

    {

      std::ifstream rs(file.c_str());

      const int size = 4;

      Word buf[size];

      if (offset != 0 && !rs.seekg(offset)) return false;

      if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false;

      seed(buf, buf + size);

      return true;

    }

    /// \brief Seding from process id and time

    ///

    /// Seding from process id and time. This function uses the

    /// current process id and the current time for initialize the

    /// random sequence.

    /// \return Currently always true.

    bool seedFromTime() {

#ifndef WIN32

      timeval tv;

      gettimeofday(&tv, 0);

      seed(getpid() + tv.tv_sec + tv.tv_usec);

#else

      FILETIME time;

      GetSystemTimeAsFileTime(&time);

      seed(GetCurrentProcessId() + time.dwHighDateTime + time.dwLowDateTime);

#endif

      return true;

    }

    /// @}

    ///\name Uniform distributions

    ///

    /// @{

    /// \brief Returns a random real number from the range [0, 1)

    ///

    /// It returns a random real number from the range [0, 1). The

    /// default Number type is \c double.

    template <typename Number>

    Number real() {

      return _random_bits::RealConversion<Number, Word>::convert(core);

    }

    double real() {

      return real<double>();

    }

    /// @}

    ///\name Uniform distributions

    ///

    /// @{

    /// \brief Returns a random real number from the range [0, 1)

    ///

    /// It returns a random double from the range [0, 1).

    double operator()() {

      return real<double>();

    }

    /// \brief Returns a random real number from the range [0, b)

    ///

    /// It returns a random real number from the range [0, b).

    }

    /// \brief Returns a random real number from the range [a, b)

    ///

    /// It returns a random real number from the range [a, b).

    }

    /// \brief Returns a random integer from a range

    ///

    /// It returns a random integer from the range {0, 1, ..., b - 1}.

    template <typename Number>

    Number integer(Number b) {

      return _random_bits::Mapping<Number, Word>::map(core, b);

    }

    /// \brief Returns a random integer from a range

    ///

    /// It returns a random integer from the range {a, a + 1, ..., b - 1}.

    template <typename Number>

    Number integer(Number a, Number b) {

      return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;

    }

    /// \brief Returns a random integer from a range

    ///

    /// It returns a random integer from the range {0, 1, ..., b - 1}.

    template <typename Number>

    Number operator[](Number b) {

      return _random_bits::Mapping<Number, Word>::map(core, b);

    }

    /// \brief Returns a random non-negative integer

    ///

    /// It returns a random non-negative integer uniformly from the

    /// whole range of the current \c Number type. The default result

    /// type of this function is <tt>unsigned int</tt>.

    template <typename Number>

    Number uinteger() {

      return _random_bits::IntConversion<Number, Word>::convert(core);

    }

    /// @}

    unsigned int uinteger() {

      return uinteger<unsigned int>();

    }

    /// \brief Returns a random integer

    ///

    /// It returns a random integer uniformly from the whole range of

    /// the current \c Number type. The default result type of this

    /// function is \c int.

    template <typename Number>

    Number integer() {

      static const int nb = std::numeric_limits<Number>::digits +

        (std::numeric_limits<Number>::is_signed ? 1 : 0);

      return _random_bits::IntConversion<Number, Word, nb>::convert(core);

    }

    int integer() {

      return integer<int>();

    }

    /// \brief Returns a random bool

    ///

    /// It returns a random bool. The generator holds a buffer for

    /// random bits. Every time when it become empty the generator makes

    /// a new random word and fill the buffer up.

    bool boolean() {

      return bool_producer.convert(core);

    }

    /// @}

    ///\name Non-uniform distributions

    ///

    ///@{

    /// \brief Returns a random bool

    ///

    /// It returns a random bool with given probability of true result.

    bool boolean(double p) {

      return operator()() < p;

    }

    /// Standard Gauss distribution

    /// Standard Gauss distribution.

    /// \note The Cartesian form of the Box-Muller

    /// transformation is used to generate a random normal distribution.

    double gauss()

    {

      double V1,V2,S;

      do {

        V1=2*real<double>()-1;

        V2=2*real<double>()-1;

        S=V1*V1+V2*V2;

      } while(S>=1);

      return std::sqrt(-2*std::log(S)/S)*V1;

    }

    /// Gauss distribution with given mean and standard deviation

    /// Gauss distribution with given mean and standard deviation.

    /// \sa gauss()

    double gauss(double mean,double std_dev)

    {

      return gauss()*std_dev+mean;

    }

    /// Exponential distribution with given mean

    /// This function generates an exponential distribution random number

    /// with mean <tt>1/lambda</tt>.

    ///

    double exponential(double lambda=1.0)

    {

      return -std::log(1.0-real<double>())/lambda;

    }

    /// Gamma distribution with given integer shape

    /// This function generates a gamma distribution random number.

    ///

    ///\param k shape parameter (<tt>k>0</tt> integer)

    double gamma(int k)

    {

      double s = 0;

      for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());

      return s;

    }

    /// Gamma distribution with given shape and scale parameter

    /// This function generates a gamma distribution random number.

    ///

    ///\param k shape parameter (<tt>k>0</tt>)

    ///\param theta scale parameter

    ///

    double gamma(double k,double theta=1.0)

    {

      double xi,nu;

      const double delta = k-std::floor(k);

      const double v0=E/(E-delta);

      do {

        double V0=1.0-real<double>();

        double V1=1.0-real<double>();

        double V2=1.0-real<double>();

        if(V2<=v0)

          {

            xi=std::pow(V1,1.0/delta);

            nu=V0*std::pow(xi,delta-1.0);

          }

        else

          {

            xi=1.0-std::log(V1);

            nu=V0*std::exp(-xi);

          }

      } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));

      return theta*(xi+gamma(int(std::floor(k))));

    }

    /// Weibull distribution

    /// This function generates a Weibull distribution random number.

    ///

    ///\param k shape parameter (<tt>k>0</tt>)

    ///\param lambda scale parameter (<tt>lambda>0</tt>)

    ///

    double weibull(double k,double lambda)

    {

      return lambda*pow(-std::log(1.0-real<double>()),1.0/k);

    }

    /// Pareto distribution

    /// This function generates a Pareto distribution random number.

    ///

    ///\param k shape parameter (<tt>k>0</tt>)

    ///\param x_min location parameter (<tt>x_min>0</tt>)

    ///

    double pareto(double k,double x_min)

    {

      return exponential(gamma(k,1.0/x_min))+x_min;

    }

    /// Poisson distribution

    /// This function generates a Poisson distribution random number with

    /// parameter \c lambda.

    ///

    /// The probability mass function of this distribusion is

    /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]

    /// \note The algorithm is taken from the book of Donald E. Knuth titled

    /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the

    /// return value.