# HG changeset patch
# User alpar
# Date 1191265041 0
# Node ID bf6d7b624d5c92608c43a4f70099ec6e8e7fbac1
# Parent  217123f59d7e1265c0de24e37fe67b030b6b0f4f
- Gamma distributon random variable.
- Test file for random.h

diff -r 217123f59d7e -r bf6d7b624d5c lemon/random.h
--- a/lemon/random.h	Mon Oct 01 18:55:58 2007 +0000
+++ b/lemon/random.h	Mon Oct 01 18:57:21 2007 +0000
@@ -748,6 +748,44 @@
       return -std::log(real<double>())/lambda;
     }
 
+    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))));
+    }
+    
+      
     ///@}
     
     ///\name Two dimensional distributions
diff -r 217123f59d7e -r bf6d7b624d5c test/Makefile.am
--- a/test/Makefile.am	Mon Oct 01 18:55:58 2007 +0000
+++ b/test/Makefile.am	Mon Oct 01 18:57:21 2007 +0000
@@ -36,6 +36,7 @@
 	test/preflow_test \
 	test/radix_sort_test \
 	test/refptr_test \
+	test/random_test \
 	test/simann_test \
 	test/suurballe_test \
 	test/test_tools_fail \
@@ -84,6 +85,7 @@
 test_preflow_test_SOURCES = test/preflow_test.cc
 test_radix_sort_test_SOURCES = test/radix_sort_test.cc
 test_refptr_test_SOURCES = test/refptr_test.cc
+test_random_test_SOURCES = test/random_test.cc
 test_simann_test_SOURCES = test/simann_test.cc
 test_suurballe_test_SOURCES = test/suurballe_test.cc
 test_test_tools_fail_SOURCES = test/test_tools_fail.cc
diff -r 217123f59d7e -r bf6d7b624d5c test/random_test.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/random_test.cc	Mon Oct 01 18:57:21 2007 +0000
@@ -0,0 +1,36 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include <lemon/random.h>
+#include "test_tools.h"
+
+///\file \brief Test cases for random.h
+///
+///\todo To be extended
+///
+
+int main()
+{
+  double a=rnd();
+  check(a<1.0&&a>0.0,"This should be in [0,1)");
+  a=rnd.gauss();
+  a=rnd.gamma(3.45,0);
+  a=rnd.gamma(4);
+  //Does gamma work with integer k?
+  a=rnd.gamma(4.0,0);
+}