--- gamma_distribution.hpp	2011-11-17 19:42:40.000000000 +0100
+++ gamma_distribution_new.hpp	2011-11-17 19:56:39.000000000 +0100
@@ -22,12 +22,15 @@
 #include <boost/limits.hpp>
 #include <boost/static_assert.hpp>
 #include <boost/random/detail/config.hpp>
-#include <boost/random/exponential_distribution.hpp>
+#include <boost/random/normal_distribution.hpp>
 
 namespace boost {
 namespace random {
 
-// The algorithm is taken from Knuth
+// The algorithm is taken from Marsaglia and Tsang, "A Simple Method
+// for generating gamma variables", ACM Transactions on Mathematical
+// Software, Vol 26, No 3 (2000), p363-372, as implemented in the
+// GNU Scientific Library, v1.15
 
 /**
  * The gamma distribution is a continuous distribution with two
@@ -112,7 +115,7 @@
      */
     explicit gamma_distribution(const result_type& alpha_arg = result_type(1.0),
                                 const result_type& beta_arg = result_type(1.0))
-      : _exp(result_type(1)), _alpha(alpha_arg), _beta(beta_arg)
+      : _normal(), _alpha(alpha_arg), _beta(beta_arg)
     {
         BOOST_ASSERT(_alpha > result_type(0));
         BOOST_ASSERT(_beta > result_type(0));
@@ -121,7 +124,7 @@
 
     /** Constructs a @c gamma_distribution from its parameters. */
     explicit gamma_distribution(const param_type& parm)
-      : _exp(result_type(1)), _alpha(parm.alpha()), _beta(parm.beta())
+      : _normal(), _alpha(parm.alpha()), _beta(parm.beta())
     {
         init();
     }
@@ -152,7 +155,7 @@
      * Effects: Subsequent uses of the distribution do not depend
      * on values produced by any engine prior to invoking reset.
      */
-    void reset() { _exp.reset(); }
+    void reset() { _normal.reset(); }
 
     /**
      * Returns a random variate distributed according to
@@ -163,45 +166,36 @@
     {
 #ifndef BOOST_NO_STDC_NAMESPACE
         // allow for Koenig lookup
-        using std::tan; using std::sqrt; using std::exp; using std::log;
-        using std::pow;
+        using std::log;  using std::pow;
 #endif
-        if(_alpha == result_type(1)) {
-            return _exp(eng) * _beta;
-        } else if(_alpha > result_type(1)) {
-            // Can we have a boost::mathconst please?
-            const result_type pi = result_type(3.14159265358979323846);
-            for(;;) {
-                result_type y = tan(pi * uniform_01<RealType>()(eng));
-                result_type x = sqrt(result_type(2)*_alpha-result_type(1))*y
-                    + _alpha-result_type(1);
-                if(x <= result_type(0))
-                    continue;
-                if(uniform_01<RealType>()(eng) >
-                    (result_type(1)+y*y) * exp((_alpha-result_type(1))
-                                               *log(x/(_alpha-result_type(1)))
-                                               - sqrt(result_type(2)*_alpha
-                                                      -result_type(1))*y))
-                    continue;
-                return x * _beta;
-            }
-        } else /* alpha < 1.0 */ {
-            for(;;) {
-                result_type u = uniform_01<RealType>()(eng);
-                result_type y = _exp(eng);
-                result_type x, q;
-                if(u < _p) {
-                    x = exp(-y/_alpha);
-                    q = _p*exp(-x);
-                } else {
-                    x = result_type(1)+y;
-                    q = _p + (result_type(1)-_p) * pow(x,_alpha-result_type(1));
-                }
-                if(u >= q)
-                    continue;
-                return x * _beta;
-            }
+        if(_alpha < result_type(1)) {
+            result_type u;
+            do u = uniform_01<RealType>()(eng);
+            while(u <= result_type(0));
+            return gamma_distribution<RealType>(_alpha + result_type(1), _beta)(eng) *
+                pow(u, result_type(1) / _alpha);
         }
+        
+        result_type v;
+        for(;;) {
+            result_type x;
+            do {
+                x = _normal(eng);
+                v = result_type(1) + _c * x;
+            } while(v <= result_type(0));
+            
+            v = v * v * v;
+            result_type u;
+            do u = uniform_01<RealType>()(eng);
+            while(u <= result_type(0));
+            
+            if (u < result_type(1) - result_type(0.0331) * x * x * x * x)
+                break;
+            
+            if (log(u) < result_type(0.5) * x * x + _d * (result_type(1) - v + log(v)))
+                break;
+        }
+        return _beta * _d * v;
     }
 
     template<class URNG>
@@ -240,7 +234,7 @@
     {
         return lhs._alpha == rhs._alpha
             && lhs._beta == rhs._beta
-            && lhs._exp == rhs._exp;
+            && lhs._normal == rhs._normal;
     }
 
     /**
@@ -269,17 +263,19 @@
     {
 #ifndef BOOST_NO_STDC_NAMESPACE
         // allow for Koenig lookup
-        using std::exp;
+        using std::sqrt;
 #endif
-        _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));
+        _d = _alpha - result_type(1) / result_type(3);
+        _c = (result_type(1) / result_type(3)) / sqrt(_d);
     }
     /// \endcond
 
-    exponential_distribution<RealType> _exp;
+    normal_distribution<RealType> _normal;
     result_type _alpha;
     result_type _beta;
     // some data precomputed from the parameters
-    result_type _p;
+    result_type _d;
+    result_type _c;
 };