boost/math/tools/roots.hpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP
   7 #define BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12 #include <boost/math/tools/complex.hpp> // test for multiprecision types.
  13
  14 #include <iostream>
  15 #include <utility>
  16 #include <boost/config/no_tr1/cmath.hpp>
  17 #include <stdexcept>
  18
  19 #include <boost/math/tools/config.hpp>
  20 #include <boost/cstdint.hpp>
  21 #include <boost/assert.hpp>
  22 #include <boost/throw_exception.hpp>
  23
  24 #ifdef BOOST_MSVC
  25 #pragma warning(push)
  26 #pragma warning(disable: 4512)
  27 #endif
  28 #include <boost/math/tools/tuple.hpp>
  29 #ifdef BOOST_MSVC
  30 #pragma warning(pop)
  31 #endif
  32
  33 #include <boost/math/special_functions/sign.hpp>
  34 #include <boost/math/special_functions/next.hpp>
  35 #include <boost/math/tools/toms748_solve.hpp>
  36 #include <boost/math/policies/error_handling.hpp>
  37
  38 namespace boost {
  39 namespace math {
  40 namespace tools {
  41
  42 namespace detail {
  43
  44 namespace dummy {
  45
  46    template<int n, class T>
  47    typename T::value_type get(const T&) BOOST_MATH_NOEXCEPT(T);
  48 }
  49
  50 template <class Tuple, class T>
  51 void unpack_tuple(const Tuple& t, T& a, T& b) BOOST_MATH_NOEXCEPT(T)
  52 {
  53    using dummy::get;
  54    // Use ADL to find the right overload for get:
  55    a = get<0>(t);
  56    b = get<1>(t);
  57 }
  58 template <class Tuple, class T>
  59 void unpack_tuple(const Tuple& t, T& a, T& b, T& c) BOOST_MATH_NOEXCEPT(T)
  60 {
  61    using dummy::get;
  62    // Use ADL to find the right overload for get:
  63    a = get<0>(t);
  64    b = get<1>(t);
  65    c = get<2>(t);
  66 }
  67
  68 template <class Tuple, class T>
  69 inline void unpack_0(const Tuple& t, T& val) BOOST_MATH_NOEXCEPT(T)
  70 {
  71    using dummy::get;
  72    // Rely on ADL to find the correct overload of get:
  73    val = get<0>(t);
  74 }
  75
  76 template <class T, class U, class V>
  77 inline void unpack_tuple(const std::pair<T, U>& p, V& a, V& b) BOOST_MATH_NOEXCEPT(T)
  78 {
  79    a = p.first;
  80    b = p.second;
  81 }
  82 template <class T, class U, class V>
  83 inline void unpack_0(const std::pair<T, U>& p, V& a) BOOST_MATH_NOEXCEPT(T)
  84 {
  85    a = p.first;
  86 }
  87
  88 template <class F, class T>
  89 void handle_zero_derivative(F f,
  90    T& last_f0,
  91    const T& f0,
  92    T& delta,
  93    T& result,
  94    T& guess,
  95    const T& min,
  96    const T& max) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
  97 {
  98    if (last_f0 == 0)
  99    {
 100       // this must be the first iteration, pretend that we had a
 101       // previous one at either min or max:
 102       if (result == min)
 103       {
 104          guess = max;
 105       }
 106       else
 107       {
 108          guess = min;
 109       }
 110       unpack_0(f(guess), last_f0);
 111       delta = guess - result;
 112    }
 113    if (sign(last_f0) * sign(f0) < 0)
 114    {
 115       // we've crossed over so move in opposite direction to last step:
 116       if (delta < 0)
 117       {
 118          delta = (result - min) / 2;
 119       }
 120       else
 121       {
 122          delta = (result - max) / 2;
 123       }
 124    }
 125    else
 126    {
 127       // move in same direction as last step:
 128       if (delta < 0)
 129       {
 130          delta = (result - max) / 2;
 131       }
 132       else
 133       {
 134          delta = (result - min) / 2;
 135       }
 136    }
 137 }
 138
 139 } // namespace
 140
 141 template <class F, class T, class Tol, class Policy>
 142 std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter, const Policy& pol) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<Policy>::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 143 {
 144    T fmin = f(min);
 145    T fmax = f(max);
 146    if (fmin == 0)
 147    {
 148       max_iter = 2;
 149       return std::make_pair(min, min);
 150    }
 151    if (fmax == 0)
 152    {
 153       max_iter = 2;
 154       return std::make_pair(max, max);
 155    }
 156
 157    //
 158    // Error checking:
 159    //
 160    static const char* function = "boost::math::tools::bisect<%1%>";
 161    if (min >= max)
 162    {
 163       return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function,
 164          "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol));
 165    }
 166    if (fmin * fmax >= 0)
 167    {
 168       return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function,
 169          "No change of sign in boost::math::tools::bisect, either there is no root to find, or there are multiple roots in the interval (f(min) = %1%).", fmin, pol));
 170    }
 171
 172    //
 173    // Three function invocations so far:
 174    //
 175    boost::uintmax_t count = max_iter;
 176    if (count < 3)
 177       count = 0;
 178    else
 179       count -= 3;
 180
 181    while (count && (0 == tol(min, max)))
 182    {
 183       T mid = (min + max) / 2;
 184       T fmid = f(mid);
 185       if ((mid == max) || (mid == min))
 186          break;
 187       if (fmid == 0)
 188       {
 189          min = max = mid;
 190          break;
 191       }
 192       else if (sign(fmid) * sign(fmin) < 0)
 193       {
 194          max = mid;
 195       }
 196       else
 197       {
 198          min = mid;
 199          fmin = fmid;
 200       }
 201       --count;
 202    }
 203
 204    max_iter -= count;
 205
 206 #ifdef BOOST_MATH_INSTRUMENT
 207    std::cout << "Bisection iteration, final count = " << max_iter << std::endl;
 208
 209    static boost::uintmax_t max_count = 0;
 210    if (max_iter > max_count)
 211    {
 212       max_count = max_iter;
 213       std::cout << "Maximum iterations: " << max_iter << std::endl;
 214    }
 215 #endif
 216
 217    return std::make_pair(min, max);
 218 }
 219
 220 template <class F, class T, class Tol>
 221 inline std::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::uintmax_t& max_iter)  BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 222 {
 223    return bisect(f, min, max, tol, max_iter, policies::policy<>());
 224 }
 225
 226 template <class F, class T, class Tol>
 227 inline std::pair<T, T> bisect(F f, T min, T max, Tol tol) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 228 {
 229    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 230    return bisect(f, min, max, tol, m, policies::policy<>());
 231 }
 232
 233
 234 template <class F, class T>
 235 T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 236 {
 237    BOOST_MATH_STD_USING
 238
 239    static const char* function = "boost::math::tools::newton_raphson_iterate<%1%>";
 240    if (min >= max)
 241    {
 242       return policies::raise_evaluation_error(function, "Range arguments in wrong order in boost::math::tools::newton_raphson_iterate(first arg=%1%)", min, boost::math::policies::policy<>());
 243    }
 244
 245    T f0(0), f1, last_f0(0);
 246    T result = guess;
 247
 248    T factor = static_cast<T>(ldexp(1.0, 1 - digits));
 249    T delta = tools::max_value<T>();
 250    T delta1 = tools::max_value<T>();
 251    T delta2 = tools::max_value<T>();
 252
 253    //
 254    // We use these to sanity check that we do actually bracket a root,
 255    // we update these to the function value when we update the endpoints
 256    // of the range.  Then, provided at some point we update both endpoints
 257    // checking that max_range_f * min_range_f <= 0 verifies there is a root
 258    // to be found somewhere.  Note that if there is no root, and we approach
 259    // a local minima, then the derivative will go to zero, and hence the next
 260    // step will jump out of bounds (or at least past the minima), so this
 261    // check *should* happen in pathological cases.
 262    //
 263    T max_range_f = 0;
 264    T min_range_f = 0;
 265
 266    boost::uintmax_t count(max_iter);
 267
 268 #ifdef BOOST_MATH_INSTRUMENT
 269    std::cout << "Newton_raphson_iterate, guess = " << guess << ", min = " << min << ", max = " << max
 270       << ", digits = " << digits << ", max_iter = " << max_iter << std::endl;
 271 #endif
 272
 273    do {
 274       last_f0 = f0;
 275       delta2 = delta1;
 276       delta1 = delta;
 277       detail::unpack_tuple(f(result), f0, f1);
 278       --count;
 279       if (0 == f0)
 280          break;
 281       if (f1 == 0)
 282       {
 283          // Oops zero derivative!!!
 284 #ifdef BOOST_MATH_INSTRUMENT
 285          std::cout << "Newton iteration, zero derivative found!" << std::endl;
 286 #endif
 287          detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
 288       }
 289       else
 290       {
 291          delta = f0 / f1;
 292       }
 293 #ifdef BOOST_MATH_INSTRUMENT
 294       std::cout << "Newton iteration " << max_iter - count << ", delta = " << delta << std::endl;
 295 #endif
 296       if (fabs(delta * 2) > fabs(delta2))
 297       {
 298          // Last two steps haven't converged.
 299          T shift = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
 300          if ((result != 0) && (fabs(shift) > fabs(result)))
 301          {
 302             delta = sign(delta) * fabs(result) * 1.1f; // Protect against huge jumps!
 303             //delta = sign(delta) * result; // Protect against huge jumps! Failed for negative result. https://github.com/boostorg/math/issues/216
 304          }
 305          else
 306             delta = shift;
 307          // reset delta1/2 so we don't take this branch next time round:
 308          delta1 = 3 * delta;
 309          delta2 = 3 * delta;
 310       }
 311       guess = result;
 312       result -= delta;
 313       if (result <= min)
 314       {
 315          delta = 0.5F * (guess - min);
 316          result = guess - delta;
 317          if ((result == min) || (result == max))
 318             break;
 319       }
 320       else if (result >= max)
 321       {
 322          delta = 0.5F * (guess - max);
 323          result = guess - delta;
 324          if ((result == min) || (result == max))
 325             break;
 326       }
 327       // Update brackets:
 328       if (delta > 0)
 329       {
 330          max = guess;
 331          max_range_f = f0;
 332       }
 333       else
 334       {
 335          min = guess;
 336          min_range_f = f0;
 337       }
 338       //
 339       // Sanity check that we bracket the root:
 340       //
 341       if (max_range_f * min_range_f > 0)
 342       {
 343          return policies::raise_evaluation_error(function, "There appears to be no root to be found in boost::math::tools::newton_raphson_iterate, perhaps we have a local minima near current best guess of %1%", guess, boost::math::policies::policy<>());
 344       }
 345    }while(count && (fabs(result * factor) < fabs(delta)));
 346
 347    max_iter -= count;
 348
 349 #ifdef BOOST_MATH_INSTRUMENT
 350    std::cout << "Newton Raphson final iteration count = " << max_iter << std::endl;
 351
 352    static boost::uintmax_t max_count = 0;
 353    if (max_iter > max_count)
 354    {
 355       max_count = max_iter;
 356       // std::cout << "Maximum iterations: " << max_iter << std::endl;
 357       // Puzzled what this tells us, so commented out for now?
 358    }
 359 #endif
 360
 361    return result;
 362 }
 363
 364 template <class F, class T>
 365 inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 366 {
 367    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 368    return newton_raphson_iterate(f, guess, min, max, digits, m);
 369 }
 370
 371 namespace detail {
 372
 373    struct halley_step
 374    {
 375       template <class T>
 376       static T step(const T& /*x*/, const T& f0, const T& f1, const T& f2) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T))
 377       {
 378          using std::fabs;
 379          T denom = 2 * f0;
 380          T num = 2 * f1 - f0 * (f2 / f1);
 381          T delta;
 382
 383          BOOST_MATH_INSTRUMENT_VARIABLE(denom);
 384          BOOST_MATH_INSTRUMENT_VARIABLE(num);
 385
 386          if ((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>()))
 387          {
 388             // possible overflow, use Newton step:
 389             delta = f0 / f1;
 390          }
 391          else
 392             delta = denom / num;
 393          return delta;
 394       }
 395    };
 396
 397    template <class F, class T>
 398    T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, boost::uintmax_t& count) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())));
 399
 400    template <class F, class T>
 401    T bracket_root_towards_max(F f, T guess, const T& f0, T& min, T& max, boost::uintmax_t& count) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 402    {
 403       using std::fabs;
 404       //
 405       // Move guess towards max until we bracket the root, updating min and max as we go:
 406       //
 407       T guess0 = guess;
 408       T multiplier = 2;
 409       T f_current = f0;
 410       if (fabs(min) < fabs(max))
 411       {
 412          while (--count && ((f_current < 0) == (f0 < 0)))
 413          {
 414             min = guess;
 415             guess *= multiplier;
 416             if (guess > max)
 417             {
 418                guess = max;
 419                f_current = -f_current;  // There must be a change of sign!
 420                break;
 421             }
 422             multiplier *= 2;
 423             unpack_0(f(guess), f_current);
 424          }
 425       }
 426       else
 427       {
 428          //
 429          // If min and max are negative we have to divide to head towards max:
 430          //
 431          while (--count && ((f_current < 0) == (f0 < 0)))
 432          {
 433             min = guess;
 434             guess /= multiplier;
 435             if (guess > max)
 436             {
 437                guess = max;
 438                f_current = -f_current;  // There must be a change of sign!
 439                break;
 440             }
 441             multiplier *= 2;
 442             unpack_0(f(guess), f_current);
 443          }
 444       }
 445
 446       if (count)
 447       {
 448          max = guess;
 449          if (multiplier > 16)
 450             return (guess0 - guess) + bracket_root_towards_min(f, guess, f_current, min, max, count);
 451       }
 452       return guess0 - (max + min) / 2;
 453    }
 454
 455    template <class F, class T>
 456    T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, boost::uintmax_t& count) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 457    {
 458       using std::fabs;
 459       //
 460       // Move guess towards min until we bracket the root, updating min and max as we go:
 461       //
 462       T guess0 = guess;
 463       T multiplier = 2;
 464       T f_current = f0;
 465
 466       if (fabs(min) < fabs(max))
 467       {
 468          while (--count && ((f_current < 0) == (f0 < 0)))
 469          {
 470             max = guess;
 471             guess /= multiplier;
 472             if (guess < min)
 473             {
 474                guess = min;
 475                f_current = -f_current;  // There must be a change of sign!
 476                break;
 477             }
 478             multiplier *= 2;
 479             unpack_0(f(guess), f_current);
 480          }
 481       }
 482       else
 483       {
 484          //
 485          // If min and max are negative we have to multiply to head towards min:
 486          //
 487          while (--count && ((f_current < 0) == (f0 < 0)))
 488          {
 489             max = guess;
 490             guess *= multiplier;
 491             if (guess < min)
 492             {
 493                guess = min;
 494                f_current = -f_current;  // There must be a change of sign!
 495                break;
 496             }
 497             multiplier *= 2;
 498             unpack_0(f(guess), f_current);
 499          }
 500       }
 501
 502       if (count)
 503       {
 504          min = guess;
 505          if (multiplier > 16)
 506             return (guess0 - guess) + bracket_root_towards_max(f, guess, f_current, min, max, count);
 507       }
 508       return guess0 - (max + min) / 2;
 509    }
 510
 511    template <class Stepper, class F, class T>
 512    T second_order_root_finder(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 513    {
 514       BOOST_MATH_STD_USING
 515
 516 #ifdef BOOST_MATH_INSTRUMENT
 517         std::cout << "Second order root iteration, guess = " << guess << ", min = " << min << ", max = " << max
 518         << ", digits = " << digits << ", max_iter = " << max_iter << std::endl;
 519 #endif
 520       static const char* function = "boost::math::tools::halley_iterate<%1%>";
 521       if (min >= max)
 522       {
 523          return policies::raise_evaluation_error(function, "Range arguments in wrong order in boost::math::tools::halley_iterate(first arg=%1%)", min, boost::math::policies::policy<>());
 524       }
 525
 526       T f0(0), f1, f2;
 527       T result = guess;
 528
 529       T factor = ldexp(static_cast<T>(1.0), 1 - digits);
 530       T delta = (std::max)(T(10000000 * guess), T(10000000));  // arbitarily large delta
 531       T last_f0 = 0;
 532       T delta1 = delta;
 533       T delta2 = delta;
 534       bool out_of_bounds_sentry = false;
 535
 536    #ifdef BOOST_MATH_INSTRUMENT
 537       std::cout << "Second order root iteration, limit = " << factor << std::endl;
 538    #endif
 539
 540       //
 541       // We use these to sanity check that we do actually bracket a root,
 542       // we update these to the function value when we update the endpoints
 543       // of the range.  Then, provided at some point we update both endpoints
 544       // checking that max_range_f * min_range_f <= 0 verifies there is a root
 545       // to be found somewhere.  Note that if there is no root, and we approach
 546       // a local minima, then the derivative will go to zero, and hence the next
 547       // step will jump out of bounds (or at least past the minima), so this
 548       // check *should* happen in pathological cases.
 549       //
 550       T max_range_f = 0;
 551       T min_range_f = 0;
 552
 553       boost::uintmax_t count(max_iter);
 554
 555       do {
 556          last_f0 = f0;
 557          delta2 = delta1;
 558          delta1 = delta;
 559          detail::unpack_tuple(f(result), f0, f1, f2);
 560          --count;
 561
 562          BOOST_MATH_INSTRUMENT_VARIABLE(f0);
 563          BOOST_MATH_INSTRUMENT_VARIABLE(f1);
 564          BOOST_MATH_INSTRUMENT_VARIABLE(f2);
 565
 566          if (0 == f0)
 567             break;
 568          if (f1 == 0)
 569          {
 570             // Oops zero derivative!!!
 571    #ifdef BOOST_MATH_INSTRUMENT
 572             std::cout << "Second order root iteration, zero derivative found!" << std::endl;
 573    #endif
 574             detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);
 575          }
 576          else
 577          {
 578             if (f2 != 0)
 579             {
 580                delta = Stepper::step(result, f0, f1, f2);
 581                if (delta * f1 / f0 < 0)
 582                {
 583                   // Oh dear, we have a problem as Newton and Halley steps
 584                   // disagree about which way we should move.  Probably
 585                   // there is cancelation error in the calculation of the
 586                   // Halley step, or else the derivatives are so small
 587                   // that their values are basically trash.  We will move
 588                   // in the direction indicated by a Newton step, but
 589                   // by no more than twice the current guess value, otherwise
 590                   // we can jump way out of bounds if we're not careful.
 591                   // See https://svn.boost.org/trac/boost/ticket/8314.
 592                   delta = f0 / f1;
 593                   if (fabs(delta) > 2 * fabs(guess))
 594                      delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess);
 595                }
 596             }
 597             else
 598                delta = f0 / f1;
 599          }
 600    #ifdef BOOST_MATH_INSTRUMENT
 601          std::cout << "Second order root iteration, delta = " << delta << std::endl;
 602    #endif
 603          T convergence = fabs(delta / delta2);
 604          if ((convergence > 0.8) && (convergence < 2))
 605          {
 606             // last two steps haven't converged.
 607             delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
 608             if ((result != 0) && (fabs(delta) > result))
 609                delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps!
 610             // reset delta2 so that this branch will *not* be taken on the
 611             // next iteration:
 612             delta2 = delta * 3;
 613             delta1 = delta * 3;
 614             BOOST_MATH_INSTRUMENT_VARIABLE(delta);
 615          }
 616          guess = result;
 617          result -= delta;
 618          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 619
 620          // check for out of bounds step:
 621          if (result < min)
 622          {
 623             T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min)))
 624                ? T(1000)
 625                : (fabs(min) < 1) && (fabs(tools::max_value<T>() * min) < fabs(result))
 626                ? ((min < 0) != (result < 0)) ? -tools::max_value<T>() : tools::max_value<T>() : T(result / min);
 627             if (fabs(diff) < 1)
 628                diff = 1 / diff;
 629             if (!out_of_bounds_sentry && (diff > 0) && (diff < 3))
 630             {
 631                // Only a small out of bounds step, lets assume that the result
 632                // is probably approximately at min:
 633                delta = 0.99f * (guess - min);
 634                result = guess - delta;
 635                out_of_bounds_sentry = true; // only take this branch once!
 636             }
 637             else
 638             {
 639                if (fabs(float_distance(min, max)) < 2)
 640                {
 641                   result = guess = (min + max) / 2;
 642                   break;
 643                }
 644                delta = bracket_root_towards_min(f, guess, f0, min, max, count);
 645                result = guess - delta;
 646                guess = min;
 647                continue;
 648             }
 649          }
 650          else if (result > max)
 651          {
 652             T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max);
 653             if (fabs(diff) < 1)
 654                diff = 1 / diff;
 655             if (!out_of_bounds_sentry && (diff > 0) && (diff < 3))
 656             {
 657                // Only a small out of bounds step, lets assume that the result
 658                // is probably approximately at min:
 659                delta = 0.99f * (guess - max);
 660                result = guess - delta;
 661                out_of_bounds_sentry = true; // only take this branch once!
 662             }
 663             else
 664             {
 665                if (fabs(float_distance(min, max)) < 2)
 666                {
 667                   result = guess = (min + max) / 2;
 668                   break;
 669                }
 670                delta = bracket_root_towards_max(f, guess, f0, min, max, count);
 671                result = guess - delta;
 672                guess = min;
 673                continue;
 674             }
 675          }
 676          // update brackets:
 677          if (delta > 0)
 678          {
 679             max = guess;
 680             max_range_f = f0;
 681          }
 682          else
 683          {
 684             min = guess;
 685             min_range_f = f0;
 686          }
 687          //
 688          // Sanity check that we bracket the root:
 689          //
 690          if (max_range_f * min_range_f > 0)
 691          {
 692             return policies::raise_evaluation_error(function, "There appears to be no root to be found in boost::math::tools::newton_raphson_iterate, perhaps we have a local minima near current best guess of %1%", guess, boost::math::policies::policy<>());
 693          }
 694       } while(count && (fabs(result * factor) < fabs(delta)));
 695
 696       max_iter -= count;
 697
 698    #ifdef BOOST_MATH_INSTRUMENT
 699       std::cout << "Second order root finder, final iteration count = " << max_iter << std::endl;
 700    #endif
 701
 702       return result;
 703    }
 704 } // T second_order_root_finder
 705
 706 template <class F, class T>
 707 T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 708 {
 709    return detail::second_order_root_finder<detail::halley_step>(f, guess, min, max, digits, max_iter);
 710 }
 711
 712 template <class F, class T>
 713 inline T halley_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 714 {
 715    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 716    return halley_iterate(f, guess, min, max, digits, m);
 717 }
 718
 719 namespace detail {
 720
 721    struct schroder_stepper
 722    {
 723       template <class T>
 724       static T step(const T& x, const T& f0, const T& f1, const T& f2) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T))
 725       {
 726          using std::fabs;
 727          T ratio = f0 / f1;
 728          T delta;
 729          if ((x != 0) && (fabs(ratio / x) < 0.1))
 730          {
 731             delta = ratio + (f2 / (2 * f1)) * ratio * ratio;
 732             // check second derivative doesn't over compensate:
 733             if (delta * ratio < 0)
 734                delta = ratio;
 735          }
 736          else
 737             delta = ratio;  // fall back to Newton iteration.
 738          return delta;
 739       }
 740    };
 741
 742 }
 743
 744 template <class F, class T>
 745 T schroder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 746 {
 747    return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);
 748 }
 749
 750 template <class F, class T>
 751 inline T schroder_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 752 {
 753    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 754    return schroder_iterate(f, guess, min, max, digits, m);
 755 }
 756 //
 757 // These two are the old spelling of this function, retained for backwards compatibity just in case:
 758 //
 759 template <class F, class T>
 760 T schroeder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 761 {
 762    return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);
 763 }
 764
 765 template <class F, class T>
 766 inline T schroeder_iterate(F f, T guess, T min, T max, int digits) BOOST_NOEXCEPT_IF(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
 767 {
 768    boost::uintmax_t m = (std::numeric_limits<boost::uintmax_t>::max)();
 769    return schroder_iterate(f, guess, min, max, digits, m);
 770 }
 771
 772 #ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS
 773 /*
 774    * Why do we set the default maximum number of iterations to the number of digits in the type?
 775    * Because for double roots, the number of digits increases linearly with the number of iterations,
 776    * so this default should recover full precision even in this somewhat pathological case.
 777    * For isolated roots, the problem is so rapidly convergent that this doesn't matter at all.
 778    */
 779 template<class Complex, class F>
 780 Complex complex_newton(F g, Complex guess, int max_iterations = std::numeric_limits<typename Complex::value_type>::digits)
 781 {
 782    typedef typename Complex::value_type Real;
 783    using std::norm;
 784    using std::abs;
 785    using std::max;
 786    // z0, z1, and z2 cannot be the same, in case we immediately need to resort to Muller's Method:
 787    Complex z0 = guess + Complex(1, 0);
 788    Complex z1 = guess + Complex(0, 1);
 789    Complex z2 = guess;
 790
 791    do {
 792       auto pair = g(z2);
 793       if (norm(pair.second) == 0)
 794       {
 795          // Muller's method. Notation follows Numerical Recipes, 9.5.2:
 796          Complex q = (z2 - z1) / (z1 - z0);
 797          auto P0 = g(z0);
 798          auto P1 = g(z1);
 799          Complex qp1 = static_cast<Complex>(1) + q;
 800          Complex A = q * (pair.first - qp1 * P1.first + q * P0.first);
 801
 802          Complex B = (static_cast<Complex>(2) * q + static_cast<Complex>(1)) * pair.first - qp1 * qp1 * P1.first + q * q * P0.first;
 803          Complex C = qp1 * pair.first;
 804          Complex rad = sqrt(B * B - static_cast<Complex>(4) * A * C);
 805          Complex denom1 = B + rad;
 806          Complex denom2 = B - rad;
 807          Complex correction = (z1 - z2) * static_cast<Complex>(2) * C;
 808          if (norm(denom1) > norm(denom2))
 809          {
 810             correction /= denom1;
 811          }
 812          else
 813          {
 814             correction /= denom2;
 815          }
 816
 817          z0 = z1;
 818          z1 = z2;
 819          z2 = z2 + correction;
 820       }
 821       else
 822       {
 823          z0 = z1;
 824          z1 = z2;
 825          z2 = z2 - (pair.first / pair.second);
 826       }
 827
 828       // See: https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root
 829       // If f' is continuous, then convergence of x_n -> x* implies f(x*) = 0.
 830       // This condition approximates this convergence condition by requiring three consecutive iterates to be clustered.
 831       Real tol = (max)(abs(z2) * std::numeric_limits<Real>::epsilon(), std::numeric_limits<Real>::epsilon());
 832       bool real_close = abs(z0.real() - z1.real()) < tol && abs(z0.real() - z2.real()) < tol && abs(z1.real() - z2.real()) < tol;
 833       bool imag_close = abs(z0.imag() - z1.imag()) < tol && abs(z0.imag() - z2.imag()) < tol && abs(z1.imag() - z2.imag()) < tol;
 834       if (real_close && imag_close)
 835       {
 836          return z2;
 837       }
 838
 839    } while (max_iterations--);
 840
 841    // The idea is that if we can get abs(f) < eps, we should, but if we go through all these iterations
 842    // and abs(f) < sqrt(eps), then roundoff error simply does not allow that we can evaluate f to < eps
 843    // This is somewhat awkward as it isn't scale invariant, but using the Daubechies coefficient example code,
 844    // I found this condition generates correct roots, whereas the scale invariant condition discussed here:
 845    // https://scicomp.stackexchange.com/questions/30597/defining-a-condition-number-and-termination-criteria-for-newtons-method
 846    // allows nonroots to be passed off as roots.
 847    auto pair = g(z2);
 848    if (abs(pair.first) < sqrt(std::numeric_limits<Real>::epsilon()))
 849    {
 850       return z2;
 851    }
 852
 853    return { std::numeric_limits<Real>::quiet_NaN(),
 854             std::numeric_limits<Real>::quiet_NaN() };
 855 }
 856 #endif
 857
 858
 859 #if !defined(BOOST_NO_CXX17_IF_CONSTEXPR)
 860 // https://stackoverflow.com/questions/48979861/numerically-stable-method-for-solving-quadratic-equations/50065711
 861 namespace detail
 862 {
 863 #if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1)
 864 float fma_workaround(float f) { return ::fmaf(f); }
 865 double fma_workaround(double f) { return ::fma(f); }
 866 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 867 long double fma_workaround(long double f) { return ::fmal(f); }
 868 #endif
 869 #endif
 870 template<class T>
 871 inline T discriminant(T const& a, T const& b, T const& c)
 872 {
 873    T w = 4 * a * c;
 874 #if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1)
 875    T e = fma_workaround(-c, 4 * a, w);
 876    T f = fma_workaround(b, b, -w);
 877 #else
 878    T e = std::fma(-c, 4 * a, w);
 879    T f = std::fma(b, b, -w);
 880 #endif
 881    return f + e;
 882 }
 883
 884 template<class T>
 885 std::pair<T, T> quadratic_roots_imp(T const& a, T const& b, T const& c)
 886 {
 887    using std::copysign;
 888    using std::sqrt;
 889    if constexpr (std::is_floating_point<T>::value)
 890    {
 891       T nan = std::numeric_limits<T>::quiet_NaN();
 892       if (a == 0)
 893       {
 894          if (b == 0 && c != 0)
 895          {
 896             return std::pair<T, T>(nan, nan);
 897          }
 898          else if (b == 0 && c == 0)
 899          {
 900             return std::pair<T, T>(0, 0);
 901          }
 902          return std::pair<T, T>(-c / b, -c / b);
 903       }
 904       if (b == 0)
 905       {
 906          T x0_sq = -c / a;
 907          if (x0_sq < 0) {
 908             return std::pair<T, T>(nan, nan);
 909          }
 910          T x0 = sqrt(x0_sq);
 911          return std::pair<T, T>(-x0, x0);
 912       }
 913       T discriminant = detail::discriminant(a, b, c);
 914       // Is there a sane way to flush very small negative values to zero?
 915       // If there is I don't know of it.
 916       if (discriminant < 0)
 917       {
 918          return std::pair<T, T>(nan, nan);
 919       }
 920       T q = -(b + copysign(sqrt(discriminant), b)) / T(2);
 921       T x0 = q / a;
 922       T x1 = c / q;
 923       if (x0 < x1)
 924       {
 925          return std::pair<T, T>(x0, x1);
 926       }
 927       return std::pair<T, T>(x1, x0);
 928    }
 929    else if constexpr (boost::math::tools::is_complex_type<T>::value)
 930    {
 931       typename T::value_type nan = std::numeric_limits<typename T::value_type>::quiet_NaN();
 932       if (a.real() == 0 && a.imag() == 0)
 933       {
 934          using std::norm;
 935          if (b.real() == 0 && b.imag() && norm(c) != 0)
 936          {
 937             return std::pair<T, T>({ nan, nan }, { nan, nan });
 938          }
 939          else if (b.real() == 0 && b.imag() && c.real() == 0 && c.imag() == 0)
 940          {
 941             return std::pair<T, T>({ 0,0 }, { 0,0 });
 942          }
 943          return std::pair<T, T>(-c / b, -c / b);
 944       }
 945       if (b.real() == 0 && b.imag() == 0)
 946       {
 947          T x0_sq = -c / a;
 948          T x0 = sqrt(x0_sq);
 949          return std::pair<T, T>(-x0, x0);
 950       }
 951       // There's no fma for complex types:
 952       T discriminant = b * b - T(4) * a * c;
 953       T q = -(b + sqrt(discriminant)) / T(2);
 954       return std::pair<T, T>(q / a, c / q);
 955    }
 956    else // Most likely the type is a boost.multiprecision.
 957    {    //There is no fma for multiprecision, and in addition it doesn't seem to be useful, so revert to the naive computation.
 958       T nan = std::numeric_limits<T>::quiet_NaN();
 959       if (a == 0)
 960       {
 961          if (b == 0 && c != 0)
 962          {
 963             return std::pair<T, T>(nan, nan);
 964          }
 965          else if (b == 0 && c == 0)
 966          {
 967             return std::pair<T, T>(0, 0);
 968          }
 969          return std::pair<T, T>(-c / b, -c / b);
 970       }
 971       if (b == 0)
 972       {
 973          T x0_sq = -c / a;
 974          if (x0_sq < 0) {
 975             return std::pair<T, T>(nan, nan);
 976          }
 977          T x0 = sqrt(x0_sq);
 978          return std::pair<T, T>(-x0, x0);
 979       }
 980       T discriminant = b * b - 4 * a * c;
 981       if (discriminant < 0)
 982       {
 983          return std::pair<T, T>(nan, nan);
 984       }
 985       T q = -(b + copysign(sqrt(discriminant), b)) / T(2);
 986       T x0 = q / a;
 987       T x1 = c / q;
 988       if (x0 < x1)
 989       {
 990          return std::pair<T, T>(x0, x1);
 991       }
 992       return std::pair<T, T>(x1, x0);
 993    }
 994 }
 995 }  // namespace detail
 996
 997 template<class T1, class T2 = T1, class T3 = T1>
 998 inline std::pair<typename tools::promote_args<T1, T2, T3>::type, typename tools::promote_args<T1, T2, T3>::type> quadratic_roots(T1 const& a, T2 const& b, T3 const& c)
 999 {
1000    typedef typename tools::promote_args<T1, T2, T3>::type value_type;
1001    return detail::quadratic_roots_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(c));
1002 }
1003
1004 #endif
1005
1006 } // namespace tools
1007 } // namespace math
1008 } // namespace boost
1009
1010 #endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP