random (class poisson_distribution<>): Add.
2006-08-14 Paolo Carlini <pcarlini@suse.de> * include/tr1/random (class poisson_distribution<>): Add. * include/tr1/random.tcc (poisson_distribution<>::operator(), operator<<(std::basic_ostream<>&, const poisson_distribution<>&), operator>>(std::basic_istream<>&, poisson_distribution<>&, poisson_distribution<>::poisson_distribution(const _RealType&)): Define. * testsuite/tr1/5_numerical_facilities/random/poisson_distribution/ requirements/typedefs.cc: New. * include/tr1/random.tcc (mersenne_twister<>::operator()): Tweak a bit for efficiency. * include/tr1/random.tcc (operator<<(std::basic_ostream<>&, const normal_distribution<>&), operator>>(std::basic_istream<>&, normal_distribution<>&)): Do not output _M_saved unnecessarily. * include/tr1/random: Trivial formatting fixes. * include/tr1/cmath: Likewise. From-SVN: r116149
This commit is contained in:
Paolo Carlini
Paolo Carlini
parent
e63d6886f4
commit
bbddd5d0c2
|
|
@ -1,3 +1,24 @@
|
|||
2006-08-14 Paolo Carlini <pcarlini@suse.de>
|
||||
|
||||
* include/tr1/random (class poisson_distribution<>): Add.
|
||||
* include/tr1/random.tcc (poisson_distribution<>::operator(),
|
||||
operator<<(std::basic_ostream<>&, const poisson_distribution<>&),
|
||||
operator>>(std::basic_istream<>&, poisson_distribution<>&,
|
||||
poisson_distribution<>::poisson_distribution(const _RealType&)):
|
||||
Define.
|
||||
* testsuite/tr1/5_numerical_facilities/random/poisson_distribution/
|
||||
requirements/typedefs.cc: New.
|
||||
|
||||
* include/tr1/random.tcc (mersenne_twister<>::operator()): Tweak
|
||||
a bit for efficiency.
|
||||
|
||||
* include/tr1/random.tcc (operator<<(std::basic_ostream<>&,
|
||||
const normal_distribution<>&), operator>>(std::basic_istream<>&,
|
||||
normal_distribution<>&)): Do not output _M_saved unnecessarily.
|
||||
|
||||
* include/tr1/random: Trivial formatting fixes.
|
||||
* include/tr1/cmath: Likewise.
|
||||
|
||||
2006-08-11 Paolo Carlini <pcarlini@suse.de>
|
||||
|
||||
* include/bits/stl_bvector.h (__fill_bvector(_Bit_iterator,
|
||||
|
|
|
|||
|
|
@ -375,7 +375,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
|
|||
std::__enable_if<typename std::tr1::__promote_2<_Tp, _Up>::__type,
|
||||
(std::__is_floating<_Tp>::__value
|
||||
|| std::__is_floating<_Up>::__value)>::__type
|
||||
atan2(_Tp __y, _Up __x)
|
||||
atan2(_Tp __y, _Up __x)
|
||||
{
|
||||
typedef typename std::tr1::__promote_2<_Tp, _Up>::__type __type;
|
||||
return std::atan2(__type(__y), __type(__x));
|
||||
|
|
|
|||
|
|
@ -44,6 +44,7 @@
|
|||
#include <iosfwd>
|
||||
#include <limits>
|
||||
#include <tr1/type_traits>
|
||||
#include <tr1/cmath>
|
||||
#include <fstream>
|
||||
|
||||
namespace std
|
||||
|
|
@ -1516,9 +1517,9 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
|
|||
/**
|
||||
* Gets the next value in the Bernoullian sequence.
|
||||
*/
|
||||
template<class UniformRandomNumberGenerator>
|
||||
template<class _UniformRandomNumberGenerator>
|
||||
result_type
|
||||
operator()(UniformRandomNumberGenerator& __urng)
|
||||
operator()(_UniformRandomNumberGenerator& __urng)
|
||||
{
|
||||
if (__urng() < _M_p)
|
||||
return true;
|
||||
|
|
@ -1585,7 +1586,6 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
|
|||
typedef _IntType result_type;
|
||||
|
||||
// constructors and member function
|
||||
|
||||
explicit
|
||||
geometric_distribution(const _RealType& __p = _RealType(0.5))
|
||||
: _M_p(__p), _M_log_p(std::log(_M_p))
|
||||
|
|
@ -1648,6 +1648,96 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
|
|||
_RealType _M_log_p;
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* @brief A discrete Poisson random number distribution.
|
||||
*
|
||||
* The formula for the poisson probability mass function is
|
||||
* @f$ p(i) = \frac{mean^i}{i!} e^{-mean} @f$ where @f$ mean @f$ is the
|
||||
* parameter of the distribution.
|
||||
*/
|
||||
template<typename _IntType = int, typename _RealType = double>
|
||||
class poisson_distribution;
|
||||
|
||||
template<typename _IntType, typename _RealType,
|
||||
typename _CharT, typename _Traits>
|
||||
std::basic_ostream<_CharT, _Traits>&
|
||||
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
||||
const poisson_distribution<_IntType, _RealType>& __x);
|
||||
|
||||
template<typename _IntType, typename _RealType,
|
||||
typename _CharT, typename _Traits>
|
||||
std::basic_istream<_CharT, _Traits>&
|
||||
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
||||
poisson_distribution<_IntType, _RealType>& __x);
|
||||
|
||||
template<typename _IntType, typename _RealType>
|
||||
class poisson_distribution
|
||||
{
|
||||
public:
|
||||
// types
|
||||
typedef _RealType input_type;
|
||||
typedef _IntType result_type;
|
||||
|
||||
// constructors and member function
|
||||
explicit
|
||||
poisson_distribution(const _RealType& __mean = _RealType(1));
|
||||
|
||||
/**
|
||||
* Gets the distribution parameter @p mean.
|
||||
*/
|
||||
_RealType
|
||||
mean() const
|
||||
{ return _M_mean; }
|
||||
|
||||
void
|
||||
reset() { }
|
||||
|
||||
template<class _UniformRandomNumberGenerator>
|
||||
result_type
|
||||
operator()(_UniformRandomNumberGenerator& __urng);
|
||||
|
||||
/**
|
||||
* Inserts a %poisson_distribution random number distribution
|
||||
* @p __x into the output stream @p __os.
|
||||
*
|
||||
* @param __os An output stream.
|
||||
* @param __x A %poisson_distribution random number distribution.
|
||||
*
|
||||
* @returns The output stream with the state of @p __x inserted or in
|
||||
* an error state.
|
||||
*/
|
||||
template<typename _IntType1, typename _RealType1,
|
||||
typename _CharT, typename _Traits>
|
||||
friend std::basic_ostream<_CharT, _Traits>&
|
||||
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
||||
const poisson_distribution<_IntType1, _RealType1>& __x);
|
||||
|
||||
/**
|
||||
* Extracts a %poisson_distribution random number distribution
|
||||
* @p __x from the input stream @p __is.
|
||||
*
|
||||
* @param __is An input stream.
|
||||
* @param __x A %poisson_distribution random number generator engine.
|
||||
*
|
||||
* @returns The input stream with @p __x extracted or in an error state.
|
||||
*/
|
||||
template<typename _IntType1, typename _RealType1,
|
||||
typename _CharT, typename _Traits>
|
||||
friend std::basic_istream<_CharT, _Traits>&
|
||||
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
||||
poisson_distribution<_IntType1, _RealType1>& __x);
|
||||
|
||||
protected:
|
||||
_RealType _M_mean;
|
||||
|
||||
_RealType _M_lm_thr;
|
||||
#if _GLIBCXX_USE_C99_MATH_TR1
|
||||
_RealType _M_lfm, _M_sm, _M_d, _M_scx4, _M_2cx, _M_c2b, _M_cb;
|
||||
#endif
|
||||
bool _M_large;
|
||||
};
|
||||
|
||||
/* @} */ // group tr1_random_distributions_discrete
|
||||
|
||||
/**
|
||||
|
|
|
|||
|
|
@ -285,13 +285,13 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
|
|||
{
|
||||
const _UIntType __upper_mask = (~_UIntType()) << __r;
|
||||
const _UIntType __lower_mask = ~__upper_mask;
|
||||
const _UIntType __fx[2] = { 0, __a };
|
||||
|
||||
for (int __k = 0; __k < (__n - __m); ++__k)
|
||||
{
|
||||
_UIntType __y = ((_M_x[__k] & __upper_mask)
|
||||
| (_M_x[__k + 1] & __lower_mask));
|
||||
_M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
|
||||
^ ((__y & 0x01) ? __a : 0));
|
||||
_M_x[__k] = _M_x[__k + __m] ^ (__y >> 1) ^ __fx[__y & 0x01];
|
||||
}
|
||||
|
||||
for (int __k = (__n - __m); __k < (__n - 1); ++__k)
|
||||
|
|
@ -299,13 +299,12 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
|
|||
_UIntType __y = ((_M_x[__k] & __upper_mask)
|
||||
| (_M_x[__k + 1] & __lower_mask));
|
||||
_M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
|
||||
^ ((__y & 0x01) ? __a : 0));
|
||||
^ __fx[__y & 0x01]);
|
||||
}
|
||||
|
||||
_UIntType __y = ((_M_x[__n - 1] & __upper_mask)
|
||||
| (_M_x[0] & __lower_mask));
|
||||
_M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
|
||||
^ ((__y & 0x01) ? __a : 0));
|
||||
_M_x[__n - 1] = _M_x[__m - 1] ^ (__y >> 1) ^ __fx[__y & 0x01];
|
||||
_M_p = 0;
|
||||
}
|
||||
|
||||
|
|
@ -655,6 +654,194 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
|
|||
}
|
||||
|
||||
|
||||
template<typename _IntType, typename _RealType>
|
||||
poisson_distribution<_IntType, _RealType>::
|
||||
poisson_distribution(const _RealType& __mean)
|
||||
: _M_mean(__mean), _M_large(false)
|
||||
{
|
||||
_GLIBCXX_DEBUG_ASSERT(_M_mean > 0.0);
|
||||
|
||||
#if _GLIBCXX_USE_C99_MATH_TR1
|
||||
if (_M_mean >= 12)
|
||||
{
|
||||
_M_large = true;
|
||||
const _RealType __m = std::floor(_M_mean);
|
||||
_M_lm_thr = std::log(_M_mean);
|
||||
_M_lfm = std::tr1::lgamma(__m + 1);
|
||||
_M_sm = std::sqrt(__m);
|
||||
|
||||
const _RealType __dx =
|
||||
std::sqrt(2 * __m
|
||||
* std::log(_RealType(40.743665431525205956834243423363677L)
|
||||
* __m));
|
||||
_M_d = std::tr1::round(std::max(_RealType(6),
|
||||
std::min(__m, __dx)));
|
||||
const _RealType __cx = 2 * (2 * __m + _M_d);
|
||||
const _RealType __cx4 = __cx / 4;
|
||||
_M_scx4 = std::sqrt(__cx4);
|
||||
_M_2cx = 2 / __cx;
|
||||
|
||||
const _RealType __pi_2 = 1.5707963267948966192313216916397514L;
|
||||
_M_c2b = std::sqrt(__pi_2 * __cx4) * std::exp(_M_2cx);
|
||||
_M_cb = __cx * std::exp(-_M_d * _M_2cx * (1 + _M_d / 2)) / _M_d;
|
||||
}
|
||||
else
|
||||
#endif
|
||||
_M_lm_thr = std::exp(-_M_mean);
|
||||
}
|
||||
|
||||
/**
|
||||
* A rejection algorithm when mean >= 12 and a simple method based
|
||||
* upon the multiplication of uniform random variates otherwise.
|
||||
* NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
|
||||
* is defined.
|
||||
*
|
||||
* Reference:
|
||||
* Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
|
||||
* New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
|
||||
*/
|
||||
template<typename _IntType, typename _RealType>
|
||||
template<class _UniformRandomNumberGenerator>
|
||||
typename poisson_distribution<_IntType, _RealType>::result_type
|
||||
poisson_distribution<_IntType, _RealType>::
|
||||
operator()(_UniformRandomNumberGenerator& __urng)
|
||||
{
|
||||
#if _GLIBCXX_USE_C99_MATH_TR1
|
||||
if (_M_large)
|
||||
{
|
||||
_RealType __x;
|
||||
|
||||
const _RealType __m = std::floor(_M_mean);
|
||||
// sqrt(mu * pi / 2)
|
||||
const _RealType __c1 = (_M_sm
|
||||
* 1.2533141373155002512078826424055226L);
|
||||
const _RealType __c2 = _M_c2b + __c1;
|
||||
const _RealType __c3 = __c2 + 1;
|
||||
const _RealType __c4 = __c3 + 1;
|
||||
// c4 + e^(1 / 78)
|
||||
const _RealType __c5 = (__c4
|
||||
+ 1.0129030479320018583185514777512983L);
|
||||
const _RealType __c = _M_cb + __c5;
|
||||
const _RealType __cx = 2 * (2 * __m + _M_d);
|
||||
|
||||
normal_distribution<_RealType> __nd;
|
||||
|
||||
bool __keepgoing = true;
|
||||
do
|
||||
{
|
||||
const _RealType __u = __c * __urng();
|
||||
const _RealType __e = -std::log(__urng());
|
||||
|
||||
_RealType __w = 0.0;
|
||||
|
||||
if (__u <= __c1)
|
||||
{
|
||||
const _RealType __n = __nd(__urng);
|
||||
const _RealType __y = -std::abs(__n) * _M_sm - 1;
|
||||
__x = std::floor(__y);
|
||||
__w = -__n * __n / 2;
|
||||
if (__x < -__m)
|
||||
continue;
|
||||
}
|
||||
else if (__u <= __c2)
|
||||
{
|
||||
const _RealType __n = __nd(__urng);
|
||||
const _RealType __y = 1 + std::abs(__n) * _M_scx4;
|
||||
__x = std::ceil(__y);
|
||||
__w = __y * (2 - __y) * _M_2cx;
|
||||
if (__x > _M_d)
|
||||
continue;
|
||||
}
|
||||
else if (__u <= __c3)
|
||||
// XXX This case not in the book, nor in the Errata...
|
||||
__x = -1;
|
||||
else if (__u <= __c4)
|
||||
__x = 0;
|
||||
else if (__u <= __c5)
|
||||
__x = 1;
|
||||
else
|
||||
{
|
||||
const _RealType __v = -std::log(__urng());
|
||||
const _RealType __y = _M_d + __v * __cx / _M_d;
|
||||
__x = std::ceil(__y);
|
||||
__w = -_M_d * _M_2cx * (1 + __y / 2);
|
||||
}
|
||||
|
||||
__keepgoing = (__w - __e - __x * _M_lm_thr
|
||||
> _M_lfm - std::tr1::lgamma(__x + __m + 1));
|
||||
|
||||
} while (__keepgoing);
|
||||
|
||||
return _IntType(std::tr1::round(__x + __m));
|
||||
}
|
||||
else
|
||||
#endif
|
||||
{
|
||||
_IntType __x = -1;
|
||||
_RealType __prod = 1.0;
|
||||
|
||||
do
|
||||
{
|
||||
__prod *= __urng();
|
||||
__x += 1;
|
||||
}
|
||||
while (__prod > _M_lm_thr);
|
||||
|
||||
return __x;
|
||||
}
|
||||
}
|
||||
|
||||
template<typename _IntType, typename _RealType,
|
||||
typename _CharT, typename _Traits>
|
||||
std::basic_ostream<_CharT, _Traits>&
|
||||
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
||||
const poisson_distribution<_IntType, _RealType>& __x)
|
||||
{
|
||||
const std::ios_base::fmtflags __flags = __os.flags();
|
||||
const _CharT __fill = __os.fill();
|
||||
const std::streamsize __precision = __os.precision();
|
||||
const _CharT __space = __os.widen(' ');
|
||||
__os.flags(std::ios_base::scientific | std::ios_base::left);
|
||||
__os.fill(__space);
|
||||
__os.precision(_Max_digits10<_RealType>::__value);
|
||||
|
||||
__os << __x._M_large << __space << __x.mean()
|
||||
<< __space << __x._M_lm_thr;
|
||||
#if _GLIBCXX_USE_C99_MATH_TR1
|
||||
if (__x._M_large)
|
||||
__os << __space << __x._M_lfm << __space << __x._M_sm
|
||||
<< __space << __x._M_d << __space << __x._M_scx4
|
||||
<< __space << __x._M_2cx << __space << __x._M_c2b
|
||||
<< __space << __x._M_cb;
|
||||
#endif
|
||||
|
||||
__os.flags(__flags);
|
||||
__os.fill(__fill);
|
||||
__os.precision(__precision);
|
||||
return __os;
|
||||
}
|
||||
|
||||
template<typename _IntType, typename _RealType,
|
||||
typename _CharT, typename _Traits>
|
||||
std::basic_istream<_CharT, _Traits>&
|
||||
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
||||
poisson_distribution<_IntType, _RealType>& __x)
|
||||
{
|
||||
const std::ios_base::fmtflags __flags = __is.flags();
|
||||
__is.flags(std::ios_base::skipws);
|
||||
|
||||
__is >> __x._M_large >> __x._M_mean >> __x._M_lm_thr;
|
||||
#if _GLIBCXX_USE_C99_MATH_TR1
|
||||
if (__x._M_large)
|
||||
__is >> __x._M_lfm >> __x._M_sm >> __x._M_d >> __x._M_scx4
|
||||
>> __x._M_2cx >> __x._M_c2b >> __x._M_cb;
|
||||
#endif
|
||||
|
||||
__is.flags(__flags);
|
||||
return __is;
|
||||
}
|
||||
|
||||
|
||||
template<typename _RealType, typename _CharT, typename _Traits>
|
||||
std::basic_ostream<_CharT, _Traits>&
|
||||
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
||||
|
|
@ -766,10 +953,11 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
|
|||
__os.fill(__space);
|
||||
__os.precision(_Max_digits10<_RealType>::__value);
|
||||
|
||||
__os << __x.mean() << __space
|
||||
<< __x.sigma() << __space
|
||||
<< __x._M_saved << __space
|
||||
<< __x._M_saved_available;
|
||||
__os << __x._M_saved_available << __space
|
||||
<< __x.mean() << __space
|
||||
<< __x.sigma();
|
||||
if (__x._M_saved_available)
|
||||
__os << __space << __x._M_saved;
|
||||
|
||||
__os.flags(__flags);
|
||||
__os.fill(__fill);
|
||||
|
|
@ -785,8 +973,10 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
|
|||
const std::ios_base::fmtflags __flags = __is.flags();
|
||||
__is.flags(std::ios_base::dec | std::ios_base::skipws);
|
||||
|
||||
__is >> __x._M_mean >> __x._M_sigma
|
||||
>> __x._M_saved >> __x._M_saved_available;
|
||||
__is >> __x._M_saved_available >> __x._M_mean
|
||||
>> __x._M_sigma;
|
||||
if (__x._M_saved_available)
|
||||
__is >> __x._M_saved;
|
||||
|
||||
__is.flags(__flags);
|
||||
return __is;
|
||||
|
|
|
|||
|
|
@ -0,0 +1,37 @@
|
|||
// { dg-do compile }
|
||||
//
|
||||
// 2006-08-13 Paolo Carlini <pcarlini@suse.de>
|
||||
//
|
||||
// Copyright (C) 2006 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 2, 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 COPYING. If not, write to the Free
|
||||
// Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
|
||||
// USA.
|
||||
|
||||
// 5.1.7.4 Class template poisson_distribution [tr.rand.dist.pois]
|
||||
// 5.1.1 [7] Table 17
|
||||
|
||||
#include <tr1/random>
|
||||
|
||||
void
|
||||
test01()
|
||||
{
|
||||
using namespace std::tr1;
|
||||
|
||||
typedef poisson_distribution<int, double> test_type;
|
||||
|
||||
typedef test_type::input_type input_type;
|
||||
typedef test_type::result_type result_type;
|
||||
}
|
||||
Loading…
Reference in New Issue