mirror/gcc
2
0
Fork 0
mirror of git://gcc.gnu.org/git/gcc.git synced 2025-02-16 12:20:11 +08:00
Code Issues Packages Projects Releases Wiki Activity

random.tcc (gamma_distribution<>::operator() (_UniformRandomNumberGenerator&, const param_type&): Redo, using the Marsaglia/Tsang algorithm.

Browse Source 
2009-06-08  Paolo Carlini  <paolo.carlini@oracle.com>

	* include/bits/random.tcc (gamma_distribution<>::operator()
	(_UniformRandomNumberGenerator&, const param_type&): Redo, using
	the Marsaglia/Tsang algorithm.
	(gamma_distribution<>::param_type::_M_initialize): Adjust.
	(operator<<(basic_ostream<>&, gamma_distribution<>),
	operator>>(basic_ostream<>&, gamma_distribution<>): Likewise.

	* include/bits/random.tcc(student_t_distribution<>::_M_gaussian):
	Remove, just use normal_distribution.
	(operator<<(basic_ostream<>&, student_t_distribution<>),
	operator>>(basic_ostream<>&, student_t_distribution<>): Adjust.
	(linear_congruential_engine<>::operator()()): Move inline.
	(lognormal_distribution<>::operator()(_UniformRandomNumberGenerator&,
	const param_type&)): Move inline, just use normal_distribution.
	(operator<<(basic_ostream<>&, lognormal_distribution<>),
	operator>>(basic_ostream<>&, lognormal_distribution<>): Adjust.
	(weibull_distribution<>::operator()(_UniformRandomNumberGenerator&,
	const param_type&)): Move here, out of line.
	(piecewise_constant_distribution<>::param_type::param_type()): Move
	inline.
	* include/bits/random.h: Adjust, minor tweaks.

From-SVN: r148276
This commit is contained in:
Paolo Carlini 2009-06-08 14:38:48 +00:00 committed by Paolo Carlini
parent 06ddd8716e
commit b01630bb3d
3 changed files with 120 additions and 217 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
libstdc++-v3/ChangeLog
View File

@ -1,3 +1,27 @@
2009-06-08 Paolo Carlini <paolo.carlini@oracle.com>
* include/bits/random.tcc (gamma_distribution<>::operator()
(_UniformRandomNumberGenerator&, const param_type&): Redo, using
the Marsaglia/Tsang algorithm.
(gamma_distribution<>::param_type::_M_initialize): Adjust.
(operator<<(basic_ostream<>&, gamma_distribution<>),
operator>>(basic_ostream<>&, gamma_distribution<>): Likewise.
* include/bits/random.tcc(student_t_distribution<>::_M_gaussian):
Remove, just use normal_distribution.
(operator<<(basic_ostream<>&, student_t_distribution<>),
operator>>(basic_ostream<>&, student_t_distribution<>): Adjust.
(linear_congruential_engine<>::operator()()): Move inline.
(lognormal_distribution<>::operator()(_UniformRandomNumberGenerator&,
const param_type&)): Move inline, just use normal_distribution.
(operator<<(basic_ostream<>&, lognormal_distribution<>),
operator>>(basic_ostream<>&, lognormal_distribution<>): Adjust.
(weibull_distribution<>::operator()(_UniformRandomNumberGenerator&,
const param_type&)): Move here, out of line.
(piecewise_constant_distribution<>::param_type::param_type()): Move
inline.
* include/bits/random.h: Adjust, minor tweaks.
2009-06-05 Benjamin Kosnik <bkoz@redhat.com>
* testsuite/29_atomics/atomic_address/cons/aggregate.cc: Remove xfail.

84
libstdc++-v3/include/bits/random.h
View File

@ -268,7 +268,11 @@ namespace std
* @brief Gets the next random number in the sequence.
*/
result_type
operator()();
operator()()
{
_M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
return _M_x;
}
/**
* @brief Compares two linear congruential random number generator
@ -1588,12 +1592,7 @@ namespace std
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
{
typedef typename _UniformRandomNumberGenerator::result_type
_UResult_type;
return _M_call(__urng, this->a(), this->b(),
typename is_integral<_UResult_type>::type());
}
{ return this->operator()(__urng, this->param()); }
/**
* Gets a uniform random number in the range @f$[0, n)@f$.
@ -1765,11 +1764,7 @@ namespace std
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
{
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
__aurng(__urng);
return (__aurng() * (this->b() - this->a())) + this->a();
}
{ return this->operator()(__urng, this->param()); }
template<typename _UniformRandomNumberGenerator>
result_type
@ -2014,12 +2009,12 @@ namespace std
explicit
lognormal_distribution(_RealType __m = _RealType(0),
_RealType __s = _RealType(1))
: _M_param(__m, __s)
: _M_param(__m, __s), _M_nd()
{ }
explicit
lognormal_distribution(const param_type& __p)
: _M_param(__p)
: _M_param(__p), _M_nd()
{ }
/**
@ -2027,7 +2022,7 @@ namespace std
*/
void
reset()
{ }
{ _M_nd.reset(); }
/**
*
@ -2077,10 +2072,13 @@ namespace std
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p);
const param_type& __p)
{ return std::exp(__p.s() * _M_nd(__urng) + __p.m()); }
private:
param_type _M_param;
normal_distribution<result_type> _M_nd;
};
/**
@ -2555,12 +2553,12 @@ namespace std
explicit
student_t_distribution(_RealType __n = _RealType(1))
: _M_param(__n)
: _M_param(__n), _M_nd()
{ }
explicit
student_t_distribution(const param_type& __p)
: _M_param(__p)
: _M_param(__p), _M_nd()
{ }
/**
@ -2568,7 +2566,7 @@ namespace std
*/
void
reset()
{ }
{ _M_nd.reset(); }
/**
*
@ -2617,12 +2615,9 @@ namespace std
const param_type& __p);
private:
template<typename _UniformRandomNumberGenerator>
result_type
_M_gaussian(_UniformRandomNumberGenerator& __urng,
const result_type __sigma);
param_type _M_param;
normal_distribution<result_type> _M_nd;
};
/**
@ -2761,14 +2756,7 @@ namespace std
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
{
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
__aurng(__urng);
if ((__aurng() - __aurng.min())
< this->p() * (__aurng.max() - __aurng.min()))
return true;
return false;
}
{ return this->operator()(__urng, this->param()); }
template<typename _UniformRandomNumberGenerator>
result_type
@ -3539,11 +3527,7 @@ namespace std
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
{
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
__aurng(__urng);
return -std::log(__aurng()) / this->lambda();
}
{ return this->operator()(__urng, this->param()); }
template<typename _UniformRandomNumberGenerator>
result_type
@ -3633,8 +3617,7 @@ namespace std
_RealType _M_alpha;
_RealType _M_beta;
// Hosts either lambda of GB or d of modified Vaduva's.
_RealType _M_l_d;
_RealType _M_malpha, _M_a2;
};
public:
@ -3645,21 +3628,20 @@ namespace std
explicit
gamma_distribution(_RealType __alpha_val = _RealType(1),
_RealType __beta_val = _RealType(1))
: _M_param(__alpha_val, __beta_val)
: _M_param(__alpha_val, __beta_val), _M_nd()
{ }
explicit
gamma_distribution(const param_type& __p)
: _M_param(__p)
: _M_param(__p), _M_nd()
{ }
/**
* @brief Resets the distribution state.
*
* Does nothing for the gamma distribution.
*/
void
reset() { }
reset()
{ _M_nd.reset(); }
/**
* @brief Returns the @f$ \alpha @f$ of the distribution.
@ -3716,6 +3698,8 @@ namespace std
private:
param_type _M_param;
normal_distribution<result_type> _M_nd;
};
/**
@ -3854,13 +3838,7 @@ namespace std
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
__aurng(__urng);
return __p.b() * std::pow(-std::log(__aurng()),
result_type(1) / __p.a());
}
const param_type& __p);
private:
param_type _M_param;
@ -4222,7 +4200,9 @@ namespace std
typedef piecewise_constant_distribution<_RealType> distribution_type;
friend class piecewise_constant_distribution<_RealType>;
param_type();
param_type()
: _M_int(), _M_den(), _M_cp()
{ _M_initialize(); }
template<typename _InputIteratorB, typename _InputIteratorW>
param_type(_InputIteratorB __bfirst,

229
libstdc++-v3/include/bits/random.tcc
View File

@ -128,19 +128,6 @@ namespace std
seed(__sum);
}
/**
* Gets the next generated value in sequence.
*/
template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
typename linear_congruential_engine<_UIntType, __a, __c, __m>::
result_type
linear_congruential_engine<_UIntType, __a, __c, __m>::
operator()()
{
_M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
return _M_x;
}
template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
typename _CharT, typename _Traits>
std::basic_ostream<_CharT, _Traits>&
@ -556,7 +543,7 @@ namespace std
{
const long double __r = static_cast<long double>(_M_b.max())
- static_cast<long double>(_M_b.min()) + 1.0L;
const result_type __m = std::log10(__r) / std::log10(2.0L);
const result_type __m = std::log(__r) / std::log(2.0L);
result_type __n, __n0, __y0, __y1, __s0, __s1;
for (size_t __i = 0; __i < 2; ++__i)
{
@ -874,17 +861,12 @@ namespace std
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{
typename gamma_distribution<>::param_type
__gamma_param(__p.k(), 1.0);
gamma_distribution<> __gamma(__gamma_param);
gamma_distribution<> __gamma(__p.k(), 1.0);
double __x = __gamma(__urng);
typename poisson_distribution<result_type>::param_type
__poisson_param(__x * __p.p() / (1.0 - __p.p()));
poisson_distribution<result_type> __poisson(__poisson_param);
result_type __m = __poisson(__urng);
return __m;
poisson_distribution<result_type> __poisson(__x * __p.p()
/ (1.0 - __p.p()));
return __poisson(__urng);
}
template<typename _IntType, typename _CharT, typename _Traits>
@ -1510,33 +1492,6 @@ namespace std
}
template<typename _RealType>
template<typename _UniformRandomNumberGenerator>
typename lognormal_distribution<_RealType>::result_type
lognormal_distribution<_RealType>::
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{
_RealType __u, __v, __r2, __normal;
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
__aurng(__urng);
do
{
// Choose x,y in uniform square (-1,-1) to (+1,+1).
__u = 2 * __aurng() - 1;
__v = 2 * __aurng() - 1;
// See if it is in the unit circle.
__r2 = __u * __u + __v * __v;
}
while (__r2 > 1 || __r2 == 0);
__normal = __u * std::sqrt(-2 * std::log(__r2) / __r2);
return std::exp(__p.s() * __normal + __p.m());
}
template<typename _RealType, typename _CharT, typename _Traits>
std::basic_ostream<_CharT, _Traits>&
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
@ -1553,7 +1508,8 @@ namespace std
__os.fill(__space);
__os.precision(std::numeric_limits<_RealType>::digits10 + 1);
__os << __x.m() << __space << __x.s();
__os << __x.m() << __space << __x.s()
<< __space << __x._M_nd;
__os.flags(__flags);
__os.fill(__fill);
@ -1573,7 +1529,7 @@ namespace std
__is.flags(__ios_base::dec | __ios_base::skipws);
_RealType __m, __s;
__is >> __m >> __s;
__is >> __m >> __s >> __x._M_nd;
__x.param(typename lognormal_distribution<_RealType>::
param_type(__m, __s));
@ -1589,9 +1545,7 @@ namespace std
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{
typename gamma_distribution<_RealType>::param_type
__gamma_param(__p.n() / 2, 1.0);
gamma_distribution<_RealType> __gamma(__gamma_param);
gamma_distribution<_RealType> __gamma(__p.n() / 2, 1.0);
return 2 * __gamma(__urng);
}
@ -1765,35 +1719,6 @@ namespace std
}
//
// This could be operator() for a Gaussian distribution.
//
template<typename _RealType>
template<typename _UniformRandomNumberGenerator>
typename student_t_distribution<_RealType>::result_type
student_t_distribution<_RealType>::
_M_gaussian(_UniformRandomNumberGenerator& __urng,
const result_type __sigma)
{
_RealType __x, __y, __r2;
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
__aurng(__urng);
do
{
// Choose x,y in uniform square (-1,-1) to (+1,+1).
__x = 2 * __aurng() - 1;
__y = 2 * __aurng() - 1;
// See if it is in the unit circle.
__r2 = __x * __x + __y * __y;
}
while (__r2 > 1 || __r2 == 0);
// Box-Muller transform.
return __sigma * __y * std::sqrt(-2 * std::log(__r2) / __r2);
}
template<typename _RealType>
template<typename _UniformRandomNumberGenerator>
typename student_t_distribution<_RealType>::result_type
@ -1803,10 +1728,8 @@ namespace std
{
if (__param.n() <= 2.0)
{
_RealType __y1 = _M_gaussian(__urng, 1.0);
typename chi_squared_distribution<_RealType>::param_type
__chisq_param(__param.n());
chi_squared_distribution<_RealType> __chisq(__chisq_param);
_RealType __y1 = _M_nd(__urng);
chi_squared_distribution<_RealType> __chisq(__param.n());
_RealType __y2 = __chisq(__urng);
return __y1 / std::sqrt(__y2 / __param.n());
@ -1814,13 +1737,12 @@ namespace std
else
{
_RealType __y1, __y2, __z;
exponential_distribution<_RealType>
__exponential(1.0 / (__param.n() / 2.0 - 1.0));
do
{
__y1 = _M_gaussian(__urng, 1.0);
typename exponential_distribution<_RealType>::param_type
__exp_param(1.0 / (__param.n() / 2.0 - 1.0));
exponential_distribution<_RealType>
__exponential(__exp_param);
__y1 = _M_nd(__urng);
__y2 = __exponential(__urng);
__z = __y1 * __y1 / (__param.n() - 2.0);
@ -1850,7 +1772,7 @@ namespace std
__os.fill(__space);
__os.precision(std::numeric_limits<_RealType>::digits10 + 1);
__os << __x.n();
__os << __x.n() << __space << __x._M_nd;
__os.flags(__flags);
__os.fill(__fill);
@ -1870,7 +1792,7 @@ namespace std
__is.flags(__ios_base::dec | __ios_base::skipws);
_RealType __n;
__is >> __n;
__is >> __n >> __x._M_nd;
__x.param(typename student_t_distribution<_RealType>::param_type(__n));
__is.flags(__flags);
@ -1883,28 +1805,16 @@ namespace std
gamma_distribution<_RealType>::param_type::
_M_initialize()
{
if (_M_alpha >= 1)
_M_l_d = std::sqrt(2 * _M_alpha - 1);
else
_M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
* (1 - _M_alpha));
_M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
_M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
}
/**
* Cheng's rejection algorithm GB for alpha >= 1 and a modification
* of Vaduva's rejection from Weibull algorithm due to Devroye for
* alpha < 1.
*
* References:
* Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral
* Shape Parameter." Applied Statistics, 26, 71-75, 1977.
*
* Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection
* and Composition Procedures." Math. Operationsforschung and Statistik,
* Series in Statistics, 8, 545-576, 1977.
*
* Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
* New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
* Marsaglia, G. and Tsang, W. W.
* "A Simple Method for Generating Gamma Variables"
* ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
*/
template<typename _RealType>
template<typename _UniformRandomNumberGenerator>
@ -1913,58 +1823,40 @@ namespace std
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __param)
{
result_type __x;
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
__aurng(__urng);
bool __reject;
const _RealType __alpha_val = __param.alpha();
const _RealType __beta_val = __param.beta();
if (__alpha_val >= 1)
result_type __u, __v, __n;
const result_type __a1 = (__param._M_malpha
- _RealType(1.0) / _RealType(3.0));
do
{
// alpha - log(4)
const result_type __b = __alpha_val
- result_type(1.3862943611198906188344642429163531L);
const result_type __c = __alpha_val + __param._M_l_d;
const result_type __1l = 1 / __param._M_l_d;
// 1 + log(9 / 2)
const result_type __k = 2.5040773967762740733732583523868748L;
do
{
const result_type __u = __aurng() / __beta_val;
const result_type __v = __aurng() / __beta_val;
const result_type __y = __1l * std::log(__v / (1 - __v));
__x = __alpha_val * std::exp(__y);
const result_type __z = __u * __v * __v;
const result_type __r = __b + __c * __y - __x;
__reject = __r < result_type(4.5) * __z - __k;
if (__reject)
__reject = __r < std::log(__z);
__n = _M_nd(__urng);
__v = result_type(1.0) + __param._M_a2 * __n;
}
while (__reject);
while (__v <= 0.0);
__v = __v * __v * __v;
__u = __aurng();
}
while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
&& (std::log(__u) > (0.5 * __n * __n + __a1
* (1.0 - __v + std::log(__v)))));
if (__param.alpha() == __param._M_malpha)
return __a1 * __v * __param.beta();
else
{
const result_type __c = 1 / __alpha_val;
do
{
const result_type __z = -std::log(__aurng() / __beta_val);
const result_type __e = -std::log(__aurng() / __beta_val);
__x = std::pow(__z, __c);
__reject = __z + __e < __param._M_l_d + __x;
}
while (__reject);
__u = __aurng();
while (__u == 0.0);
return (std::pow(__u, result_type(1.0) / __param.alpha())
* __a1 * __v * __param.beta());
}
return __beta_val * __x;
}
template<typename _RealType, typename _CharT, typename _Traits>
@ -1983,7 +1875,8 @@ namespace std
__os.fill(__space);
__os.precision(std::numeric_limits<_RealType>::digits10 + 1);
__os << __x.alpha() << __space << __x.beta();
__os << __x.alpha() << __space << __x.beta()
<< __space << __x._M_nd;
__os.flags(__flags);
__os.fill(__fill);
@ -2003,7 +1896,7 @@ namespace std
__is.flags(__ios_base::dec | __ios_base::skipws);
_RealType __alpha_val, __beta_val;
__is >> __alpha_val >> __beta_val;
__is >> __alpha_val >> __beta_val >> __x._M_nd;
__x.param(typename gamma_distribution<_RealType>::
param_type(__alpha_val, __beta_val));
@ -2012,6 +1905,19 @@ namespace std
}
template<typename _RealType>
template<typename _UniformRandomNumberGenerator>
typename weibull_distribution<_RealType>::result_type
weibull_distribution<_RealType>::
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
__aurng(__urng);
return __p.b() * std::pow(-std::log(__aurng()),
result_type(1) / __p.a());
}
template<typename _RealType, typename _CharT, typename _Traits>
std::basic_ostream<_CharT, _Traits>&
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
@ -2265,12 +2171,6 @@ namespace std
_M_cp[_M_cp.size() - 1] = 1.0;
}
template<typename _RealType>
piecewise_constant_distribution<_RealType>::param_type::
param_type()
: _M_int(), _M_den(), _M_cp()
{ _M_initialize(); }
template<typename _RealType>
template<typename _InputIteratorB, typename _InputIteratorW>
piecewise_constant_distribution<_RealType>::param_type::
@ -2315,8 +2215,7 @@ namespace std
template<typename _RealType>
template<typename _Func>
piecewise_constant_distribution<_RealType>::param_type::
param_type(size_t __nw, _RealType __xmin, _RealType __xmax,
_Func __fw)
param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
: _M_int(), _M_den(), _M_cp()
{
for (size_t __i = 0; __i <= __nw; ++__i)
@ -2713,7 +2612,7 @@ namespace std
__bits);
const long double __r = static_cast<long double>(__urng.max())
- static_cast<long double>(__urng.min()) + 1.0L;
const size_t __log2r = std::log10(__r) / std::log10(2.0L);
const size_t __log2r = std::log(__r) / std::log(2.0L);
size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
_RealType __sum = _RealType(0);
_RealType __tmp = _RealType(1);