1063 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_NEWTON_SOLVER_HPP8#define BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP9 10#ifdef _MSC_VER11#pragma once12#endif13 14#include <boost/math/tools/config.hpp>15#include <boost/math/tools/complex.hpp> // test for multiprecision types in complex Newton16#include <boost/math/tools/type_traits.hpp>17#include <boost/math/tools/cstdint.hpp>18#include <boost/math/tools/numeric_limits.hpp>19#include <boost/math/tools/tuple.hpp>20#include <boost/math/special_functions/sign.hpp>21#include <boost/math/policies/policy.hpp>22#include <boost/math/policies/error_handling.hpp>23 24#ifndef BOOST_MATH_HAS_GPU_SUPPORT25#include <boost/math/special_functions/next.hpp>26#include <boost/math/tools/toms748_solve.hpp>27#endif28 29namespace boost {30namespace math {31namespace tools {32 33namespace detail {34 35namespace dummy {36 37 template<int n, class T>38 BOOST_MATH_GPU_ENABLED typename T::value_type get(const T&) BOOST_MATH_NOEXCEPT(T);39}40 41template <class Tuple, class T>42BOOST_MATH_GPU_ENABLED void unpack_tuple(const Tuple& t, T& a, T& b) BOOST_MATH_NOEXCEPT(T)43{44 using dummy::get;45 // Use ADL to find the right overload for get:46 a = get<0>(t);47 b = get<1>(t);48}49template <class Tuple, class T>50BOOST_MATH_GPU_ENABLED void unpack_tuple(const Tuple& t, T& a, T& b, T& c) BOOST_MATH_NOEXCEPT(T)51{52 using dummy::get;53 // Use ADL to find the right overload for get:54 a = get<0>(t);55 b = get<1>(t);56 c = get<2>(t);57}58 59template <class Tuple, class T>60BOOST_MATH_GPU_ENABLED inline void unpack_0(const Tuple& t, T& val) BOOST_MATH_NOEXCEPT(T)61{62 using dummy::get;63 // Rely on ADL to find the correct overload of get:64 val = get<0>(t);65}66 67template <class T, class U, class V>68BOOST_MATH_GPU_ENABLED inline void unpack_tuple(const boost::math::pair<T, U>& p, V& a, V& b) BOOST_MATH_NOEXCEPT(T)69{70 a = p.first;71 b = p.second;72}73template <class T, class U, class V>74BOOST_MATH_GPU_ENABLED inline void unpack_0(const boost::math::pair<T, U>& p, V& a) BOOST_MATH_NOEXCEPT(T)75{76 a = p.first;77}78 79template <class F, class T>80BOOST_MATH_GPU_ENABLED void handle_zero_derivative(F f,81 T& last_f0,82 const T& f0,83 T& delta,84 T& result,85 T& guess,86 const T& min,87 const T& max) noexcept(BOOST_MATH_IS_FLOAT(T) 88 #ifndef BOOST_MATH_HAS_GPU_SUPPORT89 && noexcept(std::declval<F>()(std::declval<T>()))90 #endif91 )92{93 if (last_f0 == 0)94 {95 // this must be the first iteration, pretend that we had a96 // previous one at either min or max:97 if (result == min)98 {99 guess = max;100 }101 else102 {103 guess = min;104 }105 unpack_0(f(guess), last_f0);106 delta = guess - result;107 }108 if (sign(last_f0) * sign(f0) < 0)109 {110 // we've crossed over so move in opposite direction to last step:111 if (delta < 0)112 {113 delta = (result - min) / 2;114 }115 else116 {117 delta = (result - max) / 2;118 }119 }120 else121 {122 // move in same direction as last step:123 if (delta < 0)124 {125 delta = (result - max) / 2;126 }127 else128 {129 delta = (result - min) / 2;130 }131 }132}133 134} // namespace135 136template <class F, class T, class Tol, class Policy>137BOOST_MATH_GPU_ENABLED boost::math::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::math::uintmax_t& max_iter, const Policy& pol) noexcept(policies::is_noexcept_error_policy<Policy>::value && BOOST_MATH_IS_FLOAT(T) 138#ifndef BOOST_MATH_HAS_GPU_SUPPORT139&& noexcept(std::declval<F>()(std::declval<T>()))140#endif141)142{143 T fmin = f(min);144 T fmax = f(max);145 if (fmin == 0)146 {147 max_iter = 2;148 return boost::math::make_pair(min, min);149 }150 if (fmax == 0)151 {152 max_iter = 2;153 return boost::math::make_pair(max, max);154 }155 156 //157 // Error checking:158 //159 constexpr auto function = "boost::math::tools::bisect<%1%>";160 if (min >= max)161 {162 return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function,163 "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol));164 }165 if (fmin * fmax >= 0)166 {167 return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function,168 "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));169 }170 171 //172 // Three function invocations so far:173 //174 std::uintmax_t count = max_iter;175 if (count < 3)176 count = 0;177 else178 count -= 3;179 180 while (count && (0 == tol(min, max)))181 {182 T mid = (min + max) / 2;183 T fmid = f(mid);184 if ((mid == max) || (mid == min))185 break;186 if (fmid == 0)187 {188 min = max = mid;189 break;190 }191 else if (sign(fmid) * sign(fmin) < 0)192 {193 max = mid;194 }195 else196 {197 min = mid;198 fmin = fmid;199 }200 --count;201 }202 203 max_iter -= count;204 205#ifdef BOOST_MATH_INSTRUMENT206 std::cout << "Bisection required " << max_iter << " iterations.\n";207#endif208 209 return boost::math::make_pair(min, max);210}211 212template <class F, class T, class Tol>213BOOST_MATH_GPU_ENABLED inline boost::math::pair<T, T> bisect(F f, T min, T max, Tol tol, boost::math::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value && BOOST_MATH_IS_FLOAT(T)214#ifndef BOOST_MATH_HAS_GPU_SUPPORT215&& noexcept(std::declval<F>()(std::declval<T>()))216#endif217)218{219 return bisect(f, min, max, tol, max_iter, policies::policy<>());220}221 222template <class F, class T, class Tol>223BOOST_MATH_GPU_ENABLED inline boost::math::pair<T, T> bisect(F f, T min, T max, Tol tol) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value && BOOST_MATH_IS_FLOAT(T) 224#ifndef BOOST_MATH_HAS_GPU_SUPPORT225&& noexcept(std::declval<F>()(std::declval<T>()))226#endif227)228{229 boost::math::uintmax_t m = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();230 return bisect(f, min, max, tol, m, policies::policy<>());231}232 233 234template <class F, class T>235BOOST_MATH_GPU_ENABLED T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::math::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value && BOOST_MATH_IS_FLOAT(T)236#ifndef BOOST_MATH_HAS_GPU_SUPPORT237&& noexcept(std::declval<F>()(std::declval<T>()))238#endif239)240{241 BOOST_MATH_STD_USING242 243 constexpr auto function = "boost::math::tools::newton_raphson_iterate<%1%>";244 if (min > max)245 {246 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<>());247 }248 249 T f0(0), f1, last_f0(0);250 T result = guess;251 252 T factor = static_cast<T>(ldexp(1.0, 1 - digits));253 T delta = tools::max_value<T>();254 T delta1 = tools::max_value<T>();255 T delta2 = tools::max_value<T>();256 257 //258 // We use these to sanity check that we do actually bracket a root,259 // we update these to the function value when we update the endpoints260 // of the range. Then, provided at some point we update both endpoints261 // checking that max_range_f * min_range_f <= 0 verifies there is a root262 // to be found somewhere. Note that if there is no root, and we approach 263 // a local minima, then the derivative will go to zero, and hence the next264 // step will jump out of bounds (or at least past the minima), so this265 // check *should* happen in pathological cases.266 //267 T max_range_f = 0;268 T min_range_f = 0;269 270 boost::math::uintmax_t count(max_iter);271 272#ifdef BOOST_MATH_INSTRUMENT273 std::cout << "Newton_raphson_iterate, guess = " << guess << ", min = " << min << ", max = " << max274 << ", digits = " << digits << ", max_iter = " << max_iter << "\n";275#endif276 277 do {278 last_f0 = f0;279 delta2 = delta1;280 delta1 = delta;281 detail::unpack_tuple(f(result), f0, f1);282 --count;283 if (0 == f0)284 break;285 if (f1 == 0)286 {287 // Oops zero derivative!!!288 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);289 }290 else291 {292 delta = f0 / f1;293 }294#ifdef BOOST_MATH_INSTRUMENT295 std::cout << "Newton iteration " << max_iter - count << ", delta = " << delta << ", residual = " << f0 << "\n";296#endif297 if (fabs(delta * 2) > fabs(delta2))298 {299 // Last two steps haven't converged.300 T shift = (delta > 0) ? (result - min) / 2 : (result - max) / 2;301 if ((result != 0) && (fabs(shift) > fabs(result)))302 {303 delta = sign(delta) * fabs(result); // protect against huge jumps!304 }305 else306 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 else334 {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_INSTRUMENT350 std::cout << "Newton Raphson required " << max_iter << " iterations\n";351#endif352 353 return result;354}355 356template <class F, class T>357BOOST_MATH_GPU_ENABLED inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value && BOOST_MATH_IS_FLOAT(T)358#ifndef BOOST_MATH_HAS_GPU_SUPPORT359&& noexcept(std::declval<F>()(std::declval<T>()))360#endif361)362{363 boost::math::uintmax_t m = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();364 return newton_raphson_iterate(f, guess, min, max, digits, m);365}366 367// TODO(mborland): Disabled for now368// Recursion needs to be removed, but there is no demand at this time369#ifdef BOOST_MATH_HAS_NVRTC370}}} // Namespaces371#else372 373namespace detail {374 375 struct halley_step376 {377 template <class T>378 static T step(const T& /*x*/, const T& f0, const T& f1, const T& f2) noexcept(BOOST_MATH_IS_FLOAT(T))379 {380 using std::fabs;381 T denom = 2 * f0;382 T num = 2 * f1 - f0 * (f2 / f1);383 T delta;384 385 BOOST_MATH_INSTRUMENT_VARIABLE(denom);386 BOOST_MATH_INSTRUMENT_VARIABLE(num);387 388 if ((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>()))389 {390 // possible overflow, use Newton step:391 delta = f0 / f1;392 }393 else394 delta = denom / num;395 return delta;396 }397 };398 399 template <class F, class T>400 T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())));401 402 template <class F, class T>403 T bracket_root_towards_max(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))404 {405 using std::fabs;406 using std::ldexp;407 using std::abs;408 using std::frexp;409 if(count < 2)410 return guess - (max + min) / 2; // Not enough counts left to do anything!!411 //412 // Move guess towards max until we bracket the root, updating min and max as we go:413 //414 int e;415 frexp(max / guess, &e);416 e = abs(e);417 T guess0 = guess;418 T multiplier = e < 64 ? static_cast<T>(2) : static_cast<T>(ldexp(T(1), e / 32));419 T f_current = f0;420 if (fabs(min) < fabs(max))421 {422 while (--count && ((f_current < 0) == (f0 < 0)))423 {424 min = guess;425 guess *= multiplier;426 if (guess > max)427 {428 guess = max;429 f_current = -f_current; // There must be a change of sign!430 break;431 }432 multiplier *= e > 1024 ? 8 : 2;433 unpack_0(f(guess), f_current);434 }435 }436 else437 {438 //439 // If min and max are negative we have to divide to head towards max:440 //441 while (--count && ((f_current < 0) == (f0 < 0)))442 {443 min = guess;444 guess /= multiplier;445 if (guess > max)446 {447 guess = max;448 f_current = -f_current; // There must be a change of sign!449 break;450 }451 multiplier *= e > 1024 ? 8 : 2;452 unpack_0(f(guess), f_current);453 }454 }455 456 if (count)457 {458 max = guess;459 if (multiplier > 16)460 return (guess0 - guess) + bracket_root_towards_min(f, guess, f_current, min, max, count);461 }462 return guess0 - (max + min) / 2;463 }464 465 template <class F, class T>466 T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))467 {468 using std::fabs;469 using std::ldexp;470 using std::abs;471 using std::frexp;472 if (count < 2)473 return guess - (max + min) / 2; // Not enough counts left to do anything!!474 //475 // Move guess towards min until we bracket the root, updating min and max as we go:476 //477 int e;478 frexp(guess / min, &e);479 e = abs(e);480 T guess0 = guess;481 T multiplier = e < 64 ? static_cast<T>(2) : static_cast<T>(ldexp(T(1), e / 32));482 T f_current = f0;483 484 if (fabs(min) < fabs(max))485 {486 while (--count && ((f_current < 0) == (f0 < 0)))487 {488 max = guess;489 guess /= multiplier;490 if (guess < min)491 {492 guess = min;493 f_current = -f_current; // There must be a change of sign!494 break;495 }496 multiplier *= e > 1024 ? 8 : 2;497 unpack_0(f(guess), f_current);498 }499 }500 else501 {502 //503 // If min and max are negative we have to multiply to head towards min:504 //505 while (--count && ((f_current < 0) == (f0 < 0)))506 {507 max = guess;508 guess *= multiplier;509 if (guess < min)510 {511 guess = min;512 f_current = -f_current; // There must be a change of sign!513 break;514 }515 multiplier *= e > 1024 ? 8 : 2;516 unpack_0(f(guess), f_current);517 }518 }519 520 if (count)521 {522 min = guess;523 if (multiplier > 16)524 return (guess0 - guess) + bracket_root_towards_max(f, guess, f_current, min, max, count);525 }526 return guess0 - (max + min) / 2;527 }528 529 template <class Stepper, class F, class T>530 T second_order_root_finder(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))531 {532 BOOST_MATH_STD_USING533 534#ifdef BOOST_MATH_INSTRUMENT535 std::cout << "Second order root iteration, guess = " << guess << ", min = " << min << ", max = " << max536 << ", digits = " << digits << ", max_iter = " << max_iter << "\n";537#endif538 static const char* function = "boost::math::tools::halley_iterate<%1%>";539 if (min >= max)540 {541 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<>());542 }543 544 T f0(0), f1, f2;545 T result = guess;546 547 T factor = ldexp(static_cast<T>(1.0), 1 - digits);548 T delta = (std::max)(T(10000000 * guess), T(10000000)); // arbitrarily large delta549 T last_f0 = 0;550 T delta1 = delta;551 T delta2 = delta;552 bool out_of_bounds_sentry = false;553 554 #ifdef BOOST_MATH_INSTRUMENT555 std::cout << "Second order root iteration, limit = " << factor << "\n";556 #endif557 558 //559 // We use these to sanity check that we do actually bracket a root,560 // we update these to the function value when we update the endpoints561 // of the range. Then, provided at some point we update both endpoints562 // checking that max_range_f * min_range_f <= 0 verifies there is a root563 // to be found somewhere. Note that if there is no root, and we approach 564 // a local minima, then the derivative will go to zero, and hence the next565 // step will jump out of bounds (or at least past the minima), so this566 // check *should* happen in pathological cases.567 //568 T max_range_f = 0;569 T min_range_f = 0;570 571 std::uintmax_t count(max_iter);572 573 do {574 last_f0 = f0;575 delta2 = delta1;576 delta1 = delta;577#ifndef BOOST_MATH_NO_EXCEPTIONS578 try579#endif580 {581 detail::unpack_tuple(f(result), f0, f1, f2);582 }583#ifndef BOOST_MATH_NO_EXCEPTIONS584 catch (const std::overflow_error&)585 {586 f0 = max > 0 ? tools::max_value<T>() : -tools::min_value<T>();587 f1 = f2 = 0;588 }589#endif590 --count;591 592 BOOST_MATH_INSTRUMENT_VARIABLE(f0);593 BOOST_MATH_INSTRUMENT_VARIABLE(f1);594 BOOST_MATH_INSTRUMENT_VARIABLE(f2);595 596 if (0 == f0)597 break;598 if (f1 == 0)599 {600 // Oops zero derivative!!!601 detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max);602 }603 else604 {605 if (f2 != 0)606 {607 delta = Stepper::step(result, f0, f1, f2);608 if (delta * f1 / f0 < 0)609 {610 // Oh dear, we have a problem as Newton and Halley steps611 // disagree about which way we should move. Probably612 // there is cancelation error in the calculation of the613 // Halley step, or else the derivatives are so small614 // that their values are basically trash. We will move615 // in the direction indicated by a Newton step, but616 // by no more than twice the current guess value, otherwise617 // we can jump way out of bounds if we're not careful.618 // See https://svn.boost.org/trac/boost/ticket/8314.619 delta = f0 / f1;620 if (fabs(delta) > 2 * fabs(result))621 delta = (delta < 0 ? -1 : 1) * 2 * fabs(result);622 }623 }624 else625 delta = f0 / f1;626 }627 #ifdef BOOST_MATH_INSTRUMENT628 std::cout << "Second order root iteration, delta = " << delta << ", residual = " << f0 << "\n";629 #endif630 // We need to avoid delta/delta2 overflowing here:631 T convergence = (fabs(delta2) > 1) || (fabs(tools::max_value<T>() * delta2) > fabs(delta)) ? fabs(delta / delta2) : tools::max_value<T>();632 if ((convergence > 0.8) && (convergence < 2))633 {634 // last two steps haven't converged.635 if (fabs(min) < 1 ? fabs(1000 * min) < fabs(max) : fabs(max / min) > 1000)636 {637 if(delta > 0)638 delta = bracket_root_towards_min(f, result, f0, min, max, count);639 else640 delta = bracket_root_towards_max(f, result, f0, min, max, count);641 }642 else643 {644 delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;645 if ((result != 0) && (fabs(delta) > result))646 delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps!647 }648 // reset delta2 so that this branch will *not* be taken on the649 // next iteration:650 delta2 = delta * 3;651 delta1 = delta * 3;652 BOOST_MATH_INSTRUMENT_VARIABLE(delta);653 }654 guess = result;655 result -= delta;656 BOOST_MATH_INSTRUMENT_VARIABLE(result);657 658 // check for out of bounds step:659 if (result < min)660 {661 T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min)))662 ? T(1000)663 : (fabs(min) < 1) && (fabs(tools::max_value<T>() * min) < fabs(result))664 ? ((min < 0) != (result < 0)) ? -tools::max_value<T>() : tools::max_value<T>() : T(result / min);665 if (fabs(diff) < 1)666 diff = 1 / diff;667 if (!out_of_bounds_sentry && (diff > 0) && (diff < 3))668 {669 // Only a small out of bounds step, lets assume that the result670 // is probably approximately at min:671 delta = 0.99f * (guess - min);672 result = guess - delta;673 out_of_bounds_sentry = true; // only take this branch once!674 }675 else676 {677 if (fabs(float_distance(min, max)) < 2)678 {679 result = guess = (min + max) / 2;680 break;681 }682 delta = bracket_root_towards_min(f, guess, f0, min, max, count);683 result = guess - delta;684 if (result <= min)685 result = float_next(min);686 if (result >= max)687 result = float_prior(max);688 guess = min;689 continue;690 }691 }692 else if (result > max)693 {694 T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max);695 if (fabs(diff) < 1)696 diff = 1 / diff;697 if (!out_of_bounds_sentry && (diff > 0) && (diff < 3))698 {699 // Only a small out of bounds step, lets assume that the result700 // is probably approximately at min:701 delta = 0.99f * (guess - max);702 result = guess - delta;703 out_of_bounds_sentry = true; // only take this branch once!704 }705 else706 {707 if (fabs(float_distance(min, max)) < 2)708 {709 result = guess = (min + max) / 2;710 break;711 }712 delta = bracket_root_towards_max(f, guess, f0, min, max, count);713 result = guess - delta;714 if (result >= max)715 result = float_prior(max);716 if (result <= min)717 result = float_next(min);718 guess = min;719 continue;720 }721 }722 // update brackets:723 if (delta > 0)724 {725 max = guess;726 max_range_f = f0;727 }728 else729 {730 min = guess;731 min_range_f = f0;732 }733 //734 // Sanity check that we bracket the root:735 //736 if (max_range_f * min_range_f > 0)737 {738 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<>());739 }740 } while(count && (fabs(result * factor) < fabs(delta)));741 742 max_iter -= count;743 744 #ifdef BOOST_MATH_INSTRUMENT745 std::cout << "Second order root finder required " << max_iter << " iterations.\n";746 #endif747 748 return result;749 }750} // T second_order_root_finder751 752template <class F, class T>753T halley_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))754{755 return detail::second_order_root_finder<detail::halley_step>(f, guess, min, max, digits, max_iter);756}757 758template <class F, class T>759inline T halley_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))760{761 std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)();762 return halley_iterate(f, guess, min, max, digits, m);763}764 765namespace detail {766 767 struct schroder_stepper768 {769 template <class T>770 static T step(const T& x, const T& f0, const T& f1, const T& f2) noexcept(BOOST_MATH_IS_FLOAT(T))771 {772 using std::fabs;773 T ratio = f0 / f1;774 T delta;775 if ((x != 0) && (fabs(ratio / x) < 0.1))776 {777 delta = ratio + (f2 / (2 * f1)) * ratio * ratio;778 // check second derivative doesn't over compensate:779 if (delta * ratio < 0)780 delta = ratio;781 }782 else783 delta = ratio; // fall back to Newton iteration.784 return delta;785 }786 };787 788}789 790template <class F, class T>791T schroder_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))792{793 return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);794}795 796template <class F, class T>797inline T schroder_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))798{799 std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)();800 return schroder_iterate(f, guess, min, max, digits, m);801}802//803// These two are the old spelling of this function, retained for backwards compatibility just in case:804//805template <class F, class T>806T schroeder_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))807{808 return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter);809}810 811template <class F, class T>812inline T schroeder_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy<policies::policy<> >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))813{814 std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)();815 return schroder_iterate(f, guess, min, max, digits, m);816}817 818#ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS819/*820 * Why do we set the default maximum number of iterations to the number of digits in the type?821 * Because for double roots, the number of digits increases linearly with the number of iterations,822 * so this default should recover full precision even in this somewhat pathological case.823 * For isolated roots, the problem is so rapidly convergent that this doesn't matter at all.824 */825template<class ComplexType, class F>826ComplexType complex_newton(F g, ComplexType guess, int max_iterations = std::numeric_limits<typename ComplexType::value_type>::digits)827{828 typedef typename ComplexType::value_type Real;829 using std::norm;830 using std::abs;831 using std::max;832 // z0, z1, and z2 cannot be the same, in case we immediately need to resort to Muller's Method:833 ComplexType z0 = guess + ComplexType(1, 0);834 ComplexType z1 = guess + ComplexType(0, 1);835 ComplexType z2 = guess;836 837 do {838 auto pair = g(z2);839 if (norm(pair.second) == 0)840 {841 // Muller's method. Notation follows Numerical Recipes, 9.5.2:842 ComplexType q = (z2 - z1) / (z1 - z0);843 auto P0 = g(z0);844 auto P1 = g(z1);845 ComplexType qp1 = static_cast<ComplexType>(1) + q;846 ComplexType A = q * (pair.first - qp1 * P1.first + q * P0.first);847 848 ComplexType B = (static_cast<ComplexType>(2) * q + static_cast<ComplexType>(1)) * pair.first - qp1 * qp1 * P1.first + q * q * P0.first;849 ComplexType C = qp1 * pair.first;850 ComplexType rad = sqrt(B * B - static_cast<ComplexType>(4) * A * C);851 ComplexType denom1 = B + rad;852 ComplexType denom2 = B - rad;853 ComplexType correction = (z1 - z2) * static_cast<ComplexType>(2) * C;854 if (norm(denom1) > norm(denom2))855 {856 correction /= denom1;857 }858 else859 {860 correction /= denom2;861 }862 863 z0 = z1;864 z1 = z2;865 z2 = z2 + correction;866 }867 else868 {869 z0 = z1;870 z1 = z2;871 z2 = z2 - (pair.first / pair.second);872 }873 874 // See: https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root875 // If f' is continuous, then convergence of x_n -> x* implies f(x*) = 0.876 // This condition approximates this convergence condition by requiring three consecutive iterates to be clustered.877 Real tol = (max)(abs(z2) * std::numeric_limits<Real>::epsilon(), std::numeric_limits<Real>::epsilon());878 bool real_close = abs(z0.real() - z1.real()) < tol && abs(z0.real() - z2.real()) < tol && abs(z1.real() - z2.real()) < tol;879 bool imag_close = abs(z0.imag() - z1.imag()) < tol && abs(z0.imag() - z2.imag()) < tol && abs(z1.imag() - z2.imag()) < tol;880 if (real_close && imag_close)881 {882 return z2;883 }884 885 } while (max_iterations--);886 887 // The idea is that if we can get abs(f) < eps, we should, but if we go through all these iterations888 // and abs(f) < sqrt(eps), then roundoff error simply does not allow that we can evaluate f to < eps889 // This is somewhat awkward as it isn't scale invariant, but using the Daubechies coefficient example code,890 // I found this condition generates correct roots, whereas the scale invariant condition discussed here:891 // https://scicomp.stackexchange.com/questions/30597/defining-a-condition-number-and-termination-criteria-for-newtons-method892 // allows nonroots to be passed off as roots.893 auto pair = g(z2);894 if (abs(pair.first) < sqrt(std::numeric_limits<Real>::epsilon()))895 {896 return z2;897 }898 899 return { std::numeric_limits<Real>::quiet_NaN(),900 std::numeric_limits<Real>::quiet_NaN() };901}902#endif903 904 905#if !defined(BOOST_MATH_NO_CXX17_IF_CONSTEXPR)906// https://stackoverflow.com/questions/48979861/numerically-stable-method-for-solving-quadratic-equations/50065711907namespace detail908{909#if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1)910inline float fma_workaround(float x, float y, float z) { return ::fmaf(x, y, z); }911inline double fma_workaround(double x, double y, double z) { return ::fma(x, y, z); }912#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS913inline long double fma_workaround(long double x, long double y, long double z) { return ::fmal(x, y, z); }914#endif915#endif 916template<class T>917inline T discriminant(T const& a, T const& b, T const& c)918{919 T w = 4 * a * c;920#if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1)921 T e = fma_workaround(-c, 4 * a, w);922 T f = fma_workaround(b, b, -w);923#else924 T e = std::fma(-c, 4 * a, w);925 T f = std::fma(b, b, -w);926#endif927 return f + e;928}929 930template<class T>931std::pair<T, T> quadratic_roots_imp(T const& a, T const& b, T const& c)932{933#if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1)934 using boost::math::copysign;935#else936 using std::copysign;937#endif938 using std::sqrt;939 if constexpr (std::is_floating_point<T>::value)940 {941 T nan = std::numeric_limits<T>::quiet_NaN();942 if (a == 0)943 {944 if (b == 0 && c != 0)945 {946 return std::pair<T, T>(nan, nan);947 }948 else if (b == 0 && c == 0)949 {950 return std::pair<T, T>(0, 0);951 }952 return std::pair<T, T>(-c / b, -c / b);953 }954 if (b == 0)955 {956 T x0_sq = -c / a;957 if (x0_sq < 0) {958 return std::pair<T, T>(nan, nan);959 }960 T x0 = sqrt(x0_sq);961 return std::pair<T, T>(-x0, x0);962 }963 T discriminant = detail::discriminant(a, b, c);964 // Is there a sane way to flush very small negative values to zero?965 // If there is I don't know of it.966 if (discriminant < 0)967 {968 return std::pair<T, T>(nan, nan);969 }970 T q = -(b + copysign(sqrt(discriminant), b)) / T(2);971 T x0 = q / a;972 T x1 = c / q;973 if (x0 < x1)974 {975 return std::pair<T, T>(x0, x1);976 }977 return std::pair<T, T>(x1, x0);978 }979 else if constexpr (boost::math::tools::is_complex_type<T>::value)980 {981 typename T::value_type nan = std::numeric_limits<typename T::value_type>::quiet_NaN();982 if (a.real() == 0 && a.imag() == 0)983 {984 using std::norm;985 if (b.real() == 0 && b.imag() && norm(c) != 0)986 {987 return std::pair<T, T>({ nan, nan }, { nan, nan });988 }989 else if (b.real() == 0 && b.imag() && c.real() == 0 && c.imag() == 0)990 {991 return std::pair<T, T>({ 0,0 }, { 0,0 });992 }993 return std::pair<T, T>(-c / b, -c / b);994 }995 if (b.real() == 0 && b.imag() == 0)996 {997 T x0_sq = -c / a;998 T x0 = sqrt(x0_sq);999 return std::pair<T, T>(-x0, x0);1000 }1001 // There's no fma for complex types:1002 T discriminant = b * b - T(4) * a * c;1003 T q = -(b + sqrt(discriminant)) / T(2);1004 return std::pair<T, T>(q / a, c / q);1005 }1006 else // Most likely the type is a boost.multiprecision.1007 { //There is no fma for multiprecision, and in addition it doesn't seem to be useful, so revert to the naive computation.1008 T nan = std::numeric_limits<T>::quiet_NaN();1009 if (a == 0)1010 {1011 if (b == 0 && c != 0)1012 {1013 return std::pair<T, T>(nan, nan);1014 }1015 else if (b == 0 && c == 0)1016 {1017 return std::pair<T, T>(0, 0);1018 }1019 return std::pair<T, T>(-c / b, -c / b);1020 }1021 if (b == 0)1022 {1023 T x0_sq = -c / a;1024 if (x0_sq < 0) {1025 return std::pair<T, T>(nan, nan);1026 }1027 T x0 = sqrt(x0_sq);1028 return std::pair<T, T>(-x0, x0);1029 }1030 T discriminant = b * b - 4 * a * c;1031 if (discriminant < 0)1032 {1033 return std::pair<T, T>(nan, nan);1034 }1035 T q = -(b + copysign(sqrt(discriminant), b)) / T(2);1036 T x0 = q / a;1037 T x1 = c / q;1038 if (x0 < x1)1039 {1040 return std::pair<T, T>(x0, x1);1041 }1042 return std::pair<T, T>(x1, x0);1043 }1044}1045} // namespace detail1046 1047template<class T1, class T2 = T1, class T3 = T1>1048inline 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)1049{1050 typedef typename tools::promote_args<T1, T2, T3>::type value_type;1051 return detail::quadratic_roots_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(c));1052}1053 1054#endif1055 1056} // namespace tools1057} // namespace math1058} // namespace boost1059 1060#endif // BOOST_MATH_HAS_NVRTC1061 1062#endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP1063