random (class gamma_distribution<>): Add.

2006-07-06  Paolo Carlini  <pcarlini@suse.de>

	* include/tr1/random (class gamma_distribution<>): Add.
	* include/tr1/random.tcc (gamma_distribution<>::operator(),
	operator<<(std::basic_ostream<>&, const gamma_distribution<>&)):
	Define.
	* testsuite/tr1/5_numerical_facilities/random/gamma_distribution/
	requirements/typedefs.cc: New.

From-SVN: r115236

libstdc++-v3/ChangeLog

@@ -1,3 +1,12 @@
+2006-07-06  Paolo Carlini  <pcarlini@suse.de>
+
+	* include/tr1/random (class gamma_distribution<>): Add.
+	* include/tr1/random.tcc (gamma_distribution<>::operator(),
+	operator<<(std::basic_ostream<>&, const gamma_distribution<>&)):
+	Define.
+	* testsuite/tr1/5_numerical_facilities/random/gamma_distribution/
+	requirements/typedefs.cc: New.
+
 2006-07-06  Benjamin Kosnik  <bkoz@redhat.com>
 
 	* testsuite/util/regression/trait/assoc/trait.hpp: Format.

libstdc++-v3/include/tr1/random

@@ -1952,6 +1952,91 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
       bool        _M_saved_available;     
     };
 
+
+  /**
+   * @brief A gamma continuous distribution for random numbers.
+   *
+   * The formula for the gamma probability mass function is 
+   * @f$ p(x) = \frac{1}{\Gamma(\alpha)} x^{\alpha - 1} e^{-x} } @f$.
+   */
+  template<typename _RealType = double>
+    class gamma_distribution;
+
+  template<typename _RealType, typename _CharT, typename _Traits>
+    std::basic_ostream<_CharT, _Traits>&
+    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
+	       const gamma_distribution<_RealType>& __x);
+
+  template<typename _RealType>
+    class gamma_distribution
+    {
+    public:
+      // types
+      typedef _RealType input_type;
+      typedef _RealType result_type;
+
+    public:
+      /**
+       * Constructs a gamma distribution with parameters @f$ \alpha @f$.
+       */
+      explicit
+      gamma_distribution(const result_type& __alpha = result_type(1))
+      : _M_alpha(__alpha)
+      { 
+	_GLIBCXX_DEBUG_ASSERT(_M_alpha > 0);
+      }
+
+      /**
+       * Gets the @f$ \alpha @f$ of the distribution.
+       */
+      _RealType
+      alpha() const
+      { return _M_alpha; }
+
+      /**
+       * Resets the distribution.
+       */
+      void
+      reset() { }
+
+      template<class _UniformRandomNumberGenerator>
+        result_type
+        operator()(_UniformRandomNumberGenerator& __urng);
+
+      /**
+       * Inserts a %gamma_distribution random number distribution
+       * @p __x into the output stream @p __os.
+       *
+       * @param __os An output stream.
+       * @param __x  A %gamma_distribution random number distribution.
+       *
+       * @returns The output stream with the state of @p __x inserted or in
+       * an error state.
+       */
+      template<typename _RealType1, typename _CharT, typename _Traits>
+        friend std::basic_ostream<_CharT, _Traits>&
+        operator<<(std::basic_ostream<_CharT, _Traits>& __os,
+		   const gamma_distribution<_RealType1>& __x);
+
+      /**
+       * Extracts a %gamma_distribution random number distribution
+       * @p __x from the input stream @p __is.
+       *
+       * @param __is An input stream.
+       * @param __x  A %gamma_distribution random number generator engine.
+       *
+       * @returns The input stream with @p __x extracted or in an error state.
+       */
+      template<typename _RealType1, typename _CharT, typename _Traits>
+        friend std::basic_istream<_CharT, _Traits>&
+        operator>>(std::basic_istream<_CharT, _Traits>& __is,
+		   gamma_distribution<_RealType1>& __x)
+        { return __is >> __x._M_alpha; }
+
+    private:
+      result_type _M_alpha;
+    };
+
   /* @} */ // group tr1_random_distributions_continuous
   /* @} */ // group tr1_random_distributions
   /* @} */ // group tr1_random

libstdc++-v3/include/tr1/random.tcc

@@ -767,5 +767,95 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
       return __is;
     }
 
+
+  /**
+   * Cheng's rejection algorithm GB for alpha >= 1 and a modification
+   * of Vaduva's rejection from Weibull algorithm due to Devroye for
+   * alpha < 1.
+   *
+   * References:
+   * Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral
+   * Shape Parameter." Applied Statistics, 26, 71-75, 1977.
+   *
+   * Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection
+   * and Composition Procedures." Math. Operationsforschung and Statistik,
+   * Series in Statistics, 8, 545-576, 1977.
+   *
+   * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
+   * New York, 1986, Sect. 3.4.
+   */
+  template<typename _RealType>
+    template<class _UniformRandomNumberGenerator>
+      typename gamma_distribution<_RealType>::result_type
+      gamma_distribution<_RealType>::
+      operator()(_UniformRandomNumberGenerator& __urng)
+      {
+	result_type __x;
+
+	if (_M_alpha >= 1)
+	  {
+	    // 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);
+
+	    // 1 + log(9 / 2)
+	    const result_type __k = 2.5040773967762740733732583523868748L;
+
+	    result_type __z, __r;
+	    do
+	      {
+		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);
+
+		__z = __u * __v * __v;
+		__r = __b + __c * __y - __x;
+	      }
+	    while (__r < result_type(4.5) * __z - __k
+		   && __r < std::log(__z));
+	  }
+	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());
+
+		__x = std::pow(__z, __c);
+	      }
+	    while (__z + __e > __d + __x);
+	  }
+
+	return __x;
+      }
+
+  template<typename _RealType, typename _CharT, typename _Traits>
+    std::basic_ostream<_CharT, _Traits>&
+    operator<<(std::basic_ostream<_CharT, _Traits>& __os,
+	       const gamma_distribution<_RealType>& __x)
+    {
+      const std::ios_base::fmtflags __flags = __os.flags();
+      const _CharT __fill = __os.fill();
+      const std::streamsize __precision = __os.precision();
+      __os.flags(std::ios_base::scientific | std::ios_base::left);
+      __os.fill(__os.widen(' '));
+      __os.precision(_Private::_Max_digits10<_RealType>::__value);
+
+      __os << __x.alpha();
+
+      __os.flags(__flags);
+      __os.fill(__fill);
+      __os.precision(__precision);
+      return __os;
+    }
+
 _GLIBCXX_END_NAMESPACE
 }

libstdc++-v3/testsuite/tr1/5_numerical_facilities/random/gamma_distribution/requirements/typedefs.cc

@@ -0,0 +1,37 @@
+// { dg-do compile }
+//
+// 2006-07-06  Paolo Carlini  <pcarlini@suse.de>
+//
+// Copyright (C) 2006 Free Software Foundation, Inc.
+//
+// This file is part of the GNU ISO C++ Library.  This library is free
+// software; you can redistribute it and/or modify it under the
+// terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with this library; see the file COPYING.  If not, write to the Free
+// Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
+// USA.
+
+// 5.1.7.9 Class template gamma_distribution [tr.rand.dist.gamma]
+// 5.1.1 [7] Table 17
+
+#include <tr1/random>
+
+void
+test01() 
+{ 
+  using namespace std::tr1;
+
+  typedef gamma_distribution<double> test_type;
+
+  typedef test_type::input_type  input_type;
+  typedef test_type::result_type result_type;
+}