OpenE2K
/
gcc
2
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

random (class binomial_distribution<>): Add. 

2006-08-18  Paolo Carlini  <pcarlini@suse.de>

	* include/tr1/random (class binomial_distribution<>): Add.
	* include/tr1/random.tcc (binomial_distribution<>::operator(),
	operator<<(std::basic_ostream<>&, const binomial_distribution<>&),
	operator>>(std::basic_istream<>&, binomial_distribution<>&,
	binomial_distribution<>::_M_waiting(), binomial_distribution<>::
	_M_initialize()): Define.
	* testsuite/tr1/5_numerical_facilities/random/binomial_distribution/
	requirements/typedefs.cc: New.

	* include/tr1/random (geometric_distribution<>::
	geometric_distribution(const _RealType&)): Fix DEBUG_ASSERT
	limits.

	* include/tr1/random (poisson_distribution): Add normal_distribution
	member, adjust consistently; minor tweaks and rearrangements of the
	arithmetic.
	(operator>>(std::basic_istream<>&, poisson_distribution<>&)): Move
	out of line.
	* include/tr1/random.tcc: Adjust.

	* include/tr1/random.tcc (normal_distribution<>::operator()): Minor
	tweaks.

From-SVN: r116245
This commit is contained in:
Paolo Carlini 2006-08-18 17:15:43 +00:00 committed by Paolo Carlini
parent 7867a3f739
commit 42031254bf
4 changed files with 473 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
libstdc++-v3/ChangeLog
View File

@ -1,3 +1,28 @@
2006-08-18 Paolo Carlini <pcarlini@suse.de>
* include/tr1/random (class binomial_distribution<>): Add.
* include/tr1/random.tcc (binomial_distribution<>::operator(),
operator<<(std::basic_ostream<>&, const binomial_distribution<>&),
operator>>(std::basic_istream<>&, binomial_distribution<>&,
binomial_distribution<>::_M_waiting(), binomial_distribution<>::
_M_initialize()): Define.
* testsuite/tr1/5_numerical_facilities/random/binomial_distribution/
requirements/typedefs.cc: New.
* include/tr1/random (geometric_distribution<>::
geometric_distribution(const _RealType&)): Fix DEBUG_ASSERT
limits.
* include/tr1/random (poisson_distribution): Add normal_distribution
member, adjust consistently; minor tweaks and rearrangements of the
arithmetic.
(operator>>(std::basic_istream<>&, poisson_distribution<>&)): Move
out of line.
* include/tr1/random.tcc: Adjust.
* include/tr1/random.tcc (normal_distribution<>::operator()): Minor
tweaks.
2006-08-18 Paolo Carlini <pcarlini@suse.de>
PR libstdc++/28765

150
libstdc++-v3/include/tr1/random
View File

@ -1588,9 +1588,9 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
// constructors and member function
explicit
geometric_distribution(const _RealType& __p = _RealType(0.5))
: _M_p(__p), _M_log_p(std::log(_M_p))
: _M_p(__p), _M_log_p(std::log(__p))
{
_GLIBCXX_DEBUG_ASSERT((_M_p >= 0.0) && (_M_p <= 1.0));
_GLIBCXX_DEBUG_ASSERT((_M_p > 0.0) && (_M_p < 1.0));
}
/**
@ -1649,6 +1649,9 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
};
template<typename _RealType>
class normal_distribution;
/**
* @brief A discrete Poisson random number distribution.
*
@ -1665,6 +1668,12 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
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
{
@ -1676,7 +1685,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
// constructors and member function
explicit
poisson_distribution(const _RealType& __mean = _RealType(1))
: _M_mean(__mean)
: _M_mean(__mean), _M_nd()
{
_GLIBCXX_DEBUG_ASSERT(_M_mean > 0.0);
_M_initialize();
@ -1690,7 +1699,8 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
{ return _M_mean; }
void
reset() { }
reset()
{ _M_nd.reset(); }
template<class _UniformRandomNumberGenerator>
result_type
@ -1721,29 +1731,145 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
*
* @returns The input stream with @p __x extracted or in an error state.
*/
template<typename _CharT, typename _Traits>
template<typename _IntType1, typename _RealType1,
typename _CharT, typename _Traits>
friend std::basic_istream<_CharT, _Traits>&
operator>>(std::basic_istream<_CharT, _Traits>& __is,
poisson_distribution& __x)
{
__is >> __x._M_mean;
__x._M_initialize();
return __is;
}
poisson_distribution<_IntType1, _RealType1>& __x);
private:
void
_M_initialize();
// NB: Unused when _GLIBCXX_USE_C99_MATH_TR1 is undefined.
normal_distribution<_RealType> _M_nd;
_RealType _M_mean;
// _M_lm_thr hosts either log(mean) or the threshold of the simple
// method.
_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;
_RealType _M_lfm, _M_sm, _M_d, _M_scx, _M_1cx, _M_c2b, _M_cb;
#endif
};
/**
* @brief A discrete binomial random number distribution.
*
* The formula for the binomial probability mass function is
* @f$ p(i) = \binom{n}{i} p^i (1 - p)^{t - i} @f$ where @f$ t @f$
* and @f$ p @f$ are the parameters of the distribution.
*/
template<typename _IntType = int, typename _RealType = double>
class binomial_distribution;
template<typename _IntType, typename _RealType,
typename _CharT, typename _Traits>
std::basic_ostream<_CharT, _Traits>&
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
const binomial_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,
binomial_distribution<_IntType, _RealType>& __x);
template<typename _IntType, typename _RealType>
class binomial_distribution
{
public:
// types
typedef _RealType input_type;
typedef _IntType result_type;
// constructors and member function
explicit
binomial_distribution(_IntType __t = 1,
const _RealType& __p = _RealType(0.5))
: _M_t(__t), _M_p(__p), _M_nd()
{
_GLIBCXX_DEBUG_ASSERT((_M_t >= 0) && (_M_p >= 0.0) && (_M_p <= 1.0));
_M_initialize();
}
/**
* Gets the distribution @p t parameter.
*/
_IntType
t() const
{ return _M_t; }
/**
* Gets the distribution @p p parameter.
*/
_RealType
p() const
{ return _M_p; }
void
reset()
{ _M_nd.reset(); }
template<class _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng);
/**
* Inserts a %binomial_distribution random number distribution
* @p __x into the output stream @p __os.
*
* @param __os An output stream.
* @param __x A %binomial_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 binomial_distribution<_IntType1, _RealType1>& __x);
/**
* Extracts a %binomial_distribution random number distribution
* @p __x from the input stream @p __is.
*
* @param __is An input stream.
* @param __x A %binomial_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,
binomial_distribution<_IntType1, _RealType1>& __x);
private:
void
_M_initialize();
template<class _UniformRandomNumberGenerator>
result_type
_M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t);
// NB: Unused when _GLIBCXX_USE_C99_MATH_TR1 is undefined.
normal_distribution<_RealType> _M_nd;
_IntType _M_t;
_RealType _M_p;
_RealType _M_q;
#if _GLIBCXX_USE_C99_MATH_TR1
_RealType _M_d1, _M_d2, _M_s1, _M_s2, _M_c,
_M_a1, _M_a123, _M_s, _M_lf, _M_lp1p;
#endif
bool _M_easy;
};
/* @} */ // group tr1_random_distributions_discrete
/**

309
libstdc++-v3/include/tr1/random.tcc
View File

@ -667,21 +667,18 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
_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));
const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m
/ __pi_4));
_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 __cx = 2 * __m + _M_d;
_M_scx = std::sqrt(__cx / 2);
_M_1cx = 1 / __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;
_M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
_M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d;
}
else
#endif
@ -720,11 +717,9 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
const _RealType __e178 = 1.0129030479320018583185514777512983L;
const _RealType __c5 = __c4 + __e178;
const _RealType __c = _M_cb + __c5;
const _RealType __cx = 2 * (2 * __m + _M_d);
const _RealType __2cx = 2 * (2 * __m + _M_d);
normal_distribution<_RealType> __nd;
bool __keepgoing = true;
bool __reject = true;
do
{
const _RealType __u = __c * __urng();
@ -734,7 +729,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
if (__u <= __c1)
{
const _RealType __n = __nd(__urng);
const _RealType __n = _M_nd(__urng);
const _RealType __y = -std::abs(__n) * _M_sm - 1;
__x = std::floor(__y);
__w = -__n * __n / 2;
@ -743,10 +738,10 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
}
else if (__u <= __c2)
{
const _RealType __n = __nd(__urng);
const _RealType __y = 1 + std::abs(__n) * _M_scx4;
const _RealType __n = _M_nd(__urng);
const _RealType __y = 1 + std::abs(__n) * _M_scx;
__x = std::ceil(__y);
__w = __y * (2 - __y) * _M_2cx;
__w = __y * (2 - __y) * _M_1cx;
if (__x > _M_d)
continue;
}
@ -760,22 +755,22 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
else
{
const _RealType __v = -std::log(__urng());
const _RealType __y = _M_d + __v * __cx / _M_d;
const _RealType __y = _M_d + __v * __2cx / _M_d;
__x = std::ceil(__y);
__w = -_M_d * _M_2cx * (1 + __y / 2);
__w = -_M_d * _M_1cx * (1 + __y / 2);
}
__keepgoing = (__w - __e - __x * _M_lm_thr
> _M_lfm - std::tr1::lgamma(__x + __m + 1));
__reject = (__w - __e - __x * _M_lm_thr
> _M_lfm - std::tr1::lgamma(__x + __m + 1));
} while (__keepgoing);
} while (__reject);
return _IntType(std::tr1::round(__x + __m));
return _IntType(__x + __m + 0.5);
}
else
#endif
{
_IntType __x = -1;
_IntType __x = -1;
_RealType __prod = 1.0;
do
@ -798,11 +793,12 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
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(__os.widen(' '));
__os.fill(__space);
__os.precision(_Max_digits10<_RealType>::__value);
__os << __x.mean();
__os << __x.mean() << __space << __x._M_nd;
__os.flags(__flags);
__os.fill(__fill);
@ -810,6 +806,247 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
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_mean >> __x._M_nd;
__x._M_initialize();
__is.flags(__flags);
return __is;
}
template<typename _IntType, typename _RealType>
void
binomial_distribution<_IntType, _RealType>::
_M_initialize()
{
const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
_M_easy = true;
#if _GLIBCXX_USE_C99_MATH_TR1
if (_M_t * __p12 >= 8)
{
_M_easy = false;
const _RealType __np = std::floor(_M_t * __p12);
const _RealType __pa = __np / _M_t;
const _RealType __1p = 1 - __pa;
const _RealType __pi_4 = 0.7853981633974483096156608458198757L;
const _RealType __d1x =
std::sqrt(__np * __1p * std::log(32 * __np
/ (81 * __pi_4 * __1p)));
_M_d1 = std::tr1::round(std::max(_RealType(1), __d1x));
const _RealType __d2x =
std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
/ (__pi_4 * __pa)));
_M_d2 = std::tr1::round(std::max(_RealType(1), __d2x));
// sqrt(pi / 2)
const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
_M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
_M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
_M_c = 2 * _M_d1 / __np;
_M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
const _RealType __a12 = _M_a1 + _M_s2 * __spi_2;
const _RealType __s1s = _M_s1 * _M_s1;
_M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
* 2 * __s1s / _M_d1
* std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
const _RealType __s2s = _M_s2 * _M_s2;
_M_s = (_M_a123 + 2 * __s2s / _M_d2
* std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
_M_lf = (std::tr1::lgamma(__np + 1)
+ std::tr1::lgamma(_M_t - __np + 1));
_M_lp1p = std::log(__pa / __1p);
_M_q = -std::log(1 - (__p12 - __pa) / __1p);
}
else
#endif
_M_q = -std::log(1 - __p12);
}
template<typename _IntType, typename _RealType>
template<class _UniformRandomNumberGenerator>
typename binomial_distribution<_IntType, _RealType>::result_type
binomial_distribution<_IntType, _RealType>::
_M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t)
{
_IntType __x = 0;
_RealType __sum = 0;
do
{
const _RealType __e = -std::log(__urng());
__sum += __e / (__t - __x);
__x += 1;
}
while (__sum <= _M_q);
return __x - 1;
}
/**
* A rejection algorithm when t * p >= 8 and a simple waiting time
* method - the second in the referenced book - 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, Sect. 4 (+ Errata!).
*/
template<typename _IntType, typename _RealType>
template<class _UniformRandomNumberGenerator>
typename binomial_distribution<_IntType, _RealType>::result_type
binomial_distribution<_IntType, _RealType>::
operator()(_UniformRandomNumberGenerator& __urng)
{
result_type __ret;
const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
#if _GLIBCXX_USE_C99_MATH_TR1
if (!_M_easy)
{
_RealType __x;
const _RealType __np = std::floor(_M_t * __p12);
const _RealType __pa = __np / _M_t;
// sqrt(pi / 2)
const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
const _RealType __a1 = _M_a1;
const _RealType __a12 = __a1 + _M_s2 * __spi_2;
const _RealType __a123 = _M_a123;
const _RealType __s1s = _M_s1 * _M_s1;
const _RealType __s2s = _M_s2 * _M_s2;
bool __reject;
do
{
const _RealType __u = _M_s * __urng();
_RealType __v;
if (__u <= __a1)
{
const _RealType __n = _M_nd(__urng);
const _RealType __y = _M_s1 * std::abs(__n);
__reject = __y >= _M_d1;
if (!__reject)
{
const _RealType __e = -std::log(__urng());
__x = std::floor(__y);
__v = -__e - __n * __n / 2 + _M_c;
}
}
else if (__u <= __a12)
{
const _RealType __n = _M_nd(__urng);
const _RealType __y = _M_s2 * std::abs(__n);
__reject = __y >= _M_d2;
if (!__reject)
{
const _RealType __e = -std::log(__urng());
__x = std::floor(-__y);
__v = -__e - __n * __n / 2;
}
}
else if (__u <= __a123)
{
const _RealType __e1 = -std::log(__urng());
const _RealType __e2 = -std::log(__urng());
const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1;
__x = std::floor(__y);
__v = (-__e2 + _M_d1 * (1 / (_M_t - __np)
-__y / (2 * __s1s)));
__reject = false;
}
else
{
const _RealType __e1 = -std::log(__urng());
const _RealType __e2 = -std::log(__urng());
const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2;
__x = std::floor(-__y);
__v = -__e2 - _M_d2 * __y / (2 * __s2s);
__reject = false;
}
__reject = __reject || __x < -__np || __x > _M_t - __np;
if (!__reject)
{
const _RealType __lfx =
std::tr1::lgamma(__np + __x + 1)
+ std::tr1::lgamma(_M_t - (__np + __x) + 1);
__reject = __v > _M_lf - __lfx + __x * _M_lp1p;
}
}
while (__reject);
__x += __np + 0.5;
const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x));
__ret = _IntType(__x) + __z;
}
else
#endif
__ret = _M_waiting(__urng, _M_t);
if (__p12 != _M_p)
__ret = _M_t - __ret;
return __ret;
}
template<typename _IntType, typename _RealType,
typename _CharT, typename _Traits>
std::basic_ostream<_CharT, _Traits>&
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
const binomial_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.t() << __space << __x.p()
<< __space << __x._M_nd;
__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,
binomial_distribution<_IntType, _RealType>& __x)
{
const std::ios_base::fmtflags __flags = __is.flags();
__is.flags(std::ios_base::dec | std::ios_base::skipws);
__is >> __x._M_t >> __x._M_p >> __x._M_nd;
__x._M_initialize();
__is.flags(__flags);
return __is;
}
template<typename _RealType, typename _CharT, typename _Traits>
std::basic_ostream<_CharT, _Traits>&
@ -893,20 +1130,20 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
result_type __x, __y, __r2;
do
{
__x = result_type(2.0) * __urng() - result_type(1.0);
__y = result_type(2.0) * __urng() - result_type(1.0);
__x = result_type(2.0) * __urng() - 1.0;
__y = result_type(2.0) * __urng() - 1.0;
__r2 = __x * __x + __y * __y;
}
while (__r2 > result_type(1.0) || __r2 == result_type(0));
while (__r2 > 1.0 || __r2 == 0.0);
const result_type __mult = std::sqrt(-result_type(2.0)
* std::log(__r2) / __r2);
const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
_M_saved = __x * __mult;
_M_saved_available = true;
__ret = __y * __mult;
}
return __ret * _M_sigma + _M_mean;
__ret = __ret * _M_sigma + _M_mean;
return __ret;
}
template<typename _RealType, typename _CharT, typename _Traits>

37
libstdc++-v3/testsuite/tr1/5_numerical_facilities/random/binomial_distribution/requirements/typedefs.cc Normal file
View File

@ -0,0 +1,37 @@
// { dg-do compile }
//
// 2006-08-18 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.5 Class template binomial_distribution [tr.rand.dist.bin]
// 5.1.1 [7] Table 17
#include <tr1/random>
void
test01()
{
using namespace std::tr1;
typedef binomial_distribution<int, double> test_type;
typedef test_type::input_type input_type;
typedef test_type::result_type result_type;
}