diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp
index 574dd3021..f0212800c 100644
--- a/include/boost/math/special_functions/beta.hpp
+++ b/include/boost/math/special_functions/beta.hpp
@@ -373,7 +373,7 @@ T ibeta_power_terms(T a,
          if(a < b)
          {
             T p1 = pow(b2, b / a);
-            T l3 = a * (log(b1) + log(p1));
+            T l3 = (b1 != 0) && (p1 != 0) ? (a * (log(b1) + log(p1))) : tools::max_value<T>(); // arbitrary large value if the logs would fail!
             if((l3 < tools::log_max_value<T>())
                && (l3 > tools::log_min_value<T>()))
             {
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index 36f56d56d..ad373bd39 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -894,6 +894,7 @@ test-suite distribution_tests :
    [ run git_issue_800.cpp ../../test/build//boost_unit_test_framework  ]
    [ run git_issue_845.cpp ../../test/build//boost_unit_test_framework  ]
    [ run scipy_issue_14901.cpp ../../test/build//boost_unit_test_framework ]
+   [ run scipy_issue_15101.cpp ../../test/build//boost_unit_test_framework ]
    [ run scipy_issue_17146.cpp ../../test/build//boost_unit_test_framework  ]
    [ run scipy_issue_17388.cpp ../../test/build//boost_unit_test_framework  ]
    [ run scipy_issue_17916.cpp ../../test/build//boost_unit_test_framework  ]
diff --git a/test/scipy_issue_15101.cpp b/test/scipy_issue_15101.cpp
new file mode 100644
index 000000000..4494aa01c
--- /dev/null
+++ b/test/scipy_issue_15101.cpp
@@ -0,0 +1,33 @@
+// 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)
+
+#define BOOST_MATH_DOMAIN_ERROR_POLICY ignore_error
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
+#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false
+
+#include <cfenv>
+#include <array>
+#include <numeric>
+#include <boost/math/distributions/binomial.hpp>
+#include "math_unit_test.hpp"
+
+#pragma STDC FENV_ACCESS ON
+
+int main()
+{
+    constexpr double n {2000};
+    constexpr double p {0.9999};
+    const auto binom_dist = boost::math::binomial_distribution<double>(n, p);
+    const auto pdf1 = boost::math::pdf(binom_dist, 3);
+
+    CHECK_ULP_CLOSE(pdf1, 0, 1);
+
+    if (std::fetestexcept(FE_DIVBYZERO)) 
+    {
+        return 1;
+    }
+
+    return boost::math::test::report_errors();
+}