From a2f032055c9de2863bd889496e6454fc5f52adfe Mon Sep 17 00:00:00 2001
From: paolo <paolo@138bc75d-0d04-0410-961f-82ee72b054a4>
Date: Sun, 20 Aug 2006 16:05:05 +0000
Subject: [PATCH] 2006-08-20 Paolo Carlini <pcarlini@suse.de>
* 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.
git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@116273 138bc75d-0d04-0410-961f-82ee72b054a4
---
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 <pcarlini@suse.de>
+
+ * 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 <pcarlini@suse.de>
* 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<typename _RealType>
template<class _UniformRandomNumberGenerator>
@@ -1189,6 +1188,18 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
}
+ template<typename _RealType>
+ 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;
--
GitLab