X-Git-Url: http://git.ipfire.org/?a=blobdiff_plain;f=libstdc%2B%2B-v3%2Finclude%2Fbits%2Frandom.tcc;h=f092b5cd4cbcd29c360c6c10c0284a63b94e9fb3;hb=83ffe9cde7fe0b4deb0d1b54175fd9b19c38179c;hp=0e917b2dc39bc2bc2ea7f8ae1d9b4c5df9f3ceec;hpb=e1a029634255e159dc6817b24b6bcc6497fa400c;p=thirdparty%2Fgcc.git diff --git a/libstdc++-v3/include/bits/random.tcc b/libstdc++-v3/include/bits/random.tcc index 0e917b2dc39b..f092b5cd4cbc 100644 --- a/libstdc++-v3/include/bits/random.tcc +++ b/libstdc++-v3/include/bits/random.tcc @@ -1,6 +1,6 @@ // random number generation (out of line) -*- C++ -*- -// Copyright (C) 2009 Free Software Foundation, Inc. +// Copyright (C) 2009-2023 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 @@ -24,70 +24,90 @@ /** @file bits/random.tcc * This is an internal header file, included by other library headers. - * You should not attempt to use it directly. + * Do not attempt to use it directly. @headername{random} */ -#include -#include +#ifndef _RANDOM_TCC +#define _RANDOM_TCC 1 -namespace std +#include // std::accumulate and std::partial_sum + +namespace std _GLIBCXX_VISIBILITY(default) { - /* - * (Further) implementation-space details. - */ +_GLIBCXX_BEGIN_NAMESPACE_VERSION + + /// @cond undocumented + // (Further) implementation-space details. namespace __detail { - // General case for x = (ax + c) mod m -- use Schrage's algorithm to - // avoid integer overflow. - // - // Because a and c are compile-time integral constants the compiler - // kindly elides any unreachable paths. + // General case for x = (ax + c) mod m -- use Schrage's algorithm + // to avoid integer overflow. // // Preconditions: a > 0, m > 0. // - template - struct _Mod + // Note: only works correctly for __m % __a < __m / __a. + template + _Tp + _Mod<_Tp, __m, __a, __c, false, true>:: + __calc(_Tp __x) { - static _Tp - __calc(_Tp __x) - { - if (__a == 1) - __x %= __m; - else - { - static const _Tp __q = __m / __a; - static const _Tp __r = __m % __a; - - _Tp __t1 = __a * (__x % __q); - _Tp __t2 = __r * (__x / __q); - if (__t1 >= __t2) - __x = __t1 - __t2; - else - __x = __m - __t2 + __t1; - } + if (__a == 1) + __x %= __m; + else + { + static const _Tp __q = __m / __a; + static const _Tp __r = __m % __a; + + _Tp __t1 = __a * (__x % __q); + _Tp __t2 = __r * (__x / __q); + if (__t1 >= __t2) + __x = __t1 - __t2; + else + __x = __m - __t2 + __t1; + } - if (__c != 0) - { - const _Tp __d = __m - __x; - if (__d > __c) - __x += __c; - else - __x = __c - __d; - } - return __x; - } - }; + if (__c != 0) + { + const _Tp __d = __m - __x; + if (__d > __c) + __x += __c; + else + __x = __c - __d; + } + return __x; + } - // Special case for m == 0 -- use unsigned integer overflow as modulo - // operator. - template - struct _Mod<_Tp, __a, __c, __m, true> + template + _OutputIterator + __normalize(_InputIterator __first, _InputIterator __last, + _OutputIterator __result, const _Tp& __factor) { - static _Tp - __calc(_Tp __x) - { return __a * __x + __c; } - }; + for (; __first != __last; ++__first, ++__result) + *__result = *__first / __factor; + return __result; + } + } // namespace __detail + /// @endcond + +#if ! __cpp_inline_variables + template + constexpr _UIntType + linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier; + + template + constexpr _UIntType + linear_congruential_engine<_UIntType, __a, __c, __m>::increment; + + template + constexpr _UIntType + linear_congruential_engine<_UIntType, __a, __c, __m>::modulus; + + template + constexpr _UIntType + linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed; +#endif /** * Seeds the LCR with integral value @p __s, adjusted so that the @@ -98,48 +118,37 @@ namespace std linear_congruential_engine<_UIntType, __a, __c, __m>:: seed(result_type __s) { - if ((__detail::__mod<_UIntType, 1U, 0U, __m>(__c) == 0U) - && (__detail::__mod<_UIntType, 1U, 0U, __m>(__s) == 0U)) - _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(1U); + if ((__detail::__mod<_UIntType, __m>(__c) == 0) + && (__detail::__mod<_UIntType, __m>(__s) == 0)) + _M_x = 1; else - _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(__s); + _M_x = __detail::__mod<_UIntType, __m>(__s); } /** * Seeds the LCR engine with a value generated by @p __q. */ template - void - linear_congruential_engine<_UIntType, __a, __c, __m>:: - seed(seed_seq& __q) - { - const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits - : std::__lg(__m); - const _UIntType __k = (__k0 + 31) / 32; - _UIntType __arr[__k + 3]; - __q.generate(__arr + 0, __arr + __k + 3); - _UIntType __factor = 1U; - _UIntType __sum = 0U; - for (size_t __j = 0; __j < __k; ++__j) - { - __sum += __arr[__j + 3] * __factor; - __factor *= __detail::_Shift<_UIntType, 32>::__value; - } - seed(__sum); - } - - /** - * Gets the next generated value in sequence. - */ - template - 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 + auto + linear_congruential_engine<_UIntType, __a, __c, __m>:: + seed(_Sseq& __q) + -> _If_seed_seq<_Sseq> + { + const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits + : std::__lg(__m); + const _UIntType __k = (__k0 + 31) / 32; + uint_least32_t __arr[__k + 3]; + __q.generate(__arr + 0, __arr + __k + 3); + _UIntType __factor = 1u; + _UIntType __sum = 0u; + for (size_t __j = 0; __j < __k; ++__j) + { + __sum += __arr[__j + 3] * __factor; + __factor *= __detail::_Shift<_UIntType, 32>::__value; + } + seed(__sum); + } template @@ -148,8 +157,7 @@ namespace std const linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -169,8 +177,7 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec); @@ -181,6 +188,134 @@ namespace std return __is; } +#if ! __cpp_inline_variables + template + constexpr size_t + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::word_size; + + template + constexpr size_t + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::state_size; + + template + constexpr size_t + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::shift_size; + + template + constexpr size_t + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::mask_bits; + + template + constexpr _UIntType + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::xor_mask; + + template + constexpr size_t + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::tempering_u; + + template + constexpr _UIntType + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::tempering_d; + + template + constexpr size_t + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::tempering_s; + + template + constexpr _UIntType + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::tempering_b; + + template + constexpr size_t + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::tempering_t; + + template + constexpr _UIntType + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::tempering_c; + + template + constexpr size_t + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::tempering_l; + + template + constexpr _UIntType + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>:: + initialization_multiplier; + + template + constexpr _UIntType + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>::default_seed; +#endif template:: seed(result_type __sd) { - _M_x[0] = __detail::__mod<_UIntType, 1, 0, + _M_x[0] = __detail::__mod<_UIntType, __detail::_Shift<_UIntType, __w>::__value>(__sd); for (size_t __i = 1; __i < state_size; ++__i) @@ -200,8 +335,8 @@ namespace std _UIntType __x = _M_x[__i - 1]; __x ^= __x >> (__w - 2); __x *= __f; - __x += __i; - _M_x[__i] = __detail::__mod<_UIntType, 1, 0, + __x += __detail::__mod<_UIntType, __n>(__i); + _M_x[__i] = __detail::__mod<_UIntType, __detail::_Shift<_UIntType, __w>::__value>(__x); } _M_p = state_size; @@ -212,42 +347,99 @@ namespace std _UIntType __a, size_t __u, _UIntType __d, size_t __s, _UIntType __b, size_t __t, _UIntType __c, size_t __l, _UIntType __f> + template + auto + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>:: + seed(_Sseq& __q) + -> _If_seed_seq<_Sseq> + { + const _UIntType __upper_mask = (~_UIntType()) << __r; + const size_t __k = (__w + 31) / 32; + uint_least32_t __arr[__n * __k]; + __q.generate(__arr + 0, __arr + __n * __k); + + bool __zero = true; + for (size_t __i = 0; __i < state_size; ++__i) + { + _UIntType __factor = 1u; + _UIntType __sum = 0u; + for (size_t __j = 0; __j < __k; ++__j) + { + __sum += __arr[__k * __i + __j] * __factor; + __factor *= __detail::_Shift<_UIntType, 32>::__value; + } + _M_x[__i] = __detail::__mod<_UIntType, + __detail::_Shift<_UIntType, __w>::__value>(__sum); + + if (__zero) + { + if (__i == 0) + { + if ((_M_x[0] & __upper_mask) != 0u) + __zero = false; + } + else if (_M_x[__i] != 0u) + __zero = false; + } + } + if (__zero) + _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value; + _M_p = state_size; + } + + template void mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, - __s, __b, __t, __c, __l, __f>:: - seed(seed_seq& __q) + __s, __b, __t, __c, __l, __f>:: + _M_gen_rand(void) { const _UIntType __upper_mask = (~_UIntType()) << __r; - const size_t __k = (__w + 31) / 32; - _UIntType __arr[__k * __n]; - __q.generate(__arr + 0, __arr + __k * __n); + const _UIntType __lower_mask = ~__upper_mask; - bool __zero = true; - for (size_t __i = 0; __i < state_size; ++__i) + for (size_t __k = 0; __k < (__n - __m); ++__k) { - _UIntType __factor = 1U; - _UIntType __sum = 0U; - for (size_t __j = 0; __j < __k; ++__j) - { - __sum += __arr[__i * __k + __j] * __factor; - __factor *= __detail::_Shift<_UIntType, 32>::__value; - } - _M_x[__i] = __detail::__mod<_UIntType, 1U, 0U, - __detail::_Shift<_UIntType, __w>::__value>(__sum); - - if (__zero) - { - if (__i == 0) - { - if ((_M_x[0] & __upper_mask) != 0U) - __zero = false; - } - else if (_M_x[__i] != 0U) - __zero = false; - } + _UIntType __y = ((_M_x[__k] & __upper_mask) + | (_M_x[__k + 1] & __lower_mask)); + _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1) + ^ ((__y & 0x01) ? __a : 0)); } - if (__zero) - _M_x[0] = __detail::_Shift<_UIntType, __w - 1U>::__value; + + for (size_t __k = (__n - __m); __k < (__n - 1); ++__k) + { + _UIntType __y = ((_M_x[__k] & __upper_mask) + | (_M_x[__k + 1] & __lower_mask)); + _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1) + ^ ((__y & 0x01) ? __a : 0)); + } + + _UIntType __y = ((_M_x[__n - 1] & __upper_mask) + | (_M_x[0] & __lower_mask)); + _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1) + ^ ((__y & 0x01) ? __a : 0)); + _M_p = 0; + } + + template + void + mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, + __s, __b, __t, __c, __l, __f>:: + discard(unsigned long long __z) + { + while (__z > state_size - _M_p) + { + __z -= state_size - _M_p; + _M_gen_rand(); + } + _M_p += __z; } template= state_size) - { - const _UIntType __upper_mask = (~_UIntType()) << __r; - const _UIntType __lower_mask = ~__upper_mask; - - for (size_t __k = 0; __k < (__n - __m); ++__k) - { - _UIntType __y = ((_M_x[__k] & __upper_mask) - | (_M_x[__k + 1] & __lower_mask)); - _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1) - ^ ((__y & 0x01) ? __a : 0)); - } - - for (size_t __k = (__n - __m); __k < (__n - 1); ++__k) - { - _UIntType __y = ((_M_x[__k] & __upper_mask) - | (_M_x[__k + 1] & __lower_mask)); - _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1) - ^ ((__y & 0x01) ? __a : 0)); - } - - _UIntType __y = ((_M_x[__n - 1] & __upper_mask) - | (_M_x[0] & __lower_mask)); - _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1) - ^ ((__y & 0x01) ? __a : 0)); - _M_p = 0; - } + _M_gen_rand(); // Calculate o(x(i)). result_type __z = _M_x[_M_p++]; @@ -311,8 +478,7 @@ namespace std const mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -320,9 +486,9 @@ namespace std __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); __os.fill(__space); - for (size_t __i = 0; __i < __n - 1; ++__i) + for (size_t __i = 0; __i < __n; ++__i) __os << __x._M_x[__i] << __space; - __os << __x._M_x[__n - 1]; + __os << __x._M_p; __os.flags(__flags); __os.fill(__fill); @@ -339,45 +505,59 @@ namespace std mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); for (size_t __i = 0; __i < __n; ++__i) __is >> __x._M_x[__i]; + __is >> __x._M_p; __is.flags(__flags); return __is; } +#if ! __cpp_inline_variables + template + constexpr size_t + subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size; + + template + constexpr size_t + subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag; + + template + constexpr size_t + subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag; + + template + constexpr uint_least32_t + subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed; +#endif template void subtract_with_carry_engine<_UIntType, __w, __s, __r>:: seed(result_type __value) { - if (__value == 0) - __value = default_seed; - - std::linear_congruential_engine - __lcg(__value); + std::linear_congruential_engine + __lcg(__value == 0u ? default_seed : __value); - // I hope this is right. The "10000" tests work for the ranluxen. - const size_t __n = (word_size + 31) / 32; + const size_t __n = (__w + 31) / 32; for (size_t __i = 0; __i < long_lag; ++__i) { - _UIntType __sum = 0U; - _UIntType __factor = 1U; + _UIntType __sum = 0u; + _UIntType __factor = 1u; for (size_t __j = 0; __j < __n; ++__j) { - __sum += __detail::__mod<__detail::_UInt32Type, 1, 0, 0> + __sum += __detail::__mod::__value> (__lcg()) * __factor; __factor *= __detail::_Shift<_UIntType, 32>::__value; } - _M_x[__i] = __detail::__mod<_UIntType, 1, 0, + _M_x[__i] = __detail::__mod<_UIntType, __detail::_Shift<_UIntType, __w>::__value>(__sum); } _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; @@ -385,30 +565,31 @@ namespace std } template - void - subtract_with_carry_engine<_UIntType, __w, __s, __r>:: - seed(seed_seq& __q) - { - const size_t __n = (word_size + 31) / 32; - _UIntType __arr[long_lag + __n]; - __q.generate(__arr + 0, __arr + long_lag + __n); + template + auto + subtract_with_carry_engine<_UIntType, __w, __s, __r>:: + seed(_Sseq& __q) + -> _If_seed_seq<_Sseq> + { + const size_t __k = (__w + 31) / 32; + uint_least32_t __arr[__r * __k]; + __q.generate(__arr + 0, __arr + __r * __k); - for (size_t __i = 0; __i < long_lag; ++__i) - { - _UIntType __sum = 0U; - _UIntType __factor = 1U; - for (size_t __j = 0; __j < __n; ++__j) - { - __sum += __detail::__mod<__detail::_UInt32Type, 1, 0, 0> - (__arr[__i * __n + __j]) * __factor; - __factor *= __detail::_Shift<_UIntType, 32>::__value; - } - _M_x[__i] = __detail::__mod<_UIntType, 1, 0, - __detail::_Shift<_UIntType, __w>::__value>(__sum); - } - _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; - _M_p = 0; - } + for (size_t __i = 0; __i < long_lag; ++__i) + { + _UIntType __sum = 0u; + _UIntType __factor = 1u; + for (size_t __j = 0; __j < __k; ++__j) + { + __sum += __arr[__k * __i + __j] * __factor; + __factor *= __detail::_Shift<_UIntType, 32>::__value; + } + _M_x[__i] = __detail::__mod<_UIntType, + __detail::_Shift<_UIntType, __w>::__value>(__sum); + } + _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; + _M_p = 0; + } template typename subtract_with_carry_engine<_UIntType, __w, __s, __r>:: @@ -452,8 +633,7 @@ namespace std const subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -463,7 +643,7 @@ namespace std for (size_t __i = 0; __i < __r; ++__i) __os << __x._M_x[__i] << __space; - __os << __x._M_carry; + __os << __x._M_carry << __space << __x._M_p; __os.flags(__flags); __os.fill(__fill); @@ -476,8 +656,7 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); @@ -485,11 +664,21 @@ namespace std for (size_t __i = 0; __i < __r; ++__i) __is >> __x._M_x[__i]; __is >> __x._M_carry; + __is >> __x._M_p; __is.flags(__flags); return __is; } +#if ! __cpp_inline_variables + template + constexpr size_t + discard_block_engine<_RandomNumberEngine, __p, __r>::block_size; + + template + constexpr size_t + discard_block_engine<_RandomNumberEngine, __p, __r>::used_block; +#endif template typename discard_block_engine<_RandomNumberEngine, @@ -513,8 +702,7 @@ namespace std const discard_block_engine<_RandomNumberEngine, __p, __r>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -535,8 +723,7 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, discard_block_engine<_RandomNumberEngine, __p, __r>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); @@ -554,54 +741,116 @@ namespace std independent_bits_engine<_RandomNumberEngine, __w, _UIntType>:: operator()() { - const long double __r = static_cast(this->max()) - - static_cast(this->min()) + 1.0L; - const result_type __m = std::log10(__r) / std::log10(2.0L); - result_type __n, __n0, __y0, __y1, __s0, __s1; + typedef typename _RandomNumberEngine::result_type _Eresult_type; + const _Eresult_type __r + = (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max() + ? _M_b.max() - _M_b.min() + 1 : 0); + const unsigned __edig = std::numeric_limits<_Eresult_type>::digits; + const unsigned __m = __r ? std::__lg(__r) : __edig; + + typedef typename std::common_type<_Eresult_type, result_type>::type + __ctype; + const unsigned __cdig = std::numeric_limits<__ctype>::digits; + + unsigned __n, __n0; + __ctype __s0, __s1, __y0, __y1; + for (size_t __i = 0; __i < 2; ++__i) { __n = (__w + __m - 1) / __m + __i; __n0 = __n - __w % __n; - const result_type __w0 = __w / __n; - const result_type __w1 = __w0 + 1; - __s0 = 1UL << __w0; - __s1 = 1UL << __w1; - __y0 = __s0 * (__r / __s0); - __y1 = __s1 * (__r / __s1); - if (__r - __y0 <= __y0 / __n) + const unsigned __w0 = __w / __n; // __w0 <= __m + + __s0 = 0; + __s1 = 0; + if (__w0 < __cdig) + { + __s0 = __ctype(1) << __w0; + __s1 = __s0 << 1; + } + + __y0 = 0; + __y1 = 0; + if (__r) + { + __y0 = __s0 * (__r / __s0); + if (__s1) + __y1 = __s1 * (__r / __s1); + + if (__r - __y0 <= __y0 / __n) + break; + } + else break; } result_type __sum = 0; for (size_t __k = 0; __k < __n0; ++__k) { - result_type __u; + __ctype __u; do - __u = _M_b() - this->min(); - while (__u >= __y0); - __sum = __s0 * __sum - + __u % __s0; + __u = _M_b() - _M_b.min(); + while (__y0 && __u >= __y0); + __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u); } for (size_t __k = __n0; __k < __n; ++__k) { - result_type __u; + __ctype __u; do - __u = _M_b() - this->min(); - while (__u >= __y1); - __sum = __s1 * __sum - + __u % __s1; + __u = _M_b() - _M_b.min(); + while (__y1 && __u >= __y1); + __sum = __s1 * __sum + (__s1 ? __u % __s1 : __u); } return __sum; } +#if ! __cpp_inline_variables + template + constexpr size_t + shuffle_order_engine<_RandomNumberEngine, __k>::table_size; +#endif + + namespace __detail + { + // Determine whether an integer is representable as double. + template + constexpr bool + __representable_as_double(_Tp __x) noexcept + { + static_assert(numeric_limits<_Tp>::is_integer, ""); + static_assert(!numeric_limits<_Tp>::is_signed, ""); + // All integers <= 2^53 are representable. + return (__x <= (1ull << __DBL_MANT_DIG__)) + // Between 2^53 and 2^54 only even numbers are representable. + || (!(__x & 1) && __detail::__representable_as_double(__x >> 1)); + } + + // Determine whether x+1 is representable as double. + template + constexpr bool + __p1_representable_as_double(_Tp __x) noexcept + { + static_assert(numeric_limits<_Tp>::is_integer, ""); + static_assert(!numeric_limits<_Tp>::is_signed, ""); + return numeric_limits<_Tp>::digits < __DBL_MANT_DIG__ + || (bool(__x + 1u) // return false if x+1 wraps around to zero + && __detail::__representable_as_double(__x + 1u)); + } + } template typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type shuffle_order_engine<_RandomNumberEngine, __k>:: operator()() { - size_t __j = __k * ((_M_y - _M_b.min()) - / (_M_b.max() - _M_b.min() + 1.0L)); + constexpr result_type __range = max() - min(); + size_t __j = __k; + const result_type __y = _M_y - min(); + // Avoid using slower long double arithmetic if possible. + if _GLIBCXX17_CONSTEXPR (__detail::__p1_representable_as_double(__range)) + __j *= __y / (__range + 1.0); + else + __j *= __y / (__range + 1.0L); _M_y = _M_v[__j]; _M_v[__j] = _M_b(); @@ -614,8 +863,7 @@ namespace std operator<<(std::basic_ostream<_CharT, _Traits>& __os, const shuffle_order_engine<_RandomNumberEngine, __k>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -639,8 +887,7 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, shuffle_order_engine<_RandomNumberEngine, __k>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); @@ -655,47 +902,12 @@ namespace std } - template - template - typename uniform_int_distribution<_IntType>::result_type - uniform_int_distribution<_IntType>:: - _M_call(_UniformRandomNumberGenerator& __urng, - result_type __min, result_type __max, true_type) - { - // XXX Must be fixed to work well for *arbitrary* __urng.max(), - // __urng.min(), __max, __min. Currently works fine only in the - // most common case __urng.max() - __urng.min() >= __max - __min, - // with __urng.max() > __urng.min() >= 0. - typedef typename __gnu_cxx::__add_unsigned::__type __urntype; - typedef typename __gnu_cxx::__add_unsigned::__type - __utype; - typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype) - > sizeof(__utype)), - __urntype, __utype>::__type __uctype; - - result_type __ret; - - const __urntype __urnmin = __urng.min(); - const __urntype __urnmax = __urng.max(); - const __urntype __urnrange = __urnmax - __urnmin; - const __uctype __urange = __max - __min; - const __uctype __udenom = (__urnrange <= __urange - ? 1 : __urnrange / (__urange + 1)); - do - __ret = (__urntype(__urng()) - __urnmin) / __udenom; - while (__ret > __max - __min); - - return __ret + __min; - } - template std::basic_ostream<_CharT, _Traits>& operator<<(std::basic_ostream<_CharT, _Traits>& __os, const uniform_int_distribution<_IntType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -715,29 +927,45 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, uniform_int_distribution<_IntType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type + = typename uniform_int_distribution<_IntType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); _IntType __a, __b; - __is >> __a >> __b; - __x.param(typename uniform_int_distribution<_IntType>:: - param_type(__a, __b)); + if (__is >> __a >> __b) + __x.param(param_type(__a, __b)); __is.flags(__flags); return __is; } + template + template + void + uniform_real_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __p) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> + __aurng(__urng); + auto __range = __p.b() - __p.a(); + while (__f != __t) + *__f++ = __aurng() * __range + __p.a(); + } + template std::basic_ostream<_CharT, _Traits>& operator<<(std::basic_ostream<_CharT, _Traits>& __os, const uniform_real_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -745,7 +973,7 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); __os << __x.a() << __space << __x.b(); @@ -760,36 +988,52 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, uniform_real_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type + = typename uniform_real_distribution<_RealType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::skipws); _RealType __a, __b; - __is >> __a >> __b; - __x.param(typename uniform_real_distribution<_RealType>:: - param_type(__a, __b)); + if (__is >> __a >> __b) + __x.param(param_type(__a, __b)); __is.flags(__flags); return __is; } + template + void + std::bernoulli_distribution:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __p) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + __detail::_Adaptor<_UniformRandomNumberGenerator, double> + __aurng(__urng); + auto __limit = __p.p() * (__aurng.max() - __aurng.min()); + + while (__f != __t) + *__f++ = (__aurng() - __aurng.min()) < __limit; + } + template std::basic_ostream<_CharT, _Traits>& operator<<(std::basic_ostream<_CharT, _Traits>& __os, const bernoulli_distribution& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); const std::streamsize __precision = __os.precision(); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__os.widen(' ')); - __os.precision(std::numeric_limits::digits10 + 1); + __os.precision(std::numeric_limits::max_digits10); __os << __x.p(); @@ -814,32 +1058,63 @@ namespace std // The largest _RealType convertible to _IntType. const double __thr = std::numeric_limits<_IntType>::max() + __naf; - __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> + __detail::_Adaptor<_UniformRandomNumberGenerator, double> __aurng(__urng); double __cand; do - __cand = std::ceil(std::log(__aurng()) / __param._M_log_p); + __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p); while (__cand >= __thr); return result_type(__cand + __naf); } + template + template + void + geometric_distribution<_IntType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __param) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + // About the epsilon thing see this thread: + // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html + const double __naf = + (1 - std::numeric_limits::epsilon()) / 2; + // The largest _RealType convertible to _IntType. + const double __thr = + std::numeric_limits<_IntType>::max() + __naf; + __detail::_Adaptor<_UniformRandomNumberGenerator, double> + __aurng(__urng); + + while (__f != __t) + { + double __cand; + do + __cand = std::floor(std::log(1.0 - __aurng()) + / __param._M_log_1_p); + while (__cand >= __thr); + + *__f++ = __cand + __naf; + } + } + template std::basic_ostream<_CharT, _Traits>& operator<<(std::basic_ostream<_CharT, _Traits>& __os, const geometric_distribution<_IntType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); const std::streamsize __precision = __os.precision(); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__os.widen(' ')); - __os.precision(std::numeric_limits::digits10 + 1); + __os.precision(std::numeric_limits::max_digits10); __os << __x.p(); @@ -855,20 +1130,34 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, geometric_distribution<_IntType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type = typename geometric_distribution<_IntType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::skipws); double __p; - __is >> __p; - __x.param(typename geometric_distribution<_IntType>::param_type(__p)); + if (__is >> __p) + __x.param(param_type(__p)); __is.flags(__flags); return __is; } + // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5. + template + template + typename negative_binomial_distribution<_IntType>::result_type + negative_binomial_distribution<_IntType>:: + operator()(_UniformRandomNumberGenerator& __urng) + { + const double __y = _M_gd(__urng); + + // XXX Is the constructor too slow? + std::poisson_distribution __poisson(__y); + return __poisson(__urng); + } + template template typename negative_binomial_distribution<_IntType>::result_type @@ -876,17 +1165,55 @@ namespace std operator()(_UniformRandomNumberGenerator& __urng, const param_type& __p) { - typename gamma_distribution<>::param_type - __gamma_param(__p.k(), 1.0); - gamma_distribution<> __gamma(__gamma_param); - double __x = __gamma(__urng); + typedef typename std::gamma_distribution::param_type + param_type; + + const double __y = + _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p())); + + std::poisson_distribution __poisson(__y); + return __poisson(__urng); + } + + template + template + void + negative_binomial_distribution<_IntType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + while (__f != __t) + { + const double __y = _M_gd(__urng); - typename poisson_distribution::param_type - __poisson_param(__x * __p.p() / (1.0 - __p.p())); - poisson_distribution __poisson(__poisson_param); - result_type __m = __poisson(__urng); + // XXX Is the constructor too slow? + std::poisson_distribution __poisson(__y); + *__f++ = __poisson(__urng); + } + } - return __m; + template + template + void + negative_binomial_distribution<_IntType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __p) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + typename std::gamma_distribution::param_type + __p2(__p.k(), (1.0 - __p.p()) / __p.p()); + + while (__f != __t) + { + const double __y = _M_gd(__urng, __p2); + + std::poisson_distribution __poisson(__y); + *__f++ = __poisson(__urng); + } } template @@ -894,8 +1221,7 @@ namespace std operator<<(std::basic_ostream<_CharT, _Traits>& __os, const negative_binomial_distribution<_IntType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -903,9 +1229,10 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__os.widen(' ')); - __os.precision(std::numeric_limits::digits10 + 1); + __os.precision(std::numeric_limits::max_digits10); - __os << __x.k() << __space << __x.p(); + __os << __x.k() << __space << __x.p() + << __space << __x._M_gd; __os.flags(__flags); __os.fill(__fill); @@ -918,17 +1245,17 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, negative_binomial_distribution<_IntType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type + = typename negative_binomial_distribution<_IntType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::skipws); _IntType __k; double __p; - __is >> __k >> __p; - __x.param(typename negative_binomial_distribution<_IntType>:: - param_type(__k, __p)); + if (__is >> __k >> __p >> __x._M_gd) + __x.param(param_type(__k, __p)); __is.flags(__flags); return __is; @@ -951,7 +1278,7 @@ namespace std const double __pi_4 = 0.7853981633974483096156608458198757L; const double __dx = std::sqrt(2 * __m * std::log(32 * __m / __pi_4)); - _M_d = std::round(std::max(6.0, std::min(__m, __dx))); + _M_d = std::round(std::max(6.0, std::min(__m, __dx))); const double __cx = 2 * __m + _M_d; _M_scx = std::sqrt(__cx / 2); _M_1cx = 1 / __cx; @@ -972,7 +1299,7 @@ namespace std * is defined. * * Reference: - * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, + * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!). */ template @@ -1002,6 +1329,8 @@ namespace std const double __c2 = __param._M_c2b + __c1; const double __c3 = __c2 + 1; const double __c4 = __c3 + 1; + // 1 / 78 + const double __178 = 0.0128205128205128205128205128205128L; // e^(1 / 78) const double __e178 = 1.0129030479320018583185514777512983L; const double __c5 = __c4 + __e178; @@ -1012,7 +1341,7 @@ namespace std do { const double __u = __c * __aurng(); - const double __e = -std::log(__aurng()); + const double __e = -std::log(1.0 - __aurng()); double __w = 0.0; @@ -1041,10 +1370,14 @@ namespace std else if (__u <= __c4) __x = 0; else if (__u <= __c5) - __x = 1; + { + __x = 1; + // Only in the Errata, see libstdc++/83237. + __w = __178; + } else { - const double __v = -std::log(__aurng()); + const double __v = -std::log(1.0 - __aurng()); const double __y = __param._M_d + __v * __2cx / __param._M_d; __x = std::ceil(__y); @@ -1077,14 +1410,28 @@ namespace std } } + template + template + void + poisson_distribution<_IntType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __param) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + // We could duplicate everything from operator()... + while (__f != __t) + *__f++ = this->operator()(__urng, __param); + } + template std::basic_ostream<_CharT, _Traits>& operator<<(std::basic_ostream<_CharT, _Traits>& __os, const poisson_distribution<_IntType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -1092,7 +1439,7 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits::digits10 + 1); + __os.precision(std::numeric_limits::max_digits10); __os << __x.mean() << __space << __x._M_nd; @@ -1108,15 +1455,15 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, poisson_distribution<_IntType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type = typename poisson_distribution<_IntType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::skipws); double __mean; - __is >> __mean >> __x._M_nd; - __x.param(typename poisson_distribution<_IntType>::param_type(__mean)); + if (__is >> __mean >> __x._M_nd) + __x.param(param_type(__mean)); __is.flags(__flags); return __is; @@ -1144,11 +1491,11 @@ namespace std const double __d1x = std::sqrt(__np * __1p * std::log(32 * __np / (81 * __pi_4 * __1p))); - _M_d1 = std::round(std::max(1.0, __d1x)); + _M_d1 = std::round(std::max(1.0, __d1x)); const double __d2x = std::sqrt(__np * __1p * std::log(32 * _M_t * __1p / (__pi_4 * __pa))); - _M_d2 = std::round(std::max(1.0, __d2x)); + _M_d2 = std::round(std::max(1.0, __d2x)); // sqrt(pi / 2) const double __spi_2 = 1.2533141373155002512078826424055226L; @@ -1179,7 +1526,8 @@ namespace std template typename binomial_distribution<_IntType>::result_type binomial_distribution<_IntType>:: - _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t) + _M_waiting(_UniformRandomNumberGenerator& __urng, + _IntType __t, double __q) { _IntType __x = 0; double __sum = 0.0; @@ -1188,11 +1536,13 @@ namespace std do { - const double __e = -std::log(__aurng()); + if (__t == __x) + return __x; + const double __e = -std::log(1.0 - __aurng()); __sum += __e / (__t - __x); __x += 1; } - while (__sum <= _M_param._M_q); + while (__sum <= __q); return __x - 1; } @@ -1204,7 +1554,7 @@ namespace std * is defined. * * Reference: - * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, + * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, * New York, 1986, Ch. X, Sect. 4 (+ Errata!). */ template @@ -1216,7 +1566,7 @@ namespace std { result_type __ret; const _IntType __t = __param.t(); - const _IntType __p = __param.p(); + const double __p = __param.p(); const double __p12 = __p <= 0.5 ? __p : 1.0 - __p; __detail::_Adaptor<_UniformRandomNumberGenerator, double> __aurng(__urng); @@ -1256,7 +1606,7 @@ namespace std __reject = __y >= __param._M_d1; if (!__reject) { - const double __e = -std::log(__aurng()); + const double __e = -std::log(1.0 - __aurng()); __x = std::floor(__y); __v = -__e - __n * __n / 2 + __param._M_c; } @@ -1268,15 +1618,15 @@ namespace std __reject = __y >= __param._M_d2; if (!__reject) { - const double __e = -std::log(__aurng()); + const double __e = -std::log(1.0 - __aurng()); __x = std::floor(-__y); __v = -__e - __n * __n / 2; } } else if (__u <= __a123) { - const double __e1 = -std::log(__aurng()); - const double __e2 = -std::log(__aurng()); + const double __e1 = -std::log(1.0 - __aurng()); + const double __e2 = -std::log(1.0 - __aurng()); const double __y = __param._M_d1 + 2 * __s1s * __e1 / __param._M_d1; @@ -1287,8 +1637,8 @@ namespace std } else { - const double __e1 = -std::log(__aurng()); - const double __e2 = -std::log(__aurng()); + const double __e1 = -std::log(1.0 - __aurng()); + const double __e2 = -std::log(1.0 - __aurng()); const double __y = __param._M_d2 + 2 * __s2s * __e1 / __param._M_d2; @@ -1313,26 +1663,41 @@ namespace std __x += __np + __naf; - const _IntType __z = _M_waiting(__urng, __t - _IntType(__x)); + const _IntType __z = _M_waiting(__urng, __t - _IntType(__x), + __param._M_q); __ret = _IntType(__x) + __z; } else #endif - __ret = _M_waiting(__urng, __t); + __ret = _M_waiting(__urng, __t, __param._M_q); if (__p12 != __p) __ret = __t - __ret; return __ret; } + template + template + void + binomial_distribution<_IntType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __param) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + // We could duplicate everything from operator()... + while (__f != __t) + *__f++ = this->operator()(__urng, __param); + } + template std::basic_ostream<_CharT, _Traits>& operator<<(std::basic_ostream<_CharT, _Traits>& __os, const binomial_distribution<_IntType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -1340,7 +1705,7 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits::digits10 + 1); + __os.precision(std::numeric_limits::max_digits10); __os << __x.t() << __space << __x.p() << __space << __x._M_nd; @@ -1357,37 +1722,51 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, binomial_distribution<_IntType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type = typename binomial_distribution<_IntType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); _IntType __t; double __p; - __is >> __t >> __p >> __x._M_nd; - __x.param(typename binomial_distribution<_IntType>:: - param_type(__t, __p)); + if (__is >> __t >> __p >> __x._M_nd) + __x.param(param_type(__t, __p)); __is.flags(__flags); return __is; } + template + template + void + std::exponential_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __p) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> + __aurng(__urng); + while (__f != __t) + *__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda(); + } + template std::basic_ostream<_CharT, _Traits>& operator<<(std::basic_ostream<_CharT, _Traits>& __os, const exponential_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); const std::streamsize __precision = __os.precision(); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__os.widen(' ')); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); __os << __x.lambda(); @@ -1402,16 +1781,16 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, exponential_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type + = typename exponential_distribution<_RealType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); _RealType __lambda; - __is >> __lambda; - __x.param(typename exponential_distribution<_RealType>:: - param_type(__lambda)); + if (__is >> __lambda) + __x.param(param_type(__lambda)); __is.flags(__flags); return __is; @@ -1421,7 +1800,7 @@ namespace std /** * Polar method due to Marsaglia. * - * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag, + * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, * New York, 1986, Ch. V, Sect. 4.4. */ template @@ -1461,13 +1840,84 @@ namespace std return __ret; } + template + template + void + normal_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __param) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + + if (__f == __t) + return; + + if (_M_saved_available) + { + _M_saved_available = false; + *__f++ = _M_saved * __param.stddev() + __param.mean(); + + if (__f == __t) + return; + } + + __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> + __aurng(__urng); + + while (__f + 1 < __t) + { + result_type __x, __y, __r2; + do + { + __x = result_type(2.0) * __aurng() - 1.0; + __y = result_type(2.0) * __aurng() - 1.0; + __r2 = __x * __x + __y * __y; + } + while (__r2 > 1.0 || __r2 == 0.0); + + const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); + *__f++ = __y * __mult * __param.stddev() + __param.mean(); + *__f++ = __x * __mult * __param.stddev() + __param.mean(); + } + + if (__f != __t) + { + result_type __x, __y, __r2; + do + { + __x = result_type(2.0) * __aurng() - 1.0; + __y = result_type(2.0) * __aurng() - 1.0; + __r2 = __x * __x + __y * __y; + } + while (__r2 > 1.0 || __r2 == 0.0); + + const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); + _M_saved = __x * __mult; + _M_saved_available = true; + *__f = __y * __mult * __param.stddev() + __param.mean(); + } + } + + template + bool + operator==(const std::normal_distribution<_RealType>& __d1, + const std::normal_distribution<_RealType>& __d2) + { + if (__d1._M_param == __d2._M_param + && __d1._M_saved_available == __d2._M_saved_available) + return __d1._M_saved_available ? __d1._M_saved == __d2._M_saved : true; + else + return false; + } + template std::basic_ostream<_CharT, _Traits>& operator<<(std::basic_ostream<_CharT, _Traits>& __os, const normal_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -1475,7 +1925,7 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); __os << __x.mean() << __space << __x.stddev() << __space << __x._M_saved_available; @@ -1493,50 +1943,40 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, normal_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type = typename normal_distribution<_RealType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); double __mean, __stddev; - __is >> __mean >> __stddev - >> __x._M_saved_available; - if (__x._M_saved_available) - __is >> __x._M_saved; - __x.param(typename normal_distribution<_RealType>:: - param_type(__mean, __stddev)); - - __is.flags(__flags); - return __is; - } - - - template - template - 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; + bool __saved_avail; + if (__is >> __mean >> __stddev >> __saved_avail) + { + if (!__saved_avail || (__is >> __x._M_saved)) + { + __x._M_saved_available = __saved_avail; + __x.param(param_type(__mean, __stddev)); + } + } - // See if it is in the unit circle. - __r2 = __u * __u + __v * __v; - } - while (__r2 > 1 || __r2 == 0); + __is.flags(__flags); + return __is; + } - __normal = __u * std::sqrt(-2 * std::log(__r2) / __r2); - return std::exp(__p.s() * __normal + __p.m()); + template + template + void + lognormal_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __p) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + while (__f != __t) + *__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m()); } template @@ -1544,8 +1984,7 @@ namespace std operator<<(std::basic_ostream<_CharT, _Traits>& __os, const lognormal_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -1553,9 +1992,10 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); - __os << __x.m() << __space << __x.s(); + __os << __x.m() << __space << __x.s() + << __space << __x._M_nd; __os.flags(__flags); __os.fill(__fill); @@ -1568,33 +2008,47 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, lognormal_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type + = typename lognormal_distribution<_RealType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); _RealType __m, __s; - __is >> __m >> __s; - __x.param(typename lognormal_distribution<_RealType>:: - param_type(__m, __s)); + if (__is >> __m >> __s >> __x._M_nd) + __x.param(param_type(__m, __s)); __is.flags(__flags); return __is; } + template + template + void + std::chi_squared_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + while (__f != __t) + *__f++ = 2 * _M_gd(__urng); + } template - template - typename chi_squared_distribution<_RealType>::result_type - chi_squared_distribution<_RealType>:: - operator()(_UniformRandomNumberGenerator& __urng, - const param_type& __p) + template + void + std::chi_squared_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const typename + std::gamma_distribution::param_type& __p) { - typename gamma_distribution<_RealType>::param_type - __gamma_param(__p.n() / 2, 1.0); - gamma_distribution<_RealType> __gamma(__gamma_param); - return 2 * __gamma(__urng); + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + while (__f != __t) + *__f++ = 2 * _M_gd(__urng, __p); } template @@ -1602,8 +2056,7 @@ namespace std operator<<(std::basic_ostream<_CharT, _Traits>& __os, const chi_squared_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -1611,9 +2064,9 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); - __os << __x.n(); + __os << __x.n() << __space << __x._M_gd; __os.flags(__flags); __os.fill(__fill); @@ -1626,16 +2079,16 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, chi_squared_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type + = typename chi_squared_distribution<_RealType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); _RealType __n; - __is >> __n; - __x.param(typename chi_squared_distribution<_RealType>:: - param_type(__n)); + if (__is >> __n >> __x._M_gd) + __x.param(param_type(__n)); __is.flags(__flags); return __is; @@ -1660,13 +2113,36 @@ namespace std return __p.a() + __p.b() * std::tan(__pi * __u); } + template + template + void + cauchy_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __p) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + const _RealType __pi = 3.1415926535897932384626433832795029L; + __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> + __aurng(__urng); + while (__f != __t) + { + _RealType __u; + do + __u = __aurng(); + while (__u == 0.5); + + *__f++ = __p.a() + __p.b() * std::tan(__pi * __u); + } + } + template std::basic_ostream<_CharT, _Traits>& operator<<(std::basic_ostream<_CharT, _Traits>& __os, const cauchy_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -1674,7 +2150,7 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); __os << __x.a() << __space << __x.b(); @@ -1689,16 +2165,15 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, cauchy_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type = typename cauchy_distribution<_RealType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); _RealType __a, __b; - __is >> __a >> __b; - __x.param(typename cauchy_distribution<_RealType>:: - param_type(__a, __b)); + if (__is >> __a >> __b) + __x.param(param_type(__a, __b)); __is.flags(__flags); return __is; @@ -1706,20 +2181,35 @@ namespace std template - template - typename fisher_f_distribution<_RealType>::result_type - fisher_f_distribution<_RealType>:: - operator()(_UniformRandomNumberGenerator& __urng, - const param_type& __p) + template + void + std::fisher_f_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng) { - gamma_distribution<_RealType> __gamma; - _RealType __ym = __gamma(__urng, - typename gamma_distribution<_RealType>::param_type(__p.m() / 2, 2)); - - _RealType __yn = __gamma(__urng, - typename gamma_distribution<_RealType>::param_type(__p.n() / 2, 2)); + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + while (__f != __t) + *__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m())); + } - return (__ym * __p.n()) / (__yn * __p.m()); + template + template + void + std::fisher_f_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __p) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + typedef typename std::gamma_distribution::param_type + param_type; + param_type __p1(__p.m() / 2); + param_type __p2(__p.n() / 2); + while (__f != __t) + *__f++ = ((_M_gd_x(__urng, __p1) * n()) + / (_M_gd_y(__urng, __p2) * m())); } template @@ -1727,8 +2217,7 @@ namespace std operator<<(std::basic_ostream<_CharT, _Traits>& __os, const fisher_f_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -1736,9 +2225,10 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); - __os << __x.m() << __space << __x.n(); + __os << __x.m() << __space << __x.n() + << __space << __x._M_gd_x << __space << __x._M_gd_y; __os.flags(__flags); __os.fill(__fill); @@ -1751,89 +2241,49 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, fisher_f_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type + = typename fisher_f_distribution<_RealType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); _RealType __m, __n; - __is >> __m >> __n; - __x.param(typename fisher_f_distribution<_RealType>:: - param_type(__m, __n)); + if (__is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y) + __x.param(param_type(__m, __n)); __is.flags(__flags); return __is; } - // - // This could be operator() for a Gaussian distribution. - // template - template - typename student_t_distribution<_RealType>::result_type - student_t_distribution<_RealType>:: - _M_gaussian(_UniformRandomNumberGenerator& __urng, - const result_type __sigma) + template + void + std::student_t_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng) { - _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); + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + while (__f != __t) + *__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng)); } template - template - typename student_t_distribution<_RealType>::result_type - student_t_distribution<_RealType>:: - operator()(_UniformRandomNumberGenerator& __urng, - const param_type& __param) + template + void + std::student_t_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __p) { - 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 __y2 = __chisq(__urng); - - return __y1 / std::sqrt(__y2 / __param.n()); - } - else - { - _RealType __y1, __y2, __z; - 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); - __y2 = __exponential(__urng); - - __z = __y1 * __y1 / (__param.n() - 2.0); - } - while (1.0 - __z < 0.0 || std::exp(-__y2 - __z) > (1.0 - __z)); - - // Note that there is a typo in Knuth's formula, the line below - // is taken from the original paper of Marsaglia, Mathematics of - // Computation, 34 (1980), p 234-256 - return __y1 / std::sqrt((1.0 - 2.0 / __param.n()) * (1.0 - __z)); - } + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + typename std::gamma_distribution::param_type + __p2(__p.n() / 2, 2); + while (__f != __t) + *__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2)); } template @@ -1841,8 +2291,7 @@ namespace std operator<<(std::basic_ostream<_CharT, _Traits>& __os, const student_t_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -1850,9 +2299,9 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); - __os << __x.n(); + __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd; __os.flags(__flags); __os.fill(__fill); @@ -1865,15 +2314,16 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, student_t_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type + = typename student_t_distribution<_RealType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); _RealType __n; - __is >> __n; - __x.param(typename student_t_distribution<_RealType>::param_type(__n)); + if (__is >> __n >> __x._M_nd >> __x._M_gd) + __x.param(param_type(__n)); __is.flags(__flags); return __is; @@ -1885,28 +2335,16 @@ namespace std 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 template @@ -1915,58 +2353,106 @@ namespace std 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.0331 * __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); + __u = __aurng(); + while (__u == 0.0); + + return (std::pow(__u, result_type(1.0) / __param.alpha()) + * __a1 * __v * __param.beta()); + } + } - __x = std::pow(__z, __c); + template + template + void + gamma_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __param) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> + __aurng(__urng); - __reject = __z + __e < __param._M_l_d + __x; - } - while (__reject); - } + result_type __u, __v, __n; + const result_type __a1 = (__param._M_malpha + - _RealType(1.0) / _RealType(3.0)); - return __beta_val * __x; + if (__param.alpha() == __param._M_malpha) + while (__f != __t) + { + do + { + do + { + __n = _M_nd(__urng); + __v = result_type(1.0) + __param._M_a2 * __n; + } + while (__v <= 0.0); + + __v = __v * __v * __v; + __u = __aurng(); + } + while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n + && (std::log(__u) > (0.5 * __n * __n + __a1 + * (1.0 - __v + std::log(__v))))); + + *__f++ = __a1 * __v * __param.beta(); + } + else + while (__f != __t) + { + do + { + do + { + __n = _M_nd(__urng); + __v = result_type(1.0) + __param._M_a2 * __n; + } + while (__v <= 0.0); + + __v = __v * __v * __v; + __u = __aurng(); + } + while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n + && (std::log(__u) > (0.5 * __n * __n + __a1 + * (1.0 - __v + std::log(__v))))); + + do + __u = __aurng(); + while (__u == 0.0); + + *__f++ = (std::pow(__u, result_type(1.0) / __param.alpha()) + * __a1 * __v * __param.beta()); + } } template @@ -1974,8 +2460,7 @@ namespace std operator<<(std::basic_ostream<_CharT, _Traits>& __os, const gamma_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -1983,9 +2468,10 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); - __os << __x.alpha() << __space << __x.beta(); + __os << __x.alpha() << __space << __x.beta() + << __space << __x._M_nd; __os.flags(__flags); __os.fill(__fill); @@ -1998,29 +2484,59 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, gamma_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type = typename gamma_distribution<_RealType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); _RealType __alpha_val, __beta_val; - __is >> __alpha_val >> __beta_val; - __x.param(typename gamma_distribution<_RealType>:: - param_type(__alpha_val, __beta_val)); + if (__is >> __alpha_val >> __beta_val >> __x._M_nd) + __x.param(param_type(__alpha_val, __beta_val)); __is.flags(__flags); return __is; } + template + template + 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(result_type(1) - __aurng()), + result_type(1) / __p.a()); + } + + template + template + void + weibull_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __p) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> + __aurng(__urng); + auto __inv_a = result_type(1) / __p.a(); + + while (__f != __t) + *__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()), + __inv_a); + } + template std::basic_ostream<_CharT, _Traits>& operator<<(std::basic_ostream<_CharT, _Traits>& __os, const weibull_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -2028,7 +2544,7 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); __os << __x.a() << __space << __x.b(); @@ -2043,16 +2559,15 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, weibull_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type = typename weibull_distribution<_RealType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); _RealType __a, __b; - __is >> __a >> __b; - __x.param(typename weibull_distribution<_RealType>:: - param_type(__a, __b)); + if (__is >> __a >> __b) + __x.param(param_type(__a, __b)); __is.flags(__flags); return __is; @@ -2068,7 +2583,26 @@ namespace std { __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> __aurng(__urng); - return __p.a() - __p.b() * std::log(-std::log(__aurng())); + return __p.a() - __p.b() * std::log(-std::log(result_type(1) + - __aurng())); + } + + template + template + void + extreme_value_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __p) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> + __aurng(__urng); + + while (__f != __t) + *__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1) + - __aurng())); } template @@ -2076,8 +2610,7 @@ namespace std operator<<(std::basic_ostream<_CharT, _Traits>& __os, const extreme_value_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -2085,7 +2618,7 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); __os << __x.a() << __space << __x.b(); @@ -2100,16 +2633,16 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, extreme_value_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using param_type + = typename extreme_value_distribution<_RealType>::param_type; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); _RealType __a, __b; - __is >> __a >> __b; - __x.param(typename extreme_value_distribution<_RealType>:: - param_type(__a, __b)); + if (__is >> __a >> __b) + __x.param(param_type(__a, __b)); __is.flags(__flags); return __is; @@ -2124,34 +2657,35 @@ namespace std if (_M_prob.size() < 2) { _M_prob.clear(); - _M_prob.push_back(1.0); return; } - double __sum = std::accumulate(_M_prob.begin(), _M_prob.end(), 0.0); - // Now normalize the densities. - std::transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(), - std::bind2nd(std::divides(), __sum)); - // Accumulate partial sums. + const double __sum = std::accumulate(_M_prob.begin(), + _M_prob.end(), 0.0); + __glibcxx_assert(__sum > 0); + // Now normalize the probabilites. + __detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(), + __sum); + // Accumulate partial sums. + _M_cp.reserve(_M_prob.size()); std::partial_sum(_M_prob.begin(), _M_prob.end(), std::back_inserter(_M_cp)); - // Make sure the last cumulative probablility is one. + // Make sure the last cumulative probability is one. _M_cp[_M_cp.size() - 1] = 1.0; } template template discrete_distribution<_IntType>::param_type:: - param_type(size_t __nw, double __xmin, double __xmax, - _Func __fw) + param_type(size_t __nw, double __xmin, double __xmax, _Func __fw) : _M_prob(), _M_cp() { - for (size_t __i = 0; __i < __nw; ++__i) - { - const double __x = ((__nw - __i - 0.5) * __xmin - + (__i + 0.5) * __xmax) / __nw; - _M_prob.push_back(__fw(__x)); - } + const size_t __n = __nw == 0 ? 1 : __nw; + const double __delta = (__xmax - __xmin) / __n; + + _M_prob.reserve(__n); + for (size_t __k = 0; __k < __nw; ++__k) + _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta)); _M_initialize(); } @@ -2163,17 +2697,48 @@ namespace std operator()(_UniformRandomNumberGenerator& __urng, const param_type& __param) { - __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> + if (__param._M_cp.empty()) + return result_type(0); + + __detail::_Adaptor<_UniformRandomNumberGenerator, double> __aurng(__urng); const double __p = __aurng(); auto __pos = std::lower_bound(__param._M_cp.begin(), __param._M_cp.end(), __p); - if (__pos == __param._M_cp.end()) - return 0; - const size_t __i = __pos - __param._M_cp.begin(); - return __i; + return __pos - __param._M_cp.begin(); + } + + template + template + void + discrete_distribution<_IntType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __param) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + + if (__param._M_cp.empty()) + { + while (__f != __t) + *__f++ = result_type(0); + return; + } + + __detail::_Adaptor<_UniformRandomNumberGenerator, double> + __aurng(__urng); + + while (__f != __t) + { + const double __p = __aurng(); + auto __pos = std::lower_bound(__param._M_cp.begin(), + __param._M_cp.end(), __p); + + *__f++ = __pos - __param._M_cp.begin(); + } } template @@ -2181,8 +2746,7 @@ namespace std operator<<(std::basic_ostream<_CharT, _Traits>& __os, const discrete_distribution<_IntType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -2190,7 +2754,7 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits::digits10 + 1); + __os.precision(std::numeric_limits::max_digits10); std::vector __prob = __x.probabilities(); __os << __prob.size(); @@ -2203,31 +2767,44 @@ namespace std return __os; } +namespace __detail +{ + template + basic_istream<_CharT, _Traits>& + __extract_params(basic_istream<_CharT, _Traits>& __is, + vector<_ValT>& __vals, size_t __n) + { + __vals.reserve(__n); + while (__n--) + { + _ValT __val; + if (__is >> __val) + __vals.push_back(__val); + else + break; + } + return __is; + } +} // namespace __detail + template std::basic_istream<_CharT, _Traits>& operator>>(std::basic_istream<_CharT, _Traits>& __is, discrete_distribution<_IntType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); size_t __n; - __is >> __n; - - std::vector __prob_vec; - for (; __n != 0; --__n) + if (__is >> __n) { - double __prob; - __is >> __prob; - __prob_vec.push_back(__prob); + std::vector __prob_vec; + if (__detail::__extract_params(__is, __prob_vec, __n)) + __x.param({__prob_vec.begin(), __prob_vec.end()}); } - __x.param(typename discrete_distribution<_IntType>:: - param_type(__prob_vec.begin(), __prob_vec.end())); - __is.flags(__flags); return __is; } @@ -2238,40 +2815,33 @@ namespace std piecewise_constant_distribution<_RealType>::param_type:: _M_initialize() { - if (_M_int.size() < 2) + if (_M_int.size() < 2 + || (_M_int.size() == 2 + && _M_int[0] == _RealType(0) + && _M_int[1] == _RealType(1))) { _M_int.clear(); - _M_int.push_back(_RealType(0)); - _M_int.push_back(_RealType(1)); - _M_den.clear(); - _M_den.push_back(1.0); - return; } - double __sum = 0.0; - for (size_t __i = 0; __i < _M_den.size(); ++__i) - { - __sum += _M_den[__i] * (_M_int[__i + 1] - _M_int[__i]); - _M_cp.push_back(__sum); - } + const double __sum = std::accumulate(_M_den.begin(), + _M_den.end(), 0.0); + __glibcxx_assert(__sum > 0); - // Now normalize the densities... - std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(), - std::bind2nd(std::divides(), __sum)); - // ... and partial sums. - std::transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), - std::bind2nd(std::divides(), __sum)); - // Make sure the last cumulative probablility is one. + __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(), + __sum); + + _M_cp.reserve(_M_den.size()); + std::partial_sum(_M_den.begin(), _M_den.end(), + std::back_inserter(_M_cp)); + + // Make sure the last cumulative probability is one. _M_cp[_M_cp.size() - 1] = 1.0; - } - template - piecewise_constant_distribution<_RealType>::param_type:: - param_type() - : _M_int(), _M_den(), _M_cp() - { _M_initialize(); } + for (size_t __k = 0; __k < _M_den.size(); ++__k) + _M_den[__k] /= _M_int[__k + 1] - _M_int[__k]; + } template template @@ -2281,17 +2851,19 @@ namespace std _InputIteratorW __wbegin) : _M_int(), _M_den(), _M_cp() { - do + if (__bbegin != __bend) { - _M_int.push_back(*__bbegin); - ++__bbegin; - if (__bbegin != __bend) + for (;;) { + _M_int.push_back(*__bbegin); + ++__bbegin; + if (__bbegin == __bend) + break; + _M_den.push_back(*__wbegin); ++__wbegin; } } - while (__bbegin != __bend); _M_initialize(); } @@ -2299,17 +2871,16 @@ namespace std template template piecewise_constant_distribution<_RealType>::param_type:: - param_type(initializer_list<_RealType> __bil, _Func __fw) + param_type(initializer_list<_RealType> __bl, _Func __fw) : _M_int(), _M_den(), _M_cp() { - for (auto __biter = __bil.begin(); __biter != __bil.end(); ++__biter) + _M_int.reserve(__bl.size()); + for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) _M_int.push_back(*__biter); - for (size_t __i = 0; __i < _M_int.size() - 1; ++__i) - { - _RealType __x = 0.5 * (_M_int[__i] + _M_int[__i + 1]); - _M_den.push_back(__fw(__x)); - } + _M_den.reserve(_M_int.size() - 1); + for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) + _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k]))); _M_initialize(); } @@ -2317,22 +2888,19 @@ namespace std template template 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) - { - const _RealType __x = ((__nw - __i) * __xmin - + __i * __xmax) / __nw; - _M_int.push_back(__x); - } - for (size_t __i = 0; __i < __nw; ++__i) - { - const _RealType __x = ((__nw - __i - 0.5) * __xmin - + (__i + 0.5) * __xmax) / __nw; - _M_den.push_back(__fw(__x)); - } + const size_t __n = __nw == 0 ? 1 : __nw; + const _RealType __delta = (__xmax - __xmin) / __n; + + _M_int.reserve(__n + 1); + for (size_t __k = 0; __k <= __nw; ++__k) + _M_int.push_back(__xmin + __k * __delta); + + _M_den.reserve(__n); + for (size_t __k = 0; __k < __nw; ++__k) + _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta)); _M_initialize(); } @@ -2344,16 +2912,55 @@ namespace std operator()(_UniformRandomNumberGenerator& __urng, const param_type& __param) { - __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> + __detail::_Adaptor<_UniformRandomNumberGenerator, double> __aurng(__urng); const double __p = __aurng(); + if (__param._M_cp.empty()) + return __p; + auto __pos = std::lower_bound(__param._M_cp.begin(), __param._M_cp.end(), __p); const size_t __i = __pos - __param._M_cp.begin(); - return __param._M_int[__i] - + (__p - __param._M_cp[__i]) / __param._M_den[__i]; + const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; + + return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i]; + } + + template + template + void + piecewise_constant_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __param) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + __detail::_Adaptor<_UniformRandomNumberGenerator, double> + __aurng(__urng); + + if (__param._M_cp.empty()) + { + while (__f != __t) + *__f++ = __aurng(); + return; + } + + while (__f != __t) + { + const double __p = __aurng(); + + auto __pos = std::lower_bound(__param._M_cp.begin(), + __param._M_cp.end(), __p); + const size_t __i = __pos - __param._M_cp.begin(); + + const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; + + *__f++ = (__param._M_int[__i] + + (__p - __pref) / __param._M_den[__i]); + } } template @@ -2361,8 +2968,7 @@ namespace std operator<<(std::basic_ostream<_CharT, _Traits>& __os, const piecewise_constant_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -2370,7 +2976,7 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); std::vector<_RealType> __int = __x.intervals(); __os << __int.size() - 1; @@ -2393,34 +2999,26 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, piecewise_constant_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); size_t __n; - __is >> __n; - - std::vector<_RealType> __int_vec; - for (size_t __i = 0; __i <= __n; ++__i) + if (__is >> __n) { - _RealType __int; - __is >> __int; - __int_vec.push_back(__int); - } - - std::vector __den_vec; - for (size_t __i = 0; __i < __n; ++__i) - { - double __den; - __is >> __den; - __den_vec.push_back(__den); + std::vector<_RealType> __int_vec; + if (__detail::__extract_params(__is, __int_vec, __n + 1)) + { + std::vector __den_vec; + if (__detail::__extract_params(__is, __den_vec, __n)) + { + __x.param({ __int_vec.begin(), __int_vec.end(), + __den_vec.begin() }); + } + } } - __x.param(typename piecewise_constant_distribution<_RealType>:: - param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); - __is.flags(__flags); return __is; } @@ -2431,46 +3029,40 @@ namespace std piecewise_linear_distribution<_RealType>::param_type:: _M_initialize() { - if (_M_int.size() < 2) + if (_M_int.size() < 2 + || (_M_int.size() == 2 + && _M_int[0] == _RealType(0) + && _M_int[1] == _RealType(1) + && _M_den[0] == _M_den[1])) { _M_int.clear(); - _M_int.push_back(_RealType(0)); - _M_int.push_back(_RealType(1)); - _M_den.clear(); - _M_den.push_back(1.0); - _M_den.push_back(1.0); - return; } double __sum = 0.0; - for (size_t __i = 0; __i < _M_int.size() - 1; ++__i) + _M_cp.reserve(_M_int.size() - 1); + _M_m.reserve(_M_int.size() - 1); + for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) { - const _RealType __delta = _M_int[__i + 1] - _M_int[__i]; - __sum += 0.5 * (_M_den[__i + 1] + _M_den[__i]) * __delta; + const _RealType __delta = _M_int[__k + 1] - _M_int[__k]; + __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta; _M_cp.push_back(__sum); - _M_m.push_back((_M_den[__i + 1] - _M_den[__i]) / __delta); + _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta); } + __glibcxx_assert(__sum > 0); // Now normalize the densities... - std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(), - std::bind2nd(std::divides(),__sum)); + __detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(), + __sum); // ... and partial sums... - std::transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), - std::bind2nd(std::divides(), __sum)); + __detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum); // ... and slopes. - std::transform(_M_m.begin(), _M_m.end(), _M_m.begin(), - std::bind2nd(std::divides(), __sum)); + __detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum); + // Make sure the last cumulative probablility is one. _M_cp[_M_cp.size() - 1] = 1.0; - } - - template - piecewise_linear_distribution<_RealType>::param_type:: - param_type() - : _M_int(), _M_den(), _M_cp(), _M_m() - { _M_initialize(); } + } template template @@ -2492,10 +3084,12 @@ namespace std template template piecewise_linear_distribution<_RealType>::param_type:: - param_type(initializer_list<_RealType> __bil, _Func __fw) + param_type(initializer_list<_RealType> __bl, _Func __fw) : _M_int(), _M_den(), _M_cp(), _M_m() { - for (auto __biter = __bil.begin(); __biter != __bil.end(); ++__biter) + _M_int.reserve(__bl.size()); + _M_den.reserve(__bl.size()); + for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) { _M_int.push_back(*__biter); _M_den.push_back(__fw(*__biter)); @@ -2507,16 +3101,18 @@ namespace std template template piecewise_linear_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(), _M_m() { - for (size_t __i = 0; __i <= __nw; ++__i) + const size_t __n = __nw == 0 ? 1 : __nw; + const _RealType __delta = (__xmax - __xmin) / __n; + + _M_int.reserve(__n + 1); + _M_den.reserve(__n + 1); + for (size_t __k = 0; __k <= __nw; ++__k) { - const _RealType __x = ((__nw - __i) * __xmin - + __i * __xmax) / __nw; - _M_int.push_back(__x); - _M_den.push_back(__fw(__x)); + _M_int.push_back(__xmin + __k * __delta); + _M_den.push_back(__fw(_M_int[__k] + __delta)); } _M_initialize(); @@ -2529,31 +3125,48 @@ namespace std operator()(_UniformRandomNumberGenerator& __urng, const param_type& __param) { - result_type __x; - __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> + __detail::_Adaptor<_UniformRandomNumberGenerator, double> __aurng(__urng); const double __p = __aurng(); + if (__param._M_cp.empty()) + return __p; + auto __pos = std::lower_bound(__param._M_cp.begin(), __param._M_cp.end(), __p); const size_t __i = __pos - __param._M_cp.begin(); + + const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; + const double __a = 0.5 * __param._M_m[__i]; const double __b = __param._M_den[__i]; - const double __c = __param._M_cp[__i]; - const double __q = -0.5 * (__b -#if _GLIBCXX_USE_C99_MATH_TR1 - + std::copysign(std::sqrt(__b * __b - - 4.0 * __a * __c), __b)); -#else - + (__b < 0.0 ? -1.0 : 1.0) - * std::sqrt(__b * __b - 4.0 * __a * __c)); -#endif - const double __x0 = __param._M_int[__i]; - const double __x1 = __q / __a; - const double __x2 = __c / __q; - __x = std::max(__x0 + __x1, __x0 + __x2); + const double __cm = __p - __pref; - return __x; + _RealType __x = __param._M_int[__i]; + if (__a == 0) + __x += __cm / __b; + else + { + const double __d = __b * __b + 4.0 * __a * __cm; + __x += 0.5 * (std::sqrt(__d) - __b) / __a; + } + + return __x; + } + + template + template + void + piecewise_linear_distribution<_RealType>:: + __generate_impl(_ForwardIterator __f, _ForwardIterator __t, + _UniformRandomNumberGenerator& __urng, + const param_type& __param) + { + __glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>) + // We could duplicate everything from operator()... + while (__f != __t) + *__f++ = this->operator()(__urng, __param); } template @@ -2561,8 +3174,7 @@ namespace std operator<<(std::basic_ostream<_CharT, _Traits>& __os, const piecewise_linear_distribution<_RealType>& __x) { - typedef std::basic_ostream<_CharT, _Traits> __ostream_type; - typedef typename __ostream_type::ios_base __ios_base; + using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __os.flags(); const _CharT __fill = __os.fill(); @@ -2570,7 +3182,7 @@ namespace std const _CharT __space = __os.widen(' '); __os.flags(__ios_base::scientific | __ios_base::left); __os.fill(__space); - __os.precision(std::numeric_limits<_RealType>::digits10 + 1); + __os.precision(std::numeric_limits<_RealType>::max_digits10); std::vector<_RealType> __int = __x.intervals(); __os << __int.size() - 1; @@ -2593,52 +3205,47 @@ namespace std operator>>(std::basic_istream<_CharT, _Traits>& __is, piecewise_linear_distribution<_RealType>& __x) { - typedef std::basic_istream<_CharT, _Traits> __istream_type; - typedef typename __istream_type::ios_base __ios_base; + using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base; const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); size_t __n; - __is >> __n; - - std::vector<_RealType> __int_vec; - for (size_t __i = 0; __i <= __n; ++__i) - { - _RealType __int; - __is >> __int; - __int_vec.push_back(__int); - } - - std::vector __den_vec; - for (size_t __i = 0; __i <= __n; ++__i) + if (__is >> __n) { - double __den; - __is >> __den; - __den_vec.push_back(__den); + vector<_RealType> __int_vec; + if (__detail::__extract_params(__is, __int_vec, __n + 1)) + { + vector __den_vec; + if (__detail::__extract_params(__is, __den_vec, __n + 1)) + { + __x.param({ __int_vec.begin(), __int_vec.end(), + __den_vec.begin() }); + } + } } - - __x.param(typename piecewise_linear_distribution<_RealType>:: - param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); - __is.flags(__flags); return __is; } - template + template seed_seq::seed_seq(std::initializer_list<_IntType> __il) { + _M_v.reserve(__il.size()); for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter) - _M_v.push_back(__detail::__mod::__value>(*__iter)); } template seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end) { + if _GLIBCXX17_CONSTEXPR (__is_random_access_iter<_InputIterator>::value) + _M_v.reserve(std::distance(__begin, __end)); + for (_InputIterator __iter = __begin; __iter != __end; ++__iter) - _M_v.push_back(__detail::__mod::__value>(*__iter)); } @@ -2648,12 +3255,12 @@ namespace std _RandomAccessIterator __end) { typedef typename iterator_traits<_RandomAccessIterator>::value_type - __Type; + _Type; if (__begin == __end) return; - std::fill(__begin, __end, __Type(0x8b8b8b8bU)); + std::fill(__begin, __end, _Type(0x8b8b8b8bu)); const size_t __n = __end - __begin; const size_t __s = _M_v.size(); @@ -2664,44 +3271,72 @@ namespace std : (__n - 1) / 2; const size_t __p = (__n - __t) / 2; const size_t __q = __p + __t; - const size_t __m = std::max(__s + 1, __n); + const size_t __m = std::max(size_t(__s + 1), __n); + +#ifndef __UINT32_TYPE__ + struct _Up + { + _Up(uint_least32_t v) : _M_v(v & 0xffffffffu) { } + + operator uint_least32_t() const { return _M_v; } + + uint_least32_t _M_v; + }; + using uint32_t = _Up; +#endif - for (size_t __k = 0; __k < __m; ++__k) + // k == 0, every element in [begin,end) equals 0x8b8b8b8bu { - __Type __arg = __begin[__k % __n] - ^ __begin[(__k + __p) % __n] - ^ __begin[(__k - 1) % __n]; - __Type __r1 = __arg ^ (__arg << 27); - __r1 = __detail::__mod<__Type, 1664525U, 0U, - __detail::_Shift<__Type, 32>::__value>(__r1); - __Type __r2 = __r1; - if (__k == 0) - __r2 += __s; - else if (__k <= __s) - __r2 += __k % __n + _M_v[__k - 1]; - else - __r2 += __k % __n; - __r2 = __detail::__mod<__Type, 1U, 0U, - __detail::_Shift<__Type, 32>::__value>(__r2); - __begin[(__k + __p) % __n] += __r1; - __begin[(__k + __q) % __n] += __r2; - __begin[__k % __n] = __r2; + uint32_t __r1 = 1371501266u; + uint32_t __r2 = __r1 + __s; + __begin[__p] += __r1; + __begin[__q] = (uint32_t)__begin[__q] + __r2; + __begin[0] = __r2; + } + + for (size_t __k = 1; __k <= __s; ++__k) + { + const size_t __kn = __k % __n; + const size_t __kpn = (__k + __p) % __n; + const size_t __kqn = (__k + __q) % __n; + uint32_t __arg = (__begin[__kn] + ^ __begin[__kpn] + ^ __begin[(__k - 1) % __n]); + uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27)); + uint32_t __r2 = __r1 + (uint32_t)__kn + _M_v[__k - 1]; + __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1; + __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2; + __begin[__kn] = __r2; + } + + for (size_t __k = __s + 1; __k < __m; ++__k) + { + const size_t __kn = __k % __n; + const size_t __kpn = (__k + __p) % __n; + const size_t __kqn = (__k + __q) % __n; + uint32_t __arg = (__begin[__kn] + ^ __begin[__kpn] + ^ __begin[(__k - 1) % __n]); + uint32_t __r1 = 1664525u * (__arg ^ (__arg >> 27)); + uint32_t __r2 = __r1 + (uint32_t)__kn; + __begin[__kpn] = (uint32_t)__begin[__kpn] + __r1; + __begin[__kqn] = (uint32_t)__begin[__kqn] + __r2; + __begin[__kn] = __r2; } for (size_t __k = __m; __k < __m + __n; ++__k) { - __Type __arg = __begin[__k % __n] - + __begin[(__k + __p) % __n] - + __begin[(__k - 1) % __n]; - __Type __r3 = __arg ^ (__arg << 27); - __r3 = __detail::__mod<__Type, 1566083941U, 0U, - __detail::_Shift<__Type, 32>::__value>(__r3); - __Type __r4 = __r3 - __k % __n; - __r4 = __detail::__mod<__Type, 1U, 0U, - __detail::_Shift<__Type, 32>::__value>(__r4); - __begin[(__k + __p) % __n] ^= __r4; - __begin[(__k + __q) % __n] ^= __r3; - __begin[__k % __n] = __r4; + const size_t __kn = __k % __n; + const size_t __kpn = (__k + __p) % __n; + const size_t __kqn = (__k + __q) % __n; + uint32_t __arg = (__begin[__kn] + + __begin[__kpn] + + __begin[(__k - 1) % __n]); + uint32_t __r3 = 1566083941u * (__arg ^ (__arg >> 27)); + uint32_t __r4 = __r3 - __kn; + __begin[__kpn] ^= __r3; + __begin[__kqn] ^= __r4; + __begin[__kn] = __r4; } } @@ -2710,20 +3345,39 @@ namespace std _RealType generate_canonical(_UniformRandomNumberGenerator& __urng) { + static_assert(std::is_floating_point<_RealType>::value, + "template argument must be a floating point type"); + const size_t __b = std::min(static_cast(std::numeric_limits<_RealType>::digits), __bits); const long double __r = static_cast(__urng.max()) - static_cast(__urng.min()) + 1.0L; - const size_t __log2r = std::log10(__r) / std::log10(2.0L); - size_t __k = std::max(1UL, (__b + __log2r - 1UL) / __log2r); + const size_t __log2r = std::log(__r) / std::log(2.0L); + const size_t __m = std::max(1UL, + (__b + __log2r - 1UL) / __log2r); + _RealType __ret; _RealType __sum = _RealType(0); _RealType __tmp = _RealType(1); - for (; __k != 0; --__k) + for (size_t __k = __m; __k != 0; --__k) { __sum += _RealType(__urng() - __urng.min()) * __tmp; __tmp *= __r; } - return __sum / __tmp; + __ret = __sum / __tmp; + if (__builtin_expect(__ret >= _RealType(1), 0)) + { +#if _GLIBCXX_USE_C99_MATH_TR1 + __ret = std::nextafter(_RealType(1), _RealType(0)); +#else + __ret = _RealType(1) + - std::numeric_limits<_RealType>::epsilon() / _RealType(2); +#endif + } + return __ret; } -} + +_GLIBCXX_END_NAMESPACE_VERSION +} // namespace + +#endif