RhodeCode

  /// 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

    ///

    /// @{

    /// \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);

    }

    Number real() {

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

    }

    double real() {

      return real<double>();

    }

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

    ///

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

    template <typename Number>

    Number real(Number b) {

      return real<Number>() * b;

    }

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

    ///

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

    template <typename Number>

    Number real(Number a, Number b) {

      return real<Number>() * (b - a) + a;

    }

    /// \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).

    template <typename Number>

    Number operator()(Number b) {

      return real<Number>() * b;

    }

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

    ///

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

    template <typename Number>

    Number operator()(Number a, Number b) {

      return real<Number>() * (b - a) + a;

    }

    /// \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

    ///

    ///@{

    ///

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

    bool boolean(double p) {

      return operator()() < p;

    }

    /// \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;

    }

    /// \sa gauss()

    double gauss(double mean,double std_dev)

    {

      return gauss()*std_dev+mean;

    }

    /// Lognormal distribution

    /// Lognormal distribution. The parameters are the mean and the standard

    /// deviation of <tt>exp(X)</tt>.

    ///

    double lognormal(double n_mean,double n_std_dev)

    {

      return std::exp(gauss(n_mean,n_std_dev));

    }

    /// Lognormal distribution

    /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of

    /// the mean and the standard deviation of <tt>exp(X)</tt>.

    ///

    double lognormal(const std::pair<double,double> &params)

    {

      return std::exp(gauss(params.first,params.second));

    }

    /// Compute the lognormal parameters from mean and standard deviation

    /// This function computes the lognormal parameters from mean and

    /// standard deviation. The return value can direcly be passed to

    /// lognormal().

    std::pair<double,double> lognormalParamsFromMD(double mean,

    {

      double fr=std_dev/mean;

      fr*=fr;

      double lg=std::log(1+fr);

      return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg));

    }

    /// Lognormal distribution with given mean and standard deviation

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

    ///

    double lognormalMD(double mean,double std_dev)

    {

      return lognormal(lognormalParamsFromMD(mean,std_dev));

    }

    /// 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

    /// 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.

    int poisson(double lambda)

    {

      const double l = std::exp(-lambda);

      int k=0;

      double p = 1.0;

      do {

        k++;

        p*=real<double>();

      } while (p>=l);

      return k-1;

    }

    ///@}

    ///\name Two dimensional distributions

    ///

    ///@{

    /// Uniform distribution on the full unit circle

    /// Uniform distribution on the full unit circle.

    ///

    dim2::Point<double> disc()

    {

      double V1,V2;

      do {

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

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

      } while(V1*V1+V2*V2>=1);

      return dim2::Point<double>(V1,V2);

    }

    /// This function provides a turning symmetric two-dimensional distribution.

    /// Both coordinates are of standard normal distribution, but they are not

    /// independent.

    ///

    /// \note The coordinates are the two random variables provided by

    /// the Box-Muller method.

    dim2::Point<double> gauss2()

    {

      double V1,V2,S;

      do {

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

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

        S=V1*V1+V2*V2;

      } while(S>=1);

      double W=std::sqrt(-2*std::log(S)/S);

      return dim2::Point<double>(W*V1,W*V2);

    }

    /// A kind of two dimensional exponential distribution

    /// This function provides a turning symmetric two-dimensional distribution.

    /// The x-coordinate is of conditionally exponential distribution

    /// with the condition that x is positive and y=0. If x is negative and

    /// y=0 then, -x is of exponential distribution. The same is true for the