Re: [PATCH v4 2/2] libstdc++: Support any generator (up to 64bit) in generate_canonical [PR119739]
Jonathan Wakely <[email protected]> Tue, 4 Aug 2026 16:49:25 +0100
| Newsgroups | gmane.comp.gcc.patches,gmane.comp.gcc.libstdc++.devel |
|---|---|
| Message-ID | <CAH6eHdRDtx75aU6TGm=eMj2KfHaD9SKBCJqRbXyWzfeA0uCmgw@mail.gmail.com> |
--0000000000001fdc5106583a9ab8 Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable On Tue, 4 Aug 2026, 12:11 Tomasz Kami=C5=84ski, <[email protected]> wrote= : > This patch splits implementation of generate_canonical_any, into three > different path depending on bits generated (R), required bits of entropy > (b), and value k such that R^k > 2^d > R^(k-1): > * log2(R) >=3D d: single call of generator is sufficient (k =3D 1), > and computation can be performed in the generator value type (_UIntR); > * bit_width(R^k) <=3D 128: computation can be performed using type > __detail::_Select_uint_least_t; previous implementation is reused; > * bit_width(R^k) > 128: the computations required may overflow 128 bit; > such generators were previously unsupported, as illustrated by > r17-507-g5d17eec7d83445. > > This patch implements the support for the last branch, by observing > that while it is possible for following computation to overflow 128bit > integer: > sum_[0, k) urng() * R^i > The values of sum, such that (sum / x * 2^d) > 1, are rejected, thus > any value (sum / x) > 2^d can be rejected. As 2^d is guaranteed to fit > in 128 integer, we can use above to reject the value, before overflow > will occur. > > Firstly, we know that for l =3D k-1, the value R^l < 2^d, and > sum_[0, l) urng() * R^i <=3D R^l (because urng() is always < R), > always fit in 128bit integer. We compute that value (__lsum) using > normal loop. > > This leaves us with final iteration of multiplying __kth * R^l, where > __kth is result of invocation of urng. In that case we compute the > bitwidth (__abits) of value that can be multiplied by _R^l without > overflow 128bit integer. The __kth * R^l is performed by multiplying > the loop of __abits of __kth in loop where: > * __sumx (initially __lsum / __x) is the current value of sum / x, > * __remx (initially 2^d - __sumx) is the value that we can still, > add before overflowing (shifted with each iteration) > * __modx (initially __lsum % __x) are collected reminders. > In each iteration the lowest __kth bits are multiplied with _Rl (__val), > scaled (__valx) and checked against __rem. When fits the __sumx, __modx > and __remx are updated, and both values of __remx and __kth are shifted > right by __abits. > > If the multiplication finished without overflow (__kth is zero), we handl= e > the collected remainders (__modx), and compute the __ret value. The value= s > greater or equal 1 (e.g. due floating point approximation), are still > rejected. > > Similar approach is applied when computing x =3D R^k / 2^d. > Again while R^k will overflow 128 bits, we know that x < R (so fits > in 64bits), as R^l < 2^d < R^k, then x < R. We again compute > R * R^l / 2^d, by splitting the value of R into __abits chunks (_pR). > For result of chunks multipliations (_pRk), the bits (__xp) remaining > after division by __shift (starts from __d and reduces by __abits), > are collected into cumulative value of __x. While the lower bits > are preserved for later iteration (__pRk). > > PR libstdc++/119739 > > libstdc++-v3/ChangeLog: > > * include/bits/random.tcc (__generate_canonical_any): > Rewrite to support generator with different ranges. > * > testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon= .cc: > Add test for previously unsupported generator. > * > testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon= _eng.cc: > Add test for previously unsupported-generator ranges. > * > testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon= _eng_neg.cc: > Moved (now-supported) generator test to gencanon_eng.cc, and > removed > the file. > --- > v4: > * address typos in comments pointed out in the review > * address TODO, initializing __bwRl with actual bitwidth > (we need to handle __rand_uint128 for 32bit targets, without __int128) > * add test for 2^7 + 1 =3D 129 generator size in the gencanon_eng > * adjust gen_canon_eng numbers. > > The issues reported by LLM was real, while the bit_width(_Rl) is > guaranteed to > be less than 113, because _Rl < 2^113, the aproximation of bit_with(_R) * > l, > is not due overestimate. This addresses above. > > We should consider making __bitwitdh and other function from <bit> functi= on > available in C++11 (they are from 14 now) and providing overloads for > __rand_uint128, but this is best done as separate commit. > I don't remember why I made them C++14 only, maybe because they can't be constexpr in C++11 due to multiple statements. But that's ok, they can just be non-constexpr for C++11. > Tested on x86_64-linux locally. *gencanon* additionally tested with all > standard modes, -m32, debug and assertions and on powerpc64le. > OK for trunk? > OK, thanks > libstdc++-v3/include/bits/random.tcc | 173 +++++++++++++++--- > .../operators/gencanon.cc | 153 ++++++++++++++++ > .../operators/gencanon_eng.cc | 76 +++++++- > .../operators/gencanon_eng_neg.cc | 89 --------- > 4 files changed, 371 insertions(+), 120 deletions(-) > delete mode 100644 > libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/opera= tors/gencanon_eng_neg.cc > > diff --git a/libstdc++-v3/include/bits/random.tcc > b/libstdc++-v3/include/bits/random.tcc > index aa176db4778..b2afe411a7c 100644 > --- a/libstdc++-v3/include/bits/random.tcc > +++ b/libstdc++-v3/include/bits/random.tcc > @@ -242,7 +242,7 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION > constexpr size_t > mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d= , > __s, __b, __t, __c, __l, __f>::tempering_u; > - > + > template<typename _UIntType, > size_t __w, size_t __n, size_t __m, size_t __r, > _UIntType __a, size_t __u, _UIntType __d, size_t __s, > @@ -3807,38 +3807,155 @@ namespace __detail > constexpr unsigned __log2R > =3D sizeof(_UIntR) * __CHAR_BIT__ - __builtin_clzg(__R) - 1; > // We overstimate number of required bits, by computing > - // r such that l * log2(R) >=3D d, so: > - // R^l >=3D (2 ^ log2(R)) ^ l =3D=3D 2 ^ (log2(r) * l) >=3D 2^d > - // And then requiring l * bit_width(R) bits. > - constexpr unsigned __l =3D (__d + __log2R - 1) / __log2R; > - constexpr unsigned __bits =3D (__log2R + 1) * __l; > - using _UInt =3D typename __detail::_Select_uint_least_t<__bits>::t= ype; > - > - _GLIBCXX_GEN_CANON_CONST _UInt __rd =3D _UInt(1) << __d; > - _GLIBCXX_GEN_CANON_CONST auto __logRrd =3D __gen_canon_log(__rd, _= _R); > - _GLIBCXX_GEN_CANON_CONST unsigned __k > - =3D __logRrd.__floor_log + (__rd > __logRrd.__floor_pow); > - > - _GLIBCXX_GEN_CANON_CONST _UInt __Rk > - =3D (__k > __logRrd.__floor_log) > - ? _UInt(__logRrd.__floor_pow) * _UInt(__R) > - : _UInt(__logRrd.__floor_pow); > - _GLIBCXX_GEN_CANON_CONST _UInt __x =3D __Rk / __rd; > + // m such that m * log2(R) >=3D d, so: > + // R^m >=3D (2 ^ log2(R)) ^ m =3D=3D 2 ^ (log2(R) * m) >=3D 2^d > + // And then requiring m * bit_width(R) bits. > + constexpr unsigned __m =3D (__d + __log2R - 1) / __log2R; > + constexpr unsigned __bits =3D (__log2R + 1) * __m; > + if constexpr (__log2R >=3D __d) > + { > + // range already provide required number of bits, > + // so in this case __k =3D=3D 1, and single call to generator > + // is sufficient. Furthermore 2^d fits in _UIntR > + constexpr _UIntR __rd =3D _UIntR(1) << __d; > + constexpr _UIntR __x =3D __R >> __d; > > - while (true) > + while (true) > + { > + _UIntR __val(__urng() - _Urbg::min()); > + const _RealT __ret =3D _RealT(__val / __x) / _RealT(__rd); > + if (__ret < _RealT(1.0)) > + return __ret; > + } > + } > + // generator call produce fewer bits than __d, but R^k > rd fits > + // into 128 bits. > + else if constexpr (__bits <=3D 128) > { > - _UInt __Ri{1}; > - _UInt __sum(__urng() - _Urbg::min()); > - for (int __i =3D __k - 1; __i > 0; --__i) > + using _UInt =3D typename > __detail::_Select_uint_least_t<__bits>::type; > + > + _GLIBCXX_GEN_CANON_CONST _UInt __rd =3D _UInt(1) << __d; > + _GLIBCXX_GEN_CANON_CONST auto __logRrd =3D __gen_canon_log(__rd= , > __R); > + // __rd is power of two, and __R is not so __floor_pow is never > + // equal to __rd. > + _GLIBCXX_GEN_CANON_CONST unsigned __k =3D __logRrd.__floor_log = + 1; > + _GLIBCXX_GEN_CANON_CONST _UInt __Rk =3D __logRrd.__floor_pow * > _UInt(__R); > + // x =3D floor(R^k / 2^d) =3D R^k >> d; > + _GLIBCXX_GEN_CANON_CONST _UInt __x =3D __Rk >> __d; > + > + while (true) > { > - __Ri *=3D _UInt(__R); > - __sum +=3D _UInt(__urng() - _Urbg::min()) * __Ri; > + _UInt __Ri{1}; > + _UInt __sum(__urng() - _Urbg::min()); > + for (int __i =3D __k - 1; __i > 0; --__i) > + { > + __Ri *=3D __R; > + __sum +=3D __Ri * _UIntR(__urng() - _Urbg::min()); > + } > + const _RealT __ret =3D _RealT(__sum / __x) / _RealT(__rd); > + if (__ret < _RealT(1.0)) > + return __ret; > } > - const _RealT __ret =3D _RealT(__sum / __x) / _RealT(__rd); > - if (__ret < _RealT(1.0)) > - return __ret; > } > -#undef _GLIBCXX_GEN_CANON_CONST > + else > + { > + static_assert(__log2R < 64, "Only generators emitting up to 64 > bits are supported"); > + using _UInt =3D typename __detail::_Select_uint_least_t<128>::t= ype; > + using _UInt64 =3D typename > __detail::_Select_uint_least_t<64>::type; > + > + _GLIBCXX_GEN_CANON_CONST _UInt __rd =3D _UInt(1) << __d; > + _GLIBCXX_GEN_CANON_CONST auto __logRrd =3D __gen_canon_log(__rd= , > __R); > + // __rd is power of two, and __R is not so __floor_pow is never > + // equal to __rd, and l =3D=3D k - 1 > + _GLIBCXX_GEN_CANON_CONST unsigned __l =3D __logRrd.__floor_log; > + _GLIBCXX_GEN_CANON_CONST _UInt __Rl =3D __logRrd.__floor_pow; > + > + // __abits is maximum bit width of the value, that can > + // be multiplied by R^l without overflowing 128 bit integer > + _GLIBCXX_GEN_CANON_CONST unsigned __bwRl > +#ifndef __SIZEOF_INT128__ > + =3D __Rl._M_hi ? 128 - __builtin_clzg(__Rl._M_hi) > + : 64 - __builtin_clzg(__Rl._M_lo); > +#else > + =3D 128 - __builtin_clzg(__Rl); > +#endif > + > + // For __R close to power of two, the actual __k may be smaller > than __m, > + // and will use less than 128bits, default to two 32 bits chunk= s. > + _GLIBCXX_GEN_CANON_CONST unsigned __abits =3D __bwRl > 96 ? 128= - > __bwRl : 32; > + _GLIBCXX_GEN_CANON_CONST _UInt64 __amask =3D (_UInt64(1) << > __abits) - 1; > + > + // Compute __x =3D _Rk / __rd (_Rk >> __d), by spliting _R into > chunks > + // (_pR) of __abits, so their multiplication does not overflow > 128 bits. > + // __x fits in 64bits, becuse __R^l < __rd < __R^k, thus __x < = _R > + _UInt64 __x(0), __pR(__R); _UInt __pRk(0); > + for (unsigned __shift =3D __d; __pR; __shift -=3D __abits) > + { > + __pRk +=3D __Rl * (__pR & __amask); > + __pR >>=3D __abits; > + > + _UInt __xp(__pRk >> __shift); > + __x +=3D _UInt64(__xp); > + __pRk -=3D (__xp << __shift); > + __pRk >>=3D __abits; > + } > + > + while (true) > + { > + // Compute sum of __urng() * R^i for i in range [0, l). > + // The value is smaller than R^l which is smaller than > + // __rd so fits into 128 bit integer. > + _UInt __Ri{1}; > + _UInt __lsum(__urng() - _Urbg::min()); > + for (int __i =3D __l - 1; __i > 0; --__i) > + { > + __Ri *=3D _UInt64(__R); > + __lsum +=3D __Ri * _UInt64(__urng() - _Urbg::min()); > + } > + > + // The last iteration __urng() * R^k may overflow 128bit > + // integer. We use the fact that we reject __sum / __x >=3D > __rd, > + // and split __urng() value into chunks of __abits, so thei= r > + // products with R^l does not overflow 128bit. These produc= ts > + // are then multiplied by 2^__shift, scaled by __x and adde= d > into > + // the __sumx, until their value is smaller than remaining > value > + // __remx [(__rd - __sumx) * 2^__shift], or whole value was > + // multiplied (__kth is zero). The remainder of division by > + // __x are accumulated into __modx value, and handled later= . > + // In consequence this guarantees that __sumx never exceeds > + // __rd value, and thus does not overflow 128 bit integer. > + _UInt __sumx =3D __lsum / __x, __modx =3D __lsum % __x; > + _UInt __remx =3D __rd - __sumx; > + _UInt64 __kth(__urng() - _Urbg::min()); > + for (unsigned __shift =3D 0; __kth; __shift +=3D __abits) > + { > + const _UInt __val =3D __Rl * (__kth & __amask); > + const _UInt __valx =3D __val / __x; > + if (__valx > __remx) > + break; > + __kth >>=3D __abits; > + > + __sumx +=3D (__valx << __shift); > + __modx +=3D ((__val % __x) << __shift); > + > + __remx -=3D __valx; > + __remx >>=3D __abits; > + } > + if (__kth) // (__sum / __x) > _rd after adding __kth > iteration > + continue; > + > + // Handle accumulated remainders > + const _UInt __valx =3D __modx / __x; > + if (__valx > __rd - __sumx) // avoids overflow > + continue; > + > + __sumx +=3D __valx; > + const _RealT __ret =3D _RealT(__sumx) / _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/ope= rators/gencanon.cc > b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/ope= rators/gencanon.cc > index 1161d5603c2..775e25106d1 100644 > --- > a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/ope= rators/gencanon.cc > +++ > b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/ope= rators/gencanon.cc > @@ -249,6 +249,151 @@ test_2p31m1(const std::mt19937& rng) > } > } > > +template<std::uint64_t Max, typename Under =3D std::mt19937_64> > +struct trimmed_engine > +{ > + using result_type =3D std::uint64_t; > + > + static constexpr > + result_type min() > + { return result_type(0); } > + > + static constexpr > + result_type max() > + { return result_type(Max); } > + > + trimmed_engine() : dist(min(), max()) > + {} > + > + trimmed_engine(const Under& gen) > + : under(gen), dist(min(), max()) > + {} > + > + result_type operator()() > + { return dist(under); } > + > + void discard(std::size_t n) > + { > + for (size_t i =3D 0; i < n; ++i) > + (void)dist(under); > + } > + > + friend bool > + operator=3D=3D(const trimmed_engine& lhs, const trimmed_engine& rhs) > + { return lhs.under =3D=3D rhs.under; } > + > +private: > + Under under; > + std::uniform_int_distribution<result_type> dist; > +}; > + > +// Uses a generator that emits range of size 2^55+2. > +// To produce 128bit IEE float, we need 3 calls to above generator, > +// that can produce value that have 155 bits, and uses separate path. > +template <typename T> > +void > +test_2p55p1(const std::mt19937& rng) > +{ > + if (!std::numeric_limits<T>::is_iec559) > + return; > + > + constexpr size_t mantissa =3D std::numeric_limits<T>::digits; > + constexpr size_t call_per_elem =3D (mantissa / 55) + 1; > + const trimmed_engine<(1ULL << 55) + 1ULL, std::mt19937> lrng{rng}; > + > + int deviation =3D 0, max =3D 0, rms =3D 0, zeros =3D 0; > + int skips =3D run_generator<T, -1u, call_per_elem, 1u> > + (lrng, deviation, max, rms, zeros); > + > + switch (mantissa) > + { > + case 24: // ieee32 > + case 53: // ieee64 > + VERIFY(skips =3D=3D 0); > + VERIFY(deviation =3D=3D 7682); > + VERIFY(max =3D=3D 295); > + VERIFY(rms =3D=3D 971); > + VERIFY(zeros =3D=3D 0); > + break; > + case 64: // ieee80 > + case 106: // ibm128 > + VERIFY(skips =3D=3D 0); > + VERIFY(deviation =3D=3D 6788); > + VERIFY(max =3D=3D 248); > + VERIFY(rms =3D=3D 854); > + VERIFY(zeros =3D=3D 0); > + break; > + case 113: // ieee128 > + VERIFY(skips =3D=3D 0); > + VERIFY(deviation =3D=3D 8222); > + VERIFY(max =3D=3D 257); > + VERIFY(rms =3D=3D 999); > + VERIFY(zeros =3D=3D 0); > + break; > + default: > + VERIFY(false); > + break; > + } > +} > + > +template <typename T> > +void > +test_2p7p1(const std::mt19937& rng) > +{ > + if (!std::numeric_limits<T>::is_iec559) > + return; > + > + constexpr size_t mantissa =3D std::numeric_limits<T>::digits; > + constexpr size_t call_per_elem =3D (mantissa / 7) + 1; > + const trimmed_engine<128U, std::mt19937> lrng{rng}; > + > + int deviation =3D 0, max =3D 0, rms =3D 0, zeros =3D 0; > + int skips =3D run_generator<T, -1u, call_per_elem, 5u> > + (lrng, deviation, max, rms, zeros); > + > + switch (mantissa) > + { > + case 24: // ieee32 > + VERIFY(skips =3D=3D 31641); > + VERIFY(deviation =3D=3D 8508); > + VERIFY(max =3D=3D 294); > + VERIFY(rms =3D=3D 1077); > + VERIFY(zeros =3D=3D 0); > + break; > + case 53: // ieee64 > + VERIFY(skips =3D=3D 63971); > + VERIFY(deviation =3D=3D 7500); > + VERIFY(max =3D=3D 245); > + VERIFY(rms =3D=3D 919); > + VERIFY(zeros =3D=3D 0); > + break; > + case 64: // ieee80 > + VERIFY(skips =3D=3D 2538); > + VERIFY(deviation =3D=3D 8266); > + VERIFY(max =3D=3D 310); > + VERIFY(rms =3D=3D 1065); > + VERIFY(zeros =3D=3D 0); > + break; > + case 106: // ibm128 > + VERIFY(skips =3D=3D 6651); > + VERIFY(deviation =3D=3D 7740); > + VERIFY(max =3D=3D 282); > + VERIFY(rms =3D=3D 1000); > + VERIFY(zeros =3D=3D 0); > + break; > + case 113: // ieee128 > + VERIFY(skips =3D=3D 679); > + VERIFY(deviation =3D=3D 7450); > + VERIFY(max =3D=3D 284); > + VERIFY(rms =3D=3D 936); > + VERIFY(zeros =3D=3D 0); > + break; > + default: > + VERIFY(false); > + break; > + } > +} > + > int main() > { > std::mt19937 rng(8890); > @@ -259,20 +404,28 @@ int main() > test_2p32<float>(rng); > test_10p6<float>(rng); > test_2p31m1<float>(rng); > + test_2p55p1<float>(rng); > + test_2p7p1<float>(rng); > > test_2p32<double>(rng); > test_10p6<double>(rng); > test_2p31m1<double>(rng); > + test_2p55p1<double>(rng); > + test_2p7p1<double>(rng); > > test_2p32<long double>(rng); > test_10p6<long double>(rng); > test_2p31m1<long double>(rng); > + test_2p55p1<long double>(rng); > + test_2p7p1<long double>(rng); > > #ifndef _GLIBCXX_GENERATE_CANONICAL_STRICT > # ifdef __SIZEOF_FLOAT128__ > test_2p32<__float128>(rng); > test_10p6<__float128>(rng); > test_2p31m1<__float128>(rng); > + test_2p55p1<__float128>(rng); > + test_2p7p1<__float128>(rng); > # endif > #endif > } > diff --git > a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/ope= rators/gencanon_eng.cc > b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/ope= rators/gencanon_eng.cc > index 662e3568a30..91848468dc0 100644 > --- > a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/ope= rators/gencanon_eng.cc > +++ > b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/ope= rators/gencanon_eng.cc > @@ -10,13 +10,49 @@ test_engine() > (void)std::generate_canonical<_Real, size_t(-1)>(__engine); > } > > +template<std::uint64_t Max, typename Under =3D std::mt19937_64> > +struct trimmed_engine > +{ > + using result_type =3D std::uint64_t; > + > + static constexpr > + result_type min() > + { return result_type(0); } > + > + static constexpr > + result_type max() > + { return result_type(Max); } > + > + trimmed_engine() : dist(min(), max()) > + {} > + > + result_type operator()() > + { return dist(under); } > + > +private: > + Under under; > + std::uniform_int_distribution<result_type> dist; > +}; > + > +template<typename Real, size_t Bits> > +void > +test_bits() > +{ > + trimmed_engine<(std::uint64_t(1) << Bits) - 1> pow2_engine; > + (void)std::generate_canonical<Real, -1u>(pow2_engine); > + trimmed_engine<(std::uint64_t(1) << Bits) - 2> high_non_pow2_engine; > + (void)std::generate_canonical<Real, -1u>(high_non_pow2_engine); > + trimmed_engine<(std::uint64_t(1) << (Bits - 1))> low_non_pow2_engine; > + (void)std::generate_canonical<Real, -1u>(low_non_pow2_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_rand0>(); > test_engine<_Real, std::minstd_rand>(); > test_engine<_Real, std::mt19937>(); > test_engine<_Real, std::mt19937_64>(); > @@ -25,10 +61,44 @@ test_all_engines() > test_engine<_Real, std::ranlux24>(); > test_engine<_Real, std::ranlux48>(); > test_engine<_Real, std::knuth_b>(); > -#if __cplusplus > 202302L > +#if __cplusplus > 202302L > test_engine<_Real, std::philox4x32>(); > test_engine<_Real, std::philox4x64>(); > #endif > + > + // For 128bit floating points, generator emitting a range, which size = is > + // not power of two, but of width of B bits, such that for any N: > + // N * (B-1) < 113 (bits in ieee128) > + // (N+1) * B > 128 > + // use >128bits patch, as they would otherwise require integer with mo= re > + // than 128 bits. > + // N =3D=3D 3: B in [43, 57] > + test_bits<_Real, 43>(); > + test_bits<_Real, 57>(); > + > + // N =3D=3D 4: B in [33, 38] > + test_bits<_Real, 33>(); > + test_bits<_Real, 38>(); > + > + // N =3D=3D 5: B in [26, 29] > + test_bits<_Real, 26>(); > + test_bits<_Real, 29>(); > + > + // N =3D=3D 6: B in [22, 23] > + test_bits<_Real, 22>(); > + test_bits<_Real, 23>(); > + > + // N =3D=3D 7, B =3D=3D 19 > + test_bits<_Real, 19>(); > + > + // N =3D=3D 8, B =3D=3D 17 > + test_bits<_Real, 17>(); > + > + // N =3D=3D 9, B =3D=3D 15 > + test_bits<_Real, 15>(); > + > + // N >=3D 10 and B < 13 > + test_bits<_Real, 13>(); > } > > int main() > @@ -40,5 +110,5 @@ int main() > # ifdef __SIZEOF_FLOAT128__ > test_all_engines<__float128>(); > # endif > -#endif > +#endif > } > diff --git > a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/ope= rators/gencanon_eng_neg.cc > b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/ope= rators/gencanon_eng_neg.cc > deleted file mode 100644 > index e3b4540bf2b..00000000000 > --- > a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/ope= rators/gencanon_eng_neg.cc > +++ /dev/null > @@ -1,89 +0,0 @@ > -// { dg-do compile { target { c++11 } } } > -// { dg-require-effective-target __float128 } > -// { dg-require-effective-target base_quadfloat_support } > -// { dg-add-options __float128 } > - > -#include <random> > -#include <cstdint> > - > -template<std::uint64_t Max, typename Under =3D std::mt19937_64> > -struct trimmed_engine > -{ > - using result_type =3D std::uint64_t; > - > - static constexpr > - result_type min() > - { return result_type(0); } > - > - static constexpr > - result_type max() > - { return result_type(Max); } > - > - trimmed_engine() : dist(min(), max()) > - {} > - > - result_type operator()() > - { return dist(under); } > - > -private: > - Under under; > - std::uniform_int_distribution<result_type> dist; > -}; > - > -template<typename Real, size_t Bits> > -void > -test_non_pow2() > -{ > - trimmed_engine<(std::uint64_t(1) << Bits) - 2> non_pow2_engine; > - (void)std::generate_canonical<Real, -1u>(non_pow2_engine); > -} > - > -template<typename Real, size_t Bits> > -void > -test_pow2() > -{ > - trimmed_engine<(std::uint64_t(1) << Bits) - 1> pow2_engine; > - (void)std::generate_canonical<Real, -1u>(pow2_engine); > -} > - > -int main() > -{ > - // For 128bit floating points, generator emitting a range, which size > is > - // not power of two, but of width of B bits, such that for any N: > - // N * B < 113 (bits in ieee128) > - // (N+1) * B > 128 > - // are not supported, as they would require integer with more than 127 > bits. > - > - // N =3D=3D 3: B in [43, 57) > - test_non_pow2<__float128, 42>(); // 3 calls > - test_non_pow2<__float128, 43>(); // { dg-error "from here" } > - test_non_pow2<__float128, 56>(); // { dg-error "from here" } > - test_non_pow2<__float128, 57>(); // 2 calls > - test_pow2<__float128, 43>(); > - test_pow2<__float128, 56>(); > - > - // N =3D=3D 4: B in [33, 38) > - test_non_pow2<__float128, 32>(); // 4 calls > - test_non_pow2<__float128, 33>(); // { dg-error "from here" } > - test_non_pow2<__float128, 37>(); // { dg-error "from here" } > - test_non_pow2<__float128, 38>(); // 3 calls > - test_pow2<__float128, 33>(); > - test_pow2<__float128, 37>(); > - > - // N =3D=3D 5: B in [26, 29) > - test_non_pow2<__float128, 25>(); // 5 calls > - test_non_pow2<__float128, 26>(); // { dg-error "from here" } > - test_non_pow2<__float128, 28>(); // { dg-error "from here" } > - test_non_pow2<__float128, 29>(); // 4 calls > - test_pow2<__float128, 26>(); > - test_pow2<__float128, 28>(); > - > - // N =3D=3D 6: B =3D=3D 22 > - test_non_pow2<__float128, 21>(); // 6 calls > - test_non_pow2<__float128, 22>(); // { dg-error "from here" } > - test_non_pow2<__float128, 23>(); // 5 calls > - test_pow2<__float128, 22>(); > -} > - > -// { dg-prune-output "no type named 'type' in 'struct > std::__detail::_Select_uint_least_t" } > -// { dg-prune-output "static assertion failed: sorry, would be too much > trouble for a slow result" } > -- > 2.55.0 > > --0000000000001fdc5106583a9ab8 Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"auto"><div><br><br><div class=3D"gmail_quote gmail_quote_contai= ner"><div dir=3D"ltr" class=3D"gmail_attr">On Tue, 4 Aug 2026, 12:11 Tomasz= Kami=C5=84ski, <<a href=3D"mailto:[email protected]">tkaminsk@redhat.= com</a>> wrote:<br></div><blockquote class=3D"gmail_quote" style=3D"marg= in:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1e= x">This patch splits implementation of generate_canonical_any, into three<b= r> different path depending on bits generated (R), required bits of entropy<br= > (b), and value k such that R^k > 2^d > R^(k-1):<br> * log2(R) >=3D d: single call of generator is sufficient (k =3D 1),<br> =C2=A0 and computation can be performed in the generator value type (_UIntR= );<br> * bit_width(R^k) <=3D 128: computation can be performed using type<br> =C2=A0 __detail::_Select_uint_least_t; previous implementation is reused;<b= r> * bit_width(R^k) > 128: the computations required may overflow 128 bit;<= br> =C2=A0 such generators were previously unsupported, as illustrated by<br> =C2=A0 r17-507-g5d17eec7d83445.<br> <br> This patch implements the support for the last branch, by observing<br> that while it is possible for following computation to overflow 128bit<br> integer:<br> =C2=A0 sum_[0, k) urng() * R^i<br> The values of sum, such that (sum / x * 2^d) > 1, are rejected, thus<br> any value (sum / x) > 2^d can be rejected. As 2^d is guaranteed to fit<b= r> in 128 integer, we can use above to reject the value, before overflow<br> will occur.<br> <br> Firstly, we know that for l =3D k-1, the value R^l < 2^d, and<br> sum_[0, l) urng() * R^i <=3D R^l (because urng() is always < R),<br> always fit in 128bit integer. We compute that value (__lsum) using<br> normal loop.<br> <br> This leaves us with final iteration of multiplying __kth * R^l, where<br> __kth is result of invocation of urng. In that case we compute the<br> bitwidth (__abits) of value that can be multiplied by _R^l without<br> overflow 128bit integer. The __kth * R^l is performed by multiplying<br> the loop of __abits of __kth in loop where:<br> * __sumx (initially __lsum / __x) is the current value of sum / x,<br> * __remx (initially 2^d - __sumx) is the value that we can still,<br> =C2=A0 add before overflowing (shifted with each iteration)<br> * __modx (initially __lsum % __x) are collected reminders.<br> In each iteration the lowest __kth bits are multiplied with _Rl (__val),<br= > scaled (__valx) and checked against __rem. When fits the __sumx, __modx<br> and __remx are updated, and both values of __remx and __kth are shifted<br> right by __abits.<br> <br> If the multiplication finished without overflow (__kth is zero), we handle<= br> the collected remainders (__modx), and compute the __ret value. The values<= br> greater or equal 1 (e.g. due floating point approximation), are still<br> rejected.<br> <br> Similar approach is applied when computing x =3D R^k / 2^d.<br> Again while R^k will overflow 128 bits, we know that x < R (so fits<br> in 64bits), as R^l < 2^d < R^k, then x < R. We again compute<br> R * R^l / 2^d, by splitting the value of R into __abits chunks (_pR).<br> For result of chunks multipliations (_pRk), the bits (__xp) remaining<br> after division by __shift (starts from __d and reduces by __abits),<br> are collected into cumulative value of __x. While the lower bits<br> are preserved for later iteration (__pRk).<br> <br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 PR libstdc++/119739<br> <br> libstdc++-v3/ChangeLog:<br> <br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 * include/bits/random.tcc (__generate_canonical= _any):<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 Rewrite to support generator with different ran= ges.<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 * testsuite/26_numerics/random/uniform_real_dis= tribution/operators/gencanon.cc:<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 Add test for previously unsupported generator.<= br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 * testsuite/26_numerics/random/uniform_real_dis= tribution/operators/gencanon_eng.cc:<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 Add test for previously unsupported-generator r= anges.<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 * testsuite/26_numerics/random/uniform_real_dis= tribution/operators/gencanon_eng_neg.cc:<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 Moved (now-supported) generator test to gencano= n_eng.cc, and removed<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 the file.<br> ---<br> v4:<br> * address typos in comments pointed out in the review<br> * address TODO, initializing __bwRl with actual bitwidth<br> =C2=A0 (we need to handle __rand_uint128 for 32bit targets, without __int12= 8)<br> * add test for 2^7 + 1 =3D 129 generator size in the gencanon_eng<br> * adjust gen_canon_eng numbers.<br> <br> The issues reported by LLM was real, while the bit_width(_Rl) is guaranteed= to<br> be less than 113, because _Rl < 2^113, the aproximation of bit_with(_R) = * l,<br> is not due overestimate. This addresses above.<br> <br> We should consider making __bitwitdh and other function from <bit> fu= nction<br> available in C++11 (they are from 14 now) and providing overloads for<br> __rand_uint128, but this is best done as separate commit.<br></blockquote><= /div></div><div dir=3D"auto"><br></div><div dir=3D"auto">I don't rememb= er why I made them C++14 only, maybe because they can't be constexpr in= C++11 due to multiple statements. But that's ok, they can just be non-= constexpr for C++11.</div><div dir=3D"auto"><br></div><div dir=3D"auto"><br= ></div><div dir=3D"auto"><div class=3D"gmail_quote gmail_quote_container"><= blockquote class=3D"gmail_quote" style=3D"margin:0px 0px 0px 0.8ex;border-l= eft:1px solid rgb(204,204,204);padding-left:1ex"> <br> Tested on x86_64-linux locally. *gencanon* additionally tested with all<br> standard modes, -m32, debug and assertions and on powerpc64le.<br> OK for trunk?<br></blockquote></div></div><div dir=3D"auto"><br></div><div = dir=3D"auto">OK, thanks</div><div dir=3D"auto"><br></div><div dir=3D"auto">= <br></div><div dir=3D"auto"><div class=3D"gmail_quote gmail_quote_container= "><blockquote class=3D"gmail_quote" style=3D"margin:0px 0px 0px 0.8ex;borde= r-left:1px solid rgb(204,204,204);padding-left:1ex"> <br> =C2=A0libstdc++-v3/include/bits/random.tcc=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2= =A0 | 173 +++++++++++++++---<br> =C2=A0.../operators/gencanon.cc=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 = =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0| 153 ++++++++++++++++<br> =C2=A0.../operators/gencanon_eng.cc=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2= =A0 =C2=A0 =C2=A0 =C2=A0|=C2=A0 76 +++++++-<br> =C2=A0.../operators/gencanon_eng_neg.cc=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 = =C2=A0 =C2=A0|=C2=A0 89 ---------<br> =C2=A04 files changed, 371 insertions(+), 120 deletions(-)<br> =C2=A0delete mode 100644 libstdc++-v3/testsuite/26_numerics/random/uniform_= real_distribution/operators/gencanon_eng_neg.cc<br> <br> diff --git a/libstdc++-v3/include/bits/random.tcc b/libstdc++-v3/include/bi= ts/random.tcc<br> index aa176db4778..b2afe411a7c 100644<br> --- a/libstdc++-v3/include/bits/random.tcc<br> +++ b/libstdc++-v3/include/bits/random.tcc<br> @@ -242,7 +242,7 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION<br> =C2=A0 =C2=A0 =C2=A0constexpr size_t<br> =C2=A0 =C2=A0 =C2=A0mersenne_twister_engine<_UIntType, __w, __n, __m, __= r, __a, __u, __d,<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2= =A0 =C2=A0 =C2=A0 =C2=A0 __s, __b, __t, __c, __l, __f>::tempering_u;<br> -=C2=A0 =C2=A0<br> +<br> =C2=A0 =C2=A0template<typename _UIntType,<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0size_t __w, size_t __n, size_t __m= , size_t __r,<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UIntType __a, size_t __u, _UIntTy= pe __d, size_t __s,<br> @@ -3807,38 +3807,155 @@ namespace __detail<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0constexpr unsigned __log2R<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =3D sizeof(_UIntR) * __CHAR_BIT__ - __builtin_c= lzg(__R) - 1;<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0// We overstimate number of required bits, by co= mputing<br> -=C2=A0 =C2=A0 =C2=A0 // r such that l * log2(R) >=3D d, so:<br> -=C2=A0 =C2=A0 =C2=A0 // R^l >=3D (2 ^ log2(R)) ^ l =3D=3D 2 ^ (log2(r) = * l) >=3D 2^d<br> -=C2=A0 =C2=A0 =C2=A0 // And then requiring l * bit_width(R) bits.<br> -=C2=A0 =C2=A0 =C2=A0 constexpr unsigned __l =3D (__d + __log2R - 1) / __lo= g2R;<br> -=C2=A0 =C2=A0 =C2=A0 constexpr unsigned __bits =3D (__log2R + 1) * __l;<br= > -=C2=A0 =C2=A0 =C2=A0 using _UInt =3D typename __detail::_Select_uint_least= _t<__bits>::type;<br> -<br> -=C2=A0 =C2=A0 =C2=A0 _GLIBCXX_GEN_CANON_CONST _UInt __rd =3D _UInt(1) <= < __d;<br> -=C2=A0 =C2=A0 =C2=A0 _GLIBCXX_GEN_CANON_CONST auto __logRrd =3D __gen_cano= n_log(__rd, __R);<br> -=C2=A0 =C2=A0 =C2=A0 _GLIBCXX_GEN_CANON_CONST unsigned __k<br> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =3D __logRrd.__floor_log + (__rd > __logRrd= .__floor_pow);<br> -<br> -=C2=A0 =C2=A0 =C2=A0 _GLIBCXX_GEN_CANON_CONST _UInt __Rk<br> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =3D (__k > __logRrd.__floor_log)<br> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 ? _UInt(__logRrd.__floor_pow) * _UInt(_= _R)<br> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 : _UInt(__logRrd.__floor_pow);<br> -=C2=A0 =C2=A0 =C2=A0 _GLIBCXX_GEN_CANON_CONST _UInt __x =3D=C2=A0 __Rk / _= _rd;<br> +=C2=A0 =C2=A0 =C2=A0 // m such that m * log2(R) >=3D d, so:<br> +=C2=A0 =C2=A0 =C2=A0 // R^m >=3D (2 ^ log2(R)) ^ m =3D=3D 2 ^ (log2(R) = * m) >=3D 2^d<br> +=C2=A0 =C2=A0 =C2=A0 // And then requiring m * bit_width(R) bits.<br> +=C2=A0 =C2=A0 =C2=A0 constexpr unsigned __m =3D (__d + __log2R - 1) / __lo= g2R;<br> +=C2=A0 =C2=A0 =C2=A0 constexpr unsigned __bits =3D (__log2R + 1) * __m;<br= > +=C2=A0 =C2=A0 =C2=A0 if constexpr (__log2R >=3D __d)<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0{<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// range already provide required number= of bits,<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// so in this case __k =3D=3D 1, and sin= gle call to generator<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// is sufficient. Furthermore 2^d fits i= n _UIntR<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0constexpr _UIntR __rd =3D _UIntR(1) <= < __d;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0constexpr _UIntR __x =3D __R >> __= d;<br> <br> -=C2=A0 =C2=A0 =C2=A0 while (true)<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0while (true)<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0{<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UIntR __val(__urng() - _U= rbg::min());<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0const _RealT __ret =3D _Re= alT(__val / __x) / _RealT(__rd);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0if (__ret < _RealT(1.0)= )<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0return __ret;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0}<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0}<br> +=C2=A0 =C2=A0 =C2=A0 // generator call produce fewer bits than __d, but R^= k > rd fits<br> +=C2=A0 =C2=A0 =C2=A0 // into 128 bits.<br> +=C2=A0 =C2=A0 =C2=A0 else if constexpr (__bits <=3D 128)<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 {<br> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UInt __Ri{1};<br> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UInt __sum(__urng() - _Urbg::min());<br= > -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0for (int __i =3D __k - 1; __i > 0; --= __i)<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0using _UInt =3D typename __detail::_Sele= ct_uint_least_t<__bits>::type;<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST _UInt __rd =3D = _UInt(1) << __d;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST auto __logRrd = =3D __gen_canon_log(__rd, __R);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// __rd is power of two, and __R is not = so __floor_pow is never<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// equal to __rd.<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST unsigned __k = =3D __logRrd.__floor_log + 1;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST _UInt __Rk =3D = __logRrd.__floor_pow * _UInt(__R);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// x =3D floor(R^k / 2^d) =3D R^k >&g= t; d;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST _UInt __x =3D= =C2=A0 __Rk >> __d;<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0while (true)<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 {<br> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__Ri *=3D _UInt(__R);<br> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__sum +=3D _UInt(__urng() = - _Urbg::min()) * __Ri;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UInt __Ri{1};<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UInt __sum(__urng() - _Ur= bg::min());<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0for (int __i =3D __k - 1; = __i > 0; --__i)<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0{<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__Ri *=3D __= R;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__sum +=3D _= _Ri * _UIntR(__urng() - _Urbg::min());<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0}<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0const _RealT __ret =3D _Re= alT(__sum / __x) / _RealT(__rd);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0if (__ret < _RealT(1.0)= )<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0return __ret;<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 }<br> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0const _RealT __ret =3D _RealT(__sum / __= x) / _RealT(__rd);<br> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0if (__ret < _RealT(1.0))<br> -=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0return __ret;<br> =C2=A0 =C2=A0 =C2=A0 =C2=A0 }<br> -#undef _GLIBCXX_GEN_CANON_CONST <br> +=C2=A0 =C2=A0 =C2=A0 else<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0{<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0static_assert(__log2R < 64, "Onl= y generators emitting up to 64 bits are supported");<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0using _UInt =3D typename __detail::_Sele= ct_uint_least_t<128>::type;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0using _UInt64 =3D typename __detail::_Se= lect_uint_least_t<64>::type;<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST _UInt __rd =3D = _UInt(1) << __d;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST auto __logRrd = =3D __gen_canon_log(__rd, __R);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// __rd is power of two, and __R is not = so __floor_pow is never<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// equal to __rd, and l =3D=3D k - 1<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST unsigned __l = =3D __logRrd.__floor_log;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST _UInt __Rl =3D = __logRrd.__floor_pow;<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// __abits is maximum bit width of the v= alue, that can<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// be multiplied by R^l without overflow= ing 128 bit integer<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST unsigned __bwRl= <br> +#ifndef __SIZEOF_INT128__<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0=3D __Rl._M_hi ? 128 - __builtin_= clzg(__Rl._M_hi)<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2= =A0 =C2=A0 : 64 - __builtin_clzg(__Rl._M_lo);<br> +#else<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0=3D 128 - __builtin_clzg(__Rl);<b= r> +#endif<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// For __R close to power of two, the ac= tual __k may be smaller than __m,<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// and will use less than 128bits, defau= lt to two 32 bits chunks.<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST unsigned __abit= s =3D __bwRl > 96 ? 128 - __bwRl : 32;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_GLIBCXX_GEN_CANON_CONST _UInt64 __amask= =3D (_UInt64(1) << __abits) - 1;<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// Compute __x =3D _Rk / __rd (_Rk >&= gt; __d), by spliting _R into chunks<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// (_pR) of __abits, so their multiplica= tion does not overflow 128 bits.<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// __x fits in 64bits, becuse __R^l <= __rd < __R^k, thus __x < _R<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UInt64 __x(0), __pR(__R); _UInt __pRk(0= );<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0for (unsigned __shift =3D __d; __pR; __s= hift -=3D __abits)<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0{<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__pRk +=3D __Rl * (__pR &a= mp; __amask);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__pR >>=3D __abits;<= br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UInt __xp(__pRk >> = __shift);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__x +=3D _UInt64(__xp);<br= > +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__pRk -=3D (__xp << = __shift);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__pRk >>=3D __abits;= <br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0}<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0while (true)<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0{<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// Compute sum of __urng()= * R^i for i in range [0, l).<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// The value is smaller th= an R^l which is smaller than<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// __rd so fits into 128 b= it integer.<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UInt __Ri{1};<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UInt __lsum(__urng() - _U= rbg::min());<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0for (int __i =3D __l - 1; = __i > 0; --__i)<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0{<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__Ri *=3D _U= Int64(__R);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__lsum +=3D = __Ri * _UInt64(__urng() - _Urbg::min());<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0}<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// The last iteration __ur= ng() * R^k may overflow 128bit<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// integer. We use the fac= t that we reject __sum / __x >=3D __rd,<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// and split __urng() valu= e into chunks of __abits, so their<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// products with R^l does = not overflow 128bit. These products<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// are then multiplied by = 2^__shift, scaled by __x and added into<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// the __sumx, until their= value is smaller than remaining value<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// __remx [(__rd - __sumx)= * 2^__shift], or whole value was<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// multiplied (__kth is ze= ro). The remainder of division by<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// __x are accumulated int= o __modx value, and handled later.<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// In consequence this gua= rantees that __sumx never exceeds<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// __rd value, and thus do= es not overflow 128 bit integer.<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UInt __sumx =3D __lsum / = __x, __modx =3D __lsum % __x;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UInt __remx =3D __rd - __= sumx;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0_UInt64 __kth(__urng() - _= Urbg::min());<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0for (unsigned __shift =3D = 0; __kth; __shift +=3D __abits)<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0{<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0const _UInt = __val =3D __Rl * (__kth & __amask);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0const _UInt = __valx =3D __val / __x;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0if (__valx &= gt; __remx)<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0break= ;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__kth >&g= t;=3D __abits;<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__sumx +=3D = (__valx << __shift);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__modx +=3D = ((__val % __x) << __shift);<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__remx -=3D = __valx;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__remx >&= gt;=3D __abits;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0}<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0if (__kth) // (__sum / __x= ) > _rd after adding __kth iteration<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0continue;<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0// Handle accumulated rema= inders<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0const _UInt __valx =3D __m= odx / __x;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0if (__valx > __rd - __s= umx) // avoids overflow<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0continue;<br> +<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0__sumx +=3D __valx;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0const _RealT __ret =3D _Re= alT(__sumx) / _RealT(__rd);<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0if (__ret < _RealT(1.0)= )<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0return __ret;<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0}<br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0}<br> +#undef _GLIBCXX_GEN_CANON_CONST<br> =C2=A0 =C2=A0 =C2=A0}<br> <br> =C2=A0#if !defined(_GLIBCXX_GENERATE_CANONICAL_STRICT)<br> diff --git a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distrib= ution/operators/gencanon.cc b/libstdc++-v3/testsuite/26_numerics/random/uni= form_real_distribution/operators/gencanon.cc<br> index 1161d5603c2..775e25106d1 100644<br> --- a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/o= perators/gencanon.cc<br> +++ b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/o= perators/gencanon.cc<br> @@ -249,6 +249,151 @@ test_2p31m1(const std::mt19937& rng)<br> =C2=A0 =C2=A0}<br> =C2=A0}<br> <br> +template<std::uint64_t Max, typename Under =3D std::mt19937_64><br> +struct trimmed_engine<br> +{<br> +=C2=A0 using result_type =3D std::uint64_t;<br> +<br> +=C2=A0 static constexpr<br> +=C2=A0 result_type min()<br> +=C2=A0 { return result_type(0); }<br> +<br> +=C2=A0 static constexpr<br> +=C2=A0 result_type max()<br> +=C2=A0 { return result_type(Max); }<br> +<br> +=C2=A0 trimmed_engine() : dist(min(), max())<br> +=C2=A0 {}<br> +<br> +=C2=A0 trimmed_engine(const Under& gen)<br> +=C2=A0 : under(gen), dist(min(), max())<br> +=C2=A0 {}<br> +<br> +=C2=A0 result_type operator()()<br> +=C2=A0 { return dist(under); }<br> +<br> +=C2=A0 void discard(std::size_t n)<br> +=C2=A0 {<br> +=C2=A0 =C2=A0 for (size_t i =3D 0; i < n; ++i)<br> +=C2=A0 =C2=A0 =C2=A0 (void)dist(under);<br> +=C2=A0 }<br> +<br> +=C2=A0 friend bool<br> +=C2=A0 operator=3D=3D(const trimmed_engine& lhs, const trimmed_engine&= amp; rhs)<br> +=C2=A0 { return lhs.under =3D=3D rhs.under; }<br> +<br> +private:<br> +=C2=A0 Under under;<br> +=C2=A0 std::uniform_int_distribution<result_type> dist;<br> +};<br> +<br> +// Uses a generator that emits range of size 2^55+2.<br> +// To produce 128bit IEE float, we need 3 calls to above generator,<br> +// that can produce value that have 155 bits, and uses separate path.<br> +template <typename T><br> +void<br> +test_2p55p1(const std::mt19937& rng)<br> +{<br> +=C2=A0 if (!std::numeric_limits<T>::is_iec559)<br> +=C2=A0 =C2=A0 return;<br> +<br> +=C2=A0 constexpr size_t mantissa =3D std::numeric_limits<T>::digits;= <br> +=C2=A0 constexpr size_t call_per_elem =3D (mantissa / 55) + 1;<br> +=C2=A0 const trimmed_engine<(1ULL << 55) + 1ULL, std::mt19937>= lrng{rng};<br> +<br> +=C2=A0 int deviation =3D 0, max =3D 0, rms =3D 0, zeros =3D 0;<br> +=C2=A0 int skips =3D run_generator<T, -1u, call_per_elem, 1u><br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0(lrng, deviation, m= ax, rms, zeros);<br> +<br> +=C2=A0 switch (mantissa)<br> +=C2=A0 {<br> +=C2=A0 case 24: // ieee32<br> +=C2=A0 case 53: // ieee64<br> +=C2=A0 =C2=A0 VERIFY(skips =3D=3D 0);<br> +=C2=A0 =C2=A0 VERIFY(deviation =3D=3D 7682);<br> +=C2=A0 =C2=A0 VERIFY(max =3D=3D 295);<br> +=C2=A0 =C2=A0 VERIFY(rms =3D=3D 971);<br> +=C2=A0 =C2=A0 VERIFY(zeros =3D=3D 0);<br> +=C2=A0 =C2=A0 break;<br> +=C2=A0 case 64: // ieee80<br> +=C2=A0 case 106: // ibm128<br> +=C2=A0 =C2=A0 VERIFY(skips =3D=3D 0);<br> +=C2=A0 =C2=A0 VERIFY(deviation =3D=3D 6788);<br> +=C2=A0 =C2=A0 VERIFY(max =3D=3D 248);<br> +=C2=A0 =C2=A0 VERIFY(rms =3D=3D 854);<br> +=C2=A0 =C2=A0 VERIFY(zeros =3D=3D 0);<br> +=C2=A0 =C2=A0 break;<br> +=C2=A0 case 113: // ieee128<br> +=C2=A0 =C2=A0 VERIFY(skips =3D=3D 0);<br> +=C2=A0 =C2=A0 VERIFY(deviation =3D=3D 8222);<br> +=C2=A0 =C2=A0 VERIFY(max =3D=3D 257);<br> +=C2=A0 =C2=A0 VERIFY(rms =3D=3D 999);<br> +=C2=A0 =C2=A0 VERIFY(zeros =3D=3D 0);<br> +=C2=A0 =C2=A0 break;<br> +=C2=A0 default:<br> +=C2=A0 =C2=A0 VERIFY(false);<br> +=C2=A0 =C2=A0 break;<br> +=C2=A0 }<br> +}<br> +<br> +template <typename T><br> +void<br> +test_2p7p1(const std::mt19937& rng)<br> +{<br> +=C2=A0 if (!std::numeric_limits<T>::is_iec559)<br> +=C2=A0 =C2=A0 return;<br> +<br> +=C2=A0 constexpr size_t mantissa =3D std::numeric_limits<T>::digits;= <br> +=C2=A0 constexpr size_t call_per_elem =3D (mantissa / 7) + 1;<br> +=C2=A0 const trimmed_engine<128U, std::mt19937> lrng{rng};<br> +<br> +=C2=A0 int deviation =3D 0, max =3D 0, rms =3D 0, zeros =3D 0;<br> +=C2=A0 int skips =3D run_generator<T, -1u, call_per_elem, 5u><br> +=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0(lrng, deviation, m= ax, rms, zeros);<br> +<br> +=C2=A0 switch (mantissa)<br> +=C2=A0 {<br> +=C2=A0 case 24: // ieee32<br> +=C2=A0 =C2=A0 VERIFY(skips =3D=3D 31641);<br> +=C2=A0 =C2=A0 VERIFY(deviation =3D=3D 8508);<br> +=C2=A0 =C2=A0 VERIFY(max =3D=3D 294);<br> +=C2=A0 =C2=A0 VERIFY(rms =3D=3D 1077);<br> +=C2=A0 =C2=A0 VERIFY(zeros =3D=3D 0);<br> +=C2=A0 =C2=A0 break;<br> +=C2=A0 case 53: // ieee64<br> +=C2=A0 =C2=A0 VERIFY(skips =3D=3D 63971);<br> +=C2=A0 =C2=A0 VERIFY(deviation =3D=3D 7500);<br> +=C2=A0 =C2=A0 VERIFY(max =3D=3D 245);<br> +=C2=A0 =C2=A0 VERIFY(rms =3D=3D 919);<br> +=C2=A0 =C2=A0 VERIFY(zeros =3D=3D 0);<br> +=C2=A0 =C2=A0 break;<br> +=C2=A0 case 64: // ieee80<br> +=C2=A0 =C2=A0 VERIFY(skips =3D=3D 2538);<br> +=C2=A0 =C2=A0 VERIFY(deviation =3D=3D 8266);<br> +=C2=A0 =C2=A0 VERIFY(max =3D=3D 310);<br> +=C2=A0 =C2=A0 VERIFY(rms =3D=3D 1065);<br> +=C2=A0 =C2=A0 VERIFY(zeros =3D=3D 0);<br> +=C2=A0 =C2=A0 break;<br> +=C2=A0 case 106: // ibm128<br> +=C2=A0 =C2=A0 VERIFY(skips =3D=3D 6651);<br> +=C2=A0 =C2=A0 VERIFY(deviation =3D=3D 7740);<br> +=C2=A0 =C2=A0 VERIFY(max =3D=3D 282);<br> +=C2=A0 =C2=A0 VERIFY(rms =3D=3D 1000);<br> +=C2=A0 =C2=A0 VERIFY(zeros =3D=3D 0);<br> +=C2=A0 =C2=A0 break;<br> +=C2=A0 case 113: // ieee128<br> +=C2=A0 =C2=A0 VERIFY(skips =3D=3D 679);<br> +=C2=A0 =C2=A0 VERIFY(deviation =3D=3D 7450);<br> +=C2=A0 =C2=A0 VERIFY(max =3D=3D 284);<br> +=C2=A0 =C2=A0 VERIFY(rms =3D=3D 936);<br> +=C2=A0 =C2=A0 VERIFY(zeros =3D=3D 0);<br> +=C2=A0 =C2=A0 break;<br> +=C2=A0 default:<br> +=C2=A0 =C2=A0 VERIFY(false);<br> +=C2=A0 =C2=A0 break;<br> +=C2=A0 }<br> +}<br> +<br> =C2=A0int main()<br> =C2=A0{<br> =C2=A0 =C2=A0std::mt19937 rng(8890);<br> @@ -259,20 +404,28 @@ int main()<br> =C2=A0 =C2=A0test_2p32<float>(rng);<br> =C2=A0 =C2=A0test_10p6<float>(rng);<br> =C2=A0 =C2=A0test_2p31m1<float>(rng);<br> +=C2=A0 test_2p55p1<float>(rng);<br> +=C2=A0 test_2p7p1<float>(rng);<br> <br> =C2=A0 =C2=A0test_2p32<double>(rng);<br> =C2=A0 =C2=A0test_10p6<double>(rng);<br> =C2=A0 =C2=A0test_2p31m1<double>(rng);<br> +=C2=A0 test_2p55p1<double>(rng);<br> +=C2=A0 test_2p7p1<double>(rng);<br> <br> =C2=A0 =C2=A0test_2p32<long double>(rng);<br> =C2=A0 =C2=A0test_10p6<long double>(rng);<br> =C2=A0 =C2=A0test_2p31m1<long double>(rng);<br> +=C2=A0 test_2p55p1<long double>(rng);<br> +=C2=A0 test_2p7p1<long double>(rng);<br> <br> =C2=A0#ifndef _GLIBCXX_GENERATE_CANONICAL_STRICT<br> =C2=A0#=C2=A0 ifdef __SIZEOF_FLOAT128__<br> =C2=A0 =C2=A0test_2p32<__float128>(rng);<br> =C2=A0 =C2=A0test_10p6<__float128>(rng);<br> =C2=A0 =C2=A0test_2p31m1<__float128>(rng);<br> +=C2=A0 test_2p55p1<__float128>(rng);<br> +=C2=A0 test_2p7p1<__float128>(rng);<br> =C2=A0#=C2=A0 endif<br> =C2=A0#endif<br> =C2=A0}<br> diff --git a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distrib= ution/operators/gencanon_eng.cc b/libstdc++-v3/testsuite/26_numerics/random= /uniform_real_distribution/operators/gencanon_eng.cc<br> index 662e3568a30..91848468dc0 100644<br> --- a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/o= perators/gencanon_eng.cc<br> +++ b/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/o= perators/gencanon_eng.cc<br> @@ -10,13 +10,49 @@ test_engine()<br> =C2=A0 =C2=A0(void)std::generate_canonical<_Real, size_t(-1)>(__engin= e);<br> =C2=A0}<br> <br> +template<std::uint64_t Max, typename Under =3D std::mt19937_64><br> +struct trimmed_engine<br> +{<br> +=C2=A0 using result_type =3D std::uint64_t;<br> +<br> +=C2=A0 static constexpr<br> +=C2=A0 result_type min()<br> +=C2=A0 { return result_type(0); }<br> +<br> +=C2=A0 static constexpr<br> +=C2=A0 result_type max()<br> +=C2=A0 { return result_type(Max); }<br> +<br> +=C2=A0 trimmed_engine() : dist(min(), max())<br> +=C2=A0 {}<br> +<br> +=C2=A0 result_type operator()()<br> +=C2=A0 { return dist(under); }<br> +<br> +private:<br> +=C2=A0 Under under;<br> +=C2=A0 std::uniform_int_distribution<result_type> dist;<br> +};<br> +<br> +template<typename Real, size_t Bits><br> +void<br> +test_bits()<br> +{<br> +=C2=A0 trimmed_engine<(std::uint64_t(1) << Bits) - 1> pow2_eng= ine;<br> +=C2=A0 (void)std::generate_canonical<Real, -1u>(pow2_engine);<br> +=C2=A0 trimmed_engine<(std::uint64_t(1) << Bits) - 2> high_non= _pow2_engine;<br> +=C2=A0 (void)std::generate_canonical<Real, -1u>(high_non_pow2_engine= );<br> +=C2=A0 trimmed_engine<(std::uint64_t(1) << (Bits - 1))> low_no= n_pow2_engine;<br> +=C2=A0 (void)std::generate_canonical<Real, -1u>(low_non_pow2_engine)= ;<br> +}<br> +<br> =C2=A0template<typename _Real><br> =C2=A0void<br> =C2=A0test_all_engines()<br> =C2=A0{<br> =C2=A0 =C2=A0test_engine<_Real, std::default_random_engine>();<br> <br> -=C2=A0 test_engine<_Real, std::minstd_rand0>(); <br> +=C2=A0 test_engine<_Real, std::minstd_rand0>();<br> =C2=A0 =C2=A0test_engine<_Real, std::minstd_rand>();<br> =C2=A0 =C2=A0test_engine<_Real, std::mt19937>();<br> =C2=A0 =C2=A0test_engine<_Real, std::mt19937_64>();<br> @@ -25,10 +61,44 @@ test_all_engines()<br> =C2=A0 =C2=A0test_engine<_Real, std::ranlux24>();<br> =C2=A0 =C2=A0test_engine<_Real, std::ranlux48>();<br> =C2=A0 =C2=A0test_engine<_Real, std::knuth_b>();<br> -#if __cplusplus > 202302L <br> +#if __cplusplus > 202302L<br> =C2=A0 =C2=A0test_engine<_Real, std::philox4x32>();<br> =C2=A0 =C2=A0test_engine<_Real, std::philox4x64>();<br> =C2=A0#endif<br> +<br> +=C2=A0 // For 128bit floating points, generator emitting a range, which si= ze is<br> +=C2=A0 // not power of two, but of width of B bits, such that for any N:<b= r> +=C2=A0 // N * (B-1) < 113 (bits in ieee128)<br> +=C2=A0 // (N+1) * B > 128<br> +=C2=A0 // use >128bits patch, as they would otherwise require integer w= ith more<br> +=C2=A0 // than 128 bits.<br> +=C2=A0 // N =3D=3D 3: B in [43, 57]<br> +=C2=A0 test_bits<_Real, 43>();<br> +=C2=A0 test_bits<_Real, 57>();<br> +<br> +=C2=A0 // N =3D=3D 4: B in [33, 38]<br> +=C2=A0 test_bits<_Real, 33>();<br> +=C2=A0 test_bits<_Real, 38>();<br> +<br> +=C2=A0 // N =3D=3D 5: B in [26, 29]<br> +=C2=A0 test_bits<_Real, 26>();<br> +=C2=A0 test_bits<_Real, 29>();<br> +<br> +=C2=A0 // N =3D=3D 6: B in [22, 23]<br> +=C2=A0 test_bits<_Real, 22>();<br> +=C2=A0 test_bits<_Real, 23>();<br> +<br> +=C2=A0 // N =3D=3D 7, B =3D=3D 19<br> +=C2=A0 test_bits<_Real, 19>();<br> +<br> +=C2=A0 // N =3D=3D 8, B =3D=3D 17<br> +=C2=A0 test_bits<_Real, 17>();<br> +<br> +=C2=A0 // N =3D=3D 9, B =3D=3D 15<br> +=C2=A0 test_bits<_Real, 15>();<br> +<br> +=C2=A0 // N >=3D 10 and B < 13<br> +=C2=A0 test_bits<_Real, 13>();<br> =C2=A0}<br> <br> =C2=A0int main()<br> @@ -40,5 +110,5 @@ int main()<br> =C2=A0#=C2=A0 ifdef __SIZEOF_FLOAT128__<br> =C2=A0 =C2=A0test_all_engines<__float128>();<br> =C2=A0#=C2=A0 endif<br> -#endif <br> +#endif<br> =C2=A0}<br> diff --git a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distrib= ution/operators/gencanon_eng_neg.cc b/libstdc++-v3/testsuite/26_numerics/ra= ndom/uniform_real_distribution/operators/gencanon_eng_neg.cc<br> deleted file mode 100644<br> index e3b4540bf2b..00000000000<br> --- a/libstdc++-v3/testsuite/26_numerics/random/uniform_real_distribution/o= perators/gencanon_eng_neg.cc<br> +++ /dev/null<br> @@ -1,89 +0,0 @@<br> -// { dg-do compile { target { c++11 } } }<br> -// { dg-require-effective-target __float128 }<br> -// { dg-require-effective-target base_quadfloat_support }<br> -// { dg-add-options __float128 }<br> -<br> -#include <random><br> -#include <cstdint><br> -<br> -template<std::uint64_t Max, typename Under =3D std::mt19937_64><br> -struct trimmed_engine<br> -{<br> -=C2=A0 using result_type =3D std::uint64_t;<br> -<br> -=C2=A0 static constexpr<br> -=C2=A0 result_type min()<br> -=C2=A0 { return result_type(0); }<br> -<br> -=C2=A0 static constexpr<br> -=C2=A0 result_type max()<br> -=C2=A0 { return result_type(Max); }<br> -<br> -=C2=A0 trimmed_engine() : dist(min(), max())<br> -=C2=A0 {}<br> -<br> -=C2=A0 result_type operator()()<br> -=C2=A0 { return dist(under); }<br> -<br> -private:<br> -=C2=A0 Under under;<br> -=C2=A0 std::uniform_int_distribution<result_type> dist;<br> -};<br> -<br> -template<typename Real, size_t Bits><br> -void<br> -test_non_pow2()<br> -{<br> -=C2=A0 trimmed_engine<(std::uint64_t(1) << Bits) - 2> non_pow2= _engine;<br> -=C2=A0 (void)std::generate_canonical<Real, -1u>(non_pow2_engine);<br= > -}<br> -<br> -template<typename Real, size_t Bits><br> -void<br> -test_pow2()<br> -{<br> -=C2=A0 trimmed_engine<(std::uint64_t(1) << Bits) - 1> pow2_eng= ine;<br> -=C2=A0 (void)std::generate_canonical<Real, -1u>(pow2_engine);<br> -}<br> -<br> -int main()<br> -{<br> -=C2=A0 // For 128bit floating points, generator emitting a range, which si= ze is <br> -=C2=A0 // not power of two, but of width of B bits, such that for any N:<b= r> -=C2=A0 // N * B < 113 (bits in ieee128)<br> -=C2=A0 // (N+1) * B > 128<br> -=C2=A0 // are not supported, as they would require integer with more than = 127 bits.<br> -<br> -=C2=A0 // N =3D=3D 3: B in [43, 57)<br> -=C2=A0 test_non_pow2<__float128, 42>(); // 3 calls<br> -=C2=A0 test_non_pow2<__float128, 43>(); // { dg-error "from her= e" }<br> -=C2=A0 test_non_pow2<__float128, 56>(); // { dg-error "from her= e" }<br> -=C2=A0 test_non_pow2<__float128, 57>(); // 2 calls<br> -=C2=A0 test_pow2<__float128, 43>();<br> -=C2=A0 test_pow2<__float128, 56>();<br> -<br> -=C2=A0 // N =3D=3D 4: B in [33, 38)<br> -=C2=A0 test_non_pow2<__float128, 32>(); // 4 calls<br> -=C2=A0 test_non_pow2<__float128, 33>(); // { dg-error "from her= e" }<br> -=C2=A0 test_non_pow2<__float128, 37>(); // { dg-error "from her= e" }<br> -=C2=A0 test_non_pow2<__float128, 38>(); // 3 calls<br> -=C2=A0 test_pow2<__float128, 33>();<br> -=C2=A0 test_pow2<__float128, 37>();<br> -<br> -=C2=A0 // N =3D=3D 5: B in [26, 29)<br> -=C2=A0 test_non_pow2<__float128, 25>(); // 5 calls<br> -=C2=A0 test_non_pow2<__float128, 26>(); // { dg-error "from her= e" }<br> -=C2=A0 test_non_pow2<__float128, 28>(); // { dg-error "from her= e" }<br> -=C2=A0 test_non_pow2<__float128, 29>(); // 4 calls<br> -=C2=A0 test_pow2<__float128, 26>();<br> -=C2=A0 test_pow2<__float128, 28>();<br> -<br> -=C2=A0 // N =3D=3D 6: B =3D=3D 22<br> -=C2=A0 test_non_pow2<__float128, 21>(); // 6 calls<br> -=C2=A0 test_non_pow2<__float128, 22>(); // { dg-error "from her= e" }<br> -=C2=A0 test_non_pow2<__float128, 23>(); // 5 calls<br> -=C2=A0 test_pow2<__float128, 22>();<br> -}<br> -<br> -// { dg-prune-output "no type named 'type' in 'struct std= ::__detail::_Select_uint_least_t" }<br> -// { dg-prune-output "static assertion failed: sorry, would be too mu= ch trouble for a slow result" }<br> -- <br> 2.55.0<br> <br> </blockquote></div></div></div> --0000000000001fdc5106583a9ab8--