1.1 --- a/lemon/random.h	Mon Oct 01 18:55:58 2007 +0000
     1.2 +++ b/lemon/random.h	Mon Oct 01 18:57:21 2007 +0000
     1.3 @@ -748,6 +748,44 @@
     1.4        return -std::log(real<double>())/lambda;
     1.5      }
     1.6  
     1.7 +    double gamma(int k) 
     1.8 +    {
     1.9 +      double s = 0;
    1.10 +      for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
    1.11 +      return s;
    1.12 +    }
    1.13 +    
    1.14 +    /// Gamma distribution with given shape and scale parameter
    1.15 +
    1.16 +    /// This function generates a gamma distribution random number.
    1.17 +    /// 
    1.18 +    ///\param k shape parameter (<tt>k>0</tt>)
    1.19 +    ///\param theta scale parameter
    1.20 +    ///
    1.21 +    double gamma(double k,double theta=1.0)
    1.22 +    {
    1.23 +      double xi,nu;
    1.24 +      const double delta = k-std::floor(k);
    1.25 +      const double v0=M_E/(M_E-delta);
    1.26 +      do {
    1.27 +	double V0=1.0-real<double>();
    1.28 +	double V1=1.0-real<double>();
    1.29 +	double V2=1.0-real<double>();
    1.30 +	if(V2<=v0) 
    1.31 +	  {
    1.32 +	    xi=std::pow(V1,1.0/delta);
    1.33 +	    nu=V0*std::pow(xi,delta-1.0);
    1.34 +	  }
    1.35 +	else 
    1.36 +	  {
    1.37 +	    xi=1.0-std::log(V1);
    1.38 +	    nu=V0*std::exp(-xi);
    1.39 +	  }
    1.40 +      } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
    1.41 +      return theta*(xi-gamma(int(std::floor(k))));
    1.42 +    }
    1.43 +    
    1.44 +      
    1.45      ///@}
    1.46      
    1.47      ///\name Two dimensional distributions