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
