stc/exp/ml_stc_linux_make_v1.0/include/boost/random/gamma_distribution.hpp

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