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

diff --git a/libstdc++-v3/ChangeLog b/libstdc++-v3/ChangeLog
index 04e5977..6bf7bd7 100644 (file)
--- 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 b769bbd..5254ea6 100644 (file)
--- 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 7719e87..e809ba7 100644 (file)
--- 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;