File indexing completed on 2025-09-17 08:36:50
0001
0002
0003
0004
0005
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
0055
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
0107
0108
0109
0110
0111
0112
0113 BOOST_MATH_STD_USING
0114 T tol = tools::epsilon<T>() * 2;
0115
0116
0117
0118
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
0134
0135 T fc = f(c);
0136
0137
0138
0139 if(fc == 0)
0140 {
0141 a = c;
0142 fa = 0;
0143 d = 0;
0144 fd = 0;
0145 return;
0146 }
0147
0148
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
0171
0172
0173 BOOST_MATH_STD_USING
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
0188
0189
0190
0191
0192
0193
0194
0195 BOOST_MATH_STD_USING
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
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
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
0230 return secant_interpolate(a, b, fa, fb);
0231 }
0232
0233
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
0246
0247 for(unsigned i = 1; i <= count; ++i)
0248 {
0249
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
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
0267
0268
0269
0270
0271
0272
0273
0274
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
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 }
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
0314
0315
0316 BOOST_MATH_STD_USING
0317
0318 constexpr auto function = "boost::math::tools::toms748_solve<%1%>";
0319
0320
0321
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
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
0355 fe = e = fd = 1e5F;
0356
0357 if(fa != 0)
0358 {
0359
0360
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
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
0384 a0 = a;
0385 b0 = b;
0386
0387
0388
0389
0390
0391
0392
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
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
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
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
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
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
0464
0465
0466 if((b - a) < mu * (b0 - a0))
0467 continue;
0468
0469
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 }
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
0522
0523 T a = guess;
0524 T b = a;
0525 T fa = f(a);
0526 T fb = fa;
0527
0528
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
0538
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
0546
0547
0548
0549
0550
0551 if((max_iter - count) % step == 0)
0552 {
0553 factor *= 2;
0554 if(step > 1) step /= 2;
0555 }
0556
0557
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
0571
0572
0573 while((boost::math::sign)(fb) == (boost::math::sign)(fa))
0574 {
0575 if(fabs(a) < tools::min_value<T>())
0576 {
0577
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
0586
0587
0588
0589
0590
0591 if((max_iter - count) % step == 0)
0592 {
0593 factor *= 2;
0594 if(step > 1) step /= 2;
0595 }
0596
0597
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 }
0631 }
0632 }
0633
0634
0635 #endif
0636