diff -r 000000000000 -r e4d67989cc36 ossrv_pub/boost_apis/boost/random/gamma_distribution.hpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ossrv_pub/boost_apis/boost/random/gamma_distribution.hpp Tue Feb 02 02:01:42 2010 +0200 @@ -0,0 +1,133 @@ +/* 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 +#include +#include +#include +#include + +namespace boost { + +// Knuth +template +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::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 + 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 + friend std::basic_ostream& + operator<<(std::basic_ostream& os, const gamma_distribution& gd) + { + os << gd._alpha; + return os; + } + + template + friend std::basic_istream& + operator>>(std::basic_istream& 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 _exp; + result_type _alpha; + // some data precomputed from the parameters + result_type _p; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP