Fix for scipy issue 18302

include/boost/math/distributions/beta.hpp

@@ -396,7 +396,25 @@ namespace boost
       {
         return RealType(0);
       }
-      return ibeta_derivative(a, b, x, Policy());
+
+      RealType return_val;
+      try
+      {
+        return_val = ibeta_derivative(a, b, x, Policy());
+      }
+      catch (const std::overflow_error&)
+      {
+        // Overflows with very small x which yields a pdf approximately equal to beta
+        return_val = b;
+      }
+
+      // Bounds check
+      if (return_val > b)
+      {
+        return_val = b;
+      }
+
+      return return_val;
     } // pdf
 
     template <class RealType, class Policy>

test/Jamfile.v2

@@ -899,6 +899,7 @@ test-suite distribution_tests :
    [ run scipy_issue_17388.cpp ../../test/build//boost_unit_test_framework  ]
    [ run scipy_issue_17916.cpp ../../test/build//boost_unit_test_framework  ]
    [ run scipy_issue_17916_nct.cpp ../../test/build//boost_unit_test_framework  ]
+   [ run scipy_issue_18302.cpp ../../test/build//boost_unit_test_framework  ]
 ;
 
 test-suite mp :

test/scipy_issue_18302.cpp

@@ -0,0 +1,29 @@
+// Copyright Matt Borland, 2023
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
+//
+// See: https://github.com/scipy/scipy/issues/17916
+
+#include <boost/math/distributions/beta.hpp>
+#include <limits>
+#include "math_unit_test.hpp"
+
+int main(void)
+{
+    auto dist = boost::math::beta_distribution<double>(1, 5);
+
+    // https://www.wolframalpha.com/input?i=PDF%28beta+distribution%281%2C+5%29%2C+0%29
+    double test_pdf_spot = boost::math::pdf(dist, 0);
+    CHECK_ULP_CLOSE(test_pdf_spot, 0.0, 1);
+
+    // https://www.wolframalpha.com/input?i=PDF%28beta+distribution%281%2C+5%29%2C+1e-30%29
+    test_pdf_spot = boost::math::pdf(dist, 1e-30);
+    CHECK_ULP_CLOSE(test_pdf_spot, 5.0, 1);
+
+    // Appox equal to 5
+    test_pdf_spot = boost::math::pdf(dist, 1e-310);
+    CHECK_ULP_CLOSE(test_pdf_spot, 5.0, 1);
+
+    return boost::math::test::report_errors();
+}