* @brief Gets the next random number in the sequence.
*/
result_type
- operator()();
+ operator()()
+ {
+ _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
+ return _M_x;
+ }
/**
* @brief Compares two linear congruential random number generator
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
- {
- typedef typename _UniformRandomNumberGenerator::result_type
- _UResult_type;
- return _M_call(__urng, this->a(), this->b(),
- typename is_integral<_UResult_type>::type());
- }
+ { return this->operator()(__urng, this->param()); }
/**
* Gets a uniform random number in the range @f$[0, n)@f$.
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
- {
- __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
- __aurng(__urng);
- return (__aurng() * (this->b() - this->a())) + this->a();
- }
+ { return this->operator()(__urng, this->param()); }
template<typename _UniformRandomNumberGenerator>
result_type
explicit
lognormal_distribution(_RealType __m = _RealType(0),
_RealType __s = _RealType(1))
- : _M_param(__m, __s)
+ : _M_param(__m, __s), _M_nd()
{ }
explicit
lognormal_distribution(const param_type& __p)
- : _M_param(__p)
+ : _M_param(__p), _M_nd()
{ }
/**
*/
void
reset()
- { }
+ { _M_nd.reset(); }
/**
*
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng,
- const param_type& __p);
+ const param_type& __p)
+ { return std::exp(__p.s() * _M_nd(__urng) + __p.m()); }
private:
param_type _M_param;
+
+ normal_distribution<result_type> _M_nd;
};
/**
explicit
student_t_distribution(_RealType __n = _RealType(1))
- : _M_param(__n)
+ : _M_param(__n), _M_nd()
{ }
explicit
student_t_distribution(const param_type& __p)
- : _M_param(__p)
+ : _M_param(__p), _M_nd()
{ }
/**
*/
void
reset()
- { }
+ { _M_nd.reset(); }
/**
*
const param_type& __p);
private:
- template<typename _UniformRandomNumberGenerator>
- result_type
- _M_gaussian(_UniformRandomNumberGenerator& __urng,
- const result_type __sigma);
-
param_type _M_param;
+
+ normal_distribution<result_type> _M_nd;
};
/**
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
- {
- __detail::_Adaptor<_UniformRandomNumberGenerator, double>
- __aurng(__urng);
- if ((__aurng() - __aurng.min())
- < this->p() * (__aurng.max() - __aurng.min()))
- return true;
- return false;
- }
+ { return this->operator()(__urng, this->param()); }
template<typename _UniformRandomNumberGenerator>
result_type
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng)
- {
- __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
- __aurng(__urng);
- return -std::log(__aurng()) / this->lambda();
- }
+ { return this->operator()(__urng, this->param()); }
template<typename _UniformRandomNumberGenerator>
result_type
_RealType _M_alpha;
_RealType _M_beta;
- // Hosts either lambda of GB or d of modified Vaduva's.
- _RealType _M_l_d;
+ _RealType _M_malpha, _M_a2;
};
public:
explicit
gamma_distribution(_RealType __alpha_val = _RealType(1),
_RealType __beta_val = _RealType(1))
- : _M_param(__alpha_val, __beta_val)
+ : _M_param(__alpha_val, __beta_val), _M_nd()
{ }
explicit
gamma_distribution(const param_type& __p)
- : _M_param(__p)
+ : _M_param(__p), _M_nd()
{ }
/**
* @brief Resets the distribution state.
- *
- * Does nothing for the gamma distribution.
*/
void
- reset() { }
+ reset()
+ { _M_nd.reset(); }
/**
* @brief Returns the @f$ \alpha @f$ of the distribution.
private:
param_type _M_param;
+
+ normal_distribution<result_type> _M_nd;
};
/**
template<typename _UniformRandomNumberGenerator>
result_type
operator()(_UniformRandomNumberGenerator& __urng,
- const param_type& __p)
- {
- __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
- __aurng(__urng);
- return __p.b() * std::pow(-std::log(__aurng()),
- result_type(1) / __p.a());
- }
+ const param_type& __p);
private:
param_type _M_param;
typedef piecewise_constant_distribution<_RealType> distribution_type;
friend class piecewise_constant_distribution<_RealType>;
- param_type();
+ param_type()
+ : _M_int(), _M_den(), _M_cp()
+ { _M_initialize(); }
template<typename _InputIteratorB, typename _InputIteratorW>
param_type(_InputIteratorB __bfirst,
seed(__sum);
}
- /**
- * Gets the next generated value in sequence.
- */
- template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
- typename linear_congruential_engine<_UIntType, __a, __c, __m>::
- result_type
- linear_congruential_engine<_UIntType, __a, __c, __m>::
- operator()()
- {
- _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
- return _M_x;
- }
-
template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
typename _CharT, typename _Traits>
std::basic_ostream<_CharT, _Traits>&
{
const long double __r = static_cast<long double>(_M_b.max())
- static_cast<long double>(_M_b.min()) + 1.0L;
- const result_type __m = std::log10(__r) / std::log10(2.0L);
+ const result_type __m = std::log(__r) / std::log(2.0L);
result_type __n, __n0, __y0, __y1, __s0, __s1;
for (size_t __i = 0; __i < 2; ++__i)
{
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{
- typename gamma_distribution<>::param_type
- __gamma_param(__p.k(), 1.0);
- gamma_distribution<> __gamma(__gamma_param);
+ gamma_distribution<> __gamma(__p.k(), 1.0);
double __x = __gamma(__urng);
- typename poisson_distribution<result_type>::param_type
- __poisson_param(__x * __p.p() / (1.0 - __p.p()));
- poisson_distribution<result_type> __poisson(__poisson_param);
- result_type __m = __poisson(__urng);
-
- return __m;
+ poisson_distribution<result_type> __poisson(__x * __p.p()
+ / (1.0 - __p.p()));
+ return __poisson(__urng);
}
template<typename _IntType, typename _CharT, typename _Traits>
}
- template<typename _RealType>
- template<typename _UniformRandomNumberGenerator>
- typename lognormal_distribution<_RealType>::result_type
- lognormal_distribution<_RealType>::
- operator()(_UniformRandomNumberGenerator& __urng,
- const param_type& __p)
- {
- _RealType __u, __v, __r2, __normal;
- __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
- __aurng(__urng);
-
- do
- {
- // Choose x,y in uniform square (-1,-1) to (+1,+1).
- __u = 2 * __aurng() - 1;
- __v = 2 * __aurng() - 1;
-
- // See if it is in the unit circle.
- __r2 = __u * __u + __v * __v;
- }
- while (__r2 > 1 || __r2 == 0);
-
- __normal = __u * std::sqrt(-2 * std::log(__r2) / __r2);
-
- return std::exp(__p.s() * __normal + __p.m());
- }
-
template<typename _RealType, typename _CharT, typename _Traits>
std::basic_ostream<_CharT, _Traits>&
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
__os.fill(__space);
__os.precision(std::numeric_limits<_RealType>::digits10 + 1);
- __os << __x.m() << __space << __x.s();
+ __os << __x.m() << __space << __x.s()
+ << __space << __x._M_nd;
__os.flags(__flags);
__os.fill(__fill);
__is.flags(__ios_base::dec | __ios_base::skipws);
_RealType __m, __s;
- __is >> __m >> __s;
+ __is >> __m >> __s >> __x._M_nd;
__x.param(typename lognormal_distribution<_RealType>::
param_type(__m, __s));
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __p)
{
- typename gamma_distribution<_RealType>::param_type
- __gamma_param(__p.n() / 2, 1.0);
- gamma_distribution<_RealType> __gamma(__gamma_param);
+ gamma_distribution<_RealType> __gamma(__p.n() / 2, 1.0);
return 2 * __gamma(__urng);
}
}
- //
- // This could be operator() for a Gaussian distribution.
- //
- template<typename _RealType>
- template<typename _UniformRandomNumberGenerator>
- typename student_t_distribution<_RealType>::result_type
- student_t_distribution<_RealType>::
- _M_gaussian(_UniformRandomNumberGenerator& __urng,
- const result_type __sigma)
- {
- _RealType __x, __y, __r2;
- __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
- __aurng(__urng);
-
- do
- {
- // Choose x,y in uniform square (-1,-1) to (+1,+1).
- __x = 2 * __aurng() - 1;
- __y = 2 * __aurng() - 1;
-
- // See if it is in the unit circle.
- __r2 = __x * __x + __y * __y;
- }
- while (__r2 > 1 || __r2 == 0);
-
- // Box-Muller transform.
- return __sigma * __y * std::sqrt(-2 * std::log(__r2) / __r2);
- }
-
template<typename _RealType>
template<typename _UniformRandomNumberGenerator>
typename student_t_distribution<_RealType>::result_type
{
if (__param.n() <= 2.0)
{
- _RealType __y1 = _M_gaussian(__urng, 1.0);
- typename chi_squared_distribution<_RealType>::param_type
- __chisq_param(__param.n());
- chi_squared_distribution<_RealType> __chisq(__chisq_param);
+ _RealType __y1 = _M_nd(__urng);
+ chi_squared_distribution<_RealType> __chisq(__param.n());
_RealType __y2 = __chisq(__urng);
return __y1 / std::sqrt(__y2 / __param.n());
else
{
_RealType __y1, __y2, __z;
+ exponential_distribution<_RealType>
+ __exponential(1.0 / (__param.n() / 2.0 - 1.0));
+
do
{
- __y1 = _M_gaussian(__urng, 1.0);
- typename exponential_distribution<_RealType>::param_type
- __exp_param(1.0 / (__param.n() / 2.0 - 1.0));
- exponential_distribution<_RealType>
- __exponential(__exp_param);
+ __y1 = _M_nd(__urng);
__y2 = __exponential(__urng);
__z = __y1 * __y1 / (__param.n() - 2.0);
__os.fill(__space);
__os.precision(std::numeric_limits<_RealType>::digits10 + 1);
- __os << __x.n();
+ __os << __x.n() << __space << __x._M_nd;
__os.flags(__flags);
__os.fill(__fill);
__is.flags(__ios_base::dec | __ios_base::skipws);
_RealType __n;
- __is >> __n;
+ __is >> __n >> __x._M_nd;
__x.param(typename student_t_distribution<_RealType>::param_type(__n));
__is.flags(__flags);
gamma_distribution<_RealType>::param_type::
_M_initialize()
{
- if (_M_alpha >= 1)
- _M_l_d = std::sqrt(2 * _M_alpha - 1);
- else
- _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
- * (1 - _M_alpha));
+ _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
+
+ const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
+ _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
}
/**
- * 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, Ch. IX, Sect. 3.4 (+ Errata!).
+ * Marsaglia, G. and Tsang, W. W.
+ * "A Simple Method for Generating Gamma Variables"
+ * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
*/
template<typename _RealType>
template<typename _UniformRandomNumberGenerator>
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __param)
{
- result_type __x;
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
__aurng(__urng);
- bool __reject;
- const _RealType __alpha_val = __param.alpha();
- const _RealType __beta_val = __param.beta();
- if (__alpha_val >= 1)
- {
- // alpha - log(4)
- const result_type __b = __alpha_val
- - result_type(1.3862943611198906188344642429163531L);
- const result_type __c = __alpha_val + __param._M_l_d;
- const result_type __1l = 1 / __param._M_l_d;
-
- // 1 + log(9 / 2)
- const result_type __k = 2.5040773967762740733732583523868748L;
+ result_type __u, __v, __n;
+ const result_type __a1 = (__param._M_malpha
+ - _RealType(1.0) / _RealType(3.0));
+ do
+ {
do
{
- const result_type __u = __aurng() / __beta_val;
- const result_type __v = __aurng() / __beta_val;
-
- const result_type __y = __1l * std::log(__v / (1 - __v));
- __x = __alpha_val * std::exp(__y);
-
- const result_type __z = __u * __v * __v;
- const result_type __r = __b + __c * __y - __x;
-
- __reject = __r < result_type(4.5) * __z - __k;
- if (__reject)
- __reject = __r < std::log(__z);
+ __n = _M_nd(__urng);
+ __v = result_type(1.0) + __param._M_a2 * __n;
}
- while (__reject);
+ while (__v <= 0.0);
+
+ __v = __v * __v * __v;
+ __u = __aurng();
}
+ while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
+ && (std::log(__u) > (0.5 * __n * __n + __a1
+ * (1.0 - __v + std::log(__v)))));
+
+ if (__param.alpha() == __param._M_malpha)
+ return __a1 * __v * __param.beta();
else
{
- const result_type __c = 1 / __alpha_val;
-
do
- {
- const result_type __z = -std::log(__aurng() / __beta_val);
- const result_type __e = -std::log(__aurng() / __beta_val);
-
- __x = std::pow(__z, __c);
-
- __reject = __z + __e < __param._M_l_d + __x;
- }
- while (__reject);
+ __u = __aurng();
+ while (__u == 0.0);
+
+ return (std::pow(__u, result_type(1.0) / __param.alpha())
+ * __a1 * __v * __param.beta());
}
-
- return __beta_val * __x;
}
template<typename _RealType, typename _CharT, typename _Traits>
__os.fill(__space);
__os.precision(std::numeric_limits<_RealType>::digits10 + 1);
- __os << __x.alpha() << __space << __x.beta();
+ __os << __x.alpha() << __space << __x.beta()
+ << __space << __x._M_nd;
__os.flags(__flags);
__os.fill(__fill);
__is.flags(__ios_base::dec | __ios_base::skipws);
_RealType __alpha_val, __beta_val;
- __is >> __alpha_val >> __beta_val;
+ __is >> __alpha_val >> __beta_val >> __x._M_nd;
__x.param(typename gamma_distribution<_RealType>::
param_type(__alpha_val, __beta_val));
}
+ template<typename _RealType>
+ template<typename _UniformRandomNumberGenerator>
+ typename weibull_distribution<_RealType>::result_type
+ weibull_distribution<_RealType>::
+ operator()(_UniformRandomNumberGenerator& __urng,
+ const param_type& __p)
+ {
+ __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
+ __aurng(__urng);
+ return __p.b() * std::pow(-std::log(__aurng()),
+ result_type(1) / __p.a());
+ }
+
template<typename _RealType, typename _CharT, typename _Traits>
std::basic_ostream<_CharT, _Traits>&
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
}
template<typename _RealType>
- piecewise_constant_distribution<_RealType>::param_type::
- param_type()
- : _M_int(), _M_den(), _M_cp()
- { _M_initialize(); }
-
- template<typename _RealType>
template<typename _InputIteratorB, typename _InputIteratorW>
piecewise_constant_distribution<_RealType>::param_type::
param_type(_InputIteratorB __bbegin,
template<typename _RealType>
template<typename _Func>
piecewise_constant_distribution<_RealType>::param_type::
- param_type(size_t __nw, _RealType __xmin, _RealType __xmax,
- _Func __fw)
+ param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
: _M_int(), _M_den(), _M_cp()
{
for (size_t __i = 0; __i <= __nw; ++__i)
__bits);
const long double __r = static_cast<long double>(__urng.max())
- static_cast<long double>(__urng.min()) + 1.0L;
- const size_t __log2r = std::log10(__r) / std::log10(2.0L);
+ const size_t __log2r = std::log(__r) / std::log(2.0L);
size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
_RealType __sum = _RealType(0);
_RealType __tmp = _RealType(1);