Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-17 08:36:50

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