diff options
| -rw-r--r-- | libstdc++-v3/ChangeLog | 24 | ||||
| -rw-r--r-- | libstdc++-v3/include/bits/random.h | 84 | ||||
| -rw-r--r-- | libstdc++-v3/include/bits/random.tcc | 229 |
3 files changed, 120 insertions, 217 deletions
diff --git a/libstdc++-v3/ChangeLog b/libstdc++-v3/ChangeLog index 9708d40..9d3bf47 100644 --- a/libstdc++-v3/ChangeLog +++ b/libstdc++-v3/ChangeLog @@ -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. diff --git a/libstdc++-v3/include/bits/random.h b/libstdc++-v3/include/bits/random.h index b5ccf8a..8a21ae5 100644 --- a/libstdc++-v3/include/bits/random.h +++ b/libstdc++-v3/include/bits/random.h @@ -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, diff --git a/libstdc++-v3/include/bits/random.tcc b/libstdc++-v3/include/bits/random.tcc index b110f99..b933f6d 100644 --- a/libstdc++-v3/include/bits/random.tcc +++ b/libstdc++-v3/include/bits/random.tcc @@ -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) - { - // 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; + result_type __u, __v, __n; + const result_type __a1 = (__param._M_malpha + - _RealType(1.0) / _RealType(3.0)); + 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 @@ 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, @@ -2266,12 +2172,6 @@ namespace std } 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 @@ 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); |