Re: [v3] Fix negative_binomial_distribution

[Paolo Carlini]

... 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.carl...@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.
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())
       { }
 
       /**