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=M_E/(M_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))));

     }