[PATCH v2] libstdc++: Prevent hanging of uniform_int_distribution for non-conforming engines [PR118665]

[Tomasz Kamiński]

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 returned 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 some 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 inisde 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 from either or both size 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 front and Y values from 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  -

In consequence, of above 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.

	PR libstdc++/118665

libstdc++-v3/ChangeLog:

	* 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.
---
v2:
* expands the explanation in commit description on why the output
  remains uniform
* updates comment in code
* introduce __back_treshold variable
* applyes - __min from previous patch

OK for trunk? Backports?

 libstdc++-v3/include/bits/uniform_int_dist.h  | 13 ++++++++-
 .../operators/pr118665.cc                     | 28 +++++++++++++++++++
 2 files changed, 40 insertions(+), 1 deletion(-)
 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);
+}
-- 
2.55.0