/* boost random/gamma_distribution.hpp header file

  *

  * Copyright Jens Maurer 2002

  * Distributed under the Boost Software License, Version 1.0. (See

  * accompanying file LICENSE_1_0.txt or copy at

  * http://www.boost.org/LICENSE_1_0.txt)

  *

  * See http://www.boost.org for most recent version including documentation.

  *

  * $Id: gamma_distribution.hpp,v 1.9 2004/07/27 03:43:32 dgregor Exp $

  *

  */

 #ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP

 #define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP

 #include <cmath>

 #include <cassert>

 #include <boost/limits.hpp>

 #include <boost/static_assert.hpp>

 #include <boost/random/exponential_distribution.hpp>

 namespace boost {

 // Knuth

 template<class RealType = double>

 class gamma_distribution

 {

 public:

   typedef RealType input_type;

   typedef RealType result_type;

 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS

   BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);

 #endif

   explicit gamma_distribution(const result_type& alpha = result_type(1))

     : _exp(result_type(1)), _alpha(alpha)

   {

     assert(alpha > result_type(0));

     init();

   }

   // compiler-generated copy ctor and assignment operator are fine

   RealType alpha() const { return _alpha; }

   void reset() { _exp.reset(); }

   template<class Engine>

   result_type operator()(Engine& eng)

   {

 #ifndef BOOST_NO_STDC_NAMESPACE

     // allow for Koenig lookup

     using std::tan; using std::sqrt; using std::exp; using std::log;

     using std::pow;

 #endif

     if(_alpha == result_type(1)) {

       return _exp(eng);

     } else if(_alpha > result_type(1)) {

       // Can we have a boost::mathconst please?

       const result_type pi = result_type(3.14159265358979323846);

       for(;;) {

         result_type y = tan(pi * eng());

         result_type x = sqrt(result_type(2)*_alpha-result_type(1))*y

           + _alpha-result_type(1);

         if(x <= result_type(0))

           continue;

         if(eng() >

            (result_type(1)+y*y) * exp((_alpha-result_type(1))

                                         *log(x/(_alpha-result_type(1)))

                                         - sqrt(result_type(2)*_alpha

                                                -result_type(1))*y))

           continue;

         return x;

       }

     } else /* alpha < 1.0 */ {

       for(;;) {

         result_type u = eng();

         result_type y = _exp(eng);

         result_type x, q;

         if(u < _p) {

           x = exp(-y/_alpha);

           q = _p*exp(-x);

         } else {

           x = result_type(1)+y;

           q = _p + (result_type(1)-_p) * pow(x, _alpha-result_type(1));

         }

         if(u >= q)

           continue;

         return x;

       }

     }

   }

 #if !defined(BOOST_NO_OPERATORS_IN_NAMESPACE) && !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS)

   template<class CharT, class Traits>

   friend std::basic_ostream<CharT,Traits>&

   operator<<(std::basic_ostream<CharT,Traits>& os, const gamma_distribution& gd)

   {

     os << gd._alpha;

     return os;

   }

   template<class CharT, class Traits>

   friend std::basic_istream<CharT,Traits>&

   operator>>(std::basic_istream<CharT,Traits>& is, gamma_distribution& gd)

   {

     is >> std::ws >> gd._alpha;

     gd.init();

     return is;

   }

 #endif

 private:

   void init()

   {

 #ifndef BOOST_NO_STDC_NAMESPACE

     // allow for Koenig lookup

     using std::exp;

 #endif

     _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));

   }

   exponential_distribution<RealType> _exp;

   result_type _alpha;

   // some data precomputed from the parameters

   result_type _p;

 };

 } // namespace boost

 #endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP