From fbb4398683e093afc01a05e6250066d20ffe4482 Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Fri, 11 Nov 2022 19:36:24 +0000
Subject: [PATCH 1/2] Correct root finders in extreme cases. Fixes:
 https://github.com/boostorg/math/issues/873

---
 include/boost/math/tools/roots.hpp |  4 ++++
 test/test_ibeta_inv.hpp            | 18 ++++++++++++++++++
 2 files changed, 22 insertions(+)

diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp
index 0f2129d1f..19137481c 100644
--- a/include/boost/math/tools/roots.hpp
+++ b/include/boost/math/tools/roots.hpp
@@ -615,6 +615,8 @@ namespace detail {
                }
                delta = bracket_root_towards_min(f, guess, f0, min, max, count);
                result = guess - delta;
+               if (result <= min)
+                  result = float_next(min);
                guess = min;
                continue;
             }
@@ -641,6 +643,8 @@ namespace detail {
                }
                delta = bracket_root_towards_max(f, guess, f0, min, max, count);
                result = guess - delta;
+               if (result >= max)
+                  result = float_prior(max);
                guess = min;
                continue;
             }
diff --git a/test/test_ibeta_inv.hpp b/test/test_ibeta_inv.hpp
index 6787fdb93..5eea61dee 100644
--- a/test/test_ibeta_inv.hpp
+++ b/test/test_ibeta_inv.hpp
@@ -270,5 +270,23 @@ void test_spots(T)
          static_cast<T>(0.125),
          static_cast<T>(0.125)),
       static_cast<T>(0.99999994039535522460937500000000000000000000000L), tolerance);
+   //
+   // Bug cases, issue 873:
+   //
+   if ((std::numeric_limits<T>::max)() > static_cast<T>(1e50))
+   {
+      BOOST_CHECK_CLOSE(
+         ::boost::math::ibeta_inv(
+            static_cast<T>(1e50L),
+            static_cast<T>(10),
+            static_cast<T>(1) / static_cast<T>(10)),
+         static_cast<T>(1), tolerance);
+      BOOST_CHECK_CLOSE(
+         ::boost::math::ibetac_inv(
+            static_cast<T>(1e50L),
+            static_cast<T>(10),
+            static_cast<T>(1) / static_cast<T>(10)),
+         static_cast<T>(1), tolerance);
+   }
 }
 

From 9a9bb9ddc81001a317093c2dd76d04db18b741d6 Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Sun, 13 Nov 2022 17:08:19 +0000
Subject: [PATCH 2/2] Update root bracketing to home in faster. Special case
 for when the exponent range is so large that it basically takes "forever" to
 bracket otherwise.

---
 include/boost/math/tools/roots.hpp | 24 ++++++++++++++++++------
 1 file changed, 18 insertions(+), 6 deletions(-)

diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp
index 19137481c..cb2c89123 100644
--- a/include/boost/math/tools/roots.hpp
+++ b/include/boost/math/tools/roots.hpp
@@ -372,13 +372,19 @@ namespace detail {
    T bracket_root_towards_max(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
    {
       using std::fabs;
+      using std::ldexp;
+      using std::abs;
+      using std::frexp;
       if(count < 2)
          return guess - (max + min) / 2; // Not enough counts left to do anything!!
       //
       // Move guess towards max until we bracket the root, updating min and max as we go:
       //
+      int e;
+      frexp(max / guess, &e);
+      e = abs(e);
       T guess0 = guess;
-      T multiplier = 2;
+      T multiplier = e < 64 ? static_cast<T>(2) : static_cast<T>(ldexp(T(1), e / 32));
       T f_current = f0;
       if (fabs(min) < fabs(max))
       {
@@ -392,7 +398,7 @@ namespace detail {
                f_current = -f_current;  // There must be a change of sign!
                break;
             }
-            multiplier *= 2;
+            multiplier *= e > 1024 ? 8 : 2;
             unpack_0(f(guess), f_current);
          }
       }
@@ -411,7 +417,7 @@ namespace detail {
                f_current = -f_current;  // There must be a change of sign!
                break;
             }
-            multiplier *= 2;
+            multiplier *= e > 1024 ? 8 : 2;
             unpack_0(f(guess), f_current);
          }
       }
@@ -429,13 +435,19 @@ namespace detail {
    T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
    {
       using std::fabs;
+      using std::ldexp;
+      using std::abs;
+      using std::frexp;
       if (count < 2)
          return guess - (max + min) / 2; // Not enough counts left to do anything!!
       //
       // Move guess towards min until we bracket the root, updating min and max as we go:
       //
+      int e;
+      frexp(guess / min, &e);
+      e = abs(e);
       T guess0 = guess;
-      T multiplier = 2;
+      T multiplier = e < 64 ? static_cast<T>(2) : static_cast<T>(ldexp(T(1), e / 32));
       T f_current = f0;
 
       if (fabs(min) < fabs(max))
@@ -450,7 +462,7 @@ namespace detail {
                f_current = -f_current;  // There must be a change of sign!
                break;
             }
-            multiplier *= 2;
+            multiplier *= e > 1024 ? 8 : 2;
             unpack_0(f(guess), f_current);
          }
       }
@@ -469,7 +481,7 @@ namespace detail {
                f_current = -f_current;  // There must be a change of sign!
                break;
             }
-            multiplier *= 2;
+            multiplier *= e > 1024 ? 8 : 2;
             unpack_0(f(guess), f_current);
          }
       }