<random>: Improve binomial_distribution accuracy

stl/inc/random

@@ -2428,7 +2428,7 @@ public:
         _Ty1 _Logp;
         _Ty1 _Logp1;
 
-        _Small_poisson_distribution<_Ty> _Small;
+        _Small_poisson_distribution<_Ty> _Small; // TRANSITION, ABI: unused
     };
 
     binomial_distribution() : _Par(1, _Ty1(0.5)) {}
@@ -2505,8 +2505,26 @@ private:
 
             return _Res;
         } else if (_Par0._Mean < 1.0) {
-            // events are rare, use Poisson distribution
-            _Res = _Par0._Small(_Eng);
+            // Events are rare, use waiting time method (Luc Devroye, Non-Uniform Random Variate Generation, p. 525).
+            const _Ty1 _Rand = _NRAND(_Eng, _Ty1);
+
+            // The exit condition is log(1 - _Rand)/t < log(1-p), which is equivalent to _Rand > 1 - (1-p)^t. If
+            // we have a cheap upper bound for 1-(1-p)^t, we can exit early without having to call log. We use two
+            // such bounds, one that is tight for mean ~0 and another for mean ~1. In the first case, Bernoulli's
+            // inequality gives -1+p*t >= -(1-p)^t, so 1 - (1-p)^t <= p*t = mean. For the other bound, 1-(1-p)^t =
+            // 1-(1-p)(1-mean/t)^(t-1) <= 1-(1-p)(1-1/t)^(t-1) <= 1-(1-p)/e.
+            const _Ty1 _Ub =
+                (_STD min)(_Par0._Mean, _Ty1{3.678794411714423216e-1} * _Par0._Pp + _Ty1{6.32120558828557678e-1});
+            if (_Rand > _Ub) {
+                _Res = _Ty{0};
+            } else {
+                _Ty _Denom = _Par0._Tx;
+                _Ty1 _Sum  = _CSTD log(_Ty1{1.0} - _Rand) / _Denom;
+                while (_Sum >= _Par0._Logp1 && --_Denom != 0) {
+                    _Sum += _CSTD log(_Ty1{1.0} - _NRAND(_Eng, _Ty1)) / _Denom;
+                }
+                _Res = static_cast<_Ty>(_Par0._Tx - _Denom);
+            }
         } else { // no shortcuts
             using _Uty = make_unsigned_t<_Ty>;
             const auto _Ty1_Tx{_Float_upper_bound<_Ty1>(static_cast<_Uty>(_Par0._Tx))};

tests/libcxx/expected_results.txt

@@ -483,9 +483,6 @@ std/algorithms/robust_against_adl.pass.cpp FAIL
 # GH-1595: <bit>: bit_ceil(T(-1)) should not be a constant expression
 std/numerics/bit/bit.pow.two/bit_ceil.fail.cpp:0 FAIL
 
-# GH-1530: <random>: Poisson approximation for binomial_distribution is not accurate
-std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.PR44847.pass.cpp FAIL
-
 
 # *** CRT BUGS ***
 # We're permanently missing aligned_alloc().
@@ -669,6 +666,7 @@ std/input.output/filesystems/fs.filesystem.synopsis/file_time_type_resolution.co
 # Not yet analyzed, likely STL bugs. Assertions and other runtime failures.
 std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp FAIL
 std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp FAIL
+std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.PR44847.pass.cpp FAIL
 std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.geo/eval_param.pass.cpp FAIL
 std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.geo/eval.pass.cpp FAIL
 std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.negbin/eval_param.pass.cpp FAIL

tests/libcxx/skipped_tests.txt

@@ -483,9 +483,6 @@ algorithms\robust_against_adl.pass.cpp
 # GH-1595: <bit>: bit_ceil(T(-1)) should not be a constant expression
 numerics\bit\bit.pow.two\bit_ceil.fail.cpp
 
-# GH-1530: <random>: Poisson approximation for binomial_distribution is not accurate
-numerics\rand\rand.dis\rand.dist.bern\rand.dist.bern.bin\eval.PR44847.pass.cpp
-
 
 # *** CRT BUGS ***
 # We're permanently missing aligned_alloc().
@@ -669,6 +666,7 @@ input.output\filesystems\fs.filesystem.synopsis\file_time_type_resolution.compil
 # Not yet analyzed, likely STL bugs. Assertions and other runtime failures.
 numerics\rand\rand.dis\rand.dist.bern\rand.dist.bern.bin\eval_param.pass.cpp
 numerics\rand\rand.dis\rand.dist.bern\rand.dist.bern.bin\eval.pass.cpp
+numerics\rand\rand.dis\rand.dist.bern\rand.dist.bern.bin\eval.PR44847.pass.cpp
 numerics\rand\rand.dis\rand.dist.bern\rand.dist.bern.geo\eval_param.pass.cpp
 numerics\rand\rand.dis\rand.dist.bern\rand.dist.bern.geo\eval.pass.cpp
 numerics\rand\rand.dis\rand.dist.bern\rand.dist.bern.negbin\eval_param.pass.cpp

tests/std/test.lst

@@ -172,6 +172,7 @@ tests\GH_001103_countl_zero_correctness
 tests\GH_001105_custom_streambuf_throws
 tests\GH_001123_random_cast_out_of_range
 tests\GH_001411_core_headers
+tests\GH_001530_binomial_accuracy
 tests\GH_001541_case_sensitive_boolalpha
 tests\LWG2597_complex_branch_cut
 tests\LWG3018_shared_ptr_function

tests/std/tests/GH_001530_binomial_accuracy/env.lst

@@ -0,0 +1,4 @@
+# Copyright (c) Microsoft Corporation.
+# SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+
+RUNALL_INCLUDE ..\usual_matrix.lst

tests/std/tests/GH_001530_binomial_accuracy/test.cpp

@@ -0,0 +1,48 @@
+// Copyright (c) Microsoft Corporation.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+
+#include <cassert>
+#include <cmath>
+#include <cstddef>
+#include <random>
+#include <vector>
+
+using namespace std;
+
+template <class Generator>
+void test_binomial(const int n, const double mean, Generator& gen) {
+    const double p       = mean / n;
+    const double var     = n * p * (1.0 - p);
+    constexpr double tol = 0.01;
+    constexpr double x   = 4.0; // Standard deviation of sample variance should be less than tol / x.
+    const size_t it_max  = 2 * static_cast<size_t>(pow(var / (tol / x), 2.0));
+    binomial_distribution<> dist(n, p);
+
+    vector<size_t> counts(static_cast<size_t>(n) + 1);
+    for (size_t i = 0; i < it_max; ++i) {
+        ++counts[static_cast<size_t>(dist(gen))];
+    }
+
+    double sample_mean = 0.0;
+    for (size_t i = 1; i < counts.size(); ++i) {
+        sample_mean += static_cast<double>(i) * counts[i];
+    }
+    sample_mean /= it_max;
+
+    double sample_var = 0.0;
+    for (size_t i = 0; i < counts.size(); ++i) {
+        sample_var += counts[i] * pow(i - sample_mean, 2);
+    }
+    sample_var /= it_max - 1;
+
+    assert(abs(sample_mean / mean - 1.0) < tol);
+    assert(abs(sample_var / var - 1.0) < tol);
+}
+
+int main() {
+    mt19937 gen;
+    constexpr int n = 25;
+    test_binomial(n, 0.99, gen);
+    test_binomial(n, n - 0.99, gen);
+    return 0;
+}