From 4f5ac6ff9c505ec0f28d45a377ac7d1101694a7a Mon Sep 17 00:00:00 2001
From: KRM7 <70973547+KRM7@users.noreply.github.com>
Date: Sat, 16 Sep 2023 13:24:37 +0200
Subject: [PATCH] improved randomBinomial implementation

---
 CMakeLists.txt      |  3 +--
 src/utility/rng.hpp | 37 ++++++++++++++++++++++++++-----------
 2 files changed, 27 insertions(+), 13 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 758f993e..06a2bc5a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -101,8 +101,7 @@ else() # GNU style compiler interface
 
     # gcc specific options
     if(CMAKE_CXX_COMPILER_ID MATCHES "GNU")
-        # https://gcc.gnu.org/bugzilla//show_bug.cgi?id=71434 with gcc 12 even though it's a system header + pragma diagnostic doesn't work
-	    set(GAPP_WARN_FLAGS "${GAPP_WARN_FLAGS} -Wlogical-op -Wno-maybe-uninitialized")
+	    set(GAPP_WARN_FLAGS "${GAPP_WARN_FLAGS} -Wlogical-op")
     endif()
     # clang specific options
     if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
diff --git a/src/utility/rng.hpp b/src/utility/rng.hpp
index 04287cfa..9b7e51a5 100644
--- a/src/utility/rng.hpp
+++ b/src/utility/rng.hpp
@@ -253,7 +253,7 @@ namespace gapp::rng
 
 
     /** The global pseudo-random number generator instance used in the algorithms. */
-    inline constinit ConcurrentXoroshiro128p prng;
+    GAPP_API inline constinit ConcurrentXoroshiro128p prng;
 
 
     /** Generate a random boolean value from a uniform distribution. */
@@ -275,7 +275,7 @@ namespace gapp::rng
     template<std::floating_point RealType = double>
     RealType randomNormal(RealType mean = 0.0, RealType std_dev = 1.0);
 
-    /** Generate a random integer value from a binomial distribution with the parameters n and p. */
+    /** Generate a random integer value from an approximate binomial distribution with the parameters n and p. */
     template<std::integral IntType = int>
     IntType randomBinomial(IntType n, double p);
 
@@ -367,33 +367,44 @@ namespace gapp::rng
     }
 
     template<std::integral IntType>
-    IntType randomBinomialApprox(IntType n, double p)
+    IntType randomBinomialApproxNormal(IntType n, double p)
     {
         GAPP_ASSERT(n >= 0);
         GAPP_ASSERT(0.0 <= p && p <= 1.0);
 
         const double mean = n * p;
-        const double SD = std::sqrt(mean * (1.0 - p));
+        const double std_dev = std::sqrt(mean * (1.0 - p));
 
         const double accept_min = -0.5;
         const double accept_max = n + 0.5;
 
-        double rand = rng::randomNormal(mean, SD);
+        double rand = rng::randomNormal(mean, std_dev);
         while (!(accept_min < rand && rand < accept_max))
         {
-            rand = rng::randomNormal(mean, SD);
+            rand = rng::randomNormal(mean, std_dev);
         }
 
         return static_cast<IntType>(std::round(rand));
     }
 
+    template<std::integral IntType>
+    IntType randomBinomialApproxPoisson(IntType n, double p)
+    {
+        GAPP_ASSERT(n >= 0);
+        GAPP_ASSERT(0.0 <= p && p <= 1.0);
+
+        return std::poisson_distribution<IntType>{ n * p }(rng::prng);
+    }
+
     template<std::integral IntType>
     IntType randomBinomialExact(IntType n, double p)
     {
         GAPP_ASSERT(n >= 0);
         GAPP_ASSERT(0.0 <= p && p <= 1.0);
 
-        return std::binomial_distribution{ n, p }(rng::prng);
+        IntType res = 0;
+        while (n--) { res += (randomReal() <= p); }
+        return res;
     }
 
     template<std::integral IntType>
@@ -402,11 +413,15 @@ namespace gapp::rng
         GAPP_ASSERT(n >= 0);
         GAPP_ASSERT(0.0 <= p && p <= 1.0);
 
-        const double mean = n * p;
+        const double mean = p * n;
+        const double invp = p / n;
+
+        const bool use_poisson = (n > 50 && mean < 10.0) || (n > 5 && invp < 0.2);
+        const bool use_normal = !use_poisson && (n > 5);
 
-        return (mean >= 2.0) ?
-            rng::randomBinomialApprox(n, p) :
-            rng::randomBinomialExact(n, p);
+        if (use_poisson) return rng::randomBinomialApproxPoisson(n, p);
+        if (use_normal)  return rng::randomBinomialApproxNormal(n, p);
+        return rng::randomBinomialExact(n, p);
     }
 
     template<detail::IndexableContainer T>