From 42031254bfae5d55420bf8c5d57572c42d1b5d17 Mon Sep 17 00:00:00 2001 From: Paolo Carlini Date: Fri, 18 Aug 2006 17:15:43 +0000 Subject: [PATCH] random (class binomial_distribution<>): Add. 2006-08-18 Paolo Carlini * include/tr1/random (class binomial_distribution<>): Add. * include/tr1/random.tcc (binomial_distribution<>::operator(), operator<<(std::basic_ostream<>&, const binomial_distribution<>&), operator>>(std::basic_istream<>&, binomial_distribution<>&, binomial_distribution<>::_M_waiting(), binomial_distribution<>:: _M_initialize()): Define. * testsuite/tr1/5_numerical_facilities/random/binomial_distribution/ requirements/typedefs.cc: New. * include/tr1/random (geometric_distribution<>:: geometric_distribution(const _RealType&)): Fix DEBUG_ASSERT limits. * include/tr1/random (poisson_distribution): Add normal_distribution member, adjust consistently; minor tweaks and rearrangements of the arithmetic. (operator>>(std::basic_istream<>&, poisson_distribution<>&)): Move out of line. * include/tr1/random.tcc: Adjust. * include/tr1/random.tcc (normal_distribution<>::operator()): Minor tweaks. From-SVN: r116245 --- libstdc++-v3/ChangeLog | 25 ++ libstdc++-v3/include/tr1/random | 150 ++++++++- libstdc++-v3/include/tr1/random.tcc | 309 ++++++++++++++++-- .../requirements/typedefs.cc | 37 +++ 4 files changed, 473 insertions(+), 48 deletions(-) create mode 100644 libstdc++-v3/testsuite/tr1/5_numerical_facilities/random/binomial_distribution/requirements/typedefs.cc diff --git a/libstdc++-v3/ChangeLog b/libstdc++-v3/ChangeLog index 9d692ab132a..04e5977497a 100644 --- a/libstdc++-v3/ChangeLog +++ b/libstdc++-v3/ChangeLog @@ -1,3 +1,28 @@ +2006-08-18 Paolo Carlini + + * include/tr1/random (class binomial_distribution<>): Add. + * include/tr1/random.tcc (binomial_distribution<>::operator(), + operator<<(std::basic_ostream<>&, const binomial_distribution<>&), + operator>>(std::basic_istream<>&, binomial_distribution<>&, + binomial_distribution<>::_M_waiting(), binomial_distribution<>:: + _M_initialize()): Define. + * testsuite/tr1/5_numerical_facilities/random/binomial_distribution/ + requirements/typedefs.cc: New. + + * include/tr1/random (geometric_distribution<>:: + geometric_distribution(const _RealType&)): Fix DEBUG_ASSERT + limits. + + * include/tr1/random (poisson_distribution): Add normal_distribution + member, adjust consistently; minor tweaks and rearrangements of the + arithmetic. + (operator>>(std::basic_istream<>&, poisson_distribution<>&)): Move + out of line. + * include/tr1/random.tcc: Adjust. + + * include/tr1/random.tcc (normal_distribution<>::operator()): Minor + tweaks. + 2006-08-18 Paolo Carlini PR libstdc++/28765 diff --git a/libstdc++-v3/include/tr1/random b/libstdc++-v3/include/tr1/random index ef4390a0610..b769bbd6c74 100644 --- a/libstdc++-v3/include/tr1/random +++ b/libstdc++-v3/include/tr1/random @@ -1588,9 +1588,9 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) // constructors and member function explicit geometric_distribution(const _RealType& __p = _RealType(0.5)) - : _M_p(__p), _M_log_p(std::log(_M_p)) + : _M_p(__p), _M_log_p(std::log(__p)) { - _GLIBCXX_DEBUG_ASSERT((_M_p >= 0.0) && (_M_p <= 1.0)); + _GLIBCXX_DEBUG_ASSERT((_M_p > 0.0) && (_M_p < 1.0)); } /** @@ -1649,6 +1649,9 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) }; + template + class normal_distribution; + /** * @brief A discrete Poisson random number distribution. * @@ -1665,6 +1668,12 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) operator<<(std::basic_ostream<_CharT, _Traits>& __os, const poisson_distribution<_IntType, _RealType>& __x); + template + std::basic_istream<_CharT, _Traits>& + operator>>(std::basic_istream<_CharT, _Traits>& __is, + poisson_distribution<_IntType, _RealType>& __x); + template class poisson_distribution { @@ -1676,7 +1685,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) // constructors and member function explicit poisson_distribution(const _RealType& __mean = _RealType(1)) - : _M_mean(__mean) + : _M_mean(__mean), _M_nd() { _GLIBCXX_DEBUG_ASSERT(_M_mean > 0.0); _M_initialize(); @@ -1690,7 +1699,8 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) { return _M_mean; } void - reset() { } + reset() + { _M_nd.reset(); } template result_type @@ -1721,29 +1731,145 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) * * @returns The input stream with @p __x extracted or in an error state. */ - template + template friend std::basic_istream<_CharT, _Traits>& operator>>(std::basic_istream<_CharT, _Traits>& __is, - poisson_distribution& __x) - { - __is >> __x._M_mean; - __x._M_initialize(); - return __is; - } + poisson_distribution<_IntType1, _RealType1>& __x); private: void _M_initialize(); + // NB: Unused when _GLIBCXX_USE_C99_MATH_TR1 is undefined. + normal_distribution<_RealType> _M_nd; + _RealType _M_mean; + // _M_lm_thr hosts either log(mean) or the threshold of the simple // method. _RealType _M_lm_thr; #if _GLIBCXX_USE_C99_MATH_TR1 - _RealType _M_lfm, _M_sm, _M_d, _M_scx4, _M_2cx, _M_c2b, _M_cb; + _RealType _M_lfm, _M_sm, _M_d, _M_scx, _M_1cx, _M_c2b, _M_cb; #endif }; + + /** + * @brief A discrete binomial random number distribution. + * + * The formula for the binomial probability mass function is + * @f$ p(i) = \binom{n}{i} p^i (1 - p)^{t - i} @f$ where @f$ t @f$ + * and @f$ p @f$ are the parameters of the distribution. + */ + template + class binomial_distribution; + + template + std::basic_ostream<_CharT, _Traits>& + operator<<(std::basic_ostream<_CharT, _Traits>& __os, + const binomial_distribution<_IntType, _RealType>& __x); + + template + std::basic_istream<_CharT, _Traits>& + operator>>(std::basic_istream<_CharT, _Traits>& __is, + binomial_distribution<_IntType, _RealType>& __x); + + template + class binomial_distribution + { + public: + // types + typedef _RealType input_type; + typedef _IntType result_type; + + // constructors and member function + explicit + binomial_distribution(_IntType __t = 1, + const _RealType& __p = _RealType(0.5)) + : _M_t(__t), _M_p(__p), _M_nd() + { + _GLIBCXX_DEBUG_ASSERT((_M_t >= 0) && (_M_p >= 0.0) && (_M_p <= 1.0)); + _M_initialize(); + } + + /** + * Gets the distribution @p t parameter. + */ + _IntType + t() const + { return _M_t; } + + /** + * Gets the distribution @p p parameter. + */ + _RealType + p() const + { return _M_p; } + + void + reset() + { _M_nd.reset(); } + + template + result_type + operator()(_UniformRandomNumberGenerator& __urng); + + /** + * Inserts a %binomial_distribution random number distribution + * @p __x into the output stream @p __os. + * + * @param __os An output stream. + * @param __x A %binomial_distribution random number distribution. + * + * @returns The output stream with the state of @p __x inserted or in + * an error state. + */ + template + friend std::basic_ostream<_CharT, _Traits>& + operator<<(std::basic_ostream<_CharT, _Traits>& __os, + const binomial_distribution<_IntType1, _RealType1>& __x); + + /** + * Extracts a %binomial_distribution random number distribution + * @p __x from the input stream @p __is. + * + * @param __is An input stream. + * @param __x A %binomial_distribution random number generator engine. + * + * @returns The input stream with @p __x extracted or in an error state. + */ + template + friend std::basic_istream<_CharT, _Traits>& + operator>>(std::basic_istream<_CharT, _Traits>& __is, + binomial_distribution<_IntType1, _RealType1>& __x); + + private: + void + _M_initialize(); + + template + result_type + _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t); + + // NB: Unused when _GLIBCXX_USE_C99_MATH_TR1 is undefined. + normal_distribution<_RealType> _M_nd; + + _IntType _M_t; + _RealType _M_p; + + _RealType _M_q; +#if _GLIBCXX_USE_C99_MATH_TR1 + _RealType _M_d1, _M_d2, _M_s1, _M_s2, _M_c, + _M_a1, _M_a123, _M_s, _M_lf, _M_lp1p; +#endif + bool _M_easy; + }; + /* @} */ // group tr1_random_distributions_discrete /** diff --git a/libstdc++-v3/include/tr1/random.tcc b/libstdc++-v3/include/tr1/random.tcc index abc4d3a3df8..7719e8728f5 100644 --- a/libstdc++-v3/include/tr1/random.tcc +++ b/libstdc++-v3/include/tr1/random.tcc @@ -667,21 +667,18 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) _M_lm_thr = std::log(_M_mean); _M_lfm = std::tr1::lgamma(__m + 1); _M_sm = std::sqrt(__m); - - const _RealType __dx = - std::sqrt(2 * __m - * std::log(_RealType(40.743665431525205956834243423363677L) - * __m)); + + const _RealType __pi_4 = 0.7853981633974483096156608458198757L; + const _RealType __dx = std::sqrt(2 * __m * std::log(32 * __m + / __pi_4)); _M_d = std::tr1::round(std::max(_RealType(6), std::min(__m, __dx))); - const _RealType __cx = 2 * (2 * __m + _M_d); - const _RealType __cx4 = __cx / 4; - _M_scx4 = std::sqrt(__cx4); - _M_2cx = 2 / __cx; + const _RealType __cx = 2 * __m + _M_d; + _M_scx = std::sqrt(__cx / 2); + _M_1cx = 1 / __cx; - const _RealType __pi_2 = 1.5707963267948966192313216916397514L; - _M_c2b = std::sqrt(__pi_2 * __cx4) * std::exp(_M_2cx); - _M_cb = __cx * std::exp(-_M_d * _M_2cx * (1 + _M_d / 2)) / _M_d; + _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx); + _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) / _M_d; } else #endif @@ -720,11 +717,9 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) const _RealType __e178 = 1.0129030479320018583185514777512983L; const _RealType __c5 = __c4 + __e178; const _RealType __c = _M_cb + __c5; - const _RealType __cx = 2 * (2 * __m + _M_d); + const _RealType __2cx = 2 * (2 * __m + _M_d); - normal_distribution<_RealType> __nd; - - bool __keepgoing = true; + bool __reject = true; do { const _RealType __u = __c * __urng(); @@ -734,7 +729,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) if (__u <= __c1) { - const _RealType __n = __nd(__urng); + const _RealType __n = _M_nd(__urng); const _RealType __y = -std::abs(__n) * _M_sm - 1; __x = std::floor(__y); __w = -__n * __n / 2; @@ -743,10 +738,10 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) } else if (__u <= __c2) { - const _RealType __n = __nd(__urng); - const _RealType __y = 1 + std::abs(__n) * _M_scx4; + const _RealType __n = _M_nd(__urng); + const _RealType __y = 1 + std::abs(__n) * _M_scx; __x = std::ceil(__y); - __w = __y * (2 - __y) * _M_2cx; + __w = __y * (2 - __y) * _M_1cx; if (__x > _M_d) continue; } @@ -760,22 +755,22 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) else { const _RealType __v = -std::log(__urng()); - const _RealType __y = _M_d + __v * __cx / _M_d; + const _RealType __y = _M_d + __v * __2cx / _M_d; __x = std::ceil(__y); - __w = -_M_d * _M_2cx * (1 + __y / 2); + __w = -_M_d * _M_1cx * (1 + __y / 2); } - __keepgoing = (__w - __e - __x * _M_lm_thr - > _M_lfm - std::tr1::lgamma(__x + __m + 1)); + __reject = (__w - __e - __x * _M_lm_thr + > _M_lfm - std::tr1::lgamma(__x + __m + 1)); - } while (__keepgoing); + } while (__reject); - return _IntType(std::tr1::round(__x + __m)); + return _IntType(__x + __m + 0.5); } else #endif { - _IntType __x = -1; + _IntType __x = -1; _RealType __prod = 1.0; do @@ -798,11 +793,12 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) const std::ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); const std::streamsize __precision = __os.precision(); + const _CharT __space = __os.widen(' '); __os.flags(std::ios_base::scientific | std::ios_base::left); - __os.fill(__os.widen(' ')); + __os.fill(__space); __os.precision(_Max_digits10<_RealType>::__value); - __os << __x.mean(); + __os << __x.mean() << __space << __x._M_nd; __os.flags(__flags); __os.fill(__fill); @@ -810,6 +806,247 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) return __os; } + template + std::basic_istream<_CharT, _Traits>& + operator>>(std::basic_istream<_CharT, _Traits>& __is, + poisson_distribution<_IntType, _RealType>& __x) + { + const std::ios_base::fmtflags __flags = __is.flags(); + __is.flags(std::ios_base::skipws); + + __is >> __x._M_mean >> __x._M_nd; + __x._M_initialize(); + + __is.flags(__flags); + return __is; + } + + + template + void + binomial_distribution<_IntType, _RealType>:: + _M_initialize() + { + const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; + + _M_easy = true; + +#if _GLIBCXX_USE_C99_MATH_TR1 + if (_M_t * __p12 >= 8) + { + _M_easy = false; + const _RealType __np = std::floor(_M_t * __p12); + const _RealType __pa = __np / _M_t; + const _RealType __1p = 1 - __pa; + + const _RealType __pi_4 = 0.7853981633974483096156608458198757L; + const _RealType __d1x = + std::sqrt(__np * __1p * std::log(32 * __np + / (81 * __pi_4 * __1p))); + _M_d1 = std::tr1::round(std::max(_RealType(1), __d1x)); + const _RealType __d2x = + std::sqrt(__np * __1p * std::log(32 * _M_t * __1p + / (__pi_4 * __pa))); + _M_d2 = std::tr1::round(std::max(_RealType(1), __d2x)); + + // sqrt(pi / 2) + const _RealType __spi_2 = 1.2533141373155002512078826424055226L; + _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np)); + _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p)); + _M_c = 2 * _M_d1 / __np; + _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2; + const _RealType __a12 = _M_a1 + _M_s2 * __spi_2; + const _RealType __s1s = _M_s1 * _M_s1; + _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p)) + * 2 * __s1s / _M_d1 + * std::exp(-_M_d1 * _M_d1 / (2 * __s1s))); + const _RealType __s2s = _M_s2 * _M_s2; + _M_s = (_M_a123 + 2 * __s2s / _M_d2 + * std::exp(-_M_d2 * _M_d2 / (2 * __s2s))); + _M_lf = (std::tr1::lgamma(__np + 1) + + std::tr1::lgamma(_M_t - __np + 1)); + _M_lp1p = std::log(__pa / __1p); + + _M_q = -std::log(1 - (__p12 - __pa) / __1p); + } + else +#endif + _M_q = -std::log(1 - __p12); + } + + template + template + typename binomial_distribution<_IntType, _RealType>::result_type + binomial_distribution<_IntType, _RealType>:: + _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t) + { + _IntType __x = 0; + _RealType __sum = 0; + + do + { + const _RealType __e = -std::log(__urng()); + __sum += __e / (__t - __x); + __x += 1; + } + while (__sum <= _M_q); + + return __x - 1; + } + + /** + * A rejection algorithm when t * p >= 8 and a simple waiting time + * method - the second in the referenced book - otherwise. + * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 + * is defined. + * + * Reference: + * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, + * New York, 1986, Ch. X, Sect. 4 (+ Errata!). + */ + template + template + typename binomial_distribution<_IntType, _RealType>::result_type + binomial_distribution<_IntType, _RealType>:: + operator()(_UniformRandomNumberGenerator& __urng) + { + result_type __ret; + const _RealType __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; + +#if _GLIBCXX_USE_C99_MATH_TR1 + if (!_M_easy) + { + _RealType __x; + + const _RealType __np = std::floor(_M_t * __p12); + const _RealType __pa = __np / _M_t; + + // sqrt(pi / 2) + const _RealType __spi_2 = 1.2533141373155002512078826424055226L; + const _RealType __a1 = _M_a1; + const _RealType __a12 = __a1 + _M_s2 * __spi_2; + const _RealType __a123 = _M_a123; + const _RealType __s1s = _M_s1 * _M_s1; + const _RealType __s2s = _M_s2 * _M_s2; + + bool __reject; + do + { + const _RealType __u = _M_s * __urng(); + + _RealType __v; + + if (__u <= __a1) + { + const _RealType __n = _M_nd(__urng); + const _RealType __y = _M_s1 * std::abs(__n); + __reject = __y >= _M_d1; + if (!__reject) + { + const _RealType __e = -std::log(__urng()); + __x = std::floor(__y); + __v = -__e - __n * __n / 2 + _M_c; + } + } + else if (__u <= __a12) + { + const _RealType __n = _M_nd(__urng); + const _RealType __y = _M_s2 * std::abs(__n); + __reject = __y >= _M_d2; + if (!__reject) + { + const _RealType __e = -std::log(__urng()); + __x = std::floor(-__y); + __v = -__e - __n * __n / 2; + } + } + else if (__u <= __a123) + { + const _RealType __e1 = -std::log(__urng()); + const _RealType __e2 = -std::log(__urng()); + + const _RealType __y = _M_d1 + 2 * __s1s * __e1 / _M_d1; + __x = std::floor(__y); + __v = (-__e2 + _M_d1 * (1 / (_M_t - __np) + -__y / (2 * __s1s))); + __reject = false; + } + else + { + const _RealType __e1 = -std::log(__urng()); + const _RealType __e2 = -std::log(__urng()); + + const _RealType __y = _M_d2 + 2 * __s2s * __e1 / _M_d2; + __x = std::floor(-__y); + __v = -__e2 - _M_d2 * __y / (2 * __s2s); + __reject = false; + } + + __reject = __reject || __x < -__np || __x > _M_t - __np; + if (!__reject) + { + const _RealType __lfx = + std::tr1::lgamma(__np + __x + 1) + + std::tr1::lgamma(_M_t - (__np + __x) + 1); + __reject = __v > _M_lf - __lfx + __x * _M_lp1p; + } + } + while (__reject); + + __x += __np + 0.5; + + const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x)); + __ret = _IntType(__x) + __z; + } + else +#endif + __ret = _M_waiting(__urng, _M_t); + + if (__p12 != _M_p) + __ret = _M_t - __ret; + return __ret; + } + + template + std::basic_ostream<_CharT, _Traits>& + operator<<(std::basic_ostream<_CharT, _Traits>& __os, + const binomial_distribution<_IntType, _RealType>& __x) + { + const std::ios_base::fmtflags __flags = __os.flags(); + const _CharT __fill = __os.fill(); + const std::streamsize __precision = __os.precision(); + const _CharT __space = __os.widen(' '); + __os.flags(std::ios_base::scientific | std::ios_base::left); + __os.fill(__space); + __os.precision(_Max_digits10<_RealType>::__value); + + __os << __x.t() << __space << __x.p() + << __space << __x._M_nd; + + __os.flags(__flags); + __os.fill(__fill); + __os.precision(__precision); + return __os; + } + + template + std::basic_istream<_CharT, _Traits>& + operator>>(std::basic_istream<_CharT, _Traits>& __is, + binomial_distribution<_IntType, _RealType>& __x) + { + const std::ios_base::fmtflags __flags = __is.flags(); + __is.flags(std::ios_base::dec | std::ios_base::skipws); + + __is >> __x._M_t >> __x._M_p >> __x._M_nd; + __x._M_initialize(); + + __is.flags(__flags); + return __is; + } + template std::basic_ostream<_CharT, _Traits>& @@ -893,20 +1130,20 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1) result_type __x, __y, __r2; do { - __x = result_type(2.0) * __urng() - result_type(1.0); - __y = result_type(2.0) * __urng() - result_type(1.0); + __x = result_type(2.0) * __urng() - 1.0; + __y = result_type(2.0) * __urng() - 1.0; __r2 = __x * __x + __y * __y; } - while (__r2 > result_type(1.0) || __r2 == result_type(0)); + while (__r2 > 1.0 || __r2 == 0.0); - const result_type __mult = std::sqrt(-result_type(2.0) - * std::log(__r2) / __r2); + const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); _M_saved = __x * __mult; _M_saved_available = true; __ret = __y * __mult; } - - return __ret * _M_sigma + _M_mean; + + __ret = __ret * _M_sigma + _M_mean; + return __ret; } template diff --git a/libstdc++-v3/testsuite/tr1/5_numerical_facilities/random/binomial_distribution/requirements/typedefs.cc b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/random/binomial_distribution/requirements/typedefs.cc new file mode 100644 index 00000000000..fd2d9e47913 --- /dev/null +++ b/libstdc++-v3/testsuite/tr1/5_numerical_facilities/random/binomial_distribution/requirements/typedefs.cc @@ -0,0 +1,37 @@ +// { dg-do compile } +// +// 2006-08-18 Paolo Carlini +// +// Copyright (C) 2006 Free Software Foundation, Inc. +// +// This file is part of the GNU ISO C++ Library. This library is free +// software; you can redistribute it and/or modify it under the +// terms of the GNU General Public License as published by the +// Free Software Foundation; either version 2, or (at your option) +// any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along +// with this library; see the file COPYING. If not, write to the Free +// Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, +// USA. + +// 5.1.7.5 Class template binomial_distribution [tr.rand.dist.bin] +// 5.1.1 [7] Table 17 + +#include + +void +test01() +{ + using namespace std::tr1; + + typedef binomial_distribution test_type; + + typedef test_type::input_type input_type; + typedef test_type::result_type result_type; +}