[PATCH v2] libstdc++: Expand use of Lemire's algorithm to all generators emitting power of two range.
Jonathan Wakely
jwakely@redhat.com
Mon Sep 21 14:40:05 GMT 2026
On Mon, 21 Sep 2026 at 14:48 +0200, Tomasz Kamiński wrote:
>This patch apply's the Lemire's algorithm to all generators that emit values
"applies"
>in range which size is power of two (2^_Bits), instead only 2^32 and 2^64.
>This is achieved by masking the result of conversion of __product to lowest
>_Bits. If _Bits is equal to bitwidth of _Up type, the bitand with all set
>mask is optimized out by compilers.
>
>libstdc++-v3/ChangeLog:
>
> * include/bits/uniform_int_dist.h
> (uniform_int_distribution::operator()): Use _S_nd (Lemire's
> algorithm) if __urngrange is power of two.
> (uniform_int_distribution::_S_nd): Add _Bits parameter. Mask
> the result of conversion of __product to _Up to lowest _Bits.
> * testsuite/26_numerics/random/uniform_int_distribution/operators/values.cc:
> Add testDiscreteDist for ranlux24 and ranlux48 engines.
>---
>v2 is rebased on top of the trunk, where __high variable was removed.
OK for trunk with some typos fixed, see below.
> libstdc++-v3/include/bits/uniform_int_dist.h | 46 ++++++++++---------
> .../operators/values.cc | 24 ++++++----
> 2 files changed, 40 insertions(+), 30 deletions(-)
>
>diff --git a/libstdc++-v3/include/bits/uniform_int_dist.h b/libstdc++-v3/include/bits/uniform_int_dist.h
>index 5ee9945b61c..cdb2c3429d6 100644
>--- a/libstdc++-v3/include/bits/uniform_int_dist.h
>+++ b/libstdc++-v3/include/bits/uniform_int_dist.h
>@@ -252,7 +252,7 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
> // Lemire's nearly divisionless algorithm.
> // Returns an unbiased random number from __g downscaled to [0,__range)
> // using an unsigned type _Wp twice as wide as unsigned type _Up.
>- template<typename _Wp, typename _Urbg, typename _Up>
>+ template<size_t _Bits, typename _Wp, typename _Urbg, typename _Up>
> static _Up
> _S_nd(_Urbg& __g, _Up __range)
> {
>@@ -263,33 +263,35 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
> static_assert(_Wp_traits::__digits == (2 * _Up_traits::__digits),
> "W must be twice as wide as U");
> constexpr auto __min = _Urbg::min();
>+ constexpr _Up __mask = (_Bits < _Up_traits::__digits)
>+ ? (_Up(1) << _Bits) - 1 : ~_Up(0);
>
> // reference: Fast Random Integer Generation in an Interval
> // ACM Transactions on Modeling and Computer Simulation 29 (1), 2019
> // https://arxiv.org/abs/1805.10941
> _Wp __product = _Wp(__g() - __min) * _Wp(__range);
>- _Up __low = _Up(__product);
>+ _Up __low = _Up(__product) & __mask;
> if (__low < __range)
> {
> const _Up __threshold = -__range % __range;
> while (__low < __threshold)
> {
>- __product = _Wp(__g() - __min) * _Wp(__range);
>- __low = _Up(__product);
>+ __product = _Wp(__g() - __min) * _Wp(__range);
>+ __low = _Up(__product) & __mask;
>
> // The algorithm is modified to alternate between rejecting
> // from the beginning and end of the range. This guarantees
> // that we stop for non-uniform engines that always result
> // in values below the __threshold.
>- const _Up __back_treshold = _Up_traits::__max - __threshold;
>+ const _Up __back_treshold = __mask - __threshold;
I didn't notice when you added it, but this should be "back_threshold"
not "back_treshold", i.e. an extra 'h'
> if (__low <= __back_treshold)
> break;
>
> __product = _Wp(__g() - __min) * _Wp(__range);
>- __low = _Up(__product);
>+ __low = _Up(__product) & __mask;
> }
> }
>- return __product >> _Up_traits::__digits;
>+ return _Up(__product >> _Bits) & __mask;
> }
> };
>
>@@ -322,25 +324,25 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
>
> #pragma GCC diagnostic push
> #pragma GCC diagnostic ignored "-Wc++17-extensions" // if constexpr
>-#if defined __UINT64_TYPE__ && defined __UINT32_TYPE__
>-#if __SIZEOF_INT128__
>- if constexpr (__urngrange == __UINT64_MAX__)
>+#if defined __UINT64_TYPE__ && defined __UINT32_TYPE__ && __has_builtin(__builtin_popcountg)
>+ constexpr auto __bits = __builtin_popcountg(__urngrange);
>+ if constexpr (__detail::_Power_of_2(__urngrange + 1) && __bits <= 32)
> {
>- // __urng produces values that use exactly 64-bits,
>- // so use 128-bit integers to downscale to desired range.
>- __UINT64_TYPE__ __u64erange = __uerange;
>- __ret = __extension__ _S_nd<unsigned __int128>(__urng,
>- __u64erange);
>- }
>- else
>-#endif
>- if constexpr (__urngrange == __UINT32_MAX__)
>- {
>- // __urng produces values that use exactly 32-bits,
>+ // __urng produces values that use less or equal to 32-bits,
"no more than 32 bits,"
> // so use 64-bit integers to downscale to desired range.
> __UINT32_TYPE__ __u32erange = __uerange;
>- __ret = _S_nd<__UINT64_TYPE__>(__urng, __u32erange);
>+ __ret = _S_nd<__bits, __UINT64_TYPE__>(__urng, __u32erange);
>+ }
>+# if __SIZEOF_INT128__
>+ else if constexpr (__detail::_Power_of_2(__urngrange + 1) && __bits <= 64)
>+ {
>+ // __urng produces values that use less or equal to 64-bits,
"no more than 64 bits,"
>+ // so use 128-bit integers to downscale to desired range.
>+ __UINT64_TYPE__ __u64erange = __uerange;
>+ __ret = __extension__ _S_nd<__bits, unsigned __int128>(
>+ __urng, __u64erange);
> }
>+# endif
> else
> #endif
> {
>diff --git a/libstdc++-v3/testsuite/26_numerics/random/uniform_int_distribution/operators/values.cc b/libstdc++-v3/testsuite/26_numerics/random/uniform_int_distribution/operators/values.cc
>index fce14944d88..ebb02aa4912 100644
>--- a/libstdc++-v3/testsuite/26_numerics/random/uniform_int_distribution/operators/values.cc
>+++ b/libstdc++-v3/testsuite/26_numerics/random/uniform_int_distribution/operators/values.cc
>@@ -50,27 +50,35 @@ void test01()
> {
> using namespace __gnu_test;
>
>- std::mt19937 eng;
>+ std::mt19937 mteng;
>
> std::uniform_int_distribution<> uid1(0, 2);
>- auto buid1 = std::bind(uid1, eng);
>+ auto buid1 = std::bind(uid1, mteng);
> testDiscreteDist<ARGS>(buid1, [](int n) { return uniform_int_pdf(n, 0, 2); } );
>
> std::uniform_int_distribution<> uid2(3, 7);
>- auto buid2 = std::bind(uid2, eng);
>+ auto buid2 = std::bind(uid2, mteng);
> testDiscreteDist<ARGS>(buid2, [](int n) { return uniform_int_pdf(n, 3, 7); } );
>
> std::uniform_int_distribution<> uid3(1, 20);
>- auto buid3 = std::bind(uid3, eng);
>+ auto buid3 = std::bind(uid3, mteng);
> testDiscreteDist<ARGS>(buid3, [](int n) { return uniform_int_pdf(n, 1, 20); } );
>
>- shifted<(std::uint64_t(1) << 16)> s16e;
>- auto buid4 = std::bind(uid3, s16e);
>+ std::ranlux24 r24eng;
>+ auto buid4 = std::bind(uid3, r24eng);
> testDiscreteDist<ARGS>(buid4, [](int n) { return uniform_int_pdf(n, 1, 20); } );
>
>- shifted<(std::uint64_t(1) << 32)> s32e;
>- auto buid5 = std::bind(uid3, s32e);
>+ std::ranlux48 r48eng;
>+ auto buid5 = std::bind(uid3, r48eng);
> testDiscreteDist<ARGS>(buid5, [](int n) { return uniform_int_pdf(n, 1, 20); } );
>+
>+ shifted<(std::uint64_t(1) << 16)> s16e;
>+ auto buid6 = std::bind(uid3, s16e);
>+ testDiscreteDist<ARGS>(buid6, [](int n) { return uniform_int_pdf(n, 1, 20); } );
>+
>+ shifted<(std::uint64_t(1) << 32)> s32e;
>+ auto buid7 = std::bind(uid3, s32e);
>+ testDiscreteDist<ARGS>(buid7, [](int n) { return uniform_int_pdf(n, 1, 20); } );
> }
>
> int main()
>--
>2.55.0
>
>
More information about the Libstdc++
mailing list