This is the mail archive of the
libstdc++@gcc.gnu.org
mailing list for the libstdc++ project.
[v3] Redo gamma_distribution, other tweaks
- From: Paolo Carlini <paolo dot carlini at oracle dot com>
- To: gcc-patches at gcc dot gnu dot org
- Cc: libstdc++ at gcc dot gnu dot org
- Date: Mon, 08 Jun 2009 16:42:56 +0200
- Subject: [v3] Redo gamma_distribution, other tweaks
Hi,
as promised, this is the first stab at re-implementing
gamma_distribution + other relatively minor fixes. For now the new code
is only lightly tested (inspected some accurate histograms for 10 values
of alpha for new vs old implementation), but I'm working on the
testsuite... Tested x86_64-linux multilib, committed to mainline.
Paolo.
//////////////////////
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.
Index: include/bits/random.tcc
===================================================================
--- include/bits/random.tcc (revision 148272)
+++ include/bits/random.tcc (working copy)
@@ -128,19 +128,6 @@
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 @@
{
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 @@
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 @@
}
- 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 @@
__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 @@
__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 @@
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,48 +1719,17 @@
}
- //
- // 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
- student_t_distribution<_RealType>::
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __param)
{
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 @@
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 @@
__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 @@
__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 @@
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 @@
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)
- {
- // 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;
+ result_type __u, __v, __n;
+ const result_type __a1 = (__param._M_malpha
+ - _RealType(1.0) / _RealType(3.0));
- // 1 + log(9 / 2)
- const result_type __k = 2.5040773967762740733732583523868748L;
-
+ do
+ {
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 @@
__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 @@
__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 @@
}
+ 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,
@@ -2266,12 +2172,6 @@
}
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::
param_type(_InputIteratorB __bbegin,
@@ -2315,8 +2215,7 @@
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 @@
__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);
Index: include/bits/random.h
===================================================================
--- include/bits/random.h (revision 148272)
+++ include/bits/random.h (working copy)
@@ -268,7 +268,11 @@
* @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 @@
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 @@
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 @@
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 @@
*/
void
reset()
- { }
+ { _M_nd.reset(); }
/**
*
@@ -2077,10 +2072,13 @@
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 @@
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 @@
*/
void
reset()
- { }
+ { _M_nd.reset(); }
/**
*
@@ -2617,12 +2615,9 @@
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 @@
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 @@
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 @@
_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 @@
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 @@
private:
param_type _M_param;
+
+ normal_distribution<result_type> _M_nd;
};
/**
@@ -3854,13 +3838,7 @@
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 @@
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,