Re: [v3] Fix negative_binomial_distribution

2011-03-25 Thread 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  

* 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
 template
   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 __poisson(__y * (1.0 - p())
-/ p());
+   std::poisson_distribution __poisson(__y);
return __poisson(__urng);
   }
 
@@ -1100,10 +1099,10 @@
typedef typename std::gamma_distribution::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 __poisson(__y * (1.0 - __p.p())
-/ __p.p() );
+   std::poisson_distribution __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())
   { }
 
   /**


[v3] Fix negative_binomial_distribution

2011-03-24 Thread Paolo Carlini

Hi,

this does fix a bad thinko of mine in negative_binomial_distribution 
(the fix will certainly go in 4.6.1, unless Jakub wants it now) + I'm 
adding basic statistical tests (adapted from GSL) for all the other 
discrete distributions.


Thanks,
Paolo.

//
2011-03-24  Paolo Carlini  

* include/bits/random.h (negative_binomial_distribution<>::
negative_binomial_distribution(_IntType, double),
negative_binomial_distribution<>::
negative_binomial_distribution(const param_type&)): Fix
construction of _M_gd.
* include/bits/random.tcc (negative_binomial_distribution<>::
operator()): Fix computation, per Leger's algorithm.
* testsuite/util/testsuite_random.h (discrete_pdf,
negative_binomial_pdf, poisson_pdf, uniform_int_pdf): New.
(binomial_pdf): Swap last two parameters.
* testsuite/26_numerics/random/discrete_distribution/
operators/values.cc: New.
* testsuite/26_numerics/random/negative_binomial_distribution/
operators/values.cc: Likewise.
* testsuite/26_numerics/random/poisson_distribution/
operators/values.cc: Likewise.
* testsuite/26_numerics/random/uniform_int_distribution/
operators/values.cc: Likewise.
* testsuite/26_numerics/random/binomial_distribution/
operators/values.cc: Adjust.Index: include/bits/random.tcc
===
--- include/bits/random.tcc (revision 171401)
+++ include/bits/random.tcc (working copy)
@@ -1075,7 +1075,7 @@
   return __is;
 }
 
-
+  // This is Leger's algorithm.
   template
 template
   typename negative_binomial_distribution<_IntType>::result_type
@@ -1085,7 +1085,8 @@
const double __y = _M_gd(__urng);
 
// XXX Is the constructor too slow?
-   std::poisson_distribution __poisson(__y);
+   std::poisson_distribution __poisson(__y * (1.0 - p())
+/ p());
return __poisson(__urng);
   }
 
@@ -1099,10 +1100,10 @@
typedef typename std::gamma_distribution::param_type
  param_type;

-   const double __y =
- _M_gd(__urng, param_type(__p.k(), __p.p() / (1.0 - __p.p(;
+   const double __y = _M_gd(__urng, param_type(__p.k(), 1.0));
 
-   std::poisson_distribution __poisson(__y);
+   std::poisson_distribution __poisson(__y * (1.0 - __p.p())
+/ __p.p() );
return __poisson(__urng);
   }
 
Index: include/bits/random.h
===
--- include/bits/random.h   (revision 171401)
+++ include/bits/random.h   (working copy)
@@ -3611,8 +3611,7 @@
param_type(double __p = 0.5)
: _M_p(__p)
{
- _GLIBCXX_DEBUG_ASSERT((_M_p > 0.0)
-&& (_M_p < 1.0));
+ _GLIBCXX_DEBUG_ASSERT((_M_p > 0.0) && (_M_p < 1.0));
  _M_initialize();
}
 
@@ -3782,7 +3781,9 @@
explicit
param_type(_IntType __k = 1, double __p = 0.5)
: _M_k(__k), _M_p(__p)
-   { }
+   {
+ _GLIBCXX_DEBUG_ASSERT((_M_k > 0) && (_M_p > 0.0) && (_M_p <= 1.0));
+   }
 
_IntType
k() const
@@ -3803,12 +3804,12 @@
 
   explicit
   negative_binomial_distribution(_IntType __k = 1, double __p = 0.5)
-  : _M_param(__k, __p), _M_gd(__k, __p / (1.0 - __p))
+  : _M_param(__k, __p), _M_gd(__k, 1.0)
   { }
 
   explicit
   negative_binomial_distribution(const param_type& __p)
-  : _M_param(__p), _M_gd(__p.k(), __p.p() / (1.0 - __p.p()))
+  : _M_param(__p), _M_gd(__p.k(), 1.0)
   { }
 
   /**
Index: testsuite/26_numerics/random/uniform_int_distribution/operators/values.cc
===
--- testsuite/26_numerics/random/uniform_int_distribution/operators/values.cc   
(revision 0)
+++ testsuite/26_numerics/random/uniform_int_distribution/operators/values.cc   
(revision 0)
@@ -0,0 +1,50 @@
+// { dg-options "-std=gnu++0x" }
+// { dg-require-cstdint "" }
+//
+// Copyright (C) 2011 Free Software Foundation, Inc.
+//
+// This file is part of the GNU ISO C++ Library.  This library is free
+// software; you can redistribute it and/or modify it under the
+// terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 3, or (at your option)
+// any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with this library; see the file COPYING3.  If not see
+//