[rand.dist.pois.gamma]
git-svn-id: https://llvm.org/svn/llvm-project/libcxx/trunk@103788 91177308-0d34-0410-b5e6-96231b3b80d8
diff --git a/include/random b/include/random
index 30a9c8c..948bb9a 100644
--- a/include/random
+++ b/include/random
@@ -3348,21 +3348,22 @@
_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 @@
}
}
}
- 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>