Avoid division by 0 in beta_distribution operator()

include/boost/random/beta_distribution.hpp

@@ -101,8 +101,15 @@ class beta_distribution {
     template<class URNG>
     RealType operator()(URNG& urng) const
     {
-        RealType a = gamma_distribution<RealType>(_alpha, RealType(1.0))(urng);
-        RealType b = gamma_distribution<RealType>(_beta, RealType(1.0))(urng);
+        RealType a = 0;
+        RealType b = 0;
+
+        do
+        {
+            a = gamma_distribution<RealType>(_alpha, RealType(1.0))(urng);
+            b = gamma_distribution<RealType>(_beta, RealType(1.0))(urng);
+        } while (a + b == RealType(0));
+
         return a / (a + b);
     }
 

test/Jamfile.v2

@@ -86,6 +86,8 @@ run test_comp_xoshiro128mm.cpp ;
 run test_xoshiro128f.cpp /boost/test//boost_unit_test_framework ;
 run test_comp_xoshiro128f.cpp ;
 
+run github_issue_133.cpp ;
+
 run niederreiter_base2_validate.cpp /boost/test//boost_unit_test_framework ;
 run sobol_validate.cpp /boost/test//boost_unit_test_framework ;
 run faure_validate.cpp /boost/test//boost_unit_test_framework ;

test/github_issue_133.cpp

@@ -0,0 +1,25 @@
+// Copyright 2025 Matt Borland
+// Distributed under the Boost Software License, Version 1.0.
+// https://www.boost.org/LICENSE_1_0.txt
+//
+// See: https://github.com/boostorg/random/issues/133
+
+#include <boost/core/lightweight_test.hpp>
+#include <boost/random/beta_distribution.hpp>
+#include <random>
+#include <cmath>
+
+int main()
+{
+    constexpr double beta_param = 0.0020368700639848774;
+    boost::random::beta_distribution<double> dist(beta_param, beta_param);
+    std::mt19937_64 gen(12345);
+
+    for (int i = 0; i < 10000; ++i)
+    {
+        const double Z = dist(gen);
+        BOOST_TEST(!std::isnan(Z));
+    }
+
+    return boost::report_errors();
+}