https://gcc.gnu.org/g:72430fff7b40c9c1a156d36b3d734ebb2e850166
commit r16-6280-g72430fff7b40c9c1a156d36b3d734ebb2e850166 Author: Tomasz Kamiński <[email protected]> Date: Thu Dec 18 14:14:54 2025 +0100 libstdc++: Use smallest possible integer for __generate_cannonical_any If the span of the range R produced by uniform bit generator U passed to generate_canonical is not power of two, we need to use algorithm that requires computing power R^k that is greater than 2^d, where d is number of digits in mantissa of _RealT. Previously we have used an integer type that is has twice as many digits as d. This lead to situation that for standard engines that produced such range (like std::minstd_rand0, std::minstd_rand, std::ranlux24, ....) 256bit integer support was required for 128bit floats. However, in this cases R^4 provides more than d bits of precision, while requiring 124 bits. We overestimate the number of required bits, by computing a value l * bit_width(R) (log2(R) + 1), where l is value such that log2(R) * l >= d. As R >= 2^log2(R), then R^l >= (2^log2(R))^l == 2^(log(R) * l) >= 2^d, so k+1 >= l >= k. In consequence R^k is smaller R^l which require at most l * bit_width(R). This is an overestimate, but difference should not be higher than l bits. We replace __gen_can_pow and __gen_can_rng_calls_needed with __gen_canon_log(v, b), which computes the largest power of b that fits into v. As such a number is smaller than v, the result will always fit in it's type. Both the logarithm and the power value are returned using __gen_canon_log_res struct. libstdc++-v3/ChangeLog: * include/bits/random.h (__rand_uint128::operator>) (__rand_uint128::operator>=): Define. * include/bits/random.tcc (__generate_canonical_pow2): Adjust for use of __rand_uint128 in C++11. (__gen_can_pow, __gen_can_rng_calls_needed): Replace with __gen_canon_log. (__gen_canon_log_res, __gen_canon_log): Define. (__generate_canonical_any): Reworked how _UInt is determined. * testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon_eng.cc: New test. Reviewed-by: Jonathan Wakely <[email protected]> Signed-off-by: Tomasz Kamiński <[email protected]> Diff: --- libstdc++-v3/include/bits/random.h | 8 ++ libstdc++-v3/include/bits/random.tcc | 104 +++++++++++++++------ .../operators/gencanon_eng.cc | 39 ++++++++ 3 files changed, 121 insertions(+), 30 deletions(-) diff --git a/libstdc++-v3/include/bits/random.h b/libstdc++-v3/include/bits/random.h index 7a9ba0a7a9e9..117d63dff7f6 100644 --- a/libstdc++-v3/include/bits/random.h +++ b/libstdc++-v3/include/bits/random.h @@ -463,9 +463,17 @@ _GLIBCXX_END_INLINE_ABI_NAMESPACE(_V2) return false; } + friend _GLIBCXX14_CONSTEXPR bool + operator>(const type& __l, const type& __r) noexcept + { return __r < __l; } + friend _GLIBCXX14_CONSTEXPR bool operator<=(const type& __l, const type& __r) noexcept { return !(__r < __l); } + + friend _GLIBCXX14_CONSTEXPR bool + operator>=(const type& __l, const type& __r) noexcept + { return !(__l < __r); } #endif friend _GLIBCXX14_CONSTEXPR bool diff --git a/libstdc++-v3/include/bits/random.tcc b/libstdc++-v3/include/bits/random.tcc index 0256059ab72c..66ae14a06d59 100644 --- a/libstdc++-v3/include/bits/random.tcc +++ b/libstdc++-v3/include/bits/random.tcc @@ -3532,16 +3532,16 @@ namespace __detail const unsigned __log2_Rk_max = __k * __log2_R; const unsigned __log2_Rk = // Bits of entropy actually obtained: __log2_uint_max < __log2_Rk_max ? __log2_uint_max : __log2_Rk_max; - static_assert(__log2_Rk >= __d); // Rk = _UInt(1) << __log2_Rk; // Likely overflows, UB. - constexpr _RealT __Rk = _RealT(_UInt(1) << (__log2_Rk - 1)) * _RealT(2.0); + _GLIBCXX14_CONSTEXPR const _RealT __Rk + = _RealT(_UInt(1) << (__log2_Rk - 1)) * _RealT(2.0); #if defined(_GLIBCXX_GENERATE_CANONICAL_STRICT) const unsigned __log2_x = __log2_Rk - __d; // # of spare entropy bits. #else const unsigned __log2_x = 0; #endif - const _UInt __x = _UInt(1) << __log2_x; - constexpr _RealT __rd = __Rk / __x; + _GLIBCXX14_CONSTEXPR const _UInt __x = _UInt(1) << __log2_x; + _GLIBCXX14_CONSTEXPR const _RealT __rd = __Rk / _RealT(__x); // xrd = __x << __d; // Could overflow. while (true) @@ -3552,29 +3552,50 @@ namespace __detail __shift += __log2_R; __sum |= _UInt(__urng() - _Urbg::min()) << __shift; } - const _RealT __ret = _RealT(__sum >> __log2_x) / __rd; + const _RealT __ret = _RealT(__sum >> __log2_x) / _RealT(__rd); if (__ret < _RealT(1.0)) return __ret; } } - template <typename _UInt> - constexpr _UInt - __gen_can_pow(_UInt __base, size_t __exponent) + + template<typename _UInt> + struct __gen_canon_log_res { - _UInt __prod = __base; - while (--__exponent != 0) __prod *= __base; - return __prod; - } + unsigned __floor_log; + _UInt __floor_pow; + + constexpr __gen_canon_log_res + update(_UInt __base) const + { return {__floor_log + 1, __floor_pow * __base}; } + }; - template <typename _UInt> - constexpr unsigned - __gen_can_rng_calls_needed(_UInt __R, _UInt __rd) + + template <typename _UInt1, typename _UInt2, + typename _UComm = __conditional_t<(sizeof(_UInt2) > sizeof(_UInt1)), + _UInt2, _UInt1>> + constexpr __gen_canon_log_res<_UInt1> + __gen_canon_log(_UInt1 __val, _UInt2 __base) { - unsigned __i = 1; - for (auto __Ri = __R; __Ri < __rd; __Ri *= __R) - ++__i; - return __i; +#if __cplusplus >= 201402L + __gen_canon_log_res<_UInt1> __res{0, _UInt1(1)}; + if (_UComm(__base) > _UComm(__val)) + return __res; + + const _UInt1 __base1(__base); + do + { + __val /= __base1; + __res = __res.update(__base1); + } + while (__val >= __base1); + return __res; +#else + return (_UComm(__val) >= _UComm(__base)) + ? __gen_canon_log(__val / _UInt1(__base), _UInt1(__base)) + .update(_UInt1(__base)) + : __gen_canon_log_res<_UInt1>{0, _UInt1(1)}; +#endif } // This version must be used when the range of possible RNG results, @@ -3587,30 +3608,53 @@ namespace __detail _RealT __generate_canonical_any(_Urbg& __urng) { - using _UInt = typename __detail::_Select_uint_least_t<__d * 2>::type; - // Names below are chosen to match the description in the Standard. // Parameter d is the actual target number of bits. - const _UInt __r{2}; - const _UInt __R = _UInt(_Urbg::max() - _Urbg::min()) + 1; - const _UInt __rd = _UInt(1) << __d; - const unsigned __k = __gen_can_rng_calls_needed(__R, __rd); - const _UInt __Rk = __gen_can_pow(__R, __k); - const _UInt __x = __Rk / __rd; +#if (__cplusplus >= 201402L) || defined(__SIZEOF_INT128__) +# define _GLIBCXX_GEN_CANON_CONST constexpr +#else +# define _GLIBCXX_GEN_CANON_CONST const +#endif + + using _UIntR = typename make_unsigned<decltype(_Urbg::max())>::type; + // Cannot overflow, as _Urbg::max() - _Urbg::min() is not power of + // two minus one + constexpr _UIntR __R = _UIntR(_Urbg::max() - _Urbg::min()) + 1; + constexpr unsigned __log2R + = sizeof(_UIntR) * __CHAR_BIT__ - __builtin_clzg(__R) - 1; + // We overstimate number of required bits, by computing + // r such that l * log2(R) >= d, so: + // R^l >= (2 ^ log2(R)) ^ l == 2 ^ (log2(r) * l) >= 2^d + // And then requiring l * bit_width(R) bits. + constexpr unsigned __l = (__d + __log2R - 1) / __log2R; + constexpr unsigned __bits = (__log2R + 1) * __l; + using _UInt = typename __detail::_Select_uint_least_t<__bits>::type; + + _GLIBCXX_GEN_CANON_CONST _UInt __rd = _UInt(1) << __d; + _GLIBCXX_GEN_CANON_CONST auto __logRrd = __gen_canon_log(__rd, __R); + _GLIBCXX_GEN_CANON_CONST unsigned __k + = __logRrd.__floor_log + (__rd > __logRrd.__floor_pow); + + _GLIBCXX_GEN_CANON_CONST _UInt __Rk + = (__k > __logRrd.__floor_log) + ? _UInt(__logRrd.__floor_pow) * _UInt(__R) + : _UInt(__logRrd.__floor_pow); + _GLIBCXX_GEN_CANON_CONST _UInt __x = __Rk / __rd; while (true) { _UInt __Ri{1}; - _UInt __sum{__urng() - _Urbg::min()}; + _UInt __sum(__urng() - _Urbg::min()); for (int __i = __k - 1; __i > 0; --__i) { - __Ri *= __R; - __sum += _UInt{__urng() - _Urbg::min()} * __Ri; + __Ri *= _UInt(__R); + __sum += _UInt(__urng() - _Urbg::min()) * __Ri; } const _RealT __ret = _RealT(__sum / __x) / _RealT(__rd); if (__ret < _RealT(1.0)) return __ret; } +#undef _GLIBCXX_GEN_CANON_CONST } #if !defined(_GLIBCXX_GENERATE_CANONICAL_STRICT) diff --git a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon_eng.cc b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon_eng.cc new file mode 100644 index 000000000000..c8d5b5d263be --- /dev/null +++ b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon_eng.cc @@ -0,0 +1,39 @@ +// { dg-do compile { target { c++11 } } } + +#include <random> + +template<typename _Real, typename _URBG> +void +test_engine() +{ + _URBG __engine; + (void)std::generate_canonical<_Real, size_t(-1)>(__engine); +} + +template<typename _Real> +void +test_all_engines() +{ + test_engine<_Real, std::default_random_engine>(); + + test_engine<_Real, std::minstd_rand0>(); + test_engine<_Real, std::minstd_rand>(); + test_engine<_Real, std::mt19937>(); + test_engine<_Real, std::mt19937_64>(); + test_engine<_Real, std::ranlux24_base>(); + test_engine<_Real, std::ranlux48_base>(); + test_engine<_Real, std::ranlux24>(); + test_engine<_Real, std::ranlux48>(); + test_engine<_Real, std::knuth_b>(); +#if __cplusplus > 202302L + test_engine<_Real, std::philox4x32>(); + test_engine<_Real, std::philox4x64>(); +#endif +} + +int main() +{ + test_all_engines<float>(); + test_all_engines<double>(); + test_all_engines<long double>(); +}
