diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp
index 0f2129d1f..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);
          }
       }
@@ -615,6 +627,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 +655,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);
+   }
 }