math/tools/toms748_solve.hpp

0001 //  (C) Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_TOOLS_SOLVE_ROOT_HPP
0007 #define BOOST_MATH_TOOLS_SOLVE_ROOT_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/precision.hpp>
0014 #include <boost/math/policies/error_handling.hpp>
0015 #include <boost/math/tools/config.hpp>
0016 #include <boost/math/special_functions/sign.hpp>
0017 #include <limits>
0018 #include <utility>
0019 #include <cstdint>
0020
0021 #ifdef BOOST_MATH_LOG_ROOT_ITERATIONS
0022 #  define BOOST_MATH_LOGGER_INCLUDE <boost/math/tools/iteration_logger.hpp>
0023 #  include BOOST_MATH_LOGGER_INCLUDE
0024 #  undef BOOST_MATH_LOGGER_INCLUDE
0025 #else
0026 #  define BOOST_MATH_LOG_COUNT(count)
0027 #endif
0028
0029 namespace boost{ namespace math{ namespace tools{
0030
0031 template <class T>
0032 class eps_tolerance
0033 {
0034 public:
0035    eps_tolerance() : eps(4 * tools::epsilon<T>())
0036    {
0037
0038    }
0039    eps_tolerance(unsigned bits)
0040    {
0041       BOOST_MATH_STD_USING
0042       eps = (std::max)(T(ldexp(1.0F, 1-bits)), T(4 * tools::epsilon<T>()));
0043    }
0044    bool operator()(const T& a, const T& b)
0045    {
0046       BOOST_MATH_STD_USING
0047       return fabs(a - b) <= (eps * (std::min)(fabs(a), fabs(b)));
0048    }
0049 private:
0050    T eps;
0051 };
0052
0053 struct equal_floor
0054 {
0055    equal_floor()= default;
0056    template <class T>
0057    bool operator()(const T& a, const T& b)
0058    {
0059       BOOST_MATH_STD_USING
0060       return floor(a) == floor(b);
0061    }
0062 };
0063
0064 struct equal_ceil
0065 {
0066    equal_ceil()= default;
0067    template <class T>
0068    bool operator()(const T& a, const T& b)
0069    {
0070       BOOST_MATH_STD_USING
0071       return ceil(a) == ceil(b);
0072    }
0073 };
0074
0075 struct equal_nearest_integer
0076 {
0077    equal_nearest_integer()= default;
0078    template <class T>
0079    bool operator()(const T& a, const T& b)
0080    {
0081       BOOST_MATH_STD_USING
0082       return floor(a + 0.5f) == floor(b + 0.5f);
0083    }
0084 };
0085
0086 namespace detail{
0087
0088 template <class F, class T>
0089 void bracket(F f, T& a, T& b, T c, T& fa, T& fb, T& d, T& fd)
0090 {
0091    //
0092    // Given a point c inside the existing enclosing interval
0093    // [a, b] sets a = c if f(c) == 0, otherwise finds the new
0094    // enclosing interval: either [a, c] or [c, b] and sets
0095    // d and fd to the point that has just been removed from
0096    // the interval.  In other words d is the third best guess
0097    // to the root.
0098    //
0099    BOOST_MATH_STD_USING  // For ADL of std math functions
0100    T tol = tools::epsilon<T>() * 2;
0101    //
0102    // If the interval [a,b] is very small, or if c is too close
0103    // to one end of the interval then we need to adjust the
0104    // location of c accordingly:
0105    //
0106    if((b - a) < 2 * tol * a)
0107    {
0108       c = a + (b - a) / 2;
0109    }
0110    else if(c <= a + fabs(a) * tol)
0111    {
0112       c = a + fabs(a) * tol;
0113    }
0114    else if(c >= b - fabs(b) * tol)
0115    {
0116       c = b - fabs(b) * tol;
0117    }
0118    //
0119    // OK, lets invoke f(c):
0120    //
0121    T fc = f(c);
0122    //
0123    // if we have a zero then we have an exact solution to the root:
0124    //
0125    if(fc == 0)
0126    {
0127       a = c;
0128       fa = 0;
0129       d = 0;
0130       fd = 0;
0131       return;
0132    }
0133    //
0134    // Non-zero fc, update the interval:
0135    //
0136    if(boost::math::sign(fa) * boost::math::sign(fc) < 0)
0137    {
0138       d = b;
0139       fd = fb;
0140       b = c;
0141       fb = fc;
0142    }
0143    else
0144    {
0145       d = a;
0146       fd = fa;
0147       a = c;
0148       fa= fc;
0149    }
0150 }
0151
0152 template <class T>
0153 inline T safe_div(T num, T denom, T r)
0154 {
0155    //
0156    // return num / denom without overflow,
0157    // return r if overflow would occur.
0158    //
0159    BOOST_MATH_STD_USING  // For ADL of std math functions
0160
0161    if(fabs(denom) < 1)
0162    {
0163       if(fabs(denom * tools::max_value<T>()) <= fabs(num))
0164          return r;
0165    }
0166    return num / denom;
0167 }
0168
0169 template <class T>
0170 inline T secant_interpolate(const T& a, const T& b, const T& fa, const T& fb)
0171 {
0172    //
0173    // Performs standard secant interpolation of [a,b] given
0174    // function evaluations f(a) and f(b).  Performs a bisection
0175    // if secant interpolation would leave us very close to either
0176    // a or b.  Rationale: we only call this function when at least
0177    // one other form of interpolation has already failed, so we know
0178    // that the function is unlikely to be smooth with a root very
0179    // close to a or b.
0180    //
0181    BOOST_MATH_STD_USING  // For ADL of std math functions
0182
0183    T tol = tools::epsilon<T>() * 5;
0184    T c = a - (fa / (fb - fa)) * (b - a);
0185    if((c <= a + fabs(a) * tol) || (c >= b - fabs(b) * tol))
0186       return (a + b) / 2;
0187    return c;
0188 }
0189
0190 template <class T>
0191 T quadratic_interpolate(const T& a, const T& b, T const& d,
0192                            const T& fa, const T& fb, T const& fd,
0193                            unsigned count)
0194 {
0195    //
0196    // Performs quadratic interpolation to determine the next point,
0197    // takes count Newton steps to find the location of the
0198    // quadratic polynomial.
0199    //
0200    // Point d must lie outside of the interval [a,b], it is the third
0201    // best approximation to the root, after a and b.
0202    //
0203    // Note: this does not guarantee to find a root
0204    // inside [a, b], so we fall back to a secant step should
0205    // the result be out of range.
0206    //
0207    // Start by obtaining the coefficients of the quadratic polynomial:
0208    //
0209    T B = safe_div(T(fb - fa), T(b - a), tools::max_value<T>());
0210    T A = safe_div(T(fd - fb), T(d - b), tools::max_value<T>());
0211    A = safe_div(T(A - B), T(d - a), T(0));
0212
0213    if(A == 0)
0214    {
0215       // failure to determine coefficients, try a secant step:
0216       return secant_interpolate(a, b, fa, fb);
0217    }
0218    //
0219    // Determine the starting point of the Newton steps:
0220    //
0221    T c;
0222    if(boost::math::sign(A) * boost::math::sign(fa) > 0)
0223    {
0224       c = a;
0225    }
0226    else
0227    {
0228       c = b;
0229    }
0230    //
0231    // Take the Newton steps:
0232    //
0233    for(unsigned i = 1; i <= count; ++i)
0234    {
0235       //c -= safe_div(B * c, (B + A * (2 * c - a - b)), 1 + c - a);
0236       c -= safe_div(T(fa+(B+A*(c-b))*(c-a)), T(B + A * (2 * c - a - b)), T(1 + c - a));
0237    }
0238    if((c <= a) || (c >= b))
0239    {
0240       // Oops, failure, try a secant step:
0241       c = secant_interpolate(a, b, fa, fb);
0242    }
0243    return c;
0244 }
0245
0246 template <class T>
0247 T cubic_interpolate(const T& a, const T& b, const T& d,
0248                     const T& e, const T& fa, const T& fb,
0249                     const T& fd, const T& fe)
0250 {
0251    //
0252    // Uses inverse cubic interpolation of f(x) at points
0253    // [a,b,d,e] to obtain an approximate root of f(x).
0254    // Points d and e lie outside the interval [a,b]
0255    // and are the third and forth best approximations
0256    // to the root that we have found so far.
0257    //
0258    // Note: this does not guarantee to find a root
0259    // inside [a, b], so we fall back to quadratic
0260    // interpolation in case of an erroneous result.
0261    //
0262    BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b
0263       << " d = " << d << " e = " << e << " fa = " << fa << " fb = " << fb
0264       << " fd = " << fd << " fe = " << fe);
0265    T q11 = (d - e) * fd / (fe - fd);
0266    T q21 = (b - d) * fb / (fd - fb);
0267    T q31 = (a - b) * fa / (fb - fa);
0268    T d21 = (b - d) * fd / (fd - fb);
0269    T d31 = (a - b) * fb / (fb - fa);
0270    BOOST_MATH_INSTRUMENT_CODE(
0271       "q11 = " << q11 << " q21 = " << q21 << " q31 = " << q31
0272       << " d21 = " << d21 << " d31 = " << d31);
0273    T q22 = (d21 - q11) * fb / (fe - fb);
0274    T q32 = (d31 - q21) * fa / (fd - fa);
0275    T d32 = (d31 - q21) * fd / (fd - fa);
0276    T q33 = (d32 - q22) * fa / (fe - fa);
0277    T c = q31 + q32 + q33 + a;
0278    BOOST_MATH_INSTRUMENT_CODE(
0279       "q22 = " << q22 << " q32 = " << q32 << " d32 = " << d32
0280       << " q33 = " << q33 << " c = " << c);
0281
0282    if((c <= a) || (c >= b))
0283    {
0284       // Out of bounds step, fall back to quadratic interpolation:
0285       c = quadratic_interpolate(a, b, d, fa, fb, fd, 3);
0286    BOOST_MATH_INSTRUMENT_CODE(
0287       "Out of bounds interpolation, falling back to quadratic interpolation. c = " << c);
0288    }
0289
0290    return c;
0291 }
0292
0293 } // namespace detail
0294
0295 template <class F, class T, class Tol, class Policy>
0296 std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, std::uintmax_t& max_iter, const Policy& pol)
0297 {
0298    //
0299    // Main entry point and logic for Toms Algorithm 748
0300    // root finder.
0301    //
0302    BOOST_MATH_STD_USING  // For ADL of std math functions
0303
0304    static const char* function = "boost::math::tools::toms748_solve<%1%>";
0305
0306    //
0307    // Sanity check - are we allowed to iterate at all?
0308    //
0309    if (max_iter == 0)
0310       return std::make_pair(ax, bx);
0311
0312    std::uintmax_t count = max_iter;
0313    T a, b, fa, fb, c, u, fu, a0, b0, d, fd, e, fe;
0314    static const T mu = 0.5f;
0315
0316    // initialise a, b and fa, fb:
0317    a = ax;
0318    b = bx;
0319    if(a >= b)
0320       return boost::math::detail::pair_from_single(policies::raise_domain_error(
0321          function,
0322          "Parameters a and b out of order: a=%1%", a, pol));
0323    fa = fax;
0324    fb = fbx;
0325
0326    if(tol(a, b) || (fa == 0) || (fb == 0))
0327    {
0328       max_iter = 0;
0329       if(fa == 0)
0330          b = a;
0331       else if(fb == 0)
0332          a = b;
0333       return std::make_pair(a, b);
0334    }
0335
0336    if(boost::math::sign(fa) * boost::math::sign(fb) > 0)
0337       return boost::math::detail::pair_from_single(policies::raise_domain_error(
0338          function,
0339          "Parameters a and b do not bracket the root: a=%1%", a, pol));
0340    // dummy value for fd, e and fe:
0341    fe = e = fd = 1e5F;
0342
0343    if(fa != 0)
0344    {
0345       //
0346       // On the first step we take a secant step:
0347       //
0348       c = detail::secant_interpolate(a, b, fa, fb);
0349       detail::bracket(f, a, b, c, fa, fb, d, fd);
0350       --count;
0351       BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
0352
0353       if(count && (fa != 0) && !tol(a, b))
0354       {
0355          //
0356          // On the second step we take a quadratic interpolation:
0357          //
0358          c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2);
0359          e = d;
0360          fe = fd;
0361          detail::bracket(f, a, b, c, fa, fb, d, fd);
0362          --count;
0363          BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
0364       }
0365    }
0366
0367    while(count && (fa != 0) && !tol(a, b))
0368    {
0369       // save our brackets:
0370       a0 = a;
0371       b0 = b;
0372       //
0373       // Starting with the third step taken
0374       // we can use either quadratic or cubic interpolation.
0375       // Cubic interpolation requires that all four function values
0376       // fa, fb, fd, and fe are distinct, should that not be the case
0377       // then variable prof will get set to true, and we'll end up
0378       // taking a quadratic step instead.
0379       //
0380       T min_diff = tools::min_value<T>() * 32;
0381       bool prof = (fabs(fa - fb) < min_diff) || (fabs(fa - fd) < min_diff) || (fabs(fa - fe) < min_diff) || (fabs(fb - fd) < min_diff) || (fabs(fb - fe) < min_diff) || (fabs(fd - fe) < min_diff);
0382       if(prof)
0383       {
0384          c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2);
0385          BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!");
0386       }
0387       else
0388       {
0389          c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe);
0390       }
0391       //
0392       // re-bracket, and check for termination:
0393       //
0394       e = d;
0395       fe = fd;
0396       detail::bracket(f, a, b, c, fa, fb, d, fd);
0397       if((0 == --count) || (fa == 0) || tol(a, b))
0398          break;
0399       BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
0400       //
0401       // Now another interpolated step:
0402       //
0403       prof = (fabs(fa - fb) < min_diff) || (fabs(fa - fd) < min_diff) || (fabs(fa - fe) < min_diff) || (fabs(fb - fd) < min_diff) || (fabs(fb - fe) < min_diff) || (fabs(fd - fe) < min_diff);
0404       if(prof)
0405       {
0406          c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 3);
0407          BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!");
0408       }
0409       else
0410       {
0411          c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe);
0412       }
0413       //
0414       // Bracket again, and check termination condition, update e:
0415       //
0416       detail::bracket(f, a, b, c, fa, fb, d, fd);
0417       if((0 == --count) || (fa == 0) || tol(a, b))
0418          break;
0419       BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
0420       //
0421       // Now we take a double-length secant step:
0422       //
0423       if(fabs(fa) < fabs(fb))
0424       {
0425          u = a;
0426          fu = fa;
0427       }
0428       else
0429       {
0430          u = b;
0431          fu = fb;
0432       }
0433       c = u - 2 * (fu / (fb - fa)) * (b - a);
0434       if(fabs(c - u) > (b - a) / 2)
0435       {
0436          c = a + (b - a) / 2;
0437       }
0438       //
0439       // Bracket again, and check termination condition:
0440       //
0441       e = d;
0442       fe = fd;
0443       detail::bracket(f, a, b, c, fa, fb, d, fd);
0444       BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
0445       BOOST_MATH_INSTRUMENT_CODE(" tol = " << T((fabs(a) - fabs(b)) / fabs(a)));
0446       if((0 == --count) || (fa == 0) || tol(a, b))
0447          break;
0448       //
0449       // And finally... check to see if an additional bisection step is
0450       // to be taken, we do this if we're not converging fast enough:
0451       //
0452       if((b - a) < mu * (b0 - a0))
0453          continue;
0454       //
0455       // bracket again on a bisection:
0456       //
0457       e = d;
0458       fe = fd;
0459       detail::bracket(f, a, b, T(a + (b - a) / 2), fa, fb, d, fd);
0460       --count;
0461       BOOST_MATH_INSTRUMENT_CODE("Not converging: Taking a bisection!!!!");
0462       BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
0463    } // while loop
0464
0465    max_iter -= count;
0466    if(fa == 0)
0467    {
0468       b = a;
0469    }
0470    else if(fb == 0)
0471    {
0472       a = b;
0473    }
0474    BOOST_MATH_LOG_COUNT(max_iter)
0475    return std::make_pair(a, b);
0476 }
0477
0478 template <class F, class T, class Tol>
0479 inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, std::uintmax_t& max_iter)
0480 {
0481    return toms748_solve(f, ax, bx, fax, fbx, tol, max_iter, policies::policy<>());
0482 }
0483
0484 template <class F, class T, class Tol, class Policy>
0485 inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, std::uintmax_t& max_iter, const Policy& pol)
0486 {
0487    if (max_iter <= 2)
0488       return std::make_pair(ax, bx);
0489    max_iter -= 2;
0490    std::pair<T, T> r = toms748_solve(f, ax, bx, f(ax), f(bx), tol, max_iter, pol);
0491    max_iter += 2;
0492    return r;
0493 }
0494
0495 template <class F, class T, class Tol>
0496 inline std::pair<T, T> toms748_solve(F f, const T& ax, const T& bx, Tol tol, std::uintmax_t& max_iter)
0497 {
0498    return toms748_solve(f, ax, bx, tol, max_iter, policies::policy<>());
0499 }
0500
0501 template <class F, class T, class Tol, class Policy>
0502 std::pair<T, T> bracket_and_solve_root(F f, const T& guess, T factor, bool rising, Tol tol, std::uintmax_t& max_iter, const Policy& pol)
0503 {
0504    BOOST_MATH_STD_USING
0505    static const char* function = "boost::math::tools::bracket_and_solve_root<%1%>";
0506    //
0507    // Set up initial brackets:
0508    //
0509    T a = guess;
0510    T b = a;
0511    T fa = f(a);
0512    T fb = fa;
0513    //
0514    // Set up invocation count:
0515    //
0516    std::uintmax_t count = max_iter - 1;
0517
0518    int step = 32;
0519
0520    if((fa < 0) == (guess < 0 ? !rising : rising))
0521    {
0522       //
0523       // Zero is to the right of b, so walk upwards
0524       // until we find it:
0525       //
0526       while((boost::math::sign)(fb) == (boost::math::sign)(fa))
0527       {
0528          if(count == 0)
0529             return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol));
0530          //
0531          // Heuristic: normally it's best not to increase the step sizes as we'll just end up
0532          // with a really wide range to search for the root.  However, if the initial guess was *really*
0533          // bad then we need to speed up the search otherwise we'll take forever if we're orders of
0534          // magnitude out.  This happens most often if the guess is a small value (say 1) and the result
0535          // we're looking for is close to std::numeric_limits<T>::min().
0536          //
0537          if((max_iter - count) % step == 0)
0538          {
0539             factor *= 2;
0540             if(step > 1) step /= 2;
0541          }
0542          //
0543          // Now go ahead and move our guess by "factor":
0544          //
0545          a = b;
0546          fa = fb;
0547          b *= factor;
0548          fb = f(b);
0549          --count;
0550          BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
0551       }
0552    }
0553    else
0554    {
0555       //
0556       // Zero is to the left of a, so walk downwards
0557       // until we find it:
0558       //
0559       while((boost::math::sign)(fb) == (boost::math::sign)(fa))
0560       {
0561          if(fabs(a) < tools::min_value<T>())
0562          {
0563             // Escape route just in case the answer is zero!
0564             max_iter -= count;
0565             max_iter += 1;
0566             return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0));
0567          }
0568          if(count == 0)
0569             return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol));
0570          //
0571          // Heuristic: normally it's best not to increase the step sizes as we'll just end up
0572          // with a really wide range to search for the root.  However, if the initial guess was *really*
0573          // bad then we need to speed up the search otherwise we'll take forever if we're orders of
0574          // magnitude out.  This happens most often if the guess is a small value (say 1) and the result
0575          // we're looking for is close to std::numeric_limits<T>::min().
0576          //
0577          if((max_iter - count) % step == 0)
0578          {
0579             factor *= 2;
0580             if(step > 1) step /= 2;
0581          }
0582          //
0583          // Now go ahead and move are guess by "factor":
0584          //
0585          b = a;
0586          fb = fa;
0587          a /= factor;
0588          fa = f(a);
0589          --count;
0590          BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);
0591       }
0592    }
0593    max_iter -= count;
0594    max_iter += 1;
0595    std::pair<T, T> r = toms748_solve(
0596       f,
0597       (a < 0 ? b : a),
0598       (a < 0 ? a : b),
0599       (a < 0 ? fb : fa),
0600       (a < 0 ? fa : fb),
0601       tol,
0602       count,
0603       pol);
0604    max_iter += count;
0605    BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);
0606    BOOST_MATH_LOG_COUNT(max_iter)
0607    return r;
0608 }
0609
0610 template <class F, class T, class Tol>
0611 inline std::pair<T, T> bracket_and_solve_root(F f, const T& guess, const T& factor, bool rising, Tol tol, std::uintmax_t& max_iter)
0612 {
0613    return bracket_and_solve_root(f, guess, factor, rising, tol, max_iter, policies::policy<>());
0614 }
0615
0616 } // namespace tools
0617 } // namespace math
0618 } // namespace boost
0619
0620
0621 #endif // BOOST_MATH_TOOLS_SOLVE_ROOT_HPP
0622