[PATCH] libstdc++: Adjust generate_canonical (P0952, LWG2524) [PR119739]

[Nathan Myers]

Changes in V3:
 * Implement exactly as specified in the WP, for C++26 and up.
 * Eliminate massive overhead seen in present implementation.

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.

libstdc++-v3/ChangeLog:
	PR libstdc++/119739
	* include/bits/random.tcc
---
 libstdc++-v3/include/bits/random.tcc | 50 ++++++++++++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/libstdc++-v3/include/bits/random.tcc b/libstdc++-v3/include/bits/random.tcc
index 53ccacb2e38..941495654c0 100644
--- a/libstdc++-v3/include/bits/random.tcc
+++ b/libstdc++-v3/include/bits/random.tcc
@@ -3349,6 +3349,7 @@ namespace __detail
 	}
     }
 
+#if __cplusplus < 202300
   template<typename _RealType, size_t __bits,
 	   typename _UniformRandomNumberGenerator>
     _RealType
@@ -3385,6 +3386,55 @@ namespace __detail
 	}
       return __ret;
     }
+#else
+  template<typename _RealT, size_t __digits, typename _URBG>
+    _RealT
+    generate_canonical(_URBG& __urng)
+    {
+      static_assert(std::is_floating_point<_RealT>::value,
+	"template argument must be a floating point type");
+      static_assert(
+	requires { [](uniform_random_bit_generator auto) {}(__urng); },
+	"template argument must be an std::uniform_random_bit_generator");
+
+      const size_t __r = std::numeric_limits<_RealT>::radix;
+      const auto __rdiff = [](auto __u) constexpr -> _RealT
+	{ return _RealT(__u - _URBG::min()); };
+      constexpr static _RealT __R = __rdiff(_URBG::max()) + _RealT(1.0);
+      const size_t __reg_digits = std::numeric_limits<_RealT>::digits;
+      const size_t __d = __reg_digits < __digits ? __reg_digits : __digits;
+      const auto __pow = [](_RealT __r, auto __d) consteval -> _RealT
+	{
+	  _RealT __res = __r;
+	  while (--__d) __res *= __r;
+	  return __res;
+	};
+      constexpr static _RealT __rd = __pow(_RealT(__r), __d);
+      const int __k = [](_RealT __rd, _RealT __R) consteval
+	{
+	  unsigned __k = 1;
+	  _RealT __Ri = __R;
+	  while (__Ri < __rd)
+	    { __Ri *= __R, ++__k; }
+	  return __k;
+	} (__rd, __R);
+      constexpr static _RealT __Rk = __pow(__R, __k);
+      constexpr static  _RealT __x = [](_RealT __Rk, _RealT __rd) consteval
+	{ return std::floor(__Rk/__rd); }(__Rk, __rd);
+      constexpr static _RealT __xrd = __x * __rd;
+
+      _RealT __sum;
+      do
+	{
+	  _RealT __Ri = _RealT(1.0);
+	  __sum = __rdiff(__urng());
+	  for (int __i = 1; __i != __k; ++__i)
+	    __Ri *= __R, __sum += __rdiff(__urng()) * __Ri;
+	} while (__builtin_expect(__sum >= __xrd, false));
+      _RealT __ret = std::floor(__sum / __x) / __rd;
+      return __ret;
+    }
+#endif
 
 _GLIBCXX_END_NAMESPACE_VERSION
 } // namespace
-- 
2.50.0