From 70162dbcb8736123726ab76e5b5c3f88aa751de0 Mon Sep 17 00:00:00 2001
From: jzmaddock <john@johnmaddock.co.uk>
Date: Sat, 25 May 2019 18:18:20 +0100
Subject: [PATCH] Update second order root finding code to bracket root when
 required. When min and max differ by many orders of magnitude (as happens
 when one is zero for example), then bisection can take a very long time to
 iterate down to the root.  Instead use a bracketing strategy which doubles
 the step size with each iteration until a bracket is found, then repeat
 recursively as required until we have a reasonably small interval.  Note that
 this only kicks in when a Halley step goes out of bounds and we're therefore
 forced to thrash around looking for the root.  Fixes:
 https://github.com/boostorg/math/issues/204.

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

diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp
index e16294dc1..4cc92b717 100644
--- a/include/boost/math/tools/roots.hpp
+++ b/include/boost/math/tools/roots.hpp
@@ -32,6 +32,7 @@
 #endif
 
 #include <boost/math/special_functions/sign.hpp>
+#include <boost/math/special_functions/next.hpp>
 #include <boost/math/tools/toms748_solve.hpp>
 #include <boost/math/policies/error_handling.hpp>
 
@@ -354,6 +355,118 @@ namespace detail{
       }
    };
 
+   template <class F, class T>
+   T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, boost::uintmax_t& count);
+
+   template <class F, class T>
+   T bracket_root_towards_max(F f, T guess, const T& f0, T& min, T& max, boost::uintmax_t& count)
+   {
+      using std::fabs;
+      using boost::math::get;
+      //
+      // Move guess towards max until we bracket the root, updating min and max as we go:
+      //
+      T guess0 = guess;
+      T multiplier = 2;
+      T f_current = f0;
+      if (fabs(min) < fabs(max))
+      {
+         while (--count && ((f_current < 0) == (f0 < 0)))
+         {
+            min = guess;
+            guess *= multiplier;
+            if (guess > max)
+            {
+               guess = max;
+               break;
+            }
+            multiplier *= 2;
+            f_current = get<0>(f(guess));
+         }
+      }
+      else
+      {
+         //
+         // If min and max are negative we have to divide to head towards max:
+         //
+         while (--count && ((f_current < 0) == (f0 < 0)))
+         {
+            min = guess;
+            guess /= multiplier;
+            if (guess > max)
+            {
+               guess = max;
+               break;
+            }
+            multiplier *= 2;
+            f_current = get<0>(f(guess));
+         }
+      }
+
+      if (count)
+      {
+         max = guess;
+         if (multiplier > 16)
+            return (guess0 - guess) + bracket_root_towards_min(f, guess, f_current, min, max, count);
+      }
+      return guess0 - (max + min) / 2;
+   }
+
+   template <class F, class T>
+   T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, boost::uintmax_t& count)
+   {
+      using std::fabs;
+      using boost::math::get;
+      //
+      // Move guess towards min until we bracket the root, updating min and max as we go:
+      //
+      T guess0 = guess;
+      T multiplier = 2;
+      T f_current = f0;
+
+      if (fabs(min) < fabs(max))
+      {
+         while (--count && ((f_current < 0) == (f0 < 0)))
+         {
+            max = guess;
+            guess /= multiplier;
+            if (guess < min)
+            {
+               guess = min;
+               break;
+            }
+            multiplier *= 2;
+            f_current = get<0>(f(guess));
+         }
+      }
+      else
+      {
+         //
+         // If min and max are negative we have to multiply to head towards min:
+         //
+         while (--count && ((f_current < 0) == (f0 < 0)))
+         {
+            max = guess;
+            guess *= multiplier;
+            if (guess < min)
+            {
+               guess = min;
+               break;
+            }
+            multiplier *= 2;
+            f_current = get<0>(f(guess));
+         }
+      }
+
+      if (count)
+      {
+         min = guess;
+         if (multiplier > 16)
+            return (guess0 - guess) + bracket_root_towards_max(f, guess, f_current, min, max, count);
+      }
+      return guess0 - (max + min) / 2;
+   }
+
    template <class Stepper, class F, class T>
    T second_order_root_finder(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
    {
@@ -459,10 +572,15 @@ namespace detail{
             }
             else
             {
-               delta = (guess - min) / 2;
-               result = guess - delta;
-               if((result == min) || (result == max))
+               if (fabs(float_distance(min, max)) < 2)
+               {
+                  result = guess = (min + max) / 2;
                   break;
+               }
+               delta = bracket_root_towards_min(f, guess, f0, min, max, count);
+               result = guess - delta;
+               guess = min;
+               continue;
             }
          }
          else if(result > max)
@@ -480,10 +598,15 @@ namespace detail{
             }
             else
             {
-               delta = (guess - max) / 2;
-               result = guess - delta;
-               if((result == min) || (result == max))
+               if (fabs(float_distance(min, max)) < 2)
+               {
+                  result = guess = (min + max) / 2;
                   break;
+               }
+               delta = bracket_root_towards_max(f, guess, f0, min, max, count);
+               result = guess - delta;
+               guess = min;
+               continue;
             }
          }
          // update brackets: