Since r11-3757-g98c37d3bacbb2f the uniform_int_distribution uses Lemire's
algorithm for engines outputting power of two range. As this algorithm
(in-contrast to other) rejects low outputs from incoming generator, this
change lead to hangs (due infinite loop) when distribution was mixed
with non-uniform generator, that returns 0 constantly. Such usage has
undefined behavior according to standard, however didn't lead to hangs
prior to GCC 11 and other standard libraries.
This patch, addresses above by changing the (slower) rejection path of the
algorithm, to alternate between rejecting front/back value of the __prod % 2^N
range. This still produces an uniform output, as the range of preserved
values remain continuous and have same size.
To elaborate, we notice that the values of m (__product) are necessary
an multiplies of range size s (__range), and the shift __product >> N
distributes them along the bucket of size of 2^N. The produced distribution
is uniform if each 2^N contains the same number of multiplies of s, i.e.
we are equally likely to end inside it.
The shift itself is biased, for example given s == 5 and N == 3, we
get following buckets:
[0, 8) -> 0, 5
[8, 16) -> 10, 15
[16, 24) -> 20
[24, 32) -> 25, 30
[32, 40) -> 35
In that case, the 16 and 32 bucket (and corresponding 2 and 4 outputs)
are half less likely to occur.
To mitigate above we observe that any continuous range of size n * s
(multiply of s) contains exactly n multiplies of s, so by reducing the size
of the range from 2^N to continuous range of size that is multiply of s, we
will produce uniform distribution. This can be simply done by rejecting
2^N % s values in total from either or both sides of the range. It does not
matter from which side they are rejected.
As illustration the table bellow shows rejected value for s == 5, N == 3,
the columns corresponds to buckets and rows labels +X/-Y describe rejecting
X values from the front and Y values from the back.
0 8 16 24 32
+3/-0 0 10 - 25 -
+2/-1 0 15 - 25 -
+1/-2 0 15 - 30 -
+0/-3 5 15 - 30 -
With above changes, the algorithm finishes successfully for any generator
that returns same value constantly. While, it is possible to create an
generator that would still lead to infinite loop, by alternating between
low and high values, I believe this is inherent property on any algorithm
that rejects some of the engine's outputs.
libstdc++-v3/ChangeLog:
PR libstdc++/118665
* include/bits/uniform_int_dist.h
(uniform_int_distribution::_S_nd): Alternate between rejecting
from beginning/end of the range.
*
testsuite/26_numerics/random/uniform_int_distribution/operators/pr118665.cc:
New test.
*
testsuite/26_numerics/random/uniform_real_distribution/operators/gencanon.cc:
Updated expected values for test_2p55p1 that's depends on
uniform_int_distribution.
Reviewed-by: Jonathan Wakely <[email protected]>
Signed-off-by: Tomasz Kamiński <[email protected]>
---
v3:
* clears few remaining typos and wrong word use in commit description
* updates uniform_real_distribution/operators/gencanon.cc expected
values for test_2p55p1
Tested on x86_64-linux. Pushed to trunk.
libstdc++-v3/include/bits/uniform_int_dist.h | 13 ++++++++-
.../operators/pr118665.cc | 28 +++++++++++++++++++
.../operators/gencanon.cc | 18 ++++++------
3 files changed, 49 insertions(+), 10 deletions(-)
create mode 100644
libstdc++-v3/testsuite/26_numerics/random/uniform_int_distribution/operators/pr118665.cc
diff --git a/libstdc++-v3/include/bits/uniform_int_dist.h
b/libstdc++-v3/include/bits/uniform_int_dist.h
index 072c6351e1f..5ee9945b61c 100644
--- a/libstdc++-v3/include/bits/uniform_int_dist.h
+++ b/libstdc++-v3/include/bits/uniform_int_dist.h
@@ -271,9 +271,20 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
_Up __low = _Up(__product);
if (__low < __range)
{
- _Up __threshold = -__range % __range;
+ const _Up __threshold = -__range % __range;
while (__low < __threshold)
{
+ __product = _Wp(__g() - __min) * _Wp(__range);
+ __low = _Up(__product);
+
+ // 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;
+ if (__low <= __back_treshold)
+ break;
+
__product = _Wp(__g() - __min) * _Wp(__range);
__low = _Up(__product);
}
diff --git
a/libstdc++-v3/testsuite/26_numerics/random/uniform_int_distribution/operators/pr118665.cc
b/libstdc++-v3/testsuite/26_numerics/random/uniform_int_distribution/operators/pr118665.cc
new file mode 100644
index 00000000000..979d5fdfca9
--- /dev/null
+++
b/libstdc++-v3/testsuite/26_numerics/random/uniform_int_distribution/operators/pr118665.cc
@@ -0,0 +1,28 @@
+// { dg-do run { target c++11 } }
+// { dg-require-cstdint "" }
+
+#include <random>
+#include <testsuite_hooks.h>
+
+class Constant {
+public:
+ using result_type = unsigned;
+ static constexpr result_type min() { return 0u; }
+ static constexpr result_type max() { return ~0u; }
+
+ // always return 0
+ result_type operator()() {
+ ++calls;
+ VERIFY( calls < 1000 );
+ return 0;
+ }
+
+ unsigned calls = 0;
+};
+
+int main()
+{
+ Constant gen;
+ std::uniform_int_distribution<int> dist(0, 42);
+ dist(gen);
+}
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
index 775e25106d1..9d0161a3a6d 100644
---
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
@@ -310,24 +310,24 @@ test_2p55p1(const std::mt19937& rng)
case 24: // ieee32
case 53: // ieee64
VERIFY(skips == 0);
- VERIFY(deviation == 7682);
- VERIFY(max == 295);
- VERIFY(rms == 971);
+ VERIFY(deviation == 7636);
+ VERIFY(max == 292);
+ VERIFY(rms == 965);
VERIFY(zeros == 0);
break;
case 64: // ieee80
case 106: // ibm128
VERIFY(skips == 0);
- VERIFY(deviation == 6788);
- VERIFY(max == 248);
- VERIFY(rms == 854);
+ VERIFY(deviation == 7830);
+ VERIFY(max == 246);
+ VERIFY(rms == 963);
VERIFY(zeros == 0);
break;
case 113: // ieee128
VERIFY(skips == 0);
- VERIFY(deviation == 8222);
- VERIFY(max == 257);
- VERIFY(rms == 999);
+ VERIFY(deviation == 7930);
+ VERIFY(max == 271);
+ VERIFY(rms == 981);
VERIFY(zeros == 0);
break;
default:
--
2.55.0
