RhodeCode

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

    }

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

    {

    }  

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

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

    }

    /// A kind of two dimensional Gauss distribution