637 lines · plain
1// (C) Copyright John Maddock 2006.2// (C) Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_TOOLS_SOLVE_ROOT_HPP8#define BOOST_MATH_TOOLS_SOLVE_ROOT_HPP9 10#ifdef _MSC_VER11#pragma once12#endif13 14#include <boost/math/tools/config.hpp>15#include <boost/math/tools/precision.hpp>16#include <boost/math/tools/numeric_limits.hpp>17#include <boost/math/tools/tuple.hpp>18#include <boost/math/tools/cstdint.hpp>19#include <boost/math/policies/error_handling.hpp>20#include <boost/math/special_functions/sign.hpp>21 22#ifdef BOOST_MATH_LOG_ROOT_ITERATIONS23# define BOOST_MATH_LOGGER_INCLUDE <boost/math/tools/iteration_logger.hpp>24# include BOOST_MATH_LOGGER_INCLUDE25# undef BOOST_MATH_LOGGER_INCLUDE26#else27# define BOOST_MATH_LOG_COUNT(count)28#endif29 30namespace boost{ namespace math{ namespace tools{31 32template <class T>33class eps_tolerance34{35public:36 BOOST_MATH_GPU_ENABLED eps_tolerance() : eps(4 * tools::epsilon<T>())37 {38 39 }40 BOOST_MATH_GPU_ENABLED eps_tolerance(unsigned bits)41 {42 BOOST_MATH_STD_USING43 eps = BOOST_MATH_GPU_SAFE_MAX(T(ldexp(1.0F, 1-bits)), T(4 * tools::epsilon<T>()));44 }45 BOOST_MATH_GPU_ENABLED bool operator()(const T& a, const T& b)46 {47 BOOST_MATH_STD_USING48 return fabs(a - b) <= (eps * BOOST_MATH_GPU_SAFE_MIN(fabs(a), fabs(b)));49 }50private:51 T eps;52};53 54// CUDA warns about __host__ __device__ marker on defaulted constructor55// but the warning is benign 56#ifdef BOOST_MATH_ENABLE_CUDA57# pragma nv_diag_suppress 2001258#endif59 60struct equal_floor61{62 BOOST_MATH_GPU_ENABLED equal_floor() = default;63 64 template <class T>65 BOOST_MATH_GPU_ENABLED bool operator()(const T& a, const T& b)66 {67 BOOST_MATH_STD_USING68 return (floor(a) == floor(b)) || (fabs((b-a)/b) < boost::math::tools::epsilon<T>() * 2);69 }70};71 72struct equal_ceil73{74 BOOST_MATH_GPU_ENABLED equal_ceil() = default;75 76 template <class T>77 BOOST_MATH_GPU_ENABLED bool operator()(const T& a, const T& b)78 {79 BOOST_MATH_STD_USING80 return (ceil(a) == ceil(b)) || (fabs((b - a) / b) < boost::math::tools::epsilon<T>() * 2);81 }82};83 84struct equal_nearest_integer85{86 BOOST_MATH_GPU_ENABLED equal_nearest_integer() = default;87 88 template <class T>89 BOOST_MATH_GPU_ENABLED bool operator()(const T& a, const T& b)90 {91 BOOST_MATH_STD_USING92 return (floor(a + 0.5f) == floor(b + 0.5f)) || (fabs((b - a) / b) < boost::math::tools::epsilon<T>() * 2);93 }94};95 96#ifdef BOOST_MATH_ENABLE_CUDA97# pragma nv_diag_default 2001298#endif99 100namespace detail{101 102template <class F, class T>103BOOST_MATH_GPU_ENABLED void bracket(F f, T& a, T& b, T c, T& fa, T& fb, T& d, T& fd)104{105 //106 // Given a point c inside the existing enclosing interval107 // [a, b] sets a = c if f(c) == 0, otherwise finds the new 108 // enclosing interval: either [a, c] or [c, b] and sets109 // d and fd to the point that has just been removed from110 // the interval. In other words d is the third best guess111 // to the root.112 //113 BOOST_MATH_STD_USING // For ADL of std math functions114 T tol = tools::epsilon<T>() * 2;115 //116 // If the interval [a,b] is very small, or if c is too close 117 // to one end of the interval then we need to adjust the118 // location of c accordingly:119 //120 if((b - a) < 2 * tol * a)121 {122 c = a + (b - a) / 2;123 }124 else if(c <= a + fabs(a) * tol)125 {126 c = a + fabs(a) * tol;127 }128 else if(c >= b - fabs(b) * tol)129 {130 c = b - fabs(b) * tol;131 }132 //133 // OK, lets invoke f(c):134 //135 T fc = f(c);136 //137 // if we have a zero then we have an exact solution to the root:138 //139 if(fc == 0)140 {141 a = c;142 fa = 0;143 d = 0;144 fd = 0;145 return;146 }147 //148 // Non-zero fc, update the interval:149 //150 if(boost::math::sign(fa) * boost::math::sign(fc) < 0)151 {152 d = b;153 fd = fb;154 b = c;155 fb = fc;156 }157 else158 {159 d = a;160 fd = fa;161 a = c;162 fa= fc;163 }164}165 166template <class T>167BOOST_MATH_GPU_ENABLED inline T safe_div(T num, T denom, T r)168{169 //170 // return num / denom without overflow,171 // return r if overflow would occur.172 //173 BOOST_MATH_STD_USING // For ADL of std math functions174 175 if(fabs(denom) < 1)176 {177 if(fabs(denom * tools::max_value<T>()) <= fabs(num))178 return r;179 }180 return num / denom;181}182 183template <class T>184BOOST_MATH_GPU_ENABLED inline T secant_interpolate(const T& a, const T& b, const T& fa, const T& fb)185{186 //187 // Performs standard secant interpolation of [a,b] given188 // function evaluations f(a) and f(b). Performs a bisection189 // if secant interpolation would leave us very close to either190 // a or b. Rationale: we only call this function when at least191 // one other form of interpolation has already failed, so we know192 // that the function is unlikely to be smooth with a root very193 // close to a or b.194 //195 BOOST_MATH_STD_USING // For ADL of std math functions196 197 T tol = tools::epsilon<T>() * 5;198 T c = a - (fa / (fb - fa)) * (b - a);199 if((c <= a + fabs(a) * tol) || (c >= b - fabs(b) * tol))200 return (a + b) / 2;201 return c;202}203 204template <class T>205BOOST_MATH_GPU_ENABLED T quadratic_interpolate(const T& a, const T& b, T const& d,206 const T& fa, const T& fb, T const& fd, 207 unsigned count)208{209 //210 // Performs quadratic interpolation to determine the next point,211 // takes count Newton steps to find the location of the212 // quadratic polynomial.213 //214 // Point d must lie outside of the interval [a,b], it is the third215 // best approximation to the root, after a and b.216 //217 // Note: this does not guarantee to find a root218 // inside [a, b], so we fall back to a secant step should219 // the result be out of range.220 //221 // Start by obtaining the coefficients of the quadratic polynomial:222 //223 T B = safe_div(T(fb - fa), T(b - a), tools::max_value<T>());224 T A = safe_div(T(fd - fb), T(d - b), tools::max_value<T>());225 A = safe_div(T(A - B), T(d - a), T(0));226 227 if(A == 0)228 {229 // failure to determine coefficients, try a secant step:230 return secant_interpolate(a, b, fa, fb);231 }232 //233 // Determine the starting point of the Newton steps:234 //235 T c;236 if(boost::math::sign(A) * boost::math::sign(fa) > 0)237 {238 c = a;239 }240 else241 {242 c = b;243 }244 //245 // Take the Newton steps:246 //247 for(unsigned i = 1; i <= count; ++i)248 {249 //c -= safe_div(B * c, (B + A * (2 * c - a - b)), 1 + c - a);250 c -= safe_div(T(fa+(B+A*(c-b))*(c-a)), T(B + A * (2 * c - a - b)), T(1 + c - a));251 }252 if((c <= a) || (c >= b))253 {254 // Oops, failure, try a secant step:255 c = secant_interpolate(a, b, fa, fb);256 }257 return c;258}259 260template <class T>261BOOST_MATH_GPU_ENABLED T cubic_interpolate(const T& a, const T& b, const T& d, 262 const T& e, const T& fa, const T& fb, 263 const T& fd, const T& fe)264{265 //266 // Uses inverse cubic interpolation of f(x) at points 267 // [a,b,d,e] to obtain an approximate root of f(x).268 // Points d and e lie outside the interval [a,b]269 // and are the third and forth best approximations270 // to the root that we have found so far.271 //272 // Note: this does not guarantee to find a root273 // inside [a, b], so we fall back to quadratic274 // interpolation in case of an erroneous result.275 //276 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b277 << " d = " << d << " e = " << e << " fa = " << fa << " fb = " << fb 278 << " fd = " << fd << " fe = " << fe);279 T q11 = (d - e) * fd / (fe - fd);280 T q21 = (b - d) * fb / (fd - fb);281 T q31 = (a - b) * fa / (fb - fa);282 T d21 = (b - d) * fd / (fd - fb);283 T d31 = (a - b) * fb / (fb - fa);284 BOOST_MATH_INSTRUMENT_CODE(285 "q11 = " << q11 << " q21 = " << q21 << " q31 = " << q31286 << " d21 = " << d21 << " d31 = " << d31);287 T q22 = (d21 - q11) * fb / (fe - fb);288 T q32 = (d31 - q21) * fa / (fd - fa);289 T d32 = (d31 - q21) * fd / (fd - fa);290 T q33 = (d32 - q22) * fa / (fe - fa);291 T c = q31 + q32 + q33 + a;292 BOOST_MATH_INSTRUMENT_CODE(293 "q22 = " << q22 << " q32 = " << q32 << " d32 = " << d32294 << " q33 = " << q33 << " c = " << c);295 296 if((c <= a) || (c >= b))297 {298 // Out of bounds step, fall back to quadratic interpolation:299 c = quadratic_interpolate(a, b, d, fa, fb, fd, 3);300 BOOST_MATH_INSTRUMENT_CODE(301 "Out of bounds interpolation, falling back to quadratic interpolation. c = " << c);302 }303 304 return c;305}306 307} // namespace detail308 309template <class F, class T, class Tol, class Policy>310BOOST_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)311{312 //313 // Main entry point and logic for Toms Algorithm 748314 // root finder.315 //316 BOOST_MATH_STD_USING // For ADL of std math functions317 318 constexpr auto function = "boost::math::tools::toms748_solve<%1%>";319 320 //321 // Sanity check - are we allowed to iterate at all?322 //323 if (max_iter == 0)324 return boost::math::make_pair(ax, bx);325 326 boost::math::uintmax_t count = max_iter;327 T a, b, fa, fb, c, u, fu, a0, b0, d, fd, e, fe;328 static const T mu = 0.5f;329 330 // initialise a, b and fa, fb:331 a = ax;332 b = bx;333 if(a >= b)334 return boost::math::detail::pair_from_single(policies::raise_domain_error(335 function, 336 "Parameters a and b out of order: a=%1%", a, pol));337 fa = fax;338 fb = fbx;339 340 if(tol(a, b) || (fa == 0) || (fb == 0))341 {342 max_iter = 0;343 if(fa == 0)344 b = a;345 else if(fb == 0)346 a = b;347 return boost::math::make_pair(a, b);348 }349 350 if(boost::math::sign(fa) * boost::math::sign(fb) > 0)351 return boost::math::detail::pair_from_single(policies::raise_domain_error(352 function, 353 "Parameters a and b do not bracket the root: a=%1%", a, pol));354 // dummy value for fd, e and fe:355 fe = e = fd = 1e5F;356 357 if(fa != 0)358 {359 //360 // On the first step we take a secant step:361 //362 c = detail::secant_interpolate(a, b, fa, fb);363 detail::bracket(f, a, b, c, fa, fb, d, fd);364 --count;365 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);366 367 if(count && (fa != 0) && !tol(a, b))368 {369 //370 // On the second step we take a quadratic interpolation:371 //372 c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2);373 e = d;374 fe = fd;375 detail::bracket(f, a, b, c, fa, fb, d, fd);376 --count;377 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);378 }379 }380 381 while(count && (fa != 0) && !tol(a, b))382 {383 // save our brackets:384 a0 = a;385 b0 = b;386 //387 // Starting with the third step taken388 // we can use either quadratic or cubic interpolation.389 // Cubic interpolation requires that all four function values390 // fa, fb, fd, and fe are distinct, should that not be the case391 // then variable prof will get set to true, and we'll end up392 // taking a quadratic step instead.393 //394 T min_diff = tools::min_value<T>() * 32;395 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);396 if(prof)397 {398 c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 2);399 BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!");400 }401 else402 {403 c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe);404 }405 //406 // re-bracket, and check for termination:407 //408 e = d;409 fe = fd;410 detail::bracket(f, a, b, c, fa, fb, d, fd);411 if((0 == --count) || (fa == 0) || tol(a, b))412 break;413 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);414 //415 // Now another interpolated step:416 //417 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);418 if(prof)419 {420 c = detail::quadratic_interpolate(a, b, d, fa, fb, fd, 3);421 BOOST_MATH_INSTRUMENT_CODE("Can't take cubic step!!!!");422 }423 else424 {425 c = detail::cubic_interpolate(a, b, d, e, fa, fb, fd, fe);426 }427 //428 // Bracket again, and check termination condition, update e:429 //430 detail::bracket(f, a, b, c, fa, fb, d, fd);431 if((0 == --count) || (fa == 0) || tol(a, b))432 break;433 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);434 //435 // Now we take a double-length secant step:436 //437 if(fabs(fa) < fabs(fb))438 {439 u = a;440 fu = fa;441 }442 else443 {444 u = b;445 fu = fb;446 }447 c = u - 2 * (fu / (fb - fa)) * (b - a);448 if(fabs(c - u) > (b - a) / 2)449 {450 c = a + (b - a) / 2;451 }452 //453 // Bracket again, and check termination condition:454 //455 e = d;456 fe = fd;457 detail::bracket(f, a, b, c, fa, fb, d, fd);458 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);459 BOOST_MATH_INSTRUMENT_CODE(" tol = " << T((fabs(a) - fabs(b)) / fabs(a)));460 if((0 == --count) || (fa == 0) || tol(a, b))461 break;462 //463 // And finally... check to see if an additional bisection step is 464 // to be taken, we do this if we're not converging fast enough:465 //466 if((b - a) < mu * (b0 - a0))467 continue;468 //469 // bracket again on a bisection:470 //471 e = d;472 fe = fd;473 detail::bracket(f, a, b, T(a + (b - a) / 2), fa, fb, d, fd);474 --count;475 BOOST_MATH_INSTRUMENT_CODE("Not converging: Taking a bisection!!!!");476 BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);477 } // while loop478 479 max_iter -= count;480 if(fa == 0)481 {482 b = a;483 }484 else if(fb == 0)485 {486 a = b;487 }488 BOOST_MATH_LOG_COUNT(max_iter)489 return boost::math::make_pair(a, b);490}491 492template <class F, class T, class Tol>493BOOST_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)494{495 return toms748_solve(f, ax, bx, fax, fbx, tol, max_iter, policies::policy<>());496}497 498template <class F, class T, class Tol, class Policy>499BOOST_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)500{501 if (max_iter <= 2)502 return boost::math::make_pair(ax, bx);503 max_iter -= 2;504 boost::math::pair<T, T> r = toms748_solve(f, ax, bx, f(ax), f(bx), tol, max_iter, pol);505 max_iter += 2;506 return r;507}508 509template <class F, class T, class Tol>510BOOST_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)511{512 return toms748_solve(f, ax, bx, tol, max_iter, policies::policy<>());513}514 515template <class F, class T, class Tol, class Policy>516BOOST_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)517{518 BOOST_MATH_STD_USING519 constexpr auto function = "boost::math::tools::bracket_and_solve_root<%1%>";520 //521 // Set up initial brackets:522 //523 T a = guess;524 T b = a;525 T fa = f(a);526 T fb = fa;527 //528 // Set up invocation count:529 //530 boost::math::uintmax_t count = max_iter - 1;531 532 int step = 32;533 534 if((fa < 0) == (guess < 0 ? !rising : rising))535 {536 //537 // Zero is to the right of b, so walk upwards538 // until we find it:539 //540 while((boost::math::sign)(fb) == (boost::math::sign)(fa))541 {542 if(count == 0)543 return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol));544 //545 // Heuristic: normally it's best not to increase the step sizes as we'll just end up546 // with a really wide range to search for the root. However, if the initial guess was *really*547 // bad then we need to speed up the search otherwise we'll take forever if we're orders of548 // magnitude out. This happens most often if the guess is a small value (say 1) and the result549 // we're looking for is close to std::numeric_limits<T>::min().550 //551 if((max_iter - count) % static_cast<unsigned>(step) == 0u)552 {553 factor *= 2;554 if(step > 1) step /= 2;555 }556 //557 // Now go ahead and move our guess by "factor":558 //559 a = b;560 fa = fb;561 b *= factor;562 fb = f(b);563 --count;564 BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);565 }566 }567 else568 {569 //570 // Zero is to the left of a, so walk downwards571 // until we find it:572 //573 while((boost::math::sign)(fb) == (boost::math::sign)(fa))574 {575 if(fabs(a) < tools::min_value<T>())576 {577 // Escape route just in case the answer is zero!578 max_iter -= count;579 max_iter += 1;580 return a > 0 ? boost::math::make_pair(T(0), T(a)) : boost::math::make_pair(T(a), T(0)); 581 }582 if(count == 0)583 return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol));584 //585 // Heuristic: normally it's best not to increase the step sizes as we'll just end up586 // with a really wide range to search for the root. However, if the initial guess was *really*587 // bad then we need to speed up the search otherwise we'll take forever if we're orders of588 // magnitude out. This happens most often if the guess is a small value (say 1) and the result589 // we're looking for is close to std::numeric_limits<T>::min().590 //591 if((max_iter - count) % static_cast<unsigned>(step) == 0u)592 {593 factor *= 2;594 if(step > 1) step /= 2;595 }596 //597 // Now go ahead and move are guess by "factor":598 //599 b = a;600 fb = fa;601 a /= factor;602 fa = f(a);603 --count;604 BOOST_MATH_INSTRUMENT_CODE("a = " << a << " b = " << b << " fa = " << fa << " fb = " << fb << " count = " << count);605 }606 }607 max_iter -= count;608 max_iter += 1;609 boost::math::pair<T, T> r = toms748_solve(610 f, 611 (a < 0 ? b : a), 612 (a < 0 ? a : b), 613 (a < 0 ? fb : fa), 614 (a < 0 ? fa : fb), 615 tol, 616 count, 617 pol);618 max_iter += count;619 BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count);620 BOOST_MATH_LOG_COUNT(max_iter)621 return r;622}623 624template <class F, class T, class Tol>625BOOST_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)626{627 return bracket_and_solve_root(f, guess, factor, rising, tol, max_iter, policies::policy<>());628}629 630} // namespace tools631} // namespace math632} // namespace boost633 634 635#endif // BOOST_MATH_TOOLS_SOLVE_ROOT_HPP636 637