random (gamma_distribution<>::_M_initialize, [...]): Add.

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.

From-SVN: r116273

libstdc++-v3/ChangeLog

+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.

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

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;