Fix ibeta_inv for very small p.

include/boost/math/special_functions/detail/ibeta_inverse.hpp

@@ -118,9 +118,17 @@ T temme_method_1_ibeta_inverse(T a, T b, T z, const Policy& pol)
       x = 0.5;
    else
       x = (1 + eta * sqrt((1 + c) / eta_2)) / 2;
-
-   BOOST_MATH_ASSERT(x >= 0);
-   BOOST_MATH_ASSERT(x <= 1);
+   //
+   // These are post-conditions of the method, but the addition above
+   // may result in us being out by 1ulp either side of the boundary,
+   // so just check that we're in bounds and adjust as needed.
+   // See https://github.com/boostorg/math/issues/961
+   //
+   if (x < 0)
+      x = 0;
+   else if (x > 1)
+      x = 1;
+
    BOOST_MATH_ASSERT(eta * (x - 0.5) >= 0);
 #ifdef BOOST_INSTRUMENT
    std::cout << "Estimating x with Temme method 1: " << x << std::endl;
@@ -901,6 +909,8 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
       if(x < lower)
          x = lower;
    }
+   std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
+   std::uintmax_t max_iter_used = 0;
    //
    // Figure out how many digits to iterate towards:
    //
@@ -923,10 +933,9 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
    // Now iterate, we can use either p or q as the target here
    // depending on which is smaller:
    //
-   std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
    x = boost::math::tools::halley_iterate(
       boost::math::detail::ibeta_roots<T, Policy>(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits, max_iter);
-   policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter, pol);
+   policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter + max_iter_used, pol);
    //
    // We don't really want these asserts here, but they are useful for sanity
    // checking that we have the limits right, uncomment if you suspect bugs *only*.

include/boost/math/tools/roots.hpp

@@ -586,8 +586,8 @@ namespace detail {
                   // we can jump way out of bounds if we're not careful.
                   // See https://svn.boost.org/trac/boost/ticket/8314.
                   delta = f0 / f1;
-                  if (fabs(delta) > 2 * fabs(guess))
-                     delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess);
+                  if (fabs(delta) > 2 * fabs(result))
+                     delta = (delta < 0 ? -1 : 1) * 2 * fabs(result);
                }
             }
             else
@@ -600,9 +600,19 @@ namespace detail {
          if ((convergence > 0.8) && (convergence < 2))
          {
             // last two steps haven't converged.
-            delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
-            if ((result != 0) && (fabs(delta) > result))
-               delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps!
+            if (fabs(min) < 1 ? fabs(1000 * min) < fabs(max) : fabs(max / min) > 1000)
+            {
+               if(delta > 0)
+                  delta = bracket_root_towards_min(f, result, f0, min, max, count);
+               else
+                  delta = bracket_root_towards_max(f, result, f0, min, max, count);
+            }
+            else
+            {
+               delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
+               if ((result != 0) && (fabs(delta) > result))
+                  delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps!
+            }
             // reset delta2 so that this branch will *not* be taken on the
             // next iteration:
             delta2 = delta * 3;
@@ -641,6 +651,8 @@ namespace detail {
                result = guess - delta;
                if (result <= min)
                   result = float_next(min);
+               if (result >= max)
+                  result = float_prior(max);
                guess = min;
                continue;
             }
@@ -669,6 +681,8 @@ namespace detail {
                result = guess - delta;
                if (result >= max)
                   result = float_prior(max);
+               if (result <= min)
+                  result = float_next(min);
                guess = min;
                continue;
             }

test/Jamfile.v2

@@ -168,6 +168,7 @@ test-suite special_fun :
    [ run git_issue_705.cpp ../../test/build//boost_unit_test_framework  ]
    [ run git_issue_810.cpp ../../test/build//boost_unit_test_framework  ]
    [ run git_issue_826.cpp ../../test/build//boost_unit_test_framework  ]
+   [ run git_issue_961.cpp ]
    [ run special_functions_test.cpp ../../test/build//boost_unit_test_framework  ]
    [ run test_airy.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework  ]
    [ run test_bessel_j.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework  ]

test/git_issue_961.cpp

@@ -0,0 +1,76 @@
+//  (C) Copyright John Maddock 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 <iostream>
+#include <iomanip>
+#include <cfenv>
+#include <boost/math/special_functions/beta.hpp>
+
+using namespace std;
+using namespace boost::math;
+
+void show_fp_exception_flags()
+{
+   if (std::fetestexcept(FE_DIVBYZERO)) {
+      cout << " FE_DIVBYZERO";
+   }
+   // FE_INEXACT is common and not interesting.
+   // if (std::fetestexcept(FE_INEXACT)) {
+   //     cout << " FE_INEXACT";
+   // }
+   if (std::fetestexcept(FE_INVALID)) {
+      cout << " FE_INVALID";
+   }
+   if (std::fetestexcept(FE_OVERFLOW)) {
+      cout << " FE_OVERFLOW";
+   }
+   if (std::fetestexcept(FE_UNDERFLOW)) {
+      cout << " FE_UNDERFLOW";
+   }
+   cout << endl;
+}
+
+template <class Policy>
+int test()
+{
+   double a = 14.208308325339239;
+   double b = a;
+
+   double p = 6.4898872103239473e-300;  // Throws exception: Assertion `x >= 0' failed.
+   // double p = 7.8e-307;  // No flags set, returns 8.57354094063444939e-23
+   // double p = 7.7e-307;  // FE_UNDERFLOW set, returns 0.0
+
+   while (p > (std::numeric_limits<double>::min)())
+   {
+      std::feclearexcept(FE_ALL_EXCEPT);
+
+      try {
+
+         double x = ibeta_inv(a, b, p, Policy());
+
+         show_fp_exception_flags();
+
+         std::cout << std::scientific << std::setw(24)
+            << std::setprecision(17) << x << std::endl;
+      }
+      catch (const std::exception& e)
+      {
+         std::cout << e.what() << std::endl;
+         return 1;
+      }
+
+      p /= 1.25;
+   }
+
+   return 0;
+}
+
+int main(int argc, char* argv[])
+{
+   using namespace boost::math::policies;
+   if (test<policy<>>())
+      return 1;
+   return test<policy<promote_double<false>>>();
+}