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

     2.1 --- a/test/Makefile.am	Mon Oct 01 18:55:58 2007 +0000
     2.2 +++ b/test/Makefile.am	Mon Oct 01 18:57:21 2007 +0000
     2.3 @@ -36,6 +36,7 @@
     2.4  	test/preflow_test \
     2.5  	test/radix_sort_test \
     2.6  	test/refptr_test \
     2.7 +	test/random_test \
     2.8  	test/simann_test \
     2.9  	test/suurballe_test \
    2.10  	test/test_tools_fail \
    2.11 @@ -84,6 +85,7 @@
    2.12  test_preflow_test_SOURCES = test/preflow_test.cc
    2.13  test_radix_sort_test_SOURCES = test/radix_sort_test.cc
    2.14  test_refptr_test_SOURCES = test/refptr_test.cc
    2.15 +test_random_test_SOURCES = test/random_test.cc
    2.16  test_simann_test_SOURCES = test/simann_test.cc
    2.17  test_suurballe_test_SOURCES = test/suurballe_test.cc
    2.18  test_test_tools_fail_SOURCES = test/test_tools_fail.cc

     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/test/random_test.cc	Mon Oct 01 18:57:21 2007 +0000
     3.3 @@ -0,0 +1,36 @@
     3.4 +/* -*- C++ -*-
     3.5 + *
     3.6 + * This file is a part of LEMON, a generic C++ optimization library
     3.7 + *
     3.8 + * Copyright (C) 2003-2007
     3.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    3.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    3.11 + *
    3.12 + * Permission to use, modify and distribute this software is granted
    3.13 + * provided that this copyright notice appears in all copies. For
    3.14 + * precise terms see the accompanying LICENSE file.
    3.15 + *
    3.16 + * This software is provided "AS IS" with no warranty of any kind,
    3.17 + * express or implied, and with no claim as to its suitability for any
    3.18 + * purpose.
    3.19 + *
    3.20 + */
    3.21 +
    3.22 +#include <lemon/random.h>
    3.23 +#include "test_tools.h"
    3.24 +
    3.25 +///\file \brief Test cases for random.h
    3.26 +///
    3.27 +///\todo To be extended
    3.28 +///
    3.29 +
    3.30 +int main()
    3.31 +{
    3.32 +  double a=rnd();
    3.33 +  check(a<1.0&&a>0.0,"This should be in [0,1)");
    3.34 +  a=rnd.gauss();
    3.35 +  a=rnd.gamma(3.45,0);
    3.36 +  a=rnd.gamma(4);
    3.37 +  //Does gamma work with integer k?
    3.38 +  a=rnd.gamma(4.0,0);
    3.39 +}