Changes in v5: * Static-assert movable RNG object correctly. * Add a more comprehensive test gencanon.cc
Changes in v4: * Static-assert that arg is floating-point, coercible from bigger unsigned. * Static-assert that arg satisfies uniform_random_bit_generator, movable. * #include uniform_int_dist.h for concept uniform_random_bit_generator * Coerce floating consts from unsigned literals, matching other usage. * Remove redundant type decorations static and -> . * Use __builtin_floor for broader applicability than std::floor. * Use __builtin_pow in place of immediate lambda. * Rename, rearrange local objects for clarity. * Add redundant "if constexpr (k > 1)" for clarity. * #if-guard new impl against Clang, which lacks constexpr floor() etc. Changes in v3: * Implement exactly as specified in the WP, for C++26 and up. Implement P0952R2 "A new specification for std::generate_canonical". It has us start over if the naive generation of a float in 0..1 comes up exactly equal to 1, which occurs too rarely for the change to measurably affect performance. A test for the case already appears in 64351.cc. This change amounts to a different resolution for PR64351 and LWG 2524. A new test, gencanon.cc, is more thorough. libstdc++-v3/ChangeLog: PR libstdc++/119739 * include/bits/random.tcc: Add generate_canonical impl for C++26. * include/std/random: Include uniform_int_dist.h for concept. * testsuite/26_numerics/random/uniform_real_distribution/operators/64351.cc: Adapt test for both pre- and post- C++26. * testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon.cc: Test for float and double from 32-bit RNG, including 1.0 do-over. --- libstdc++-v3/include/bits/random.tcc | 48 +++++++++++ libstdc++-v3/include/std/random | 1 + .../operators/64351.cc | 8 +- .../operators/gencanon.cc | 85 +++++++++++++++++++ 4 files changed, 140 insertions(+), 2 deletions(-) create mode 100644 libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon.cc diff --git a/libstdc++-v3/include/bits/random.tcc b/libstdc++-v3/include/bits/random.tcc index 53ccacb2e38..65f457a6659 100644 --- a/libstdc++-v3/include/bits/random.tcc +++ b/libstdc++-v3/include/bits/random.tcc @@ -3349,6 +3349,9 @@ namespace __detail } } +#if __cplusplus < 202400 || defined(__clang__) + // Clang lacks constexpr __builtin_floor and __builtin_pow. + template<typename _RealType, size_t __bits, typename _UniformRandomNumberGenerator> _RealType @@ -3385,6 +3388,51 @@ namespace __detail } return __ret; } +#else + template<typename _RealT, size_t __digits, typename _URBG> + _RealT + generate_canonical(_URBG& __urng) + { + static_assert(requires { + [](uniform_random_bit_generator auto) {}( + static_cast<std::remove_cv_t<_URBG>&&>(__urng)); }, + "argument must meet uniform_random_bit_generator requirements"); + static_assert( + std::is_floating_point_v<_RealT> && requires { _RealT(_URBG::max()); }, + "template argument must be floating point, coercible from unsigned"); + + // Names are chosen to match the description in the Standard. + const size_t __r = std::numeric_limits<_RealT>::radix; + const auto __sample = [](auto __high) constexpr + { return _RealT(__high - _URBG::min()); }; + constexpr _RealT __R = _RealT(1) + __sample(_URBG::max()); + const size_t __d_max = std::numeric_limits<_RealT>::digits; + const size_t __d = __digits < __d_max ? __digits : __d_max; + constexpr _RealT __rd = __builtin_pow(_RealT(__r), __d); + const unsigned __k = [](_RealT __R, _RealT __rd) consteval + { + unsigned __i = 1; + for (auto __Ri = __R; __Ri < __rd; __Ri *= __R) + ++__i; + return __i; + } (__R, __rd); + constexpr _RealT __Rk = __builtin_pow(__R, __k); + constexpr _RealT __x = __builtin_floor(__Rk / __rd); + constexpr _RealT __xrd = __x * __rd; + + _RealT __sum; + do + { + _RealT __Ri = _RealT(1); + __sum = __sample(__urng()); + if constexpr (__k > 1) + for (int __i = 1; __i != __k; ++__i) + __Ri *= __R, __sum += __sample(__urng()) * __Ri; + } while (__builtin_expect(__sum >= __xrd, false)); + _RealT __ret = __builtin_floor(__sum / __x) / __rd; + return __ret; + } +#endif _GLIBCXX_END_NAMESPACE_VERSION } // namespace diff --git a/libstdc++-v3/include/std/random b/libstdc++-v3/include/std/random index 0e058a04bd9..86d4f0a8380 100644 --- a/libstdc++-v3/include/std/random +++ b/libstdc++-v3/include/std/random @@ -49,6 +49,7 @@ #include <type_traits> #include <bits/random.h> #include <bits/opt_random.h> +#include <bits/uniform_int_dist.h> #include <bits/random.tcc> #endif // C++11 diff --git a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/operators/64351.cc b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/operators/64351.cc index 3c0cb8bbd49..6212a32a364 100644 --- a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/operators/64351.cc +++ b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/operators/64351.cc @@ -45,16 +45,20 @@ test02() rng.seed(sequence); rng.discard(12 * 629143); std::mt19937 rng2{rng}; + int mismatch = 0; for (int i = 0; i < 20; ++i) { float n = std::generate_canonical<float, std::numeric_limits<float>::digits>(rng); VERIFY( n != 1.f ); - // PR libstdc++/80137 rng2.discard(1); - VERIFY( rng == rng2 ); } + // PR libstdc++/80137 +#if __cplusplus >= 202400 + rng2.discard(1); // account for a 1.0 generated and discarded. +#endif + VERIFY( rng == rng2 ); } int diff --git a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon.cc b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon.cc new file mode 100644 index 00000000000..ba781315e0d --- /dev/null +++ b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon.cc @@ -0,0 +1,85 @@ +// { dg-do run { target { c++26 && { ! simulator } } } } + +#include <random> +#include <limits> +#include <type_traits> +#include <cmath> +#include <testsuite_hooks.h> + +// Verify P0952R9 implementation requiring a second round-trip +// if first yields exactly 1. In this test, the RNG delivering +// 32 bits per call is seeded such that this occurs once on the +// sixth iteration for float, and not at all for double. +// However, each double iteration requires two calls to the RNG. + +template <typename T, typename RNG> +void +test01(RNG& rng, int& deviation, int& max, int& rms, int& zero, int& skips) +{ + const auto size = 1'000'000, buckets = 100; + std::array<number, buckets> histo{}; + int zero = 0, skips = 0; + for (auto i = 0; i != size; ++i) { + T sample = std::generate_canonical<T, 100>(rng); + VERIFY(sample >= 0.0); + VERIFY(sample < 1.0) + if (sample == 0.0) { + ++zero; + } + auto bucket = static_cast<int>(std::floor(sample * buckets)); + ++histo[bucket]; + rng2.discard(1); + if (rng != rng2) { + ++skips; + rng2.discard(1); + VERIFY(rng == rng2) + } + } + int devsquare = 0; + for (int i = 0; i < buckets; ++i) { + const auto expected = size / buckets; + auto count = histo[i]; + auto diff = count - expected; + if (diff < 0) diff = -diff; + deviation += diff; + devsquare += diff * diff; + if (diff > max) max = diff; + } + rms = std::sqrt(devsquare); +} + +int +main() +{ + std::mt19937 rng(8890); + std::seed_seq sequence{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}; + rng.seed(sequence); + rng.discard(12 * 629143); + + { + int deviation{}, max{}, rms{}, zero{}, skips{}; + auto rng2{rng}; + test01<float>(rng2, deviation, max, rms, zero, skips); + + if (std::numeric_limits<float>::is_iec559) { + VERIFY(deviation == 7032); + VERIFY(max == 276); + VERIFY(rms == 906); + VERIFY(zero == 0); + } + VERIFY(skips == 1); + } + { + auto rng3{rng} + test01<double>(rng2, deviation, max, rms, zero, skips); + + if (std::numeric_limits<double>::is_iec559) { + VERIFY(deviation == 7650); + VERIFY(max == 259); + VERIFY(rms == 975); + VERIFY(zero == 0); + } + VERIFY(skips == 1'000'000); + } + return 0; +} -- 2.50.0