From 96ddac7425dd56d6852f0ad53bf2283d9976aae8 Mon Sep 17 00:00:00 2001 From: Paolo Carlini Date: Sun, 20 Aug 2006 16:05:05 +0000 Subject: [PATCH] random (gamma_distribution<>::_M_initialize, [...]): Add. 2006-08-20 Paolo Carlini * include/tr1/random (gamma_distribution<>::_M_initialize, gamma_distribution<>::_M_l_d): Add. (gamma_distribution<>::gamma_distribution(const result_type&), operator>>(std::basic_istream<>&, gamma_distribution&)): Use it. include/tr1/random.tcc (gamma_distribution<>::_M_initialize): Define. (gamma_distribution<>::operator()): Adjust. * include/tr1/random (geometric_distribution<>::_M_initialize): Add. (geometric_distribution<>::geometric_distribution(const _RealType&), operator>>(std::basic_istream<>&, geometric_distribution&)): Use it. From-SVN: r116273 --- libstdc++-v3/ChangeLog | 14 +++++++++ libstdc++-v3/include/tr1/random | 25 ++++++++++++--- libstdc++-v3/include/tr1/random.tcc | 48 ++++++++++++++++++----------- 3 files changed, 64 insertions(+), 23 deletions(-) diff --git a/libstdc++-v3/ChangeLog b/libstdc++-v3/ChangeLog index 04e5977497a6..6bf7bd738845 100644 --- a/libstdc++-v3/ChangeLog +++ b/libstdc++-v3/ChangeLog @@ -1,3 +1,17 @@ +2006-08-20 Paolo Carlini + + * include/tr1/random (gamma_distribution<>::_M_initialize, + gamma_distribution<>::_M_l_d): Add. + (gamma_distribution<>::gamma_distribution(const result_type&), + operator>>(std::basic_istream<>&, gamma_distribution&)): Use it. + include/tr1/random.tcc (gamma_distribution<>::_M_initialize): + Define. + (gamma_distribution<>::operator()): Adjust. + + * include/tr1/random (geometric_distribution<>::_M_initialize): Add. + (geometric_distribution<>::geometric_distribution(const _RealType&), + operator>>(std::basic_istream<>&, geometric_distribution&)): Use it. + 2006-08-18 Paolo Carlini * include/tr1/random (class binomial_distribution<>): Add. diff --git a/libstdc++-v3/include/tr1/random b/libstdc++-v3/include/tr1/random index b769bbd6c747..5254ea6e8b33 100644 --- a/libstdc++-v3/include/tr1/random +++ b/libstdc++-v3/include/tr1/random @@ -1588,9 +1588,10 @@ _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(__p)) + : _M_p(__p) { _GLIBCXX_DEBUG_ASSERT((_M_p > 0.0) && (_M_p < 1.0)); + _M_initialize(); } /** @@ -1639,11 +1640,15 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) geometric_distribution& __x) { __is >> __x._M_p; - __x._M_log_p = std::log(__x._M_p); + __x._M_initialize(); return __is; } private: + void + _M_initialize() + { _M_log_p = std::log(_M_p); } + _RealType _M_p; _RealType _M_log_p; }; @@ -1746,8 +1751,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) _RealType _M_mean; - // _M_lm_thr hosts either log(mean) or the threshold of the simple - // method. + // 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_scx, _M_1cx, _M_c2b, _M_cb; @@ -2204,6 +2208,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) : _M_alpha(__alpha_val) { _GLIBCXX_DEBUG_ASSERT(_M_alpha > 0); + _M_initialize(); } /** @@ -2251,10 +2256,20 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) friend std::basic_istream<_CharT, _Traits>& operator>>(std::basic_istream<_CharT, _Traits>& __is, gamma_distribution& __x) - { return __is >> __x._M_alpha; } + { + __is >> __x._M_alpha; + __x._M_initialize(); + return __is; + } private: + void + _M_initialize(); + result_type _M_alpha; + + // Hosts either lambda of GB or d of modified Vaduva's. + result_type _M_l_d; }; /* @} */ // group tr1_random_distributions_continuous diff --git a/libstdc++-v3/include/tr1/random.tcc b/libstdc++-v3/include/tr1/random.tcc index 7719e8728f58..e809ba73ab6a 100644 --- a/libstdc++-v3/include/tr1/random.tcc +++ b/libstdc++-v3/include/tr1/random.tcc @@ -1106,11 +1106,10 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) /** - * Classic Box-Muller method. + * Polar method due to Marsaglia. * - * Reference: - * Box, G. E. P. and Muller, M. E. "A Note on the Generation of - * Random Normal Deviates." Ann. Math. Stat. 29, 610-611, 1958. + * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, + * New York, 1986, Ch. V, Sect. 4.4. */ template template @@ -1189,6 +1188,18 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) } + template + void + gamma_distribution<_RealType>:: + _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)); + } + /** * Cheng's rejection algorithm GB for alpha >= 1 and a modification * of Vaduva's rejection from Weibull algorithm due to Devroye for @@ -1213,19 +1224,18 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) { result_type __x; + bool __reject; if (_M_alpha >= 1) { // alpha - log(4) const result_type __b = _M_alpha - result_type(1.3862943611198906188344642429163531L); - const result_type __l = std::sqrt(2 * _M_alpha - 1); - const result_type __c = _M_alpha + __l; - const result_type __1l = 1 / __l; + const result_type __c = _M_alpha + _M_l_d; + const result_type __1l = 1 / _M_l_d; // 1 + log(9 / 2) const result_type __k = 2.5040773967762740733732583523868748L; - result_type __z, __r; do { const result_type __u = __urng(); @@ -1234,27 +1244,29 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) const result_type __y = __1l * std::log(__v / (1 - __v)); __x = _M_alpha * std::exp(__y); - __z = __u * __v * __v; - __r = __b + __c * __y - __x; + 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); } - while (__r < result_type(4.5) * __z - __k - && __r < std::log(__z)); + while (__reject); } else { const result_type __c = 1 / _M_alpha; - const result_type __d = - std::pow(_M_alpha, _M_alpha / (1 - _M_alpha)) * (1 - _M_alpha); - result_type __z, __e; do { - __z = -std::log(__urng()); - __e = -std::log(__urng()); + const result_type __z = -std::log(__urng()); + const result_type __e = -std::log(__urng()); __x = std::pow(__z, __c); + + __reject = __z + __e < _M_l_d + __x; } - while (__z + __e < __d + __x); + while (__reject); } return __x;