diff mbox

[v3] Fix negative_binomial_distribution

Message ID 4D8C69FC.30902@oracle.com
State New
Headers show

Commit Message

Paolo Carlini March 25, 2011, 10:10 a.m. UTC 
... tweaking the fix like this makes for slightly faster repeated calls 
(spares a division) and also makes clearer that we had a plain typo p / 
(1 - p) for (1 - p) / p, grrr.. Double checked Devroye in the meanwhile.

Committed to mainline, will be in 4.6.1.

Paolo.

///////////////
2011-03-25  Paolo Carlini  <paolo.carlini@oracle.com>

	* include/bits/random.h (negative_binomial_distribution<>::
	negative_binomial_distribution(_IntType, double),
	negative_binomial_distribution<>::
	negative_binomial_distribution(const param_type&)): Tweak
	construction of _M_gd.
	* include/bits/random.tcc (negative_binomial_distribution<>::
	operator()): Adjust.
diff mbox

Patch

Index: include/bits/random.tcc
===================================================================
--- include/bits/random.tcc	(revision 171411)
+++ include/bits/random.tcc	(working copy)
@@ -1075,7 +1075,7 @@ 
       return __is;
     }
 
-  // This is Leger's algorithm.
+  // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
   template<typename _IntType>
     template<typename _UniformRandomNumberGenerator>
       typename negative_binomial_distribution<_IntType>::result_type
@@ -1085,8 +1085,7 @@ 
 	const double __y = _M_gd(__urng);
 
 	// XXX Is the constructor too slow?
-	std::poisson_distribution<result_type> __poisson(__y * (1.0 - p())
-							 / p());
+	std::poisson_distribution<result_type> __poisson(__y);
 	return __poisson(__urng);
       }
 
@@ -1100,10 +1099,10 @@ 
 	typedef typename std::gamma_distribution<result_type>::param_type
 	  param_type;
 	
-	const double __y = _M_gd(__urng, param_type(__p.k(), 1.0));
+	const double __y =
+	  _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
 
-	std::poisson_distribution<result_type> __poisson(__y * (1.0 - __p.p())
-							 / __p.p() );
+	std::poisson_distribution<result_type> __poisson(__y);
 	return __poisson(__urng);
       }
 
Index: include/bits/random.h
===================================================================
--- include/bits/random.h	(revision 171411)
+++ include/bits/random.h	(working copy)
@@ -3804,12 +3804,12 @@ 
 
       explicit
       negative_binomial_distribution(_IntType __k = 1, double __p = 0.5)
-      : _M_param(__k, __p), _M_gd(__k, 1.0)
+      : _M_param(__k, __p), _M_gd(__k, (1.0 - __p) / __p)
       { }
 
       explicit
       negative_binomial_distribution(const param_type& __p)
-      : _M_param(__p), _M_gd(__p.k(), 1.0)
+      : _M_param(__p), _M_gd(__p.k(), (1.0 - __p.p()) / __p.p())
       { }
 
       /**