| // (C) Copyright John Maddock 2006. |
| // (C) Copyright Matt Borland 2024. |
| // Use, modification and distribution are subject to the |
| // Boost Software License, Version 1.0. (See accompanying file |
| // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
| |
| #ifndef BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP |
| #define BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP |
| |
| #ifdef _MSC_VER |
| #pragma once |
| #endif |
| |
| #include <boost/math/tools/config.hpp> |
| #include <boost/math/tools/complex.hpp> // test for multiprecision types in complex Newton |
| #include <boost/math/tools/type_traits.hpp> |
| #include <boost/math/tools/cstdint.hpp> |
| #include <boost/math/tools/numeric_limits.hpp> |
| #include <boost/math/tools/tuple.hpp> |
| #include <boost/math/special_functions/sign.hpp> |
| #include <boost/math/policies/policy.hpp> |
| #include <boost/math/policies/error_handling.hpp> |
| |
| #ifndef BOOST_MATH_HAS_GPU_SUPPORT |
| #include <boost/math/special_functions/next.hpp> |
| #include <boost/math/tools/toms748_solve.hpp> |
| #endif |
| |
| namespace boost { |
| namespace math { |
| namespace tools { |
| |
| namespace detail { |
| |
| namespace dummy { |
| |
| template<int n, class T> |
| BOOST_MATH_GPU_ENABLED typename T::value_type get(const T&) BOOST_MATH_NOEXCEPT(T); |
| } |
| |
| template <class Tuple, class T> |
| BOOST_MATH_GPU_ENABLED void unpack_tuple(const Tuple& t, T& a, T& b) BOOST_MATH_NOEXCEPT(T) |
| { |
| using dummy::get; |
| // Use ADL to find the right overload for get: |
| a = get<0>(t); |
| b = get<1>(t); |
| } |
| template <class Tuple, class T> |
| BOOST_MATH_GPU_ENABLED void unpack_tuple(const Tuple& t, T& a, T& b, T& c) BOOST_MATH_NOEXCEPT(T) |
| { |
| using dummy::get; |
| // Use ADL to find the right overload for get: |
| a = get<0>(t); |
| b = get<1>(t); |
| c = get<2>(t); |
| } |
| |
| template <class Tuple, class T> |
| BOOST_MATH_GPU_ENABLED inline void unpack_0(const Tuple& t, T& val) BOOST_MATH_NOEXCEPT(T) |
| { |
| using dummy::get; |
| // Rely on ADL to find the correct overload of get: |
| val = get<0>(t); |
| } |
| |
| template <class T, class U, class V> |
| BOOST_MATH_GPU_ENABLED inline void unpack_tuple(const boost::math::pair<T, U>& p, V& a, V& b) BOOST_MATH_NOEXCEPT(T) |
| { |
| a = p.first; |
| b = p.second; |
| } |
| template <class T, class U, class V> |
| BOOST_MATH_GPU_ENABLED inline void unpack_0(const boost::math::pair<T, U>& p, V& a) BOOST_MATH_NOEXCEPT(T) |
| { |
| a = p.first; |
| } |
| |
| template <class F, class T> |
| BOOST_MATH_GPU_ENABLED void handle_zero_derivative(F f, |
| T& last_f0, |
| const T& f0, |
| T& delta, |
| T& result, |
| T& guess, |
| const T& min, |
| const T& max) noexcept(BOOST_MATH_IS_FLOAT(T) |
| #ifndef BOOST_MATH_HAS_GPU_SUPPORT |
| && noexcept(std::declval<F>()(std::declval<T>())) |
| #endif |
| ) |
| { |
| if (last_f0 == 0) |
| { |
| // this must be the first iteration, pretend that we had a |
| // previous one at either min or max: |
| if (result == min) |
| { |
| guess = max; |
| } |
| else |
| { |
| guess = min; |
| } |
| unpack_0(f(guess), last_f0); |
| delta = guess - result; |
| } |
| if (sign(last_f0) * sign(f0) < 0) |
| { |
| // we've crossed over so move in opposite direction to last step: |
| if (delta < 0) |
| { |
| delta = (result - min) / 2; |
| } |
| else |
| { |
| delta = (result - max) / 2; |
| } |
| } |
| else |
| { |
| // move in same direction as last step: |
| if (delta < 0) |
| { |
| delta = (result - max) / 2; |
| } |
| else |
| { |
| delta = (result - min) / 2; |
| } |
| } |
| } |
| |
| } // namespace |
| |
| template <class F, class T, class Tol, class Policy> |
| BOOST_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) |
| #ifndef BOOST_MATH_HAS_GPU_SUPPORT |
| && noexcept(std::declval<F>()(std::declval<T>())) |
| #endif |
| ) |
| { |
| T fmin = f(min); |
| T fmax = f(max); |
| if (fmin == 0) |
| { |
| max_iter = 2; |
| return boost::math::make_pair(min, min); |
| } |
| if (fmax == 0) |
| { |
| max_iter = 2; |
| return boost::math::make_pair(max, max); |
| } |
| |
| // |
| // Error checking: |
| // |
| constexpr auto function = "boost::math::tools::bisect<%1%>"; |
| if (min >= max) |
| { |
| return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, |
| "Arguments in wrong order in boost::math::tools::bisect (first arg=%1%)", min, pol)); |
| } |
| if (fmin * fmax >= 0) |
| { |
| return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, |
| "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)); |
| } |
| |
| // |
| // Three function invocations so far: |
| // |
| std::uintmax_t count = max_iter; |
| if (count < 3) |
| count = 0; |
| else |
| count -= 3; |
| |
| while (count && (0 == tol(min, max))) |
| { |
| T mid = (min + max) / 2; |
| T fmid = f(mid); |
| if ((mid == max) || (mid == min)) |
| break; |
| if (fmid == 0) |
| { |
| min = max = mid; |
| break; |
| } |
| else if (sign(fmid) * sign(fmin) < 0) |
| { |
| max = mid; |
| } |
| else |
| { |
| min = mid; |
| fmin = fmid; |
| } |
| --count; |
| } |
| |
| max_iter -= count; |
| |
| #ifdef BOOST_MATH_INSTRUMENT |
| std::cout << "Bisection required " << max_iter << " iterations.\n"; |
| #endif |
| |
| return boost::math::make_pair(min, max); |
| } |
| |
| template <class F, class T, class Tol> |
| BOOST_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) |
| #ifndef BOOST_MATH_HAS_GPU_SUPPORT |
| && noexcept(std::declval<F>()(std::declval<T>())) |
| #endif |
| ) |
| { |
| return bisect(f, min, max, tol, max_iter, policies::policy<>()); |
| } |
| |
| template <class F, class T, class Tol> |
| BOOST_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) |
| #ifndef BOOST_MATH_HAS_GPU_SUPPORT |
| && noexcept(std::declval<F>()(std::declval<T>())) |
| #endif |
| ) |
| { |
| boost::math::uintmax_t m = (boost::math::numeric_limits<boost::math::uintmax_t>::max)(); |
| return bisect(f, min, max, tol, m, policies::policy<>()); |
| } |
| |
| |
| template <class F, class T> |
| BOOST_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) |
| #ifndef BOOST_MATH_HAS_GPU_SUPPORT |
| && noexcept(std::declval<F>()(std::declval<T>())) |
| #endif |
| ) |
| { |
| BOOST_MATH_STD_USING |
| |
| constexpr auto function = "boost::math::tools::newton_raphson_iterate<%1%>"; |
| if (min > max) |
| { |
| 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<>()); |
| } |
| |
| T f0(0), f1, last_f0(0); |
| T result = guess; |
| |
| T factor = static_cast<T>(ldexp(1.0, 1 - digits)); |
| T delta = tools::max_value<T>(); |
| T delta1 = tools::max_value<T>(); |
| T delta2 = tools::max_value<T>(); |
| |
| // |
| // We use these to sanity check that we do actually bracket a root, |
| // we update these to the function value when we update the endpoints |
| // of the range. Then, provided at some point we update both endpoints |
| // checking that max_range_f * min_range_f <= 0 verifies there is a root |
| // to be found somewhere. Note that if there is no root, and we approach |
| // a local minima, then the derivative will go to zero, and hence the next |
| // step will jump out of bounds (or at least past the minima), so this |
| // check *should* happen in pathological cases. |
| // |
| T max_range_f = 0; |
| T min_range_f = 0; |
| |
| boost::math::uintmax_t count(max_iter); |
| |
| #ifdef BOOST_MATH_INSTRUMENT |
| std::cout << "Newton_raphson_iterate, guess = " << guess << ", min = " << min << ", max = " << max |
| << ", digits = " << digits << ", max_iter = " << max_iter << "\n"; |
| #endif |
| |
| do { |
| last_f0 = f0; |
| delta2 = delta1; |
| delta1 = delta; |
| detail::unpack_tuple(f(result), f0, f1); |
| --count; |
| if (0 == f0) |
| break; |
| if (f1 == 0) |
| { |
| // Oops zero derivative!!! |
| detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max); |
| } |
| else |
| { |
| delta = f0 / f1; |
| } |
| #ifdef BOOST_MATH_INSTRUMENT |
| std::cout << "Newton iteration " << max_iter - count << ", delta = " << delta << ", residual = " << f0 << "\n"; |
| #endif |
| if (fabs(delta * 2) > fabs(delta2)) |
| { |
| // Last two steps haven't converged. |
| T shift = (delta > 0) ? (result - min) / 2 : (result - max) / 2; |
| if ((result != 0) && (fabs(shift) > fabs(result))) |
| { |
| delta = sign(delta) * fabs(result); // protect against huge jumps! |
| } |
| else |
| delta = shift; |
| // reset delta1/2 so we don't take this branch next time round: |
| delta1 = 3 * delta; |
| delta2 = 3 * delta; |
| } |
| guess = result; |
| result -= delta; |
| if (result <= min) |
| { |
| delta = 0.5F * (guess - min); |
| result = guess - delta; |
| if ((result == min) || (result == max)) |
| break; |
| } |
| else if (result >= max) |
| { |
| delta = 0.5F * (guess - max); |
| result = guess - delta; |
| if ((result == min) || (result == max)) |
| break; |
| } |
| // Update brackets: |
| if (delta > 0) |
| { |
| max = guess; |
| max_range_f = f0; |
| } |
| else |
| { |
| min = guess; |
| min_range_f = f0; |
| } |
| // |
| // Sanity check that we bracket the root: |
| // |
| if (max_range_f * min_range_f > 0) |
| { |
| 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<>()); |
| } |
| }while(count && (fabs(result * factor) < fabs(delta))); |
| |
| max_iter -= count; |
| |
| #ifdef BOOST_MATH_INSTRUMENT |
| std::cout << "Newton Raphson required " << max_iter << " iterations\n"; |
| #endif |
| |
| return result; |
| } |
| |
| template <class F, class T> |
| BOOST_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) |
| #ifndef BOOST_MATH_HAS_GPU_SUPPORT |
| && noexcept(std::declval<F>()(std::declval<T>())) |
| #endif |
| ) |
| { |
| boost::math::uintmax_t m = (boost::math::numeric_limits<boost::math::uintmax_t>::max)(); |
| return newton_raphson_iterate(f, guess, min, max, digits, m); |
| } |
| |
| // TODO(mborland): Disabled for now |
| // Recursion needs to be removed, but there is no demand at this time |
| #ifdef BOOST_MATH_HAS_NVRTC |
| }}} // Namespaces |
| #else |
| |
| namespace detail { |
| |
| struct halley_step |
| { |
| template <class T> |
| static T step(const T& /*x*/, const T& f0, const T& f1, const T& f2) noexcept(BOOST_MATH_IS_FLOAT(T)) |
| { |
| using std::fabs; |
| T denom = 2 * f0; |
| T num = 2 * f1 - f0 * (f2 / f1); |
| T delta; |
| |
| BOOST_MATH_INSTRUMENT_VARIABLE(denom); |
| BOOST_MATH_INSTRUMENT_VARIABLE(num); |
| |
| if ((fabs(num) < 1) && (fabs(denom) >= fabs(num) * tools::max_value<T>())) |
| { |
| // possible overflow, use Newton step: |
| delta = f0 / f1; |
| } |
| else |
| delta = denom / num; |
| return delta; |
| } |
| }; |
| |
| template <class F, class T> |
| 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>()))); |
| |
| template <class F, class T> |
| 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>()))) |
| { |
| using std::fabs; |
| using std::ldexp; |
| using std::abs; |
| using std::frexp; |
| if(count < 2) |
| return guess - (max + min) / 2; // Not enough counts left to do anything!! |
| // |
| // Move guess towards max until we bracket the root, updating min and max as we go: |
| // |
| int e; |
| frexp(max / guess, &e); |
| e = abs(e); |
| T guess0 = guess; |
| T multiplier = e < 64 ? static_cast<T>(2) : static_cast<T>(ldexp(T(1), e / 32)); |
| T f_current = f0; |
| if (fabs(min) < fabs(max)) |
| { |
| while (--count && ((f_current < 0) == (f0 < 0))) |
| { |
| min = guess; |
| guess *= multiplier; |
| if (guess > max) |
| { |
| guess = max; |
| f_current = -f_current; // There must be a change of sign! |
| break; |
| } |
| multiplier *= e > 1024 ? 8 : 2; |
| unpack_0(f(guess), f_current); |
| } |
| } |
| else |
| { |
| // |
| // If min and max are negative we have to divide to head towards max: |
| // |
| while (--count && ((f_current < 0) == (f0 < 0))) |
| { |
| min = guess; |
| guess /= multiplier; |
| if (guess > max) |
| { |
| guess = max; |
| f_current = -f_current; // There must be a change of sign! |
| break; |
| } |
| multiplier *= e > 1024 ? 8 : 2; |
| unpack_0(f(guess), f_current); |
| } |
| } |
| |
| if (count) |
| { |
| max = guess; |
| if (multiplier > 16) |
| return (guess0 - guess) + bracket_root_towards_min(f, guess, f_current, min, max, count); |
| } |
| return guess0 - (max + min) / 2; |
| } |
| |
| template <class F, class T> |
| 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>()))) |
| { |
| using std::fabs; |
| using std::ldexp; |
| using std::abs; |
| using std::frexp; |
| if (count < 2) |
| return guess - (max + min) / 2; // Not enough counts left to do anything!! |
| // |
| // Move guess towards min until we bracket the root, updating min and max as we go: |
| // |
| int e; |
| frexp(guess / min, &e); |
| e = abs(e); |
| T guess0 = guess; |
| T multiplier = e < 64 ? static_cast<T>(2) : static_cast<T>(ldexp(T(1), e / 32)); |
| T f_current = f0; |
| |
| if (fabs(min) < fabs(max)) |
| { |
| while (--count && ((f_current < 0) == (f0 < 0))) |
| { |
| max = guess; |
| guess /= multiplier; |
| if (guess < min) |
| { |
| guess = min; |
| f_current = -f_current; // There must be a change of sign! |
| break; |
| } |
| multiplier *= e > 1024 ? 8 : 2; |
| unpack_0(f(guess), f_current); |
| } |
| } |
| else |
| { |
| // |
| // If min and max are negative we have to multiply to head towards min: |
| // |
| while (--count && ((f_current < 0) == (f0 < 0))) |
| { |
| max = guess; |
| guess *= multiplier; |
| if (guess < min) |
| { |
| guess = min; |
| f_current = -f_current; // There must be a change of sign! |
| break; |
| } |
| multiplier *= e > 1024 ? 8 : 2; |
| unpack_0(f(guess), f_current); |
| } |
| } |
| |
| if (count) |
| { |
| min = guess; |
| if (multiplier > 16) |
| return (guess0 - guess) + bracket_root_towards_max(f, guess, f_current, min, max, count); |
| } |
| return guess0 - (max + min) / 2; |
| } |
| |
| template <class Stepper, class F, class T> |
| 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>()))) |
| { |
| BOOST_MATH_STD_USING |
| |
| #ifdef BOOST_MATH_INSTRUMENT |
| std::cout << "Second order root iteration, guess = " << guess << ", min = " << min << ", max = " << max |
| << ", digits = " << digits << ", max_iter = " << max_iter << "\n"; |
| #endif |
| static const char* function = "boost::math::tools::halley_iterate<%1%>"; |
| if (min >= max) |
| { |
| 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<>()); |
| } |
| |
| T f0(0), f1, f2; |
| T result = guess; |
| |
| T factor = ldexp(static_cast<T>(1.0), 1 - digits); |
| T delta = (std::max)(T(10000000 * guess), T(10000000)); // arbitrarily large delta |
| T last_f0 = 0; |
| T delta1 = delta; |
| T delta2 = delta; |
| bool out_of_bounds_sentry = false; |
| |
| #ifdef BOOST_MATH_INSTRUMENT |
| std::cout << "Second order root iteration, limit = " << factor << "\n"; |
| #endif |
| |
| // |
| // We use these to sanity check that we do actually bracket a root, |
| // we update these to the function value when we update the endpoints |
| // of the range. Then, provided at some point we update both endpoints |
| // checking that max_range_f * min_range_f <= 0 verifies there is a root |
| // to be found somewhere. Note that if there is no root, and we approach |
| // a local minima, then the derivative will go to zero, and hence the next |
| // step will jump out of bounds (or at least past the minima), so this |
| // check *should* happen in pathological cases. |
| // |
| T max_range_f = 0; |
| T min_range_f = 0; |
| |
| std::uintmax_t count(max_iter); |
| |
| do { |
| last_f0 = f0; |
| delta2 = delta1; |
| delta1 = delta; |
| #ifndef BOOST_MATH_NO_EXCEPTIONS |
| try |
| #endif |
| { |
| detail::unpack_tuple(f(result), f0, f1, f2); |
| } |
| #ifndef BOOST_MATH_NO_EXCEPTIONS |
| catch (const std::overflow_error&) |
| { |
| f0 = max > 0 ? tools::max_value<T>() : -tools::min_value<T>(); |
| f1 = f2 = 0; |
| } |
| #endif |
| --count; |
| |
| BOOST_MATH_INSTRUMENT_VARIABLE(f0); |
| BOOST_MATH_INSTRUMENT_VARIABLE(f1); |
| BOOST_MATH_INSTRUMENT_VARIABLE(f2); |
| |
| if (0 == f0) |
| break; |
| if (f1 == 0) |
| { |
| // Oops zero derivative!!! |
| detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max); |
| } |
| else |
| { |
| if (f2 != 0) |
| { |
| delta = Stepper::step(result, f0, f1, f2); |
| if (delta * f1 / f0 < 0) |
| { |
| // Oh dear, we have a problem as Newton and Halley steps |
| // disagree about which way we should move. Probably |
| // there is cancelation error in the calculation of the |
| // Halley step, or else the derivatives are so small |
| // that their values are basically trash. We will move |
| // in the direction indicated by a Newton step, but |
| // by no more than twice the current guess value, otherwise |
| // we can jump way out of bounds if we're not careful. |
| // See https://svn.boost.org/trac/boost/ticket/8314. |
| delta = f0 / f1; |
| if (fabs(delta) > 2 * fabs(result)) |
| delta = (delta < 0 ? -1 : 1) * 2 * fabs(result); |
| } |
| } |
| else |
| delta = f0 / f1; |
| } |
| #ifdef BOOST_MATH_INSTRUMENT |
| std::cout << "Second order root iteration, delta = " << delta << ", residual = " << f0 << "\n"; |
| #endif |
| // We need to avoid delta/delta2 overflowing here: |
| T convergence = (fabs(delta2) > 1) || (fabs(tools::max_value<T>() * delta2) > fabs(delta)) ? fabs(delta / delta2) : tools::max_value<T>(); |
| if ((convergence > 0.8) && (convergence < 2)) |
| { |
| // last two steps haven't converged. |
| if (fabs(min) < 1 ? fabs(1000 * min) < fabs(max) : fabs(max / min) > 1000) |
| { |
| if(delta > 0) |
| delta = bracket_root_towards_min(f, result, f0, min, max, count); |
| else |
| delta = bracket_root_towards_max(f, result, f0, min, max, count); |
| } |
| else |
| { |
| delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; |
| if ((result != 0) && (fabs(delta) > result)) |
| delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps! |
| } |
| // reset delta2 so that this branch will *not* be taken on the |
| // next iteration: |
| delta2 = delta * 3; |
| delta1 = delta * 3; |
| BOOST_MATH_INSTRUMENT_VARIABLE(delta); |
| } |
| guess = result; |
| result -= delta; |
| BOOST_MATH_INSTRUMENT_VARIABLE(result); |
| |
| // check for out of bounds step: |
| if (result < min) |
| { |
| T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) |
| ? T(1000) |
| : (fabs(min) < 1) && (fabs(tools::max_value<T>() * min) < fabs(result)) |
| ? ((min < 0) != (result < 0)) ? -tools::max_value<T>() : tools::max_value<T>() : T(result / min); |
| if (fabs(diff) < 1) |
| diff = 1 / diff; |
| if (!out_of_bounds_sentry && (diff > 0) && (diff < 3)) |
| { |
| // Only a small out of bounds step, lets assume that the result |
| // is probably approximately at min: |
| delta = 0.99f * (guess - min); |
| result = guess - delta; |
| out_of_bounds_sentry = true; // only take this branch once! |
| } |
| else |
| { |
| if (fabs(float_distance(min, max)) < 2) |
| { |
| result = guess = (min + max) / 2; |
| break; |
| } |
| delta = bracket_root_towards_min(f, guess, f0, min, max, count); |
| result = guess - delta; |
| if (result <= min) |
| result = float_next(min); |
| if (result >= max) |
| result = float_prior(max); |
| guess = min; |
| continue; |
| } |
| } |
| else if (result > max) |
| { |
| T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max); |
| if (fabs(diff) < 1) |
| diff = 1 / diff; |
| if (!out_of_bounds_sentry && (diff > 0) && (diff < 3)) |
| { |
| // Only a small out of bounds step, lets assume that the result |
| // is probably approximately at min: |
| delta = 0.99f * (guess - max); |
| result = guess - delta; |
| out_of_bounds_sentry = true; // only take this branch once! |
| } |
| else |
| { |
| if (fabs(float_distance(min, max)) < 2) |
| { |
| result = guess = (min + max) / 2; |
| break; |
| } |
| delta = bracket_root_towards_max(f, guess, f0, min, max, count); |
| result = guess - delta; |
| if (result >= max) |
| result = float_prior(max); |
| if (result <= min) |
| result = float_next(min); |
| guess = min; |
| continue; |
| } |
| } |
| // update brackets: |
| if (delta > 0) |
| { |
| max = guess; |
| max_range_f = f0; |
| } |
| else |
| { |
| min = guess; |
| min_range_f = f0; |
| } |
| // |
| // Sanity check that we bracket the root: |
| // |
| if (max_range_f * min_range_f > 0) |
| { |
| 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<>()); |
| } |
| } while(count && (fabs(result * factor) < fabs(delta))); |
| |
| max_iter -= count; |
| |
| #ifdef BOOST_MATH_INSTRUMENT |
| std::cout << "Second order root finder required " << max_iter << " iterations.\n"; |
| #endif |
| |
| return result; |
| } |
| } // T second_order_root_finder |
| |
| template <class F, class T> |
| T 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>()))) |
| { |
| return detail::second_order_root_finder<detail::halley_step>(f, guess, min, max, digits, max_iter); |
| } |
| |
| template <class F, class T> |
| inline 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>()))) |
| { |
| std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)(); |
| return halley_iterate(f, guess, min, max, digits, m); |
| } |
| |
| namespace detail { |
| |
| struct schroder_stepper |
| { |
| template <class T> |
| static T step(const T& x, const T& f0, const T& f1, const T& f2) noexcept(BOOST_MATH_IS_FLOAT(T)) |
| { |
| using std::fabs; |
| T ratio = f0 / f1; |
| T delta; |
| if ((x != 0) && (fabs(ratio / x) < 0.1)) |
| { |
| delta = ratio + (f2 / (2 * f1)) * ratio * ratio; |
| // check second derivative doesn't over compensate: |
| if (delta * ratio < 0) |
| delta = ratio; |
| } |
| else |
| delta = ratio; // fall back to Newton iteration. |
| return delta; |
| } |
| }; |
| |
| } |
| |
| template <class F, class T> |
| T 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>()))) |
| { |
| return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter); |
| } |
| |
| template <class F, class T> |
| inline 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>()))) |
| { |
| std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)(); |
| return schroder_iterate(f, guess, min, max, digits, m); |
| } |
| // |
| // These two are the old spelling of this function, retained for backwards compatibility just in case: |
| // |
| template <class F, class T> |
| T 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>()))) |
| { |
| return detail::second_order_root_finder<detail::schroder_stepper>(f, guess, min, max, digits, max_iter); |
| } |
| |
| template <class F, class T> |
| inline 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>()))) |
| { |
| std::uintmax_t m = (std::numeric_limits<std::uintmax_t>::max)(); |
| return schroder_iterate(f, guess, min, max, digits, m); |
| } |
| |
| #ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS |
| /* |
| * Why do we set the default maximum number of iterations to the number of digits in the type? |
| * Because for double roots, the number of digits increases linearly with the number of iterations, |
| * so this default should recover full precision even in this somewhat pathological case. |
| * For isolated roots, the problem is so rapidly convergent that this doesn't matter at all. |
| */ |
| template<class ComplexType, class F> |
| ComplexType complex_newton(F g, ComplexType guess, int max_iterations = std::numeric_limits<typename ComplexType::value_type>::digits) |
| { |
| typedef typename ComplexType::value_type Real; |
| using std::norm; |
| using std::abs; |
| using std::max; |
| // z0, z1, and z2 cannot be the same, in case we immediately need to resort to Muller's Method: |
| ComplexType z0 = guess + ComplexType(1, 0); |
| ComplexType z1 = guess + ComplexType(0, 1); |
| ComplexType z2 = guess; |
| |
| do { |
| auto pair = g(z2); |
| if (norm(pair.second) == 0) |
| { |
| // Muller's method. Notation follows Numerical Recipes, 9.5.2: |
| ComplexType q = (z2 - z1) / (z1 - z0); |
| auto P0 = g(z0); |
| auto P1 = g(z1); |
| ComplexType qp1 = static_cast<ComplexType>(1) + q; |
| ComplexType A = q * (pair.first - qp1 * P1.first + q * P0.first); |
| |
| ComplexType B = (static_cast<ComplexType>(2) * q + static_cast<ComplexType>(1)) * pair.first - qp1 * qp1 * P1.first + q * q * P0.first; |
| ComplexType C = qp1 * pair.first; |
| ComplexType rad = sqrt(B * B - static_cast<ComplexType>(4) * A * C); |
| ComplexType denom1 = B + rad; |
| ComplexType denom2 = B - rad; |
| ComplexType correction = (z1 - z2) * static_cast<ComplexType>(2) * C; |
| if (norm(denom1) > norm(denom2)) |
| { |
| correction /= denom1; |
| } |
| else |
| { |
| correction /= denom2; |
| } |
| |
| z0 = z1; |
| z1 = z2; |
| z2 = z2 + correction; |
| } |
| else |
| { |
| z0 = z1; |
| z1 = z2; |
| z2 = z2 - (pair.first / pair.second); |
| } |
| |
| // See: https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root |
| // If f' is continuous, then convergence of x_n -> x* implies f(x*) = 0. |
| // This condition approximates this convergence condition by requiring three consecutive iterates to be clustered. |
| Real tol = (max)(abs(z2) * std::numeric_limits<Real>::epsilon(), std::numeric_limits<Real>::epsilon()); |
| bool real_close = abs(z0.real() - z1.real()) < tol && abs(z0.real() - z2.real()) < tol && abs(z1.real() - z2.real()) < tol; |
| bool imag_close = abs(z0.imag() - z1.imag()) < tol && abs(z0.imag() - z2.imag()) < tol && abs(z1.imag() - z2.imag()) < tol; |
| if (real_close && imag_close) |
| { |
| return z2; |
| } |
| |
| } while (max_iterations--); |
| |
| // The idea is that if we can get abs(f) < eps, we should, but if we go through all these iterations |
| // and abs(f) < sqrt(eps), then roundoff error simply does not allow that we can evaluate f to < eps |
| // This is somewhat awkward as it isn't scale invariant, but using the Daubechies coefficient example code, |
| // I found this condition generates correct roots, whereas the scale invariant condition discussed here: |
| // https://scicomp.stackexchange.com/questions/30597/defining-a-condition-number-and-termination-criteria-for-newtons-method |
| // allows nonroots to be passed off as roots. |
| auto pair = g(z2); |
| if (abs(pair.first) < sqrt(std::numeric_limits<Real>::epsilon())) |
| { |
| return z2; |
| } |
| |
| return { std::numeric_limits<Real>::quiet_NaN(), |
| std::numeric_limits<Real>::quiet_NaN() }; |
| } |
| #endif |
| |
| |
| #if !defined(BOOST_MATH_NO_CXX17_IF_CONSTEXPR) |
| // https://stackoverflow.com/questions/48979861/numerically-stable-method-for-solving-quadratic-equations/50065711 |
| namespace detail |
| { |
| #if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1) |
| inline float fma_workaround(float x, float y, float z) { return ::fmaf(x, y, z); } |
| inline double fma_workaround(double x, double y, double z) { return ::fma(x, y, z); } |
| #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
| inline long double fma_workaround(long double x, long double y, long double z) { return ::fmal(x, y, z); } |
| #endif |
| #endif |
| template<class T> |
| inline T discriminant(T const& a, T const& b, T const& c) |
| { |
| T w = 4 * a * c; |
| #if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1) |
| T e = fma_workaround(-c, 4 * a, w); |
| T f = fma_workaround(b, b, -w); |
| #else |
| T e = std::fma(-c, 4 * a, w); |
| T f = std::fma(b, b, -w); |
| #endif |
| return f + e; |
| } |
| |
| template<class T> |
| std::pair<T, T> quadratic_roots_imp(T const& a, T const& b, T const& c) |
| { |
| #if defined(BOOST_GNU_STDLIB) && !defined(_GLIBCXX_USE_C99_MATH_TR1) |
| using boost::math::copysign; |
| #else |
| using std::copysign; |
| #endif |
| using std::sqrt; |
| if constexpr (std::is_floating_point<T>::value) |
| { |
| T nan = std::numeric_limits<T>::quiet_NaN(); |
| if (a == 0) |
| { |
| if (b == 0 && c != 0) |
| { |
| return std::pair<T, T>(nan, nan); |
| } |
| else if (b == 0 && c == 0) |
| { |
| return std::pair<T, T>(0, 0); |
| } |
| return std::pair<T, T>(-c / b, -c / b); |
| } |
| if (b == 0) |
| { |
| T x0_sq = -c / a; |
| if (x0_sq < 0) { |
| return std::pair<T, T>(nan, nan); |
| } |
| T x0 = sqrt(x0_sq); |
| return std::pair<T, T>(-x0, x0); |
| } |
| T discriminant = detail::discriminant(a, b, c); |
| // Is there a sane way to flush very small negative values to zero? |
| // If there is I don't know of it. |
| if (discriminant < 0) |
| { |
| return std::pair<T, T>(nan, nan); |
| } |
| T q = -(b + copysign(sqrt(discriminant), b)) / T(2); |
| T x0 = q / a; |
| T x1 = c / q; |
| if (x0 < x1) |
| { |
| return std::pair<T, T>(x0, x1); |
| } |
| return std::pair<T, T>(x1, x0); |
| } |
| else if constexpr (boost::math::tools::is_complex_type<T>::value) |
| { |
| typename T::value_type nan = std::numeric_limits<typename T::value_type>::quiet_NaN(); |
| if (a.real() == 0 && a.imag() == 0) |
| { |
| using std::norm; |
| if (b.real() == 0 && b.imag() && norm(c) != 0) |
| { |
| return std::pair<T, T>({ nan, nan }, { nan, nan }); |
| } |
| else if (b.real() == 0 && b.imag() && c.real() == 0 && c.imag() == 0) |
| { |
| return std::pair<T, T>({ 0,0 }, { 0,0 }); |
| } |
| return std::pair<T, T>(-c / b, -c / b); |
| } |
| if (b.real() == 0 && b.imag() == 0) |
| { |
| T x0_sq = -c / a; |
| T x0 = sqrt(x0_sq); |
| return std::pair<T, T>(-x0, x0); |
| } |
| // There's no fma for complex types: |
| T discriminant = b * b - T(4) * a * c; |
| T q = -(b + sqrt(discriminant)) / T(2); |
| return std::pair<T, T>(q / a, c / q); |
| } |
| else // Most likely the type is a boost.multiprecision. |
| { //There is no fma for multiprecision, and in addition it doesn't seem to be useful, so revert to the naive computation. |
| T nan = std::numeric_limits<T>::quiet_NaN(); |
| if (a == 0) |
| { |
| if (b == 0 && c != 0) |
| { |
| return std::pair<T, T>(nan, nan); |
| } |
| else if (b == 0 && c == 0) |
| { |
| return std::pair<T, T>(0, 0); |
| } |
| return std::pair<T, T>(-c / b, -c / b); |
| } |
| if (b == 0) |
| { |
| T x0_sq = -c / a; |
| if (x0_sq < 0) { |
| return std::pair<T, T>(nan, nan); |
| } |
| T x0 = sqrt(x0_sq); |
| return std::pair<T, T>(-x0, x0); |
| } |
| T discriminant = b * b - 4 * a * c; |
| if (discriminant < 0) |
| { |
| return std::pair<T, T>(nan, nan); |
| } |
| T q = -(b + copysign(sqrt(discriminant), b)) / T(2); |
| T x0 = q / a; |
| T x1 = c / q; |
| if (x0 < x1) |
| { |
| return std::pair<T, T>(x0, x1); |
| } |
| return std::pair<T, T>(x1, x0); |
| } |
| } |
| } // namespace detail |
| |
| template<class T1, class T2 = T1, class T3 = T1> |
| inline 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) |
| { |
| typedef typename tools::promote_args<T1, T2, T3>::type value_type; |
| return detail::quadratic_roots_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(c)); |
| } |
| |
| #endif |
| |
| } // namespace tools |
| } // namespace math |
| } // namespace boost |
| |
| #endif // BOOST_MATH_HAS_NVRTC |
| |
| #endif // BOOST_MATH_TOOLS_NEWTON_SOLVER_HPP |
