random.h (negative_binomial_distribution<>:: negative_binomial_distribution(_IntType, double), [...]): Fix construction of _M_gd.
2011-03-24 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&)): 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. From-SVN: r171411
This commit is contained in:
parent
67295642aa
commit
113b21bdaf
9 changed files with 313 additions and 23 deletions
|
@ -1,3 +1,26 @@
|
|||
2011-03-24 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&)): 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.
|
||||
|
||||
2011-03-24 Rainer Orth <ro@CeBiTec.Uni-Bielefeld.DE>
|
||||
|
||||
* config/abi/post/solaris2.8/baseline_symbols.txt: Regenerate.
|
||||
|
|
|
@ -3611,8 +3611,7 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
|
|||
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 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
|
|||
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 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
|
|||
|
||||
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)
|
||||
{ }
|
||||
|
||||
/**
|
||||
|
|
|
@ -1075,7 +1075,7 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
|
|||
return __is;
|
||||
}
|
||||
|
||||
|
||||
// This is Leger's algorithm.
|
||||
template<typename _IntType>
|
||||
template<typename _UniformRandomNumberGenerator>
|
||||
typename negative_binomial_distribution<_IntType>::result_type
|
||||
|
@ -1085,7 +1085,8 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
|
|||
const double __y = _M_gd(__urng);
|
||||
|
||||
// XXX Is the constructor too slow?
|
||||
std::poisson_distribution<result_type> __poisson(__y);
|
||||
std::poisson_distribution<result_type> __poisson(__y * (1.0 - p())
|
||||
/ p());
|
||||
return __poisson(__urng);
|
||||
}
|
||||
|
||||
|
@ -1099,10 +1100,10 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
|
|||
typedef typename std::gamma_distribution<result_type>::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<result_type> __poisson(__y);
|
||||
std::poisson_distribution<result_type> __poisson(__y * (1.0 - __p.p())
|
||||
/ __p.p() );
|
||||
return __poisson(__urng);
|
||||
}
|
||||
|
||||
|
|
|
@ -33,16 +33,16 @@ void test01()
|
|||
|
||||
std::binomial_distribution<> bd1(5, 0.3);
|
||||
auto bbd1 = std::bind(bd1, eng);
|
||||
testDiscreteDist(bbd1, [](int n) { return binomial_pdf(n, 0.3, 5); } );
|
||||
testDiscreteDist(bbd1, [](int n) { return binomial_pdf(n, 5, 0.3); } );
|
||||
|
||||
std::binomial_distribution<> bd2(55, 0.3);
|
||||
auto bbd2 = std::bind(bd2, eng);
|
||||
testDiscreteDist(bbd2, [](int n) { return binomial_pdf(n, 0.3, 55); } );
|
||||
testDiscreteDist(bbd2, [](int n) { return binomial_pdf(n, 55, 0.3); } );
|
||||
|
||||
// libstdc++/48114
|
||||
std::binomial_distribution<> bd3(10, 0.75);
|
||||
auto bbd3 = std::bind(bd3, eng);
|
||||
testDiscreteDist(bbd3, [](int n) { return binomial_pdf(n, 0.75, 10); } );
|
||||
testDiscreteDist(bbd3, [](int n) { return binomial_pdf(n, 10, 0.75); } );
|
||||
}
|
||||
|
||||
int main()
|
||||
|
|
|
@ -0,0 +1,52 @@
|
|||
// { 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
|
||||
// <http://www.gnu.org/licenses/>.
|
||||
|
||||
// 26.5.8.6.1 Class template discrete_distribution [rand.dist.samp.discrete]
|
||||
|
||||
#include <random>
|
||||
#include <functional>
|
||||
#include <testsuite_random.h>
|
||||
|
||||
void test01()
|
||||
{
|
||||
using namespace __gnu_test;
|
||||
|
||||
std::mt19937 eng;
|
||||
|
||||
std::discrete_distribution<> dd1({ });
|
||||
auto bdd1 = std::bind(dd1, eng);
|
||||
testDiscreteDist(bdd1, [](int n) { return discrete_pdf(n, { }); } );
|
||||
|
||||
std::discrete_distribution<> dd2({ 1.0, 3.0, 2.0});
|
||||
auto bdd2 = std::bind(dd2, eng);
|
||||
testDiscreteDist(bdd2, [](int n)
|
||||
{ return discrete_pdf(n, { 1.0, 3.0, 2.0}); } );
|
||||
|
||||
std::discrete_distribution<> dd3({ 2.0, 2.0, 1.0, 0.0, 4.0});
|
||||
auto bdd3 = std::bind(dd3, eng);
|
||||
testDiscreteDist(bdd3, [](int n)
|
||||
{ return discrete_pdf(n, { 2.0, 2.0, 1.0, 0.0, 4.0}); } );
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
test01();
|
||||
return 0;
|
||||
}
|
|
@ -0,0 +1,55 @@
|
|||
// { dg-options "-std=gnu++0x" }
|
||||
// { dg-require-cstdint "" }
|
||||
// { dg-require-cmath "" }
|
||||
//
|
||||
// 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
|
||||
// <http://www.gnu.org/licenses/>.
|
||||
|
||||
// 26.5.8.3.4 Class template negative_binomial_distribution
|
||||
// [rand.dist.bern.negbin]
|
||||
|
||||
#include <random>
|
||||
#include <functional>
|
||||
#include <testsuite_random.h>
|
||||
|
||||
void test01()
|
||||
{
|
||||
using namespace __gnu_test;
|
||||
|
||||
std::mt19937 eng;
|
||||
|
||||
std::negative_binomial_distribution<> nbd1(5, 0.3);
|
||||
auto bnbd1 = std::bind(nbd1, eng);
|
||||
testDiscreteDist(bnbd1, [](int n)
|
||||
{ return negative_binomial_pdf(n, 5, 0.3); } );
|
||||
|
||||
std::negative_binomial_distribution<> nbd2(55, 0.3);
|
||||
auto bnbd2 = std::bind(nbd2, eng);
|
||||
testDiscreteDist(bnbd2, [](int n)
|
||||
{ return negative_binomial_pdf(n, 55, 0.3); } );
|
||||
|
||||
std::negative_binomial_distribution<> nbd3(10, 0.75);
|
||||
auto bnbd3 = std::bind(nbd3, eng);
|
||||
testDiscreteDist(bnbd3, [](int n)
|
||||
{ return negative_binomial_pdf(n, 10, 0.75); } );
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
test01();
|
||||
return 0;
|
||||
}
|
|
@ -0,0 +1,51 @@
|
|||
// { dg-options "-std=gnu++0x" }
|
||||
// { dg-require-cstdint "" }
|
||||
// { dg-require-cmath "" }
|
||||
//
|
||||
// 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
|
||||
// <http://www.gnu.org/licenses/>.
|
||||
|
||||
// 26.5.8.4.1 Class template poisson_distribution [rand.dist.pois.poisson]
|
||||
|
||||
#include <random>
|
||||
#include <functional>
|
||||
#include <testsuite_random.h>
|
||||
|
||||
void test01()
|
||||
{
|
||||
using namespace __gnu_test;
|
||||
|
||||
std::mt19937 eng;
|
||||
|
||||
std::poisson_distribution<> pd1(3.0);
|
||||
auto bpd1 = std::bind(pd1, eng);
|
||||
testDiscreteDist(bpd1, [](int n) { return poisson_pdf(n, 3.0); } );
|
||||
|
||||
std::poisson_distribution<> pd2(15.0);
|
||||
auto bpd2 = std::bind(pd2, eng);
|
||||
testDiscreteDist(bpd2, [](int n) { return poisson_pdf(n, 15.0); } );
|
||||
|
||||
std::poisson_distribution<> pd3(30.0);
|
||||
auto bpd3 = std::bind(pd3, eng);
|
||||
testDiscreteDist(bpd3, [](int n) { return poisson_pdf(n, 30.0); } );
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
test01();
|
||||
return 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
|
||||
// <http://www.gnu.org/licenses/>.
|
||||
|
||||
// 26.5.8.2.1 Class template uniform_int_distribution [rand.dist.uni.int]
|
||||
|
||||
#include <random>
|
||||
#include <functional>
|
||||
#include <testsuite_random.h>
|
||||
|
||||
void test01()
|
||||
{
|
||||
using namespace __gnu_test;
|
||||
|
||||
std::mt19937 eng;
|
||||
|
||||
std::uniform_int_distribution<> uid1(0, 2);
|
||||
auto buid1 = std::bind(uid1, eng);
|
||||
testDiscreteDist(buid1, [](int n) { return uniform_int_pdf(n, 0, 2); } );
|
||||
|
||||
std::uniform_int_distribution<> uid2(3, 7);
|
||||
auto buid2 = std::bind(uid2, eng);
|
||||
testDiscreteDist(buid2, [](int n) { return uniform_int_pdf(n, 3, 7); } );
|
||||
|
||||
std::uniform_int_distribution<> uid3(1, 20);
|
||||
auto buid3 = std::bind(uid3, eng);
|
||||
testDiscreteDist(buid3, [](int n) { return uniform_int_pdf(n, 1, 20); } );
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
test01();
|
||||
return 0;
|
||||
}
|
|
@ -25,6 +25,7 @@
|
|||
#define _GLIBCXX_TESTSUITE_RANDOM_H
|
||||
|
||||
#include <cmath>
|
||||
#include <initializer_list>
|
||||
#include <testsuite_hooks.h>
|
||||
|
||||
namespace __gnu_test
|
||||
|
@ -79,27 +80,27 @@ namespace __gnu_test
|
|||
else if (k == 1)
|
||||
return p;
|
||||
else
|
||||
return 0;
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
#ifdef _GLIBCXX_USE_C99_MATH_TR1
|
||||
inline double
|
||||
binomial_pdf(int k, double p, int n)
|
||||
binomial_pdf(int k, int n, double p)
|
||||
{
|
||||
if (k < 0 || k > n)
|
||||
return 0;
|
||||
return 0.0;
|
||||
else
|
||||
{
|
||||
double q;
|
||||
|
||||
if (p == 0)
|
||||
q = (k == 0) ? 1 : 0;
|
||||
else if (p == 1)
|
||||
q = (k == n) ? 1 : 0;
|
||||
if (p == 0.0)
|
||||
q = (k == 0) ? 1.0 : 0.0;
|
||||
else if (p == 1.0)
|
||||
q = (k == n) ? 1.0 : 0.0;
|
||||
else
|
||||
{
|
||||
double ln_Cnk = (std::lgamma(n + 1) - std::lgamma(k + 1)
|
||||
- std::lgamma(n - k + 1));
|
||||
double ln_Cnk = (std::lgamma(n + 1.0) - std::lgamma(k + 1.0)
|
||||
- std::lgamma(n - k + 1.0));
|
||||
q = ln_Cnk + k * std::log(p) + (n - k) * std::log1p(-p);
|
||||
q = std::exp(q);
|
||||
}
|
||||
|
@ -109,16 +110,72 @@ namespace __gnu_test
|
|||
}
|
||||
#endif
|
||||
|
||||
inline double
|
||||
discrete_pdf(int k, std::initializer_list<double> wl)
|
||||
{
|
||||
if (!wl.size())
|
||||
wl = { 1.0 };
|
||||
|
||||
if (k < 0 || k >= wl.size())
|
||||
return 0.0;
|
||||
else
|
||||
{
|
||||
double sum = 0.0;
|
||||
for (auto it = wl.begin(); it != wl.end(); ++it)
|
||||
sum += *it;
|
||||
return wl.begin()[k] / sum;
|
||||
}
|
||||
}
|
||||
|
||||
inline double
|
||||
geometric_pdf(int k, double p)
|
||||
{
|
||||
if (k < 0)
|
||||
return 0;
|
||||
return 0.0;
|
||||
else if (k == 0)
|
||||
return p;
|
||||
else
|
||||
return p * std::pow(1 - p, k);
|
||||
}
|
||||
|
||||
#ifdef _GLIBCXX_USE_C99_MATH_TR1
|
||||
inline double
|
||||
negative_binomial_pdf(int k, int n, double p)
|
||||
{
|
||||
if (k < 0)
|
||||
return 0.0;
|
||||
else
|
||||
{
|
||||
double f = std::lgamma(k + (double)n);
|
||||
double a = std::lgamma(n);
|
||||
double b = std::lgamma(k + 1.0);
|
||||
|
||||
return std::exp(f - a - b) * std::pow(p, n) * std::pow(1 - p, k);
|
||||
}
|
||||
}
|
||||
|
||||
inline double
|
||||
poisson_pdf(int k, double mu)
|
||||
{
|
||||
if (k < 0)
|
||||
return 0.0;
|
||||
else
|
||||
{
|
||||
double lf = std::lgamma(k + 1.0);
|
||||
return std::exp(std::log(mu) * k - lf - mu);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
inline double
|
||||
uniform_int_pdf(int k, int a, int b)
|
||||
{
|
||||
if (k < 0 || k < a || k > b)
|
||||
return 0.0;
|
||||
else
|
||||
return 1.0 / (b - a + 1.0);
|
||||
}
|
||||
|
||||
} // namespace __gnu_test
|
||||
|
||||
#endif // #ifndef _GLIBCXX_TESTSUITE_RANDOM_H
|
||||
|
|
Loading…
Add table
Reference in a new issue