boost/random/gamma_distribution.hpp

   1 /* boost random/gamma_distribution.hpp header file
   2  *
   3  * Copyright Jens Maurer 2002
   4  * Distributed under the Boost Software License, Version 1.0. (See
   5  * accompanying file LICENSE_1_0.txt or copy at
   6  * http://www.boost.org/LICENSE_1_0.txt)
   7  *
   8  * See http://www.boost.org for most recent version including documentation.
   9  *
  10  * $Id$
  11  *
  12  */
  13
  14 #ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
  15 #define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
  16
  17 #include <boost/config/no_tr1/cmath.hpp>
  18 #include <cassert>
  19 #include <boost/limits.hpp>
  20 #include <boost/static_assert.hpp>
  21 #include <boost/random/detail/config.hpp>
  22 #include <boost/random/exponential_distribution.hpp>
  23
  24 namespace boost {
  25
  26 // Knuth
  27 template<class RealType = double>
  28 class gamma_distribution
  29 {
  30 public:
  31   typedef RealType input_type;
  32   typedef RealType result_type;
  33
  34 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
  35   BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
  36 #endif
  37
  38   explicit gamma_distribution(const result_type& alpha_arg = result_type(1))
  39     : _exp(result_type(1)), _alpha(alpha_arg)
  40   {
  41     assert(_alpha > result_type(0));
  42     init();
  43   }
  44
  45   // compiler-generated copy ctor and assignment operator are fine
  46
  47   RealType alpha() const { return _alpha; }
  48
  49   void reset() { _exp.reset(); }
  50
  51   template<class Engine>
  52   result_type operator()(Engine& eng)
  53   {
  54 #ifndef BOOST_NO_STDC_NAMESPACE
  55     // allow for Koenig lookup
  56     using std::tan; using std::sqrt; using std::exp; using std::log;
  57     using std::pow;
  58 #endif
  59     if(_alpha == result_type(1)) {
  60       return _exp(eng);
  61     } else if(_alpha > result_type(1)) {
  62       // Can we have a boost::mathconst please?
  63       const result_type pi = result_type(3.14159265358979323846);
  64       for(;;) {
  65         result_type y = tan(pi * eng());
  66         result_type x = sqrt(result_type(2)*_alpha-result_type(1))*y
  67           + _alpha-result_type(1);
  68         if(x <= result_type(0))
  69           continue;
  70         if(eng() >
  71            (result_type(1)+y*y) * exp((_alpha-result_type(1))
  72                                         *log(x/(_alpha-result_type(1)))
  73                                         - sqrt(result_type(2)*_alpha
  74                                                -result_type(1))*y))
  75           continue;
  76         return x;
  77       }
  78     } else /* alpha < 1.0 */ {
  79       for(;;) {
  80         result_type u = eng();
  81         result_type y = _exp(eng);
  82         result_type x, q;
  83         if(u < _p) {
  84           x = exp(-y/_alpha);
  85           q = _p*exp(-x);
  86         } else {
  87           x = result_type(1)+y;
  88           q = _p + (result_type(1)-_p) * pow(x, _alpha-result_type(1));
  89         }
  90         if(u >= q)
  91           continue;
  92         return x;
  93       }
  94     }
  95   }
  96
  97 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
  98   template<class CharT, class Traits>
  99   friend std::basic_ostream<CharT,Traits>&
 100   operator<<(std::basic_ostream<CharT,Traits>& os, const gamma_distribution& gd)
 101   {
 102     os << gd._alpha;
 103     return os;
 104   }
 105
 106   template<class CharT, class Traits>
 107   friend std::basic_istream<CharT,Traits>&
 108   operator>>(std::basic_istream<CharT,Traits>& is, gamma_distribution& gd)
 109   {
 110     is >> std::ws >> gd._alpha;
 111     gd.init();
 112     return is;
 113   }
 114 #endif
 115
 116 private:
 117   void init()
 118   {
 119 #ifndef BOOST_NO_STDC_NAMESPACE
 120     // allow for Koenig lookup
 121     using std::exp;
 122 #endif
 123     _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));
 124   }
 125
 126   exponential_distribution<RealType> _exp;
 127   result_type _alpha;
 128   // some data precomputed from the parameters
 129   result_type _p;
 130 };
 131
 132 } // namespace boost
 133
 134 #endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP