random.tcc (gamma_distribution<>::operator()): Fixes from the Errata of Devroye's...

2006-08-10  Paolo Carlini  <pcarlini@suse.de>

	* include/tr1/random.tcc (gamma_distribution<>::operator()): Fixes
	from the Errata of Devroye's book.

From-SVN: r116061

diff --git a/libstdc++-v3/ChangeLog b/libstdc++-v3/ChangeLog
index 2753d050afce520bdc9eee8623b896e791adf253..866147f8e367bd0644db1f8b17d63ea26cea360d 100644 (file)
--- a/libstdc++-v3/ChangeLog
+++ b/libstdc++-v3/ChangeLog
@@ -1,3 +1,8 @@
+2006-08-10  Paolo Carlini  <pcarlini@suse.de>
+
+       * include/tr1/random.tcc (gamma_distribution<>::operator()): Fixes
+       from the Errata of Devroye's book.
+
 2006-08-10  Paolo Carlini  <pcarlini@suse.de>
 
        * include/bits/stl_bvector.h (_Bit_iterator_base::_M_incr(ptrdiff_t)):

diff --git a/libstdc++-v3/include/tr1/random.tcc b/libstdc++-v3/include/tr1/random.tcc
index 5ce415b1829c6cd58ca4a0c7ac1f51788b0b836b..17732e3d68443a277d303a61df02a62366c26c17 100644 (file)
--- a/libstdc++-v3/include/tr1/random.tcc
+++ b/libstdc++-v3/include/tr1/random.tcc
@@ -807,7 +807,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
    * Series in Statistics, 8, 545-576, 1977.
    *
    * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
-   * New York, 1986, Sect. 3.4.
+   * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
    */
   template<typename _RealType>
     template<class _UniformRandomNumberGenerator>
@@ -822,7 +822,9 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
            // alpha - log(4)
            const result_type __b = _M_alpha
              - result_type(1.3862943611198906188344642429163531L);
-           const result_type __c = _M_alpha + std::sqrt(2 * _M_alpha - 1);
+           const result_type __l = std::sqrt(2 * _M_alpha - 1);
+           const result_type __c = _M_alpha + __l;
+           const result_type __1l = 1 / __l;
 
            // 1 + log(9 / 2)
            const result_type __k = 2.5040773967762740733732583523868748L;
@@ -833,8 +835,8 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
                const result_type __u = __urng();
                const result_type __v = __urng();
 
-               const result_type __y = _M_alpha * std::log(__v / (1 - __v));
-               __x = _M_alpha * std::exp(__v);
+               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;
@@ -856,7 +858,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
 
                __x = std::pow(__z, __c);
              }
-           while (__z + __e > __d + __x);
+           while (__z + __e < __d + __x);
          }
 
        return __x;