[rand.dist.pois.gamma]

git-svn-id: https://llvm.org/svn/llvm-project/libcxx/trunk@103788 91177308-0d34-0410-b5e6-96231b3b80d8
This commit is contained in:
Howard Hinnant
2010-05-14 18:43:10 +00:00
parent d0e811a37d
commit f417abe683
19 changed files with 755 additions and 8 deletions

View File

@@ -3348,21 +3348,22 @@ template<class _URNG>
_RealType
gamma_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p)
{
result_type __a = __p_.alpha();
result_type __a = __p.alpha();
uniform_real_distribution<result_type> __gen(0, 1);
exponential_distribution<result_type> __egen;
result_type __x;
if (__a == 1)
return exponential_distribution<result_type>(1/__p_.beta())(__g);
__x = __egen(__g);
else if (__a > 1)
{
const result_type __b = __a - 1;
const result_type __c = 3 * __a - result_type(0.75);
uniform_real_distribution<result_type> __gen(0, 1);
result_type __x;
while (true)
{
const result_type __u = __gen(__g);
const result_type __v = __gen(__g);
const result_type __w = __u * (1 - __u);
if (__w =! 0)
if (__w != 0)
{
const result_type __y = _STD::sqrt(__c / __w) *
(__u - result_type(0.5));
@@ -3377,10 +3378,29 @@ gamma_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p)
}
}
}
return __x * __p_.beta();
}
// else __a < 1
return 0; // temp!!!
else // __a < 1
{
while (true)
{
const result_type __u = __gen(__g);
const result_type __es = __egen(__g);
if (__u <= 1 - __a)
{
__x = _STD::pow(__u, 1 / __a);
if (__x <= __es)
break;
}
else
{
const result_type __e = -_STD::log((1-__u)/__a);
__x = _STD::pow(1 - __a + __a * __e, 1 / __a);
if (__x <= __e + __es)
break;
}
}
}
return __x * __p.beta();
}
template <class _CharT, class _Traits, class _RT>