diff options
| author | Howard Hinnant <hhinnant@apple.com> | 2010-05-14 22:43:10 +0400 |
|---|---|---|
| committer | Howard Hinnant <hhinnant@apple.com> | 2010-05-14 22:43:10 +0400 |
| commit | 7070922ff8954467121a62bd2c2a6e6e0a0bfa86 (patch) | |
| tree | c535413976e4a50b27cd48fdf33d1b8182d60bba /libcxx/include/random | |
| parent | bdb1b0d6cb9688b70491265566f76c1e401fa376 (diff) | |
[rand.dist.pois.gamma]
llvm-svn: 103788
Diffstat (limited to 'libcxx/include/random')
| -rw-r--r-- | libcxx/include/random | 36 |
1 files changed, 28 insertions, 8 deletions
diff --git a/libcxx/include/random b/libcxx/include/random index 30a9c8ca376b..948bb9aa9b37 100644 --- a/libcxx/include/random +++ b/libcxx/include/random @@ -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> |