From cb068993561cb143da401d5cf4ad52efc4b6ea57 Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Sat, 7 Sep 2024 17:36:39 +0100
Subject: [PATCH] Catch evaluation_error's in temme_method_2_ibeta_inverse
 (#1175)

Catch evaluation_error's in temme_method_2_ibeta_inverse
fixes https://github.com/boostorg/math/issues/1169.
Add test case.
---
 .../detail/ibeta_inverse.hpp                  | 21 +++++++++++++---
 test/Jamfile.v2                               |  1 +
 test/git_issue_1175.cpp                       | 25 +++++++++++++++++++
 3 files changed, 44 insertions(+), 3 deletions(-)
 create mode 100644 test/git_issue_1175.cpp

diff --git a/include/boost/math/special_functions/detail/ibeta_inverse.hpp b/include/boost/math/special_functions/detail/ibeta_inverse.hpp
index 90f6e9070..6f222cf77 100644
--- a/include/boost/math/special_functions/detail/ibeta_inverse.hpp
+++ b/include/boost/math/special_functions/detail/ibeta_inverse.hpp
@@ -304,9 +304,23 @@ BOOST_MATH_GPU_ENABLED T temme_method_2_ibeta_inverse(T /*a*/, T /*b*/, T z, T r
    //
    // And iterate:
    //
-   x = tools::newton_raphson_iterate(
-      temme_root_finder<T>(-lu, alpha), x, lower, upper, policies::digits<T, Policy>() / 2);
-
+#ifndef BOOST_MATH_NO_EXCEPTIONS
+   try {
+#endif
+      x = tools::newton_raphson_iterate(
+         temme_root_finder<T>(-lu, alpha), x, lower, upper, policies::digits<T, Policy>() / 2);
+#ifndef BOOST_MATH_NO_EXCEPTIONS
+   }
+   catch (const boost::math::evaluation_error&)
+   {
+      // Due to numerical instability we may have cases where no root is found when
+      // in fact we should just touch the origin.  We simply ignore the error here
+      // and return our best guess for x so far...
+      // Maybe we should special case the symmetrical parameter case, but it's not clear 
+      // whether that is the only situation when problems can occur.
+      // See https://github.com/boostorg/math/issues/1169
+   }
+#endif
    return x;
 }
 //
@@ -321,6 +335,7 @@ BOOST_MATH_GPU_ENABLED T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const
 {
    BOOST_MATH_STD_USING // ADL of std names
 
+
    //
    // Begin by getting an initial approximation for the quantity
    // eta from the dominant part of the incomplete beta:
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index e6814fd86..e19957777 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -192,6 +192,7 @@ test-suite special_fun :
    [ run git_issue_184.cpp ]
    [ run git_issue_1137.cpp ]
    [ run git_issue_1139.cpp ]
+   [ run git_issue_1175.cpp ]
    [ run special_functions_test.cpp /boost/test//boost_unit_test_framework  ]
    [ run test_airy.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework  ]
    [ run test_bessel_j.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework  ]
diff --git a/test/git_issue_1175.cpp b/test/git_issue_1175.cpp
new file mode 100644
index 000000000..9770acf53
--- /dev/null
+++ b/test/git_issue_1175.cpp
@@ -0,0 +1,25 @@
+//  (C) Copyright Matt Borland 2023.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include "math_unit_test.hpp"
+#include <iostream>
+#include <boost/math/distributions/beta.hpp>
+
+using namespace std;
+using boost::math::beta_distribution;
+
+int main(int argc, char* argv[])
+{
+   double a = 5.0;
+   double b = 5.0;
+   double p = 0.5;
+
+   beta_distribution<> dist(a, b);
+   double x = quantile(dist, p);
+
+   CHECK_ULP_CLOSE(x, 0.5, 2);
+
+   return boost::math::test::report_errors();
+}