From d30f4c50ae2696ecbd485db0d4a66a09e4e9e08f Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 27 Jun 2023 17:09:24 -0400 Subject: [PATCH 01/18] basic structure of solution --- include/boost/math/tools/roots.hpp | 274 +++++++++++++++++++++++++++++ test/test_root_iterations.cpp | 24 ++- 2 files changed, 294 insertions(+), 4 deletions(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index 9fd7bee488..c236709b68 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -24,6 +24,10 @@ #include #include +#include +#include +#include + namespace boost { namespace math { namespace tools { @@ -213,6 +217,275 @@ inline std::pair bisect(F f, T min, T max, Tol tol) noexcept(policies::is_ } + +#if 1 + + +#if 0 + // Newton's method says that: + // dx = f(x)/f'(x) + // x_next = x - dx + // + // Consider f(x) = x, ... f'(x) = 1, so f(x) / f'(x) = x + // x = 1/1, dx = 1/1, x_next = 1/1 - 1/1 = 0 + // x = 0, dx = 0, x_next = 0 - 0 = 0 + // + // Consider f(x) = x^2,... f'(x) = 2x, so f(x) / f'(x) = x/2 + // x = 1/1, dx = 1/2, x_next = 1/1 - 1/2 = 1/2 + // x = 1/2, dx = 1/4, x_next = 1/2 - 1/4 = 1/4 + // This is a double root + // + // Consider f(x) = x^3,... f'(x) = 3x^2, so f(x) / f'(x) = x/3 + // x = 1/1, dx = 1/3, x_next = 1/1 - 1/3 = 2/3 + // x = 2/3, dx = 2/9, x_next = 2/3 - 2/9 = 4/9 + // + // Consider f(x) = x^n, ... f'(x) = n * x^(n-1), so f(x) / f'(x) = x/n + // x = 1/1, dx = 1/n, x_next = 1/1 - 1/n = (n-1)/n + // x = (n-1)/n, dx = (n-1)/n^2, x_next = (n-1)/n - (n-1)/n^2 = (n-1)^2/n^2 + // + // Ratio of two consecutive dx: + // f(x) = x, ratio = 0 + // f(x) = x^2, ratio = 1/2 + // f(x) = x^3, ratio = 2/3 + // + // Ratio of two consecutive dx: + // dx_1 = 1/n + // dx_2 = (n-1)/n^2 + // dx_2 / dx_1 = (n-1) / n + // + + if (is_last_newton_reg) { + if (sign(dx) == sign(dx_prev)) { + // Estimate multiplicity n + T ratio = fabs(dx) / fabs(dx_prev); + T n = 1 / (1 - ratio); + + if (1 < n) { + T dx_want = dx * n; + T x_want = x - dx_want; + if (xl < x_want && x_want < xh) { + is_last_newton_reg = false; + if (is_print) { + std::cout << "TRIGGER -- "; + std::cout << "n: " << n; + std::cout << ", dx: " << dx << ", dx_prev: " << dx_prev; + std::cout << std::endl; + } + dx = dx_want; + x = x_want; + } + } + } + } else { + is_last_newton_reg = true; + } +#endif + + + + + + + + + + + + + + + + + + + + +// = = = = = = = = = = = = = = = = = = = = = = = = = +// = = = = = = = = = = = = = = = = = = = = = = = = = +// = = = = = = = = = = = = = = = = = = = = = = = = = +template +T newton_bracket(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) +{ + BOOST_MATH_STD_USING + static const char* function = "boost::math::tools::newton_bracket<%1%>"; + T factor = static_cast(ldexp(1.0, 1 - digits)); + + std::uintmax_t count(max_iter); + + T f0(0), f1; + T f0_l, f0_h; + detail::unpack_tuple(f(xl), f0_l, f1); + detail::unpack_tuple(f(xh), f0_h, f1); + // count -= 2; + + if (f0_l == 0.0) { // Left bound is solution + return xl; + } else if (f0_h == 0.0) { // Right bound is solution + return xh; + } else if (0 < f0_l * f0_h) { // Root not bracketed + return policies::raise_evaluation_error(function, "Does not bracket boost::math::tools::newton_raphson_iterate %1%", xl, boost::math::policies::policy<>()); + } + + // Set to middle if out of range + if (x < std::min(xl, xh) || std::max(xl, xh) < x) { + x = (xl + xh) / 2; + } + + // Flip l and h if required + if (0 < f0_l) { + std::swap(xl, xh); + std::swap(f0_l, f0_h); + } + + T dx = fabs(xh - xl); + T dx_p = dx * 3; + T dx_pp = dx_p * 3; + + const bool is_print = false; + // const bool is_print = true; + if (is_print) std::cout << "=========================================================" << std::endl; + + do { + if (is_print) std::cout << std::endl; + + // Evaluate + detail::unpack_tuple(f(x), f0, f1); + count--; + + // Exit success + if (fabs(dx) < fabs(x * factor)) { + max_iter -= count; + return x; + } + + // Update brackets + if (f0 < 0.0) { + xl = x; + f0_l = f0; + } else { + xh = x; + f0_h = f0; + } + + // Calculate dx + dx_pp = dx_p; + dx_p = dx; + dx = f0 / f1; + + // Last two steps haven't converged. + if (fabs(dx_pp) < fabs(dx * 2)) { + if (is_print) std::cout << "hack -- "; + T shift = (0 < dx) ? (x - xl) / 2 : (x - xh) / 2; + T dx_des; + if ((x != 0) && (fabs(shift) > fabs(x))) { + dx_des = copysign(x, dx); // Protect against huge jumps! + } else { + dx_des = shift; + } + dx = dx_des; + dx_p = 3 * dx; + dx_pp = 3 * dx; + } + + T x_newton = x - dx; // Calculate x_newton + + const bool is_too_L = x_newton < xl; + const bool is_too_H = x_newton > xh; + if (is_too_L || is_too_H ) { // Out of bounds + if (is_print) std::cout << "bisect -- " << is_too_L << is_too_H << ", "; + dx = (xh - xl) / 2; + x = xl + dx; + } else { // Newton step + if (is_print) std::cout << "newton -- "; + if (x == x_newton || x_newton == xl || x_newton == xh) { + max_iter -= count; + return x; + } + x = x_newton; + } + + if (is_print) { + // std::cout << std::setprecision(21) << std::fixed; + std::cout << "x: " << x << ", f0: " << f0 << ", f1: " << f1 << ", f0/f1: " << f0 / f1 << std::endl; + std::cout << "xl: " << xl << ", f0_l: " << f0_l << std::endl; + std::cout << "xh: " << xh << ", f0_h: " << f0_h << std::endl; + std::cout << "dx: " << dx << ", " << "xh - xl: " << fabs(xh - xl) << ", "; + std::cout << std::endl; + } + } while (count); + + max_iter -= count; + return x; +} + + +// NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN +// NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN +// NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN +template +T newton_raphson_iterate(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) +{ + BOOST_MATH_STD_USING + + static const char* function = "boost::math::tools::newton_raphson_iterate<%1%>"; + if (xl > xh) { + return policies::raise_evaluation_error(function, "Range arguments in wrong order in boost::math::tools::newton_raphson_iterate(first arg=%1%)", xl, boost::math::policies::policy<>()); + } + + T f0(0), f1, last_f0(0); + T x_prev = x; + + T delta = 0.0; + + T f0_max; + T f0_min; + detail::unpack_tuple(f(xh), f0_max, f1); + detail::unpack_tuple(f(xl), f0_min, f1); + + std::uintmax_t count(max_iter); + + // delta1 = delta; + detail::unpack_tuple(f(x), f0, f1); + --count; + + // Exit success + if (0 == f0) { + max_iter -= count; + return x; + // Oops zero derivative!!! + } else if (f1 == 0) { + detail::handle_zero_derivative(f, last_f0, f0, delta, x, x_prev, xl, xh); + // Normal case + } else { + delta = f0 / f1; + } + + // Update x + x_prev = x; + x -= delta; + + if (delta < 0) { + if (0 < f0 * f0_max) { + // TODO: bracket other side if possible + 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%", x_prev, boost::math::policies::policy<>()); + } else { + return newton_bracket(f, x, x_prev, xh, digits, max_iter); + } + } else { + if (0 < f0 * f0_min) { + // TODO: bracket other side if possible + 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%", x_prev, boost::math::policies::policy<>()); + } else { + return newton_bracket(f, x, xl, x_prev, digits, max_iter); + } + } +} + + + + + +#else template T newton_raphson_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { @@ -331,6 +604,7 @@ T newton_raphson_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t& return result; } +#endif template inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) diff --git a/test/test_root_iterations.cpp b/test/test_root_iterations.cpp index fada5d9bb0..d5c5c80bbf 100644 --- a/test/test_root_iterations.cpp +++ b/test/test_root_iterations.cpp @@ -169,7 +169,7 @@ BOOST_AUTO_TEST_CASE( test_main ) // Newton next: iters = 1000; dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, guess / 2, result * 10, newton_limits, iters); - BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits::epsilon() * 2); + BOOST_CHECK_CLOSE_FRACTION(dr, result, ldexp(1.0, 1 - newton_limits)); BOOST_CHECK_LE(iters, 12); // Halley next: iters = 1000; @@ -192,7 +192,7 @@ BOOST_AUTO_TEST_CASE( test_main ) // Newton next: iters = 1000; dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters); - BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits::epsilon() * 2); + BOOST_CHECK_CLOSE_FRACTION(dr, result, ldexp(1.0, 1 - newton_limits)); BOOST_CHECK_LE(iters, 12); // Halley next: iters = 1000; @@ -212,11 +212,16 @@ BOOST_AUTO_TEST_CASE( test_main ) r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance(), iters); BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits::epsilon() * 4); BOOST_CHECK_LE(iters, 12); + + // Newton next: iters = 1000; dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters); - BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits::epsilon() * 2); + BOOST_CHECK_CLOSE_FRACTION(dr, result, ldexp(1.0, 1 - newton_limits)); BOOST_CHECK_LE(iters, 5); + // BOOST_CHECK_LE(iters, 8); + + // Halley next: iters = 1000; dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); @@ -235,11 +240,18 @@ BOOST_AUTO_TEST_CASE( test_main ) r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance(), iters); BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits::epsilon() * 4); BOOST_CHECK_LE(iters, 12); + + // ******************************************************************** + // ******************************************************************** + // ******************************************************************** // Newton next: iters = 1000; dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters); - BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits::epsilon() * 2); + BOOST_CHECK_CLOSE_FRACTION(dr, result, ldexp(1.0, 1 - newton_limits)); BOOST_CHECK_LE(iters, 5); + // BOOST_CHECK_LE(iters, 8); + + // Halley next: iters = 1000; dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); @@ -264,6 +276,8 @@ BOOST_AUTO_TEST_CASE( test_main ) # include "ibeta_small_data.ipp" + +#if 1 for (unsigned i = 0; i < ibeta_small_data.size(); ++i) { // @@ -318,5 +332,7 @@ BOOST_AUTO_TEST_CASE( test_main ) } } +#endif + } From 2cb0f125954bde6954060aa28bc8ee3d1adcc0c2 Mon Sep 17 00:00:00 2001 From: Ryan Date: Wed, 28 Jun 2023 16:08:29 -0400 Subject: [PATCH 02/18] created helper class for bracket --- include/boost/math/tools/roots.hpp | 503 ++++++++--------------------- 1 file changed, 143 insertions(+), 360 deletions(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index c236709b68..3b69de2175 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -216,395 +216,178 @@ inline std::pair bisect(F f, T min, T max, Tol tol) noexcept(policies::is_ return bisect(f, min, max, tol, m, policies::policy<>()); } +namespace detail { + // Stores x, f(x), and f'(x) + template + class CacheOrder1 { + public: + template + CacheOrder1(F f, T x) : x_(x) { detail::unpack_tuple(f(x_), f0_, f1_); } + + T x() const { return x_; } + T f0() const { return f0_; } // f(x) + T f1() const { return f1_; } // f'(x) + + private: + T x_; + T f0_; // f(x) + T f1_; // f'(x) + }; + // Calculates the step with Newton's method + template + class StepNewton { + public: + T calc_dx(const CacheOrder1& cache) { + return cache.f0() / cache.f1(); + } + }; -#if 1 - + // Enforces bracketing when root finding with Newton's method + template + class BracketHelper { + public: + static T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { + BracketHelper helper(f, x, xl, xh, digits, max_iter); + T x_sol = helper.solve_impl(); + max_iter -= helper.count(); // Make max_iter the number of iterations used. + return x_sol; + } -#if 0 - // Newton's method says that: - // dx = f(x)/f'(x) - // x_next = x - dx - // - // Consider f(x) = x, ... f'(x) = 1, so f(x) / f'(x) = x - // x = 1/1, dx = 1/1, x_next = 1/1 - 1/1 = 0 - // x = 0, dx = 0, x_next = 0 - 0 = 0 - // - // Consider f(x) = x^2,... f'(x) = 2x, so f(x) / f'(x) = x/2 - // x = 1/1, dx = 1/2, x_next = 1/1 - 1/2 = 1/2 - // x = 1/2, dx = 1/4, x_next = 1/2 - 1/4 = 1/4 - // This is a double root - // - // Consider f(x) = x^3,... f'(x) = 3x^2, so f(x) / f'(x) = x/3 - // x = 1/1, dx = 1/3, x_next = 1/1 - 1/3 = 2/3 - // x = 2/3, dx = 2/9, x_next = 2/3 - 2/9 = 4/9 - // - // Consider f(x) = x^n, ... f'(x) = n * x^(n-1), so f(x) / f'(x) = x/n - // x = 1/1, dx = 1/n, x_next = 1/1 - 1/n = (n-1)/n - // x = (n-1)/n, dx = (n-1)/n^2, x_next = (n-1)/n - (n-1)/n^2 = (n-1)^2/n^2 - // - // Ratio of two consecutive dx: - // f(x) = x, ratio = 0 - // f(x) = x^2, ratio = 1/2 - // f(x) = x^3, ratio = 2/3 - // - // Ratio of two consecutive dx: - // dx_1 = 1/n - // dx_2 = (n-1)/n^2 - // dx_2 / dx_1 = (n-1) / n - // + private: + // Stores the last two Newton step sizes + class StepSizeHistory { + public: + StepSizeHistory() { reset(); } - if (is_last_newton_reg) { - if (sign(dx) == sign(dx_prev)) { - // Estimate multiplicity n - T ratio = fabs(dx) / fabs(dx_prev); - T n = 1 / (1 - ratio); - - if (1 < n) { - T dx_want = dx * n; - T x_want = x - dx_want; - if (xl < x_want && x_want < xh) { - is_last_newton_reg = false; - if (is_print) { - std::cout << "TRIGGER -- "; - std::cout << "n: " << n; - std::cout << ", dx: " << dx << ", dx_prev: " << dx_prev; - std::cout << std::endl; - } - dx = dx_want; - x = x_want; - } - } - } - } else { - is_last_newton_reg = true; + void reset() { + dx_p_ = (std::numeric_limits::max)(); + dx_pp_ = (std::numeric_limits::max)(); } -#endif - - - - - - - - - - - - - - - - - + void push(T input) { + dx_pp_ = dx_p_; + dx_p_ = input; + } + T dx_p() const { return dx_p_; } + T dx_pp() const { return dx_pp_; } + + private: + T dx_p_; + T dx_pp_; + }; + + // Same as the inputs to newton_raphson_iterate + BracketHelper(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) + : f_(f) + , cache_x_(f, x) + , cache_xl_(f, xl) + , cache_xh_(f, xh) + , count_(max_iter) + , factor_(static_cast(ldexp(1.0, 1 - digits))) // factor_ = 1.0 * 2^(1-digits) + { + if (xh < xl) { + static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; + policies::raise_evaluation_error(function, "Range arguments in wrong order in boost::math::tools::detail::BracketHelper(first arg=%1%)", xl, boost::math::policies::policy<>()); + } + } -// = = = = = = = = = = = = = = = = = = = = = = = = = -// = = = = = = = = = = = = = = = = = = = = = = = = = -// = = = = = = = = = = = = = = = = = = = = = = = = = -template -T newton_bracket(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) -{ - BOOST_MATH_STD_USING - static const char* function = "boost::math::tools::newton_bracket<%1%>"; - T factor = static_cast(ldexp(1.0, 1 - digits)); - - std::uintmax_t count(max_iter); - - T f0(0), f1; - T f0_l, f0_h; - detail::unpack_tuple(f(xl), f0_l, f1); - detail::unpack_tuple(f(xh), f0_h, f1); - // count -= 2; - - if (f0_l == 0.0) { // Left bound is solution - return xl; - } else if (f0_h == 0.0) { // Right bound is solution - return xh; - } else if (0 < f0_l * f0_h) { // Root not bracketed - return policies::raise_evaluation_error(function, "Does not bracket boost::math::tools::newton_raphson_iterate %1%", xl, boost::math::policies::policy<>()); - } - - // Set to middle if out of range - if (x < std::min(xl, xh) || std::max(xl, xh) < x) { - x = (xl + xh) / 2; - } - - // Flip l and h if required - if (0 < f0_l) { - std::swap(xl, xh); - std::swap(f0_l, f0_h); - } - - T dx = fabs(xh - xl); - T dx_p = dx * 3; - T dx_pp = dx_p * 3; - - const bool is_print = false; - // const bool is_print = true; - if (is_print) std::cout << "=========================================================" << std::endl; - - do { - if (is_print) std::cout << std::endl; - - // Evaluate - detail::unpack_tuple(f(x), f0, f1); - count--; + T solve_impl() { + if (0 == f0()) return x(); // Solved at initial x - // Exit success - if (fabs(dx) < fabs(x * factor)) { - max_iter -= count; - return x; - } + // Calculate direction towards root + const T dx = Step().calc_dx(cache_x_); + const T x_next = x() - dx; // x_next calculated with step dx - // Update brackets - if (f0 < 0.0) { - xl = x; - f0_l = f0; - } else { - xh = x; - f0_h = f0; - } + const bool can_bracket_l = f0() * f0_l() <= 0; // Sign flip between x and xl + const bool can_bracket_h = f0() * f0_h() <= 0; // Sign flip between x and xh + const bool want_to_go_l = 0 <= dx; // Because of negative sign in update formula x_next = x() - dx - // Calculate dx - dx_pp = dx_p; - dx_p = dx; - dx = f0 / f1; - - // Last two steps haven't converged. - if (fabs(dx_pp) < fabs(dx * 2)) { - if (is_print) std::cout << "hack -- "; - T shift = (0 < dx) ? (x - xl) / 2 : (x - xh) / 2; - T dx_des; - if ((x != 0) && (fabs(shift) > fabs(x))) { - dx_des = copysign(x, dx); // Protect against huge jumps! + if (can_bracket_l && (want_to_go_l || !can_bracket_h)) { + set_cache_H_to_cache_x(); // Evaluated point x becomes new high bound + return solve_bracket(x_next); + } else if (can_bracket_h) { + set_cache_L_to_cache_x(); // Evaluated point x becomes new low bound + return solve_bracket(x_next); } else { - dx_des = shift; - } - dx = dx_des; - dx_p = 3 * dx; - dx_pp = 3 * dx; - } - - T x_newton = x - dx; // Calculate x_newton - - const bool is_too_L = x_newton < xl; - const bool is_too_H = x_newton > xh; - if (is_too_L || is_too_H ) { // Out of bounds - if (is_print) std::cout << "bisect -- " << is_too_L << is_too_H << ", "; - dx = (xh - xl) / 2; - x = xl + dx; - } else { // Newton step - if (is_print) std::cout << "newton -- "; - if (x == x_newton || x_newton == xl || x_newton == xh) { - max_iter -= count; - return x; + static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; + return policies::raise_evaluation_error(function, "There appears to be no root to be found in boost::math::tools::detail::BracketHelper, perhaps we have a local minima near current best guess of %1%", x(), boost::math::policies::policy<>()); } - x = x_newton; } - if (is_print) { - // std::cout << std::setprecision(21) << std::fixed; - std::cout << "x: " << x << ", f0: " << f0 << ", f1: " << f1 << ", f0/f1: " << f0 / f1 << std::endl; - std::cout << "xl: " << xl << ", f0_l: " << f0_l << std::endl; - std::cout << "xh: " << xh << ", f0_h: " << f0_h << std::endl; - std::cout << "dx: " << dx << ", " << "xh - xl: " << fabs(xh - xl) << ", "; - std::cout << std::endl; - } - } while (count); - - max_iter -= count; - return x; -} + T solve_bracket(T x_next) { + BOOST_MATH_STD_USING + static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; + if (f0_l() == 0.0) { // Low bound is root + return xl(); + } else if (f0_h() == 0.0) { // High bound is root + return xh(); + } else if (0 < f0_l() * f0_h()) { // Bounds do not bracket root + return policies::raise_evaluation_error(function, "Does not bracket boost::math::tools::detail::BracketHelper %1%", xl(), boost::math::policies::policy<>()); + } + + // Orient cache_l and cache_h so that f0_l is negative and f0_h is positive + if (0 < f0_l()) std::swap(cache_xl_, cache_xh_); -// NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN -// NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN -// NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN -template -T newton_raphson_iterate(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) -{ - BOOST_MATH_STD_USING + do { + count_--; - static const char* function = "boost::math::tools::newton_raphson_iterate<%1%>"; - if (xl > xh) { - return policies::raise_evaluation_error(function, "Range arguments in wrong order in boost::math::tools::newton_raphson_iterate(first arg=%1%)", xl, boost::math::policies::policy<>()); - } + // Calculate Newton step + const T dx_newton = Step().calc_dx(cache_x_); + x_next = x() - dx_newton; // Checked against bracket bounds below - T f0(0), f1, last_f0(0); - T x_prev = x; + if (!is_x_between_l_and_h(x_next) || + fabs(dx_pp()) < fabs(dx_newton * 4)) { // Bisection halves step size each iteration + x_next = (xl() + xh()) / 2; // Midpoint + dx_history_.reset(); // Invalitade step size history + } - T delta = 0.0; - - T f0_max; - T f0_min; - detail::unpack_tuple(f(xh), f0_max, f1); - detail::unpack_tuple(f(xl), f0_min, f1); + const T dx = x_next - x(); + cache_x_ = detail::CacheOrder1(f_, x_next); // Update cache_x_ with x_next - std::uintmax_t count(max_iter); + if (fabs(dx) < fabs(x() * factor_)) { break; } // Tolerance met + if (x_next == xl() || x_next == xh() || dx == 0) { break; } + + // Update bracket + (f0() < 0.0) ? set_cache_L_to_cache_x() : set_cache_H_to_cache_x(); - // delta1 = delta; - detail::unpack_tuple(f(x), f0, f1); - --count; + dx_history_.push(dx); // Store step size + } while (count_); - // Exit success - if (0 == f0) { - max_iter -= count; - return x; - // Oops zero derivative!!! - } else if (f1 == 0) { - detail::handle_zero_derivative(f, last_f0, f0, delta, x, x_prev, xl, xh); - // Normal case - } else { - delta = f0 / f1; - } - - // Update x - x_prev = x; - x -= delta; - - if (delta < 0) { - if (0 < f0 * f0_max) { - // TODO: bracket other side if possible - 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%", x_prev, boost::math::policies::policy<>()); - } else { - return newton_bracket(f, x, x_prev, xh, digits, max_iter); + return x(); } - } else { - if (0 < f0 * f0_min) { - // TODO: bracket other side if possible - 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%", x_prev, boost::math::policies::policy<>()); - } else { - return newton_bracket(f, x, xl, x_prev, digits, max_iter); - } - } -} - - - + bool is_x_between_l_and_h(T x) const { return (std::min(xl(), xh()) <= x && x <= std::max(xl(), xh())); } + + void set_cache_L_to_cache_x() { cache_xl_ = cache_x_; } + void set_cache_H_to_cache_x() { cache_xh_ = cache_x_; } + + T x() const { return cache_x_.x(); } + T f0() const { return cache_x_.f0(); } + T xl() const { return cache_xl_.x(); } + T f0_l() const { return cache_xl_.f0(); } + T xh() const { return cache_xh_.x(); } + T f0_h() const { return cache_xh_.f0(); } + T dx_pp() const { return dx_history_.dx_pp(); } + std::uintmax_t count() const { return count_; } + + F f_; + detail::CacheOrder1 cache_x_; + detail::CacheOrder1 cache_xl_; + detail::CacheOrder1 cache_xh_; + std::uintmax_t count_; + T factor_; + StepSizeHistory dx_history_; + }; +} // namespace detail -#else template -T newton_raphson_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) +T newton_raphson_iterate(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { - BOOST_MATH_STD_USING - - static const char* 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(ldexp(1.0, 1 - digits)); - T delta = tools::max_value(); - T delta1 = tools::max_value(); - T delta2 = tools::max_value(); - - // - // 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); - -#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) * 1.1f; // Protect against huge jumps! - //delta = sign(delta) * result; // Protect against huge jumps! Failed for negative result. https://github.com/boostorg/math/issues/216 - } - 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; + return detail::BracketHelper>::solve(f, x, xl, xh, digits, max_iter); } -#endif template inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) From 7bb9a80bc2ea4d0b82d32b8bf950493654c8c743 Mon Sep 17 00:00:00 2001 From: Ryan Date: Wed, 28 Jun 2023 16:09:08 -0400 Subject: [PATCH 03/18] added unit test new behavior --- test/test_root_iterations.cpp | 5 ----- test/test_roots.cpp | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 34 insertions(+), 5 deletions(-) diff --git a/test/test_root_iterations.cpp b/test/test_root_iterations.cpp index d5c5c80bbf..f998b6d724 100644 --- a/test/test_root_iterations.cpp +++ b/test/test_root_iterations.cpp @@ -276,8 +276,6 @@ BOOST_AUTO_TEST_CASE( test_main ) # include "ibeta_small_data.ipp" - -#if 1 for (unsigned i = 0; i < ibeta_small_data.size(); ++i) { // @@ -331,8 +329,5 @@ BOOST_AUTO_TEST_CASE( test_main ) BOOST_CHECK_LE(iters, 10); } } - -#endif - } diff --git a/test/test_roots.cpp b/test/test_roots.cpp index dda14c2e6f..7dcfa1201a 100644 --- a/test/test_roots.cpp +++ b/test/test_roots.cpp @@ -437,6 +437,38 @@ void test_beta(T, const char* /* name */) test_inverses(ibeta_large_data); } +void test_newton_bracket() { + int newton_limits = static_cast(std::numeric_limits::digits * 0.6); + const auto tol = ldexp(1, 1 - newton_limits); + + // Function from #808 + const auto quartic_pos = [](const double x) { + const auto y = (x * x - 1) * (x * x + 0.1); + const auto dy = 2 * x * ((x * x - 1) + (x * x + 0.1)); + return std::make_pair(y, dy); + }; + + // Negation of the above + const auto quartic_neg = [&](const double x) { + const auto p = quartic_pos(x); + return std::make_pair(-p.first, -p.second); + }; + + const auto fn_test_bracket = [&](const auto fn_quartic) { + // Initial search direction has no root + BOOST_CHECK_CLOSE(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.5, 0.1, 1.1, newton_limits), tol); + // Initial search direction has root + BOOST_CHECK_CLOSE(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.8, 0.1, 1.1, newton_limits), tol); + // Initial search direction has root + BOOST_CHECK_CLOSE(-1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.5, -1.1, 1.1, newton_limits), tol); + // Initial search direction has root + BOOST_CHECK_CLOSE(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.8, -1.1, 1.1, newton_limits), tol); + }; + + fn_test_bracket(quartic_pos); + fn_test_bracket(quartic_neg); +} + #if !defined(BOOST_NO_CXX11_AUTO_DECLARATIONS) && !defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) && !defined(BOOST_NO_CXX11_LAMBDAS) template void test_complex_newton() @@ -654,6 +686,8 @@ BOOST_AUTO_TEST_CASE( test_main ) test_beta(0.1, "double"); + test_newton_bracket(); + #if !defined(BOOST_NO_CXX11_AUTO_DECLARATIONS) && !defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) && !defined(BOOST_NO_CXX11_LAMBDAS) test_complex_newton>(); test_complex_newton>(); From 4da9504347bc05910eb0d6603a483a54f9440ac2 Mon Sep 17 00:00:00 2001 From: Ryan Date: Thu, 6 Jul 2023 19:07:13 -0400 Subject: [PATCH 04/18] passes unit tests locally, but needs refactor --- include/boost/math/special_functions/beta.hpp | 17 +- include/boost/math/tools/roots.hpp | 2451 ++++++++++++++++- test/test_root_iterations.cpp | 16 +- test/test_roots.cpp | 10 +- 4 files changed, 2433 insertions(+), 61 deletions(-) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 77ca35dcd8..3caa736b89 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -1021,18 +1021,15 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de BOOST_MATH_ASSERT((p_derivative == 0) || normalised); if(!(boost::math::isfinite)(a)) - return policies::raise_domain_error(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol); + return policies::raise_domain_error(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol); if(!(boost::math::isfinite)(b)) - return policies::raise_domain_error(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol); - if(!(boost::math::isfinite)(x)) + return policies::raise_domain_error(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol); + if (!(0 <= x && x <= 1)) return policies::raise_domain_error(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol); if(p_derivative) *p_derivative = -1; // value not set. - if((x < 0) || (x > 1)) - return policies::raise_domain_error(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol); - if(normalised) { if(a < 0) @@ -1422,18 +1419,16 @@ T ibeta_derivative_imp(T a, T b, T x, const Policy& pol) // start with the usual error checks: // if (!(boost::math::isfinite)(a)) - return policies::raise_domain_error(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol); + return policies::raise_domain_error(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol); if (!(boost::math::isfinite)(b)) - return policies::raise_domain_error(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol); - if (!(boost::math::isfinite)(x)) + return policies::raise_domain_error(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol); + if (!(0 <= x && x <= 1)) return policies::raise_domain_error(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol); if(a <= 0) return policies::raise_domain_error(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol); if(b <= 0) return policies::raise_domain_error(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol); - if((x < 0) || (x > 1)) - return policies::raise_domain_error(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol); // // Now the corner cases: // diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index 3b69de2175..b70befbb4f 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -24,9 +24,8 @@ #include #include +#include #include -#include -#include namespace boost { namespace math { @@ -216,45 +215,2099 @@ inline std::pair bisect(F f, T min, T max, Tol tol) noexcept(policies::is_ return bisect(f, min, max, tol, m, policies::policy<>()); } + +#if 0 + + // + template + class Optional { + public: + Optional() : has_data_(false) {} + Optional(C c) : has_data_(true), c_(c) {} + + bool has_data() const { return has_data_; } + const C& c() const { + assert(has_data_); + return c_; + } + + void replace(C c) { c_ = c; } + + private: + bool has_data_; + C c_; + }; + + // Enum for bracketing states: case_1, case_2U, case_2B + enum class BracketState { case_1, case_2U, case_2B, case_2I }; + + // Enum for solve outputs + enum class OutputResult { success, failure }; + + // Solve output type + template + using SolveOutput = std::pair; + + + + + + +#endif + + + + + + + + + + + + + namespace detail { - // Stores x, f(x), and f'(x) + + +#if 1 + + // constexpr bool DEBUG_PRINT = false; + constexpr bool DEBUG_PRINT = true; + + + // Stores the last two Newton step sizes template + class StepSizeHistory { + public: + StepSizeHistory() { reset(); } + + void reset() { + // NOTE: for T = boost::math::concepts::real_concept, both std::numeric_limits::infinity() + // and std::numeric_limits::max() are 0 + dx_p_ = static_cast(1) / 0; + dx_pp_ = static_cast(1) / 0; + } + + void push(T input) { + dx_pp_ = dx_p_; + dx_p_ = input; + } + + T dx_p() const { return dx_p_; } + T dx_pp() const { return dx_pp_; } + + private: + T dx_p_; + T dx_pp_; + }; + + // Stores x, f(x), and f'(x) + template class CacheOrder1 { public: template - CacheOrder1(F f, T x) : x_(x) { detail::unpack_tuple(f(x_), f0_, f1_); } + CacheOrder1(F f, T x) : x_(x), has_f0_f1_(true) { detail::unpack_tuple(f(x_), f0_, f1_); } + CacheOrder1(T x) : x_(x) , has_f0_f1_(false) {} + + T calc_dx() const { return Step().calc_dx(*this); } + + bool is_l() const { return 0 <= calc_dx(); } + bool is_h() const { return !is_l(); } + + bool is_B(const CacheOrder1& z) const { return !is_same_sign(z); } + bool is_I(const CacheOrder1& z) const { return !is_B(z) && !is_U(z); } + bool is_U(const CacheOrder1& z) const { + if (is_ordered(z)) { + return is_U_ordered(*this, z); + } else { + return is_U_ordered(z, *this); + } + } + std::pair, CacheOrder1> sort(const CacheOrder1& z) const { + if (is_ordered(z)) { + return std::make_pair(*this, z); + } else { + return std::make_pair(z, *this); + } + } + bool is_ordered(const CacheOrder1& z) const { return x() < z.x(); } + bool is_U_ordered(const CacheOrder1& l, const CacheOrder1& h) const { + const bool is_pos = 0 < l.f0(); + + T l_slope = l.f1(); + T h_slope = h.f1(); + + if (!is_pos) { + l_slope *= -1; + h_slope *= -1; + } + + return (l_slope < 0) && (0 < h_slope); + } + bool is_smaller(const CacheOrder1& z) const { + using std::fabs; + return fabs(f0_) < fabs(z.f0()); + } + + bool is_same_sign(const CacheOrder1& z) const { return sign(f0_) == sign(z.f0()); } + bool is_same_slope(const CacheOrder1& z) const { return sign(f1_) == sign(z.f1()); } + + void print() const { + if (!DEBUG_PRINT) return; + + std::cout << "x: " << x_; + if (has_f0_f1_) { + assert(has_f0_f1_); + std::cout << ", f0: " << f0_ << ", f1: " << f1_; + } + std::cout << std::endl; + } + + bool is_x_only() const { return !has_f0_f1_; } + bool is_x_f0_f1() const { return has_f0_f1_; } + + bool is_dx_small(const T& factor) const { + // |dx| ≤ |x| * factor + // |f0| / |f1| ≤ |x| * factor + // |f0| ≤ |x| * |f1| * factor + using std::fabs; + + const bool bool_eval = fabs(f0_) <= fabs(x_ * f1_) * factor; + const bool bool_eval_2 = fabs(f0_) <= fabs(x_ * f1_ * factor); + + if (DEBUG_PRINT) { + std::cout.precision(200); + std::cout << "fabs(f0_): " << fabs(f0_) << std::endl; + std::cout << "fabs(x_ * f1_): " << fabs(x_ * f1_) << std::endl; + std::cout << "factor: " << factor << std::endl; + std::cout << "bool_eval; " << bool_eval << std::endl; + std::cout << "bool_eval_2; " << bool_eval_2 << std::endl; + std::cout << "fabs(f0_) == 0 " << (fabs(f0_) == 0) << std::endl; + std::cout << "fabs(x_ * f1_) * factor == 0 " << (fabs(x_ * f1_) * factor == 0) << std::endl; + + // Print data type T to std::cout + std::cout << "T: " << typeid(T).name() << std::endl; + std::cout << std::endl; + } + + + return bool_eval; + } + + T x() const { return x_; } + T f0() const { assert(has_f0_f1_); return f0_; } // f(x) + T f1() const { assert(has_f0_f1_); return f1_; } // f'(x) + + private: + T x_; + T f0_; // f(x) + T f1_; // f'(x) + bool has_f0_f1_; + }; + + // + // + template + class NEXT { + public: + NEXT static FromX(T x) { + return NEXT(x); + } + + NEXT static TakeStep(const CacheOrder1& cx) { + const T dx = cx.calc_dx(); + return NEXT(cx.x() + dx, dx); + } + + NEXT static From_x_dx(T x, T dx) { + return NEXT(x, dx); + } + + T x() const { return x_; } + T dx() const { return dx_; } + + void print() const { + if (!DEBUG_PRINT) return; + std::cout << "X_DX_NEXT -- x: " << x_; + std::cout << ", dx: " << dx_; + std::cout << std::endl; + } + + private: + NEXT(T x) : x_(x), dx_(static_cast(1) / 0) {} + NEXT(T x, T dx) : x_(x), dx_(dx) {} + + T x_; + T dx_; + }; + + // Calculates the step with Newton's method + template + class StepNewton { + public: + T calc_dx(const CacheOrder1& z) const { return -z.f0() / z.f1(); } + }; + + // Enum for status flags + enum class StatusFlag { normal, dx_small, out_of_iterations, + case_1, case_2U, case_2B, case_2I, + success, failure, stall, f0_is_zero, case_0 + }; + + + +//////////////////////////////////////////////////////////////////////////////////////// + + + // + // + template + class ContextF { + public: + ContextF(F f, T xl, T xh, int digits, std::uintmax_t max_iter) + : f_(f) + , factor_(static_cast(ldexp(1.0, 1 - digits))) + , max_iter_(max_iter) + , count_(max_iter) + , l_orig_(xl) + , h_orig_(xh) + { + max_iter_ = std::min(max_iter_, static_cast(400)); + count_ = std::min(count_, static_cast(400)); + } + + CacheOrder1 operator()(const NEXT& next) { + const T x = next.x(); + + if (count_ <= 0) { + static const char* function = "boost::math::tools::detail::ContextF<%1%>"; + policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::ContextF, last closest guess was %1%", x, boost::math::policies::policy<>()); + } + + if (DEBUG_PRINT) std::cout << "EVAL IS: " << count_ << ", at: " << x << std::endl; + + --count_; + return CacheOrder1(f_, x); + } + + T dx_p() const { return dx_history_.dx_p(); } + T dx_pp() const { return dx_history_.dx_pp(); } + std::uintmax_t iters_used() const { return max_iter_ - count_; } + + T factor() const { return factor_; } + std::uintmax_t max_iter() const { return max_iter_; } + std::uintmax_t count() const { return count_; } + StepSizeHistory& dx_history() { return dx_history_; } + + private: + F f_; + T factor_; + std::uintmax_t max_iter_; + std::uintmax_t count_; + StepSizeHistory dx_history_; + T l_orig_; + T h_orig_; + }; + + + +//////////////////////////////////////////////////////////////////////////////////////// + + // + // + template + class BoundState { + public: + BoundState(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) + : cache_l_(xl) + , cache_h_(xh) + , flag_(StatusFlag::case_1) + , context_(f, xl, xh, digits, max_iter) + { + const auto next = NEXT::FromX(x); + const auto cache_x_temp = context_(next); + update_cache_sign_dx(cache_x_temp); + } + + T calc_dx_true(const CacheOrder1& cx) const { + return cx.x() - cache_x().x(); + } + + // CacheOrder1 operator()(T x) { return context_(x); } + CacheOrder1 operator()(const NEXT& next) { return context_(next); } + + void print_dx_history() const { + if (!DEBUG_PRINT) return; + // std::cout << "dx_p: " << dx_p() << std::endl; + // std::cout << "dx_pp: " << dx_pp() << std::endl; + } + + void update_cache_l(const CacheOrder1& cx) { + if (DEBUG_PRINT) std::cout << "set dx_true to: " << calc_dx_true(cx) << std::endl; + if (DEBUG_PRINT) std::cout << "dx_newton: " << cx.calc_dx() << std::endl; + dx_history().push(calc_dx_true(cx)); + print_dx_history(); + + if (DEBUG_PRINT) std::cout << "UPDATE_CACHE_L to: " << cx.x() << std::endl; + + cache_l_ = cx; + is_last_eval_l_ = true; + } + void update_cache_h(const CacheOrder1& cx) { + if (DEBUG_PRINT) std::cout << "set dx_true to: " << calc_dx_true(cx) << std::endl; + if (DEBUG_PRINT) std::cout << "dx_newton: " << cx.calc_dx() << std::endl; + dx_history().push(calc_dx_true(cx)); + print_dx_history(); + + if (DEBUG_PRINT) std::cout << "UPDATE_CACHE_H to: " << cx.x() << std::endl; + + cache_h_ = cx; + is_last_eval_l_ = false; + } + + void update_cache_sign_f0(const CacheOrder1& cx) { + if (cx.is_same_sign(cache_l_)) { + update_cache_l(cx); + } else { + update_cache_h(cx); + } + } + + void update_cache_sign_dx(const CacheOrder1& cx) { + if (cx.is_l()) { + update_cache_l(cx); + } else { + update_cache_h(cx); + } + } + + bool is_converge_too_slow(const T& dx) const { + using std::fabs; + return fabs(dx_pp()) < fabs(dx * 4); + } + + T midpoint() const { return (cache_l_.x() + cache_h_.x()) / 2; } + bool is_inbounds(T x) const { return (l_smart() <= x) && (x <= h_smart()); } + + T l_smart() const { return cache_l_.x(); } + T h_smart() const { return cache_h_.x(); } + + bool is_dx_small(const NEXT& next) const { + using std::fabs; + + const T side_1 = fabs(next.dx()); + const T side_2 = fabs(next.x()) * factor(); + + // if (DEBUG_PRINT) std::cout << "side_1: " << side_1 << std::endl; + // if (DEBUG_PRINT) std::cout << "side_2: " << side_2 << std::endl; + + return side_1 <= side_2; + + // return fabs(dx_p()) <= fabs(cache_x().x()) * factor(); + } + + void set_flag_if_dx_small(const NEXT& next) { + if (cache_x().f0() == 0) { + set_flag(StatusFlag::f0_is_zero); + } else if (is_dx_small(next)) { + set_flag(StatusFlag::dx_small); + } + } + void set_flag_if_count_zero() { + if (count() <= 0) { + set_flag(StatusFlag::out_of_iterations); + if (DEBUG_PRINT) std::cout << "OUT OF ITERATIONS: count: " << count() << std::endl; + } + } + void set_flag_if_stall(const NEXT& next) { + if (next.dx() == 0) { set_flag(StatusFlag::stall); } + } + + void materialize_l() { + if (is_last_eval_l_) { return; } + cache_l_ = (*this)(l_smart()); + is_last_eval_l_ = true; + } + void materialize_h() { + if (!is_last_eval_l_) { return; } + cache_h_ = (*this)(h_smart()); + is_last_eval_l_ = false; + } + + const CacheOrder1& cache_x() const { + if (is_last_eval_l_) { return cache_l_; } else { return cache_h_; } + } + + int num_x_only() const { return cache_l_.is_x_only() + cache_h_.is_x_only(); } + + T get_limit() const { + assert(num_x_only() == 1); + return is_last_eval_l_ ? cache_h_.x() : cache_l_.x(); + } + + void set_state_2() { + assert(num_x_only() == 0); + + const auto& cache_l = cache_l_; + const auto& cache_h = cache_h_; + + if (cache_l.is_B(cache_h)) { + flag_ = StatusFlag::case_2B; + } else if (cache_l.is_U(cache_h)) { + flag_ = StatusFlag::case_2U; + } else { + assert(cache_l.is_I(cache_h)); + flag_ = StatusFlag::case_2I; + } + } + + void print() const { + if (!DEBUG_PRINT) return; + std::cout << "is_last_eval_l_: " << is_last_eval_l_ << std::endl; + std::cout << "l: "; cache_l_.print(); + std::cout << "h: "; cache_h_.print(); + std::cout << "flag: "; + print_root_status_flag(); + } + + void print_root_status_flag() const { + if (!DEBUG_PRINT) return; + switch (flag()) { + case StatusFlag::normal: std::cout << "normal" << std::endl; break; + case StatusFlag::dx_small: std::cout << "dx_small" << std::endl; break; + case StatusFlag::out_of_iterations: std::cout << "out_of_iterations" << std::endl; break; + case StatusFlag::case_1: std::cout << "case_1" << std::endl; break; + case StatusFlag::case_2U: std::cout << "case_2U" << std::endl; break; + case StatusFlag::case_2B: std::cout << "case_2B" << std::endl; break; + case StatusFlag::case_2I: std::cout << "case_2I" << std::endl; break; + case StatusFlag::success: std::cout << "success" << std::endl; break; + case StatusFlag::failure: std::cout << "failure" << std::endl; break; + case StatusFlag::stall: std::cout << "stall" << std::endl; break; + case StatusFlag::f0_is_zero: std::cout << "f0_is_zero" << std::endl; break; + case StatusFlag::case_0: std::cout << "case_0" << std::endl; break; + } + } + + void set_flag(StatusFlag flag) { flag_ = flag; } + + bool is_B() const { return cache_l_.is_B(cache_h_); } + bool is_U() const { return cache_l_.is_U(cache_h_); } + bool is_I() const { return cache_l_.is_I(cache_h_); } + + StepSizeHistory& dx_history() { return context_.dx_history(); } + std::uintmax_t count() const { return context_.count(); } + T factor() const { return context_.factor(); } + std::uintmax_t max_iter() const { return context_.max_iter(); } + T dx_p() const { return context_.dx_p(); } + T dx_pp() const { return context_.dx_pp(); } + std::uintmax_t iters_used() const { return context_.iters_used(); } + + NEXT create_next_for_x(T x) { + const T dx = x - cache_x().x(); + return NEXT::From_x_dx(x, dx); + } + + const CacheOrder1& cache_l() const { return cache_l_; } + const CacheOrder1& cache_h() const { return cache_h_; } + StatusFlag flag() const { return flag_; } + + void set_x_next_value(T x) { x_next_value_ = x; } + T x_next_value() const { return x_next_value_; } + + private: + bool is_last_eval_l_; + CacheOrder1 cache_l_; + CacheOrder1 cache_h_; + StatusFlag flag_; + ContextF context_; + T x_next_value_; + }; + +//////////////////////////////////////////////////////////////////////////////////////// + + template + class Bound2B; + + template + class Bound2U; + + template + class Bound1; + +//////////////////////////////////////////////////////////////////////////////////////// + + // + // + template + class BoundFn { + public: + void do_debug_printout(BoundState& b) { + b.print(); + if (DEBUG_PRINT) std::cout << std::endl; + } + + void solve_next(BoundState& b) { + do_debug_printout(b); + + while (is_bound_still_valid(b)) { + const auto next = calc_next(b); + if (DEBUG_PRINT) next.print(); + + b.set_x_next_value(next.x()); + + if (is_stop(b, next)) { + if (DEBUG_PRINT) std::cout << "GOT STOP SIGNAL" << std::endl; + break; + } + + const auto cx = b(next); + this->update_cache(b, cx); + do_debug_printout(b); + } + if (DEBUG_PRINT) std::cout << "DONE WITH LOOP" << std::endl; + } + + NEXT calc_next(BoundState& b) { + // Attempt Newton step + const auto& cx = b.cache_x(); + const auto next = NEXT::TakeStep(cx); + + const bool c1 = !(b.is_inbounds(next.x())); + const bool c2 = b.is_converge_too_slow(next.dx()); + + if (c1) { + if (DEBUG_PRINT) std::cout << "NOT INBOUNDS x_next: " << next.x() << std::endl; + } + if (c2) { + if (DEBUG_PRINT) std::cout << "CONVERGE TOO SLOW dx_next: " << next.dx() << std::endl; + } + + // if (!(b.is_inbounds(next.x())) || b.is_converge_too_slow(next.dx())) { + if (c1 || c2) { + b.dx_history().reset(); + if (DEBUG_PRINT) std::cout << "TOOK BISECTION wanted: " << next.x() << ", dx: " << next.dx() << std::endl; + return this->calc_next_bisection(b); + } else { + if (DEBUG_PRINT) std::cout << "TOOK NEWTON -------------" << std::endl; + return next; + } + } + + NEXT calc_next_bisection_regular(BoundState& b) { + const auto x_mid = b.midpoint(); + if (x_mid == b.l_smart() || x_mid == b.h_smart()) { + b.set_flag(StatusFlag::dx_small); + } + // const auto next = NEXT::FromX(x_mid); + const auto next = b.create_next_for_x(x_mid); + return next; + } + + virtual bool is_stop(BoundState& b, const NEXT next) = 0; + virtual NEXT calc_next_bisection(BoundState& b) = 0; + virtual void update_cache(BoundState& b, const CacheOrder1& cx) = 0; + virtual void print_class() const = 0; + // virtual void is_bound_still_valid() const = 0; + virtual bool is_bound_still_valid(const BoundState& b) const = 0; + }; + + // + // + template + class Bound1 : public BoundFn { + public: + bool is_stop(BoundState& b, const NEXT next) override { + // if (b.num_x_only() != 2) { + b.set_flag_if_dx_small(next); + b.set_flag_if_count_zero(); + b.set_flag_if_stall(next); + // } + return b.flag() != StatusFlag::case_1; + } + + NEXT calc_next_bisection(BoundState& b) override { + const T x = b.cache_x().x(); + const T x_limit = b.get_limit(); + const T dx_limit = x_limit - x; + const T dx_8 = b.cache_x().calc_dx() * 1000; + + using std::fabs; + if (fabs(dx_limit) < fabs(dx_8)) { + return b.create_next_for_x(x_limit); + // b.set_x_dx_next_for_x(x_limit); + // b.materialize_l(); + // b.materialize_h(); + // b.set_state_2(); + } else { + // update_cache(b, b(x + dx_8)); + // b.set_x_dx_next_for_x(x + dx_8); + return b.create_next_for_x(x + dx_8); + } + } + + void update_cache(BoundState& b, const CacheOrder1& cx) override { + b.update_cache_sign_dx(cx); + if (b.num_x_only() == 0) { + b.set_state_2(); + } else { + if (b.get_limit() == cx.x()) { + b.set_flag(StatusFlag::case_2I); + } + } + } + + void print_class() const override { + if (!DEBUG_PRINT) return; + std::cout << "Bound1" << std::endl; } + + bool is_bound_still_valid(const BoundState& b) const override { return b.flag() == StatusFlag::case_1; } + }; + + // + // + template + class Bound2B : public BoundFn { + public: + bool is_stop(BoundState& b, const NEXT next) override { + b.set_flag_if_dx_small(next); + b.set_flag_if_count_zero(); + b.set_flag_if_stall(next); + return b.flag() != StatusFlag::case_2B; + } + + NEXT calc_next_bisection(BoundState& b) override { + return this->calc_next_bisection_regular(b); + } + + void update_cache(BoundState& b, const CacheOrder1& cx) override { + b.update_cache_sign_f0(cx); + } + + void print_class() const override { + if (!DEBUG_PRINT) return; + std::cout << "Bound2B" << std::endl; } + + bool is_bound_still_valid(const BoundState& b) const override { return b.flag() == StatusFlag::case_2B; } + }; + + // + // + template + class Bound2U : public BoundFn { + public: + bool is_stop(BoundState& b, const NEXT next) override { + b.set_flag_if_dx_small(next); + b.set_flag_if_count_zero(); + return b.flag() != StatusFlag::case_2U; + } + + NEXT calc_next_bisection(BoundState& b) override { + return this->calc_next_bisection_regular(b); + } + + void update_cache(BoundState& b, const CacheOrder1& cx) override { + if (is_dx_bigger_than_prev(b, cx)) { + b.set_flag(StatusFlag::case_2I); + } + b.update_cache_sign_dx(cx); + } + + bool is_dx_bigger_than_prev(BoundState& b, const CacheOrder1& cx) const { + using std::fabs; + const auto& cache_same_side = cx.is_l() ? b.cache_l() : b.cache_h(); + const T dx_cx = cx.calc_dx(); + const T dx_same_side = cache_same_side.calc_dx(); + const bool is_dx_bigger = fabs(dx_same_side) < fabs(dx_cx); + return is_dx_bigger; + } + + void print_class() const override { + if (!DEBUG_PRINT) return; + std::cout << "Bound2U" << std::endl; } + + bool is_bound_still_valid(const BoundState& b) const override { return b.flag() == StatusFlag::case_2U; } + }; + + +//////////////// +//////////////// +//////////////// + + + // template + // std::pair solve(BoundState& state) { + + // state.is_finished(); + // } + + template + std::pair solve(BoundState& state) { + + if (DEBUG_PRINT) { + std::cout << std::endl; + std::cout << "=========================================" << std::endl; + std::cout << "Count: " << state.count() << std::endl; + std::cout << "STATUS IS: "; + state.print_root_status_flag(); + std::cout << "-----------------------------------------" << std::endl; + std::cout << std::endl; + } + + if (state.flag() == StatusFlag::case_1) { + auto solver_case_1 = Bound1(); + solver_case_1.solve_next(state); + return solve(state); + + } else if (state.flag() == StatusFlag::case_2B) { + auto solver_case_2b = Bound2B(); + solver_case_2b.solve_next(state); + return solve(state); + + } else if (state.flag() == StatusFlag::case_2U) { + auto solver_case_2u = Bound2U(); + solver_case_2u.solve_next(state); + return solve(state); + + } else if (state.flag() == StatusFlag::case_2I) { + return std::make_pair(false, state.x_next_value()); + + } else if (state.flag() == StatusFlag::stall) { + return std::make_pair(true, state.x_next_value()); + + } else if (state.flag() == StatusFlag::dx_small) { + return std::make_pair(true, state.x_next_value()); + + } else if (state.flag() == StatusFlag::f0_is_zero) { + return std::make_pair(true, state.x_next_value()); + + } else if (state.flag() == StatusFlag::out_of_iterations) { + if (DEBUG_PRINT) std::cout << "ERROR OUT OF ITERATIONS" << std::endl; + assert(false); + + } else { + + // std::cout << "NOT HANDLED: " << static_cast(state.flag()) << std::endl; + if (DEBUG_PRINT) std::cout << "NOT HANDLED: " << std::endl; + state.print_root_status_flag(); + assert(false); + } + + assert(false); + return std::make_pair(false, state.x_next_value()); + } + + + +///////////////////////// + + template + T solve_2(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { + // Create BoundState + BoundState state(f, x, xl, xh, digits, max_iter); + + + } + +///////////////////////// + + + + template + T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { + if (DEBUG_PRINT) std::cout << "================= newton_raphson_iterate ==================== with max_iter: " << max_iter << std::endl; + + if (DEBUG_PRINT) { + std::cout << "l_orig: " << xl << std::endl; + std::cout << "x_orig: " << x << std::endl; + std::cout << "h_orig: " << xh << std::endl; + } + + + + // Create BoundState +#if 1 + BoundState state(f, x, xl, xh, digits, max_iter); + +#else + BoundState state(f, xl, xh, digits, max_iter); + + // Calc NEXT + const auto next = NEXT(x); + + // Calc CacheX + const auto cache_x_orig = state(next); + + + + + // // Set up problem + // state.update_cache_sign_dx(cache_x_orig); + + // // auto solver_case_1 = Bound1(); + // Bound1().solve(state); +#endif + + + + + + + + + + + +////////////////////////////////////////////////////////////////////////// + + + // Get CacheX + const auto cache_x_orig = state.cache_x(); + + // Solve problem + auto p = solve(state); + + +////////////////////////////////////////////////////////////////////////// + // Define lambdas + auto fn_form_bracket_l = [&](){ + const auto next = NEXT::FromX(xl); + const auto cache_l_orig = state(next); + state.update_cache_l(cache_l_orig); + state.update_cache_h(cache_x_orig); + }; + auto fn_form_bracket_h = [&](){ + const auto next = NEXT::FromX(xh); + const auto cache_h_orig = state(next); + state.update_cache_l(cache_x_orig); + state.update_cache_h(cache_h_orig); + }; + auto fn_return = [&](const BoundState& s, const std::pair p) { + s.print_root_status_flag(); + assert(p.first); + max_iter = s.iters_used(); + return p.second; + }; +////////////////////////////////////////////////////////////////////////// + + // Do the natural + if (p.first) { + if (DEBUG_PRINT) std::cout << "succuss first time" << std::endl; + max_iter = state.iters_used(); + return fn_return(state, p); + } + if (DEBUG_PRINT) std::cout << "FAILURE first time" << std::endl << std::endl; + + + cache_x_orig.is_l() ? fn_form_bracket_h() : fn_form_bracket_l(); + state.set_state_2(); + if (state.is_B()) { + if (DEBUG_PRINT) std::cout << "BRACKETED PRIORITY 1" << std::endl; + p = solve(state); + return fn_return(state, p); + } else { + if (DEBUG_PRINT) std::cout << "COULD NOT BRACKET PRIORITY 1" << std::endl << std::endl; + } + + cache_x_orig.is_l() ? fn_form_bracket_l() : fn_form_bracket_h(); + state.set_state_2(); + if (state.is_B()) { + if (DEBUG_PRINT) std::cout << "BRACKETED PRIORITY 2" << std::endl; + p = solve(state); + return fn_return(state, p); + } else { + if (DEBUG_PRINT) std::cout << "COULD NOT BRACKET PRIORITY 2" << std::endl << std::endl; + } + + static const char* function = "boost::math::tools::detail::solve<%3%>"; + const T x_last = state.x_next_value(); + return policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::solve, last closest guess was %1%", x_last, boost::math::policies::policy<>()); + }; + +// from compile_test/gauss_kronrod_concept_test.cpp:11: +// ../../../libs/math/test/test_ibeta_inv.hpp(64): last checkpoint +// ../../../bin.v2/libs/math/test/test_bessel_airy_zeros.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/test_bessel_airy_zeros.o... +// test_inverse_gaussian +// cat "../../../bin.v2/libs/math/test/test_ibeta_inv_double.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/test_ibeta_inv_double.output" + +// test_hypergeometric_laplace_transform + + + +// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/test_vector_barycentric_rational.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/test_vector_barycentric_rational.o... +// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/cstdfloat_concept_check_1.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/compile_test/cstdfloat_concept_check_1.o... +// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/std_real_concept_check.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/compile_test/std_real_concept_check.o... +// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/std_real_concept_check_80.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/std_real_concept_check_80.o... +// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/std_real_concept_check_32.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/std_real_concept_check_32.o... +// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/std_real_concept_check_128.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/std_real_concept_check_128.o... +// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/cstdfloat_concept_check_2.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/compile_test/cstdfloat_concept_check_2.o... +// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/cstdfloat_concept_check_3.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/compile_test/cstdfloat_concept_check_3.o... +// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/cstdfloat_concept_check_4.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/compile_test/cstdfloat_concept_check_4.o... +// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/std_real_concept_check_64.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/std_real_concept_check_64.o... + + + +// ../../../b2 --report_level=detailed --show_progress=false | grep -E "test case|failed" + + + + +// ### ??? +// +// ### CONFIRMED TO BE REAL TESTS ### +// test_legendre +// test_ibeta_inv_double +// test_laguerre +// test_bessel_airy_zeros +// test_hermite +// +// testing.capture-output ../../../bin.v2/libs/math/example/legendre_stieltjes_example.test/gcc-9/debug/threading-multi/visibility-hidden/legendre_stieltjes_example.run + + + + + + +// Start here +#if 0 +//////////////////////////////////////////////////////////////////////////////////////// + + // Forward declare ContextF + template + class ContextF; + + + // + // + // + // + // + template + class BoundsState { + public: + virtual ~BoundsState() = default; + + void set_b(BoundState& b) { b_ = &b; } + + void do_step_eval_then_x_next(ContextF& helper_f) { + // Attempt Newton step + const auto cx = cache_x(); + const T dx_newton = cx.calc_dx(); + const T x_next = cx.x() + dx_newton; + + // Need to do bisection + if (!is_inbounds(x_next) || is_converge_too_slow(dx_newton)) { + calc_next_bisection(); + return; + } + + // Newton's step + update_cache(helper_f(x_next)); + return; + } + + bool is_converge_too_slow(const T& dx) const { + return fabs(b_->dx_pp()) < fabs(dx * 4); + } + + virtual const CacheOrder1& cache_x() const = 0; + virtual NEXT calc_next_bisection() = 0; + virtual void update_cache(const CacheOrder1& cache_x) = 0; + + bool is_inbounds(T x) const { return (l_smart() <= x) && (x <= h_smart()); } + virtual T l_smart() const = 0; + virtual T h_smart() const = 0; + + + protected: + ContextF* b_; + }; + + + + // + template + class ContextF { + public: + ContextF(F f, int digits, std::uintmax_t max_iter) + : f_(f) + , factor_(static_cast(ldexp(1.0, 1 - digits))) + , count_(max_iter) + , max_iter_(max_iter) {} + + ~ContextF() { + delete bounds_; + } + + CacheOrder1 operator()(T x) { + if (0 == count_) { + static const char* function = "boost::math::tools::detail::ContextF<%1%>"; + policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::ContextF, last closest guess was %1%", x, boost::math::policies::policy<>()); + } + --count_; + return CacheOrder1(f_, x); + } + + void transition_to_state(BoundsState* bounds) { + if (this->bounds_ != nullptr) + delete this->bounds_; + this->bounds_ = bounds; + this->bounds_->set_b(*this); + } + + SolveOutput solve(BoundsState* bounds) { + transition_to_state(bounds); + return solve(); + } + + T dx_pp() const { return dx_history_.dx_pp(); } + int iters_used() const { return max_iter_ - count_; } + + private: + SolveOutput solve() { + BOOST_MATH_STD_USING + std::cout << std::scientific << std::setprecision(6) << std::endl; + + while (count_) { + bounds_->do_step_eval_then_x_next(*this); + } + assert(false); + return std::make_pair(OutputResult::failure, 0); + } + +#if 0 + std::pair solve(Bounds& bounds, CacheOrder1 cache_x) { + BOOST_MATH_STD_USING + + std::cout << std::scientific << std::setprecision(6) << std::endl; + + T dx_true; + bool is_prev_newton = false; + + T x_next; + while (count_) { + // Calculate Newton step + const T dx_newton = cache_x.calc_dx(); + x_next = cache_x.x() + dx_newton; // Checked against bracket bounds below + + + const bool c1 = !bounds.is_inbounds(x_next); + if (c1 || fabs(dx_pp()) < fabs(dx_newton * 4)) { // Bisection halves step size each iteration + const auto p_bisect = bounds.handle_bisection(*this, cache_x); + if (bounds.is_case_2I()) { std::cout << "is_case_2I" << std::endl; return std::make_pair(false, 0); } + + const auto cache_x_next = p_bisect.second; + std::cout << "STEP: BISECTION: " << " -- wanted: " << x_next << " -- got: " << cache_x_next.x() << " -- bounds: " << c1 << " -- converge: " << (fabs(dx_pp()) < fabs(dx_newton * 4)) << std::endl; + x_next = cache_x_next.x(); + dx_history_.reset(); // Invalidate step size history + + dx_true = x_next - cache_x.x(); + + cache_x = cache_x_next; + is_prev_newton = false; + } else { + std::cout << "STEP: NEWTON" << std::endl; + dx_true = x_next - cache_x.x(); + cache_x = (*this)(x_next); // Update cache_x_ with x_next + is_prev_newton = true; + } + + std::cout << "AFTER STEP-- x_next: " << x_next << ", dx_true: " << dx_true << std::endl; + bounds.print(); + + assert(bounds.is_inbounds(x_next)); + + // std::cout << " UPDATED cache_x: " << std::endl; + // std::cout << " "; + // cache_x.print(); + // std::cout << " count: " << count_ << std::endl; + + // Determine if successful + const bool is_dx_small = fabs(cache_x.calc_dx()) < fabs(cache_x.x() * factor_); + if (is_dx_small) { + if (bounds.is_case_2B() || bounds.is_case_1()) { + std::cout << "is_case_2B && is_dx_small" << std::endl; + std::cout << "dx_true: " << fabs(dx_true) << ", cache_x.x(): " << fabs(cache_x.x()) << ", factor_: " << factor_ << std::endl; + break; // Tolerance met + } else { + bounds.print_case(); + std::cout << "dx_true: " << fabs(dx_true) << ", cache_x.x(): " << fabs(cache_x.x()) << ", factor_: " << factor_ << std::endl; + std::cout << "FAILURE dx_small" << std::endl; + return std::make_pair(false, 0); // Did not bracket root + } + } + + if (cache_x.f0() == 0) { + std::cout << "SUCCESS f0 = 0" << std::endl; + break; // Function is zero + } + + const bool is_numerical = (x_next == bounds.l_smart() || x_next == bounds.h_smart() || dx_true == 0); + if (is_numerical) { + if (bounds.is_case_2B()) { + std::cout << "is_case_2B && is_numerical" << std::endl; + break; // Tolerance met + } else { + std::cout << "x_next == bounds.l_smart(): " << (x_next == bounds.l_smart()) << std::endl; + std::cout << "x_next == bounds.h_smart(): " << (x_next == bounds.h_smart()) << std::endl; + std::cout << "dx_true == 0: " << (dx_true == 0) << std::endl; + std::cout << "FAILURE numerical" << std::endl; + return std::make_pair(false, 0); // Did not bracket root + } + } + + + if (bounds.is_case_2I()) { + std::cout << "FAILURE is_case_2I" << std::endl; + return std::make_pair(false, 0); + } + + bounds.update_state(cache_x); // Update bracket state + + if (bounds.is_case_2I()) { + std::cout << "FAILURE is_case_2I" << std::endl; + return std::make_pair(false, 0); + } + + std::cout << std::endl; + dx_history_.push(dx_newton); // Store step size + } + + return std::make_pair(true, cache_x.x()); + } +#endif + + private: + F f_; + T factor_; + std::uintmax_t count_; + std::uintmax_t max_iter_; + StepSizeHistory dx_history_; + BoundsState* bounds_; + }; + + + + template + class Bounds2; // : public BoundsState; + + template + class Bounds2B; // : public BoundsState; + + template + class Bounds2U; // : public BoundsState; + + template + class Bounds2I; // : public BoundsState; + + // + // + // + // + // + template + class Bounds2 : public BoundsState { + public: + static Bounds2* Identify(CacheOrder1 cache_a, CacheOrder1 cache_b) { + const auto cache_sorted = cache_a.sort(cache_b); + const auto cache_l = cache_sorted.first; + const auto cache_h = cache_sorted.second; + + if (cache_l.is_B(cache_h)) { + return new Bounds2B(cache_l, cache_h); + } else if (cache_l.is_U(cache_h)) { + return new Bounds2U(cache_l, cache_h); + } else { + return new Bounds2I(cache_l, cache_h); + } + } + + NEXT calc_next_bisection() override { + const auto midpoint = (cache_l_.x() + cache_h_.x()) / 2; + const auto cache_midpoint = helper_f(midpoint); + update_cache(cache_midpoint); + } + + void update_cache_sign_f0(const CacheOrder1& cache_x) { + if (cache_x.is_same_sign(cache_l_)) { + cache_l_ = cache_x; + } else { + cache_h_ = cache_x; + } + } + + void update_cache_sign_dx(const CacheOrder1& cache_x) { + if (cache_x.is_same_slope(cache_l_)) { // TODO: make better + cache_l_ = cache_x; + } else { + cache_h_ = cache_x; + } + } + + T l_smart() const override { return cache_l_.x(); } + T h_smart() const override { return cache_h_.x(); } + + const detail::CacheOrder1& cache_x() const override { + return (is_last_eval_l_) ? cache_l_ : cache_h_; + } + + bool is_B() const { return cache_l_.is_B(cache_h_); } + bool is_U() const { return cache_l_.is_U(cache_h_); } + bool is_I() const { return cache_l_.is_I(cache_h_); } + + const CacheOrder1& cache_l() { return cache_l_; } + const CacheOrder1& cache_h() { return cache_h_; } + + private: + bool is_last_eval_l_; + CacheOrder1 cache_l_; + CacheOrder1 cache_h_; + }; + + // + // + // + // + // + template + class Bounds2B : public Bounds2 { + public: + Bounds2B(const CacheOrder1& cache_l, const CacheOrder1& cache_h) + : Bounds2(cache_l, cache_h) + { + assert(this->is_B()); + } + + void update_cache(const CacheOrder1& cache_x) override { + update_cache_sign_f0(cache_x); + } + }; + + // + // + // + // + // + template + class Bounds2U : public Bounds2 { + public: + Bounds2U(const CacheOrder1& cache_l, const CacheOrder1& cache_h) + : Bounds2(cache_l, cache_h) + { + assert(this->is_U()); + } + + void update_cache(const CacheOrder1& cache_x) override { + update_cache_sign_dx(cache_x); + if (this->is_B()) { + this->b_->transition_to_state( + new Bounds2B(this->cache_l_, this->cache_h_) + ); + } + } + }; + + + // + // + // + // + // + template + class Bounds2I : public Bounds2 { + public: + Bounds2I(const CacheOrder1& cache_l, const CacheOrder1& cache_h) + : Bounds2(cache_l, cache_h) + { + assert(this->is_U()); + } + + void update_cache(const CacheOrder1& cache_x) override { + assert(false); + } + }; + + + // + // + // + // + // + template + class Bounds1 : public BoundsState { + public: + Bounds1(T l, const detail::CacheOrder1& cache_x, T h) + : cache_x_(cache_x) + { + limit_ = cache_x_.is_l() ? h : l; + } + + NEXT calc_next_bisection() override { + const auto cache_limit = (*this->b_)(limit_); + auto* bounds_new = Bounds2::Identify(cache_limit, cache_x_); + this->b_->transition_to_state(bounds_new); + } + + T l_smart() const override { return cache_x_.is_l() ? cache_x_.x() : limit_; } + T h_smart() const override { return cache_x_.is_h() ? cache_x_.x() : limit_; } + + void update_cache(const CacheOrder1& cache_x) override { cache_x_ = cache_x; } + + const CacheOrder1& cache_x() const override { return cache_x_; } + + private: + CacheOrder1 cache_x_; + T limit_; + }; + + + +/////////////////////////////////////////////////////////////////////// + + + template + T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { + // Create ContextF + ContextF b(f, digits, max_iter); + + // Create cache_x + const auto cache_x = b(x); + if (0 == cache_x.f0()) return x; // Check if root is already found + + // Create Bounds + auto bounds = Bounds1(xl, cache_x, xh); + + // Solve + const auto p = b.solve(&bounds); + + return p.second; + +#if 0 + if (p.first) { + std::cout << "succuss first time" << std::endl; + max_iter = b.iters_used(); + return p.second; + } else { + std::cout << "FAILURE first time" << std::endl << std::endl; + + bounds = Bounds(xl, cache_x, xh); + cache_x.is_l() ? bounds.eval_h_orig(b) : bounds.eval_l_orig(b); + if (bounds.is_case_2B()) { + std::cout << "BRACKETED PRIORITY 1" << std::endl; + return b.solve(bounds, cache_x).second; + } else { + std::cout << "COULD NOT BRACKET PRIORITY 1" << std::endl << std::endl; + } + + bounds = Bounds(xl, cache_x, xh); + cache_x.is_l() ? bounds.eval_l_orig(b) : bounds.eval_h_orig(b); + if (bounds.is_case_2B()) { + std::cout << "BRACKETED PRIORITY 2" << std::endl; + return b.solve(bounds, cache_x).second; + } else { + std::cout << "l: " << bounds.l_smart() << ", h: " << bounds.h_smart() << std::endl; + std::cout << "f0_l: " << bounds.cache_l().f0() << ", f0_h: " << bounds.cache_h().f0() << std::endl; + std::cout << "COULD NOT BRACKET PRIORITY 2" << std::endl << std::endl; + } + } +#endif + + static const char* function = "boost::math::tools::detail::ContextF<%1%>"; + policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::ContextF, last closest guess was %1%", x, boost::math::policies::policy<>()); + return 0; + }; + + + + + + +#if 0 + + + // + // + // + // + // + template + class Bounds2 { + public: + Bounds2(T l, const detail::CacheOrder1& cache_x, T h) + : var_l_(LimitOrig(l)), var_h_(LimitOrig(h)) + { + if (cache_x.is_l()) { + is_last_eval_l_ = true; + var_l_ = cache_x; + } else { + is_last_eval_l_ = false; + var_h_ = cache_x; + } + } + + template + void do_step_eval_then_x_next(H& helper_f) { + const auto cx = cache_x(); + const T x_next = cx.x() + cx.calc_dx(); + + if (!is_inbounds(x_next) || is_converge_too_slow()) { + return calc_next_bisection(helper_f); + } + + // Newton's step + update_cache(helper_f(x_next)); + } + + template + NEXT calc_next_bisection(H& helper_f) { + if (is_case_1()) { + materialize(helper_f); + return; + } + const T x_mid = (x_l() + x_h()) / 2; + update_cache_2(helper_f(x_mid)); + return; + } + + void update_cache(const detail::CacheOrder1& cache_x) { + if (is_case_1()) { + update_cache_1(cache_x); + } else { + update_cache_2(cache_x); + } + } + + void update_cache_1(const detail::CacheOrder1& cache_x) { + update_cache_using_dx(cache_x); + if (size() == 2) { + determine_state_2(); + } + } + + void update_cache_2(const detail::CacheOrder1& cache_x) { + if (is_case_2B()) { + update_cache_using_sign(cache_x); + } else { + assert(is_case_2U()); + update_cache_2U(cache_x); + } + } + + void update_cache_2U(const detail::CacheOrder1& cache_x) { + const T dx_l_orig = cache_l().calc_dx(); + const T dx_h_orig = cache_h().calc_dx(); + + update_cache_using_dx(cache_x); + + if (cache_l().is_B(cache_h())) { // Now bracket root! + state_ = BracketState::case_2B; + } else { + const T dx_l = cache_l().calc_dx(); + const T dx_h = cache_h().calc_dx(); + + if (dx_l_orig < dx_l || dx_h_orig < dx_h) { + state_ = BracketState::case_2I; + } + } + } + + template + void materialize(H& helper_f) { + // If var_l_ holds a LimitOrig, then we need to evaluate it + if (var_l_.is_limit_orig()) { + set_l(helper_f(var_l_.x())); + } + // If var_h_ holds a LimitOrig, then we need to evaluate it + if (var_h_.is_limit_orig()) { + set_h(helper_f(var_h_.x())); + } + determine_state_2(); + } + + void determine_state_2() { + assert(size() == 2); + if (cache_l().is_B(cache_h())) { + state_ = BracketState::case_2B; + } else if (cache_l().is_U(cache_h())) { + state_ = BracketState::case_2U; + } else { + state_ = BracketState::case_2I; + } + } + + void update_state_using_dx(const CacheOrder1& cache_x) { + if (cache_x.is_l()) { + var_l_ = cache_x; + } else { + var_h_ = cache_x; + } + } + + void update_state_using_sign(const CacheOrder1& cache_x) { + if (has_l()) { + if (cache_x.is_same_sign(cache_l())) { + var_l_ = cache_x; + } else { + var_h_ = cache_x; + } + } else { + if (cache_x.is_same_sign(cache_h())) { + var_h_ = cache_x; + } else { + var_l_ = cache_x; + } + } + } + +//////////////////////////////////////////////////////////////////////////////// + + detail::CacheOrder1 cache_x() const { return is_last_eval_l_ ? val_l_ : var_h_; } + + int size() const { return var_l_.is_cache() + var_h_.is_cache(); } + + bool is_inbounds(T x) const { return (l() <= x) && (x <= h()); } + + void set_l(const detail::CacheOrder1& cache_x) { + is_last_eval_l_ = true; + var_l_ = cache_x; + } + void set_h(const detail::CacheOrder1& cache_x) { + is_last_eval_l_ = false; + var_h_ = cache_x; + } + + T l() const { return var_l_.x(); } + T h() const { return var_h_.x(); } + + private: + BracketState state_; + bool is_last_eval_l_; + booth::variant, detail::CacheOrder1> var_l_; + booth::variant, detail::CacheOrder1> var_h_; + StepSizeHistory dx_history_; + }; + + + + + + + + + + + + + + + + + + + // + template + class Bounds { + public: + Bounds(T l, detail::CacheOrder1 cache_x, T h) + : state_(BracketState::case_1) + , l_orig_(l) + , h_orig_(h) + { + update_state_using_dx(cache_x); + } + + void update_state_using_dx(const CacheOrder1& cache_x) { + if (cache_x.is_l()) { + opt_cache_l_ = cache_x; + } else { + opt_cache_h_ = cache_x; + } + } + + void update_state_using_sign(const CacheOrder1& cache_x) { + if (has_l()) { + if (cache_x.is_same_sign(cache_l())) { + opt_cache_l_ = cache_x; + } else { + opt_cache_h_ = cache_x; + } + } else { + if (cache_x.is_same_sign(cache_h())) { + opt_cache_h_ = cache_x; + } else { + opt_cache_l_ = cache_x; + } + } + } + + void determine_state_exhaustive() { + if (cache_l().is_B(cache_h())) { + state_ = BracketState::case_2B; + } else if (cache_l().is_U(cache_h())) { + state_ = BracketState::case_2U; + } else { + state_ = BracketState::case_2I; + } + } + + bool update_state(const CacheOrder1& cache_x) { + if (is_case_1()) { + const bool c1 = cache_x.is_l() && has_h(); + const bool c2 = cache_x.is_h() && has_l(); + update_state_using_dx(cache_x); + if (c1 || c2) { + determine_state_exhaustive(); + } + } else if (is_case_2B()) { + update_state_using_sign(cache_x); + if (!cache_l().is_B(cache_h())) { + std::cout << std::endl; + std::cout << "SOMETHING IS WRONG WRONG WRONG" << std::endl; + print(); + cache_x.print(); + assert(false); + } + } else { + assert(is_case_2U()); + update_state_2U(cache_x); + } + return true; + } + + void update_state_2U(const CacheOrder1& cache_x) { + const T dx_l_orig = cache_l().calc_dx(); + const T dx_h_orig = cache_h().calc_dx(); + + update_state_using_dx(cache_x); + + if (cache_l().is_B(cache_h())) { // Now bracket root! + state_ = BracketState::case_2B; + } else { + const T dx_l = cache_l().calc_dx(); + const T dx_h = cache_h().calc_dx(); + + if (dx_l_orig < dx_l || dx_h_orig < dx_h) { + state_ = BracketState::case_2I; + } + } + } + + static T midpoint_smart(T l, T h) { + BOOST_MATH_STD_USING + + const T t_min = std::numeric_limits::min(); + const T t_max = std::numeric_limits::max(); + const T t_lowest = std::numeric_limits::lowest(); + const T t_inf = std::numeric_limits::infinity(); + + const T l_abs = fabs(l); + const T h_abs = fabs(h); + if (std::max(l_abs, h_abs) < 100 || // Numbers are small + (h - l) < std::min(l_abs, h_abs)) { // Numbers are close + return (l + h) / 2; + } + + // Want both numbers to have the same sign + if (sign(l) * sign(h) == -1) { return 0; } + + // If l is -Inf make it the lowest float + if (l == -t_inf) { + l = t_lowest; + } else if (h == t_inf) { + h = t_max; + } + + // If l is 0 make it the smallest float + if (l == 0) { + l = t_min; + } else if (h == 0) { + h = -t_min, l; + } + + return copysign(sqrt(l * h), l); + } + + template + std::pair> handle_bisection(H& helper_f, const CacheOrder1& cache_x) { + if (is_case_1()) { + const T x_next = cache_x.is_l() ? h_orig_ : l_orig_; + + CacheOrder1 cache_x_next; + + if (cache_x.is_l()) { + cache_x_next = helper_f(x_next); + opt_cache_h_ = cache_x_next; + } else { + cache_x_next = helper_f(x_next); + opt_cache_l_ = cache_x_next; + } + + determine_state_exhaustive(); + print_case(); + } + + // const T x_mid = (l_smart() + h_smart()) / 2; + const T x_mid = midpoint_smart(l_smart(), h_smart()); + + std::cout << "x_mid: " << x_mid << std::endl; + + const auto cache_x_next = helper_f(x_mid); + + std::cout << "post eval" << std::endl; + + return std::make_pair(true, cache_x_next); + } + + template + void eval_l_orig(H& helper_f) { + const auto cache_x = helper_f(l_orig_); + push_front(cache_x); + determine_state_exhaustive(); + } + template + void eval_h_orig(H& helper_f) { + const auto cache_x = helper_f(h_orig_); + push_back(cache_x); + determine_state_exhaustive(); + } + + void push_front(const CacheOrder1& cache_x) { + assert(size() == 1); + if (has_l()) { opt_cache_h_ = opt_cache_l_; } + opt_cache_l_ = cache_x; + } + void push_back(const CacheOrder1& cache_x) { + assert(size() == 1); + if (has_h()) { opt_cache_l_ = opt_cache_h_; } + opt_cache_h_ = cache_x; + } + + int size() const { return has_l() + has_h(); } + bool has_l() const { return opt_cache_l_.has_data(); } + bool has_h() const { return opt_cache_h_.has_data(); } + + const CacheOrder1& cache_l() const { return opt_cache_l_.c(); } + const CacheOrder1& cache_h() const { return opt_cache_h_.c(); } + + T l_smart() const { return has_l() ? cache_l().x() : l_orig_; } + T h_smart() const { return has_h() ? cache_h().x() : h_orig_; } + + bool is_inbounds(T x) const { return (l_smart() <= x) && (x <= h_smart()); } + + bool is_case_1() const { return state_ == BracketState::case_1; } + bool is_case_2U() const { return state_ == BracketState::case_2U; } + bool is_case_2B() const { return state_ == BracketState::case_2B; } + bool is_case_2I() const { return state_ == BracketState::case_2I; } + + std::string get_case_string() const { + if (is_case_1()) { return "case_1"; } + if (is_case_2U()) { return "case_2U"; } + if (is_case_2B()) { return "case_2B"; } + assert(is_case_2I()); + return "case_2I"; + } + + void print_case() const { std::cout << get_case_string() << std::endl; } + + void print() const { + std::cout << "BOUNDS: " << get_case_string() << std::endl; + std::cout << " l_orig: " << l_orig_ << std::endl; + if (has_l()) { cache_l().print(); } + if (has_h()) { cache_h().print(); } + std::cout << " h_orig: " << h_orig_ << std::endl; + } + + private: + BracketState state_; + T l_orig_; + detail::Optional> opt_cache_l_; + detail::Optional> opt_cache_h_; + T h_orig_; + }; + + template + T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { + // Create ContextF + ContextF helper_f(f, digits, max_iter); + + // Create cache_x + const auto cache_x = helper_f(x); + if (0 == cache_x.f0()) return x; // Check if root is already found + + // Create Bounds2 + auto bounds = Bounds1(xl, cache_x, xh); + + // Do Newton loop + const auto p = helper_f.solve(bounds, cache_x); + + if (p.first) { + std::cout << "succuss first time" << std::endl; + max_iter = helper_f.iters_used(); + return p.second; + } else { + std::cout << "FAILURE first time" << std::endl << std::endl; + + bounds = Bounds2(xl, cache_x, xh); + cache_x.is_l() ? bounds.eval_h_orig(helper_f) : bounds.eval_l_orig(helper_f); + if (bounds.is_case_2B()) { + std::cout << "BRACKETED PRIORITY 1" << std::endl; + return helper_f.solve(bounds, cache_x).second; + } else { + std::cout << "COULD NOT BRACKET PRIORITY 1" << std::endl << std::endl; + } + + bounds = Bounds2(xl, cache_x, xh); + cache_x.is_l() ? bounds.eval_l_orig(helper_f) : bounds.eval_h_orig(helper_f); + if (bounds.is_case_2B()) { + std::cout << "BRACKETED PRIORITY 2" << std::endl; + return helper_f.solve(bounds, cache_x).second; + } else { + std::cout << "l: " << bounds.l_smart() << ", h: " << bounds.h_smart() << std::endl; + std::cout << "f0_l: " << bounds.cache_l().f0() << ", f0_h: " << bounds.cache_h().f0() << std::endl; + std::cout << "COULD NOT BRACKET PRIORITY 2" << std::endl << std::endl; + } + } + + static const char* function = "boost::math::tools::detail::ContextF<%1%>"; + policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::ContextF, last closest guess was %1%", x, boost::math::policies::policy<>()); + return 0; + }; + +#if 0 + template + T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { + // Create ContextF + ContextF helper_f(f, digits, max_iter); + + // Create cache_x + const auto cache_x = helper_f(x); + if (0 == cache_x.f0()) return x; // Check if root is already found + + // Create Bounds + auto bounds = Bounds(xl, cache_x, xh); + + // Do Newton loop + const auto p = helper_f.solve(bounds, cache_x); + + if (p.first) { + std::cout << "succuss first time" << std::endl; + max_iter = helper_f.iters_used(); + return p.second; + } else { + std::cout << "FAILURE first time" << std::endl << std::endl; + + bounds = Bounds(xl, cache_x, xh); + cache_x.is_l() ? bounds.eval_h_orig(helper_f) : bounds.eval_l_orig(helper_f); + if (bounds.is_case_2B()) { + std::cout << "BRACKETED PRIORITY 1" << std::endl; + return helper_f.solve(bounds, cache_x).second; + } else { + std::cout << "COULD NOT BRACKET PRIORITY 1" << std::endl << std::endl; + } + + bounds = Bounds(xl, cache_x, xh); + cache_x.is_l() ? bounds.eval_l_orig(helper_f) : bounds.eval_h_orig(helper_f); + if (bounds.is_case_2B()) { + std::cout << "BRACKETED PRIORITY 2" << std::endl; + return helper_f.solve(bounds, cache_x).second; + } else { + std::cout << "l: " << bounds.l_smart() << ", h: " << bounds.h_smart() << std::endl; + std::cout << "f0_l: " << bounds.cache_l().f0() << ", f0_h: " << bounds.cache_h().f0() << std::endl; + std::cout << "COULD NOT BRACKET PRIORITY 2" << std::endl << std::endl; + } + } + + static const char* function = "boost::math::tools::detail::ContextF<%1%>"; + policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::ContextF, last closest guess was %1%", x, boost::math::policies::policy<>()); + return 0; + }; +#endif +#endif +#endif + +} // namespace detail + - T x() const { return x_; } - T f0() const { return f0_; } // f(x) - T f1() const { return f1_; } // f'(x) - private: - T x_; - T f0_; // f(x) - T f1_; // f'(x) - }; - // Calculates the step with Newton's method - template - class StepNewton { - public: - T calc_dx(const CacheOrder1& cache) { - return cache.f0() / cache.f1(); +template +T newton_raphson_iterate(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) +{ + return detail::solve>(f, x, xl, xh, digits, max_iter); +} + +template +inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) +{ + std::uintmax_t m = (std::numeric_limits::max)(); + return newton_raphson_iterate(f, guess, min, max, digits, m); +} + + + +#if 1 + +template +T newton_raphson_iterate_classic(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) +{ + BOOST_MATH_STD_USING + + static const char* function = "boost::math::tools::newton_raphson_iterate_classic<%1%>"; + if (min > max) + { + return policies::raise_evaluation_error(function, "Range arguments in wrong order in boost::math::tools::newton_raphson_iterate_classic(first arg=%1%)", min, boost::math::policies::policy<>()); + } + + T f0(0), f1, last_f0(0); + T result = guess; + + T factor = static_cast(ldexp(1.0, 1 - digits)); + T delta = tools::max_value(); + T delta1 = tools::max_value(); + T delta2 = tools::max_value(); + + // + // 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); + +#ifdef BOOST_MATH_INSTRUMENT + std::cout << "newton_raphson_iterate_classic, guess = " << guess << ", min = " << min << ", max = " << max + << ", digits = " << digits << ", max_iter = " << max_iter << "\n"; +#endif + + do { + last_f0 = f0; + delta2 = delta1; + delta1 = delta; + std::cout << std::endl; + std::cout << "EVAL IS: " << count << ", at: " << result << std::endl; + detail::unpack_tuple(f(result), f0, f1); + --count; + + std::cout << "result: " << result << ", f0: " << f0 << ", f1: " << f1; + + 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; + std::cout << ", delta: " << delta << std::endl; + } +#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) * 1.1f; // Protect against huge jumps! + //delta = sign(delta) * result; // Protect against huge jumps! Failed for negative result. https://github.com/boostorg/math/issues/216 + } + 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_classic, perhaps we have a local minima near current best guess of %1%", guess, boost::math::policies::policy<>()); + } + }while(count && (fabs(result * factor) < fabs(delta))); + + std::cout << "delta: " << delta << std::endl; + std::cout << "fabs(result * factor)" << fabs(result * factor) << std::endl; + std::cout << "fabs(delta) " << fabs(delta) << std::endl; + + max_iter -= count; + +#ifdef BOOST_MATH_INSTRUMENT + std::cout << "Newton Raphson required " << max_iter << " iterations\n"; +#endif + + return result; +} + +template +inline T newton_raphson_iterate_classic(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) +{ + std::uintmax_t m = (std::numeric_limits::max)(); + return newton_raphson_iterate_classic(f, guess, min, max, digits, m); +} + + +#endif + + + + + + + + + + + + + +#else + + + +///////////////////////////////// + // Enforces bracketing when root finding with Newton's method template class BracketHelper { - public: - static T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { - BracketHelper helper(f, x, xl, xh, digits, max_iter); - T x_sol = helper.solve_impl(); - max_iter -= helper.count(); // Make max_iter the number of iterations used. - return x_sol; - } - - private: // Stores the last two Newton step sizes class StepSizeHistory { public: @@ -278,6 +2331,138 @@ namespace detail { T dx_pp_; }; + // + + + + + static Bounds create_bracket_candidate(ContextF& helper_f, const detail::CacheOrder1& cache_x, T bound) { + if (is_can_eval_bound(bound)) { + return Bounds(cache_x, helper_f(bound)); + } else { + return rollout_for_bracket_candidate(helper_f, cache_x, bound); + } + } + + static + + public: + static T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { + // // Create ContextF + // ContextF helper_f(f, digits, max_iter); + + // // Create cache_x + // const auto cache_x = helper_f(x); + + // // Check if root is already found + // if (0 == cache_x.f0()) return x; + + // // Get preferred bound + // const T dx = Step().calc_dx(cache_x); + // const T bound_pref = (0 <= dx) ? xl : xh; + + // // Check if bound is valid + // const auto bc_pref = create_bracket_candidate(helper_f, cache_x, bound_pref); + + + + + const bool is_can_eval_bound = is_can_eval_bound(bound); + if (is_can_eval_bound) { + const auto cache_b = detail::CacheOrder1(f, bound); + const auto bc = Bounds(cache_x, cache_b); + if (bc.may_contain_root()) { + return solve_bc(f, bc, digits, max_iter, cache_x); + } else { + return std::make_pair(false, 0); + } + } + + + if (0 <= dx) { // Want to go lower + const auto p_l = solve_l(f, x, xl, xh, digits, max_iter, cache_x); + } else { + + } + } + + static std::pain solve_bound_prefered(F f, BranchCandidate bc, int digits, std::uintmax_t& max_iter) { + if (bc.can_bracket()) { + return solve_bracket(f, bc, digits, max_iter, cache_x); + } else if (bc.can_brent()) { + return solve_brent(f, bc, digits, max_iter, cache_x); + } else { + return std::make_pair(false, 0); + } + } + + static std::pair solve_l(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter, const detail::CacheOrder1& cache_x) { + const bool is_can_eval_bound(xl); + if (is_can_eval_bound) { + const auto bc = BranchCandidate(cache_x, detail::CacheOrder1(f, xl)); + if (bc.may_contain_root()) { + return solve_bc(f, x, xl, xh, digits, max_iter, cache_x); + } else { + return std::make_pair(false, 0); + } + } + + if (!is_can_eval_bound) { + // Find lower bound + } + + + } + + + + + // std::cout << "================= solve ====================" << std::endl; + + // const auto cache_debug = detail::CacheOrder1(f, x); + // std::cout << "x: " << cache_debug.x() << ", f0: " << cache_debug.f0() << ", f1: " << cache_debug.f1() << std::endl; + // // x: 0.56763, f0: -0.00218081, f1: 0.0721978 + + // const auto cache_debug_l = detail::CacheOrder1(f, xl); + // std::cout << "xl: " << cache_debug_l.x() << ", f0: " << cache_debug_l.f0() << ", f1: " << cache_debug_l.f1() << std::endl; + // // xl: 0, f0: -0.999001, f1: 0 + + // const auto cache_debug_y = detail::CacheOrder1(f, 1.0); + // std::cout << "y: " << cache_debug_y.x() << ", f0: " << cache_debug_y.f0() << ", f1: " << cache_debug_y.f1() << std::endl; + // // y: 1, f0: 0.000998974, f1: 1.64892e-07 + + // const auto cache_debug_z = detail::CacheOrder1(f, 2.0); + // std::cout << "z: " << cache_debug_z.x() << ", f0: " << cache_debug_z.f0() << ", f1: " << cache_debug_z.f1() << std::endl; + // // z: 2, f0: 0.000998974,flip f1: 2.88776e-33 + + + // std::cout << "xh: " << xh << std::endl; + // // xh: 3.40282e+38 + + + // Get maximum value of T + // const T x_limit = std::numeric_limits::max(); + + // std::cout << "x_limit: " << x_limit << std::endl; + // const bool bb = x_limit == xh; + // std::cout << bb << std::endl; + + // assert(false); + + const auto cache_debug_h = detail::CacheOrder1(f, xh); + std::cout << "xh: " << cache_debug_h.x() << ", f0: " << cache_debug_h.f0() << ", f1: " << cache_debug_h.f1() << std::endl; + // Did not finish + + + BracketHelper helper(f, x, xl, xh, digits, max_iter); + + T x_sol = helper.solve_impl(); + max_iter -= helper.count(); // Make max_iter the number of iterations used. + std::cout << "max_iter: " << max_iter << std::endl; + return x_sol; + } + + private: // Same as the inputs to newton_raphson_iterate BracketHelper(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) : f_(f) @@ -291,8 +2476,41 @@ namespace detail { static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; policies::raise_evaluation_error(function, "Range arguments in wrong order in boost::math::tools::detail::BracketHelper(first arg=%1%)", xl, boost::math::policies::policy<>()); } + + std::cout << "================= factor_: " << factor_ << "====================" << std::endl; + } + + + T calc_dx_next_default() const { return Step().calc_dx(cache_x_); } + T calc_x_next_default() const { return x() - calc_dx_next_default(); } + + bool can_bracket_l() const { return f0() * f0_l() <= 0; } + bool can_bracket_h() const { return f0() * f0_h() <= 0; } + + T solve_bracket_l() { + set_cache_H_to_cache_x(); // Evaluated point x becomes new high bound + return solve_bracket(); + } + T solve_bracket_h() { + set_cache_L_to_cache_x(); // Evaluated point x becomes new low bound + return solve_bracket(); + } + + bool want_to_go_l() const { return 0 <= calc_dx_next_default(); } // Because of negative sign in update formula x_next = x() - dx + + // + // + // + // + // + T solve_impl() { + if (0 == f0()) return x(); // Solved at initial x + + } + +#if 0 T solve_impl() { if (0 == f0()) return x(); // Solved at initial x @@ -315,22 +2533,152 @@ namespace detail { return policies::raise_evaluation_error(function, "There appears to be no root to be found in boost::math::tools::detail::BracketHelper, perhaps we have a local minima near current best guess of %1%", x(), boost::math::policies::policy<>()); } } +#else + +//////////////////////////////////// + + T solve_impl() { + if (0 == f0()) return x(); // Solved at initial x + + std::cout << "SOLVE_IMPL" << std::endl; + + // Bracket root in direction of dx if possible + if (want_to_go_l() && can_bracket_l()) { + std::cout << "DIRECT --> SOLVE_BRACKET_L" << std::endl; + return solve_bracket_l(); + } else if (can_bracket_h()) { + std::cout << "DIRECT --> SOLVE_BRACKET_H" << std::endl; + return solve_bracket_h(); + } + + return solve_classic_with_fallback(); + } + + T solve_classic_with_fallback() { + std::cout << "Doing classic with fallback" << std::endl; + + // Copy caches + const auto cache_x_copy = cache_x_; + const auto cache_xh_copy = cache_xh_; + const auto cache_xl_copy = cache_xl_; + + bool is_success_classic = false; + const T root = solve_classic(&is_success_classic); + + if (is_success_classic) { + return root; + } else { + // Restore caches + cache_x_ = cache_x_copy; + cache_xh_ = cache_xh_copy; + cache_xl_ = cache_xl_copy; + + // Try bracketing + if (can_bracket_l()) { + return solve_bracket_l(); + } else if (can_bracket_h()) { + return solve_bracket_h(); + } + } + + static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; + return policies::raise_evaluation_error(function, "There appears to be no root to be found in boost::math::tools::detail::BracketHelper, perhaps we have a local minima near current best guess of %1%", x(), boost::math::policies::policy<>()); + } + + // + // The classic solver is used when the root is not bracketed. If the classic solver detects + // a sign flip it passes control to the bracket solver. + // + T solve_classic(bool* is_success) { + BOOST_MATH_STD_USING + static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; + + const bool same_sign_lh = 1 == sign(f0_l()) * sign(f0_h()); + const bool same_sign_lx = 1 == sign(f0_l()) * sign(f0()); + + if (!same_sign_lh || !same_sign_lx) { + *is_success = false; + // 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<>()); + return policies::raise_evaluation_error(function, "Not same sign boost::math::tools::detail::BracketHelper %1%", xl(), boost::math::policies::policy<>()); + } + + const bool is_negate = (1 == sign(f0_l())) ? false : true; + + const auto fn_unoriented = [&](T x) { + T brent_0; + T brent_1; + detail::unpack_tuple(f_(x), brent_0, brent_1); + std::cout << " BRENT eval at: " << x << ", f(x): " << brent_0 << std::endl; + return brent_0; + }; + + const auto fn_brent = [&](T x) { + T brent_0 = fn_unoriented(x); + return is_negate ? -brent_0 : brent_0; + }; + + + for (double i = 0.0; i <= 1.0; i += 0.0001) { + // std::cout << "fn_brent(" << i << "): "; // << fn_brent(i) << std::endl; + fn_brent(i); + } + + std::cout << std::endl; + + + const auto pair_x_fx = brent_find_minima(fn_brent, xl(), xh(), 1024); + const auto x_min = pair_x_fx.first; + const auto fx_min = fn_unoriented(x_min); + + std::cout << "x_min: " << x_min << std::endl; + std::cout << "fx_min: " << fx_min << std::endl; + + if (1 == sign(f0_l()) * sign(fx_min)) { + std::cout << "CLASSIC could not find root bracket" << std::endl; + std::cout << "x: " << x() << ", f0: " << f0() << ", f1: " << cache_x_.f1() << std::endl; + std::cout << "xl: " << xl() << ", f0_l: " << f0_l() << ", f1_l: " << f1_l() << std::endl; + std::cout << "xh: " << xh() << ", f0_h: " << f0_h() << ", f1_h: " << f1_h() << std::endl; + *is_success = false; + return x_min; // Return doesn't matter + } else { + *is_success = true; + std::cout << "CLASSIC found root bracket" << std::endl; + cache_x_ = detail::CacheOrder1(f_, x_min); // Update cache_x_ with x_next + + + + return (sign(fx_min) == 1) ? solve_bracket_h() : solve_bracket_l(); + } + } + +#endif - T solve_bracket(T x_next) { + T solve_bracket() { BOOST_MATH_STD_USING static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; + std::cout << "SOLVE_BRACKET" + << "RealType: " << typeid(T).name() + << " -- count: " << count_ << "====================" << std::endl; + + if (f0_l() == 0.0) { // Low bound is root + std::cout << "Low bound is root" << std::endl; return xl(); } else if (f0_h() == 0.0) { // High bound is root + std::cout << "High bound is root" << std::endl; return xh(); } else if (0 < f0_l() * f0_h()) { // Bounds do not bracket root + std::cout << "Bounds do not bracket root" << std::endl; return policies::raise_evaluation_error(function, "Does not bracket boost::math::tools::detail::BracketHelper %1%", xl(), boost::math::policies::policy<>()); } // Orient cache_l and cache_h so that f0_l is negative and f0_h is positive if (0 < f0_l()) std::swap(cache_xl_, cache_xh_); + T x_next; do { count_--; @@ -340,20 +2688,50 @@ namespace detail { if (!is_x_between_l_and_h(x_next) || fabs(dx_pp()) < fabs(dx_newton * 4)) { // Bisection halves step size each iteration - x_next = (xl() + xh()) / 2; // Midpoint - dx_history_.reset(); // Invalitade step size history + + // if xl and xh have different signs set x_next to 0 + if (-1 == sign(xl()) * sign(xh())) { + x_next = 0; + } else { + const T x_S = (fabs(xl()) < fabs(xh())) ? xl() : xh(); + const T x_B = (fabs(xl()) < fabs(xh())) ? xh() : xl(); + const T SB_mag = fabs(xl() - xh()); + + if (x_S < SB_mag && 1 < fabs(x_S)) { // Spans orders of magnitude + x_next = copysign(sqrt(fabs(x_S)) * sqrt(fabs(x_B)), x_S); + } else { + x_next = (xl() + xh()) / 2; // Midpoint + } + } + + dx_history_.reset(); // Invalidate step size history } const T dx = x_next - x(); cache_x_ = detail::CacheOrder1(f_, x_next); // Update cache_x_ with x_next - if (fabs(dx) < fabs(x() * factor_)) { break; } // Tolerance met - if (x_next == xl() || x_next == xh() || dx == 0) { break; } + std::cout << "x: " << x() << ", dx: " << dx << ", dx_newton: " << dx_newton + << ", f1: " << cache_x_.f1() << ", f0: " << cache_x_.f0() + << ", xl: " << xl() << ", xh: " << xh() + << std::endl; + + if (fabs(dx) < fabs(x() * factor_)) { + std::cout << "EXIT: tolerance met" << std::endl; + break; } // Tolerance met + if (x_next == xl() || x_next == xh() || dx == 0) { + std::cout << "EXIT: numerical reason" << std::endl; + break; } + // Update bracket (f0() < 0.0) ? set_cache_L_to_cache_x() : set_cache_H_to_cache_x(); dx_history_.push(dx); // Store step size + + if (count_ == 0) { + std::cout << "EXIT: max_iter reached" << std::endl; + break; + } } while (count_); return x(); @@ -368,8 +2746,10 @@ namespace detail { T f0() const { return cache_x_.f0(); } T xl() const { return cache_xl_.x(); } T f0_l() const { return cache_xl_.f0(); } + T f1_l() const { return cache_xl_.f1(); } T xh() const { return cache_xh_.x(); } T f0_h() const { return cache_xh_.f0(); } + T f1_h() const { return cache_xh_.f1(); } T dx_pp() const { return dx_history_.dx_pp(); } std::uintmax_t count() const { return count_; } @@ -386,6 +2766,7 @@ namespace detail { template T newton_raphson_iterate(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { + std::cout << "================= newton_raphson_iterate ====================" << std::endl; return detail::BracketHelper>::solve(f, x, xl, xh, digits, max_iter); } @@ -396,6 +2777,10 @@ inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) noexcept return newton_raphson_iterate(f, guess, min, max, digits, m); } +#endif + + + namespace detail { struct halley_step diff --git a/test/test_root_iterations.cpp b/test/test_root_iterations.cpp index f998b6d724..9f2b4d213d 100644 --- a/test/test_root_iterations.cpp +++ b/test/test_root_iterations.cpp @@ -212,16 +212,11 @@ BOOST_AUTO_TEST_CASE( test_main ) r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance(), iters); BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits::epsilon() * 4); BOOST_CHECK_LE(iters, 12); - - // Newton next: iters = 1000; dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters); - BOOST_CHECK_CLOSE_FRACTION(dr, result, ldexp(1.0, 1 - newton_limits)); + BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits::epsilon() * 2); BOOST_CHECK_LE(iters, 5); - // BOOST_CHECK_LE(iters, 8); - - // Halley next: iters = 1000; dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); @@ -240,18 +235,11 @@ BOOST_AUTO_TEST_CASE( test_main ) r = boost::math::tools::bracket_and_solve_root(cbrt_functor_noderiv(arg), guess, 2.0, true, boost::math::tools::eps_tolerance(), iters); BOOST_CHECK_CLOSE_FRACTION((r.first + r.second) / 2, result, std::numeric_limits::epsilon() * 4); BOOST_CHECK_LE(iters, 12); - - // ******************************************************************** - // ******************************************************************** - // ******************************************************************** // Newton next: iters = 1000; dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters); BOOST_CHECK_CLOSE_FRACTION(dr, result, ldexp(1.0, 1 - newton_limits)); BOOST_CHECK_LE(iters, 5); - // BOOST_CHECK_LE(iters, 8); - - // Halley next: iters = 1000; dr = boost::math::tools::halley_iterate(cbrt_functor_2deriv(arg), guess, result / 10, result * 10, newton_limits, iters); @@ -274,6 +262,7 @@ BOOST_AUTO_TEST_CASE( test_main ) #endif #define T double + # include "ibeta_small_data.ipp" for (unsigned i = 0; i < ibeta_small_data.size(); ++i) @@ -329,5 +318,6 @@ BOOST_AUTO_TEST_CASE( test_main ) BOOST_CHECK_LE(iters, 10); } } + } diff --git a/test/test_roots.cpp b/test/test_roots.cpp index 7dcfa1201a..351f40d1a2 100644 --- a/test/test_roots.cpp +++ b/test/test_roots.cpp @@ -456,13 +456,13 @@ void test_newton_bracket() { const auto fn_test_bracket = [&](const auto fn_quartic) { // Initial search direction has no root - BOOST_CHECK_CLOSE(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.5, 0.1, 1.1, newton_limits), tol); + BOOST_CHECK_CLOSE_FRACTION(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.5, 0.1, 1.1, newton_limits), tol); // Initial search direction has root - BOOST_CHECK_CLOSE(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.8, 0.1, 1.1, newton_limits), tol); + BOOST_CHECK_CLOSE_FRACTION(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.8, 0.1, 1.1, newton_limits), tol); // Initial search direction has root - BOOST_CHECK_CLOSE(-1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.5, -1.1, 1.1, newton_limits), tol); + BOOST_CHECK_CLOSE_FRACTION(-1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.5, -1.1, 1.1, newton_limits), tol); // Initial search direction has root - BOOST_CHECK_CLOSE(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.8, -1.1, 1.1, newton_limits), tol); + BOOST_CHECK_CLOSE_FRACTION(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.8, -1.1, 1.1, newton_limits), tol); }; fn_test_bracket(quartic_pos); @@ -671,7 +671,9 @@ void test_failures() BOOST_CHECK_THROW(boost::math::tools::newton_raphson_iterate([](double x) { return std::make_pair(x * x + 1, 2 * x); }, 10.0, -12.0, 12.0, 52), boost::math::evaluation_error); BOOST_CHECK_THROW(boost::math::tools::newton_raphson_iterate([](double x) { return std::make_pair(x * x + 1, 2 * x); }, -10.0, -12.0, 12.0, 52), boost::math::evaluation_error); // There is a root, but a bad guess takes us into a local minima: +#if 0 BOOST_CHECK_THROW(boost::math::tools::newton_raphson_iterate([](double x) { return std::make_pair(boost::math::pow<6>(x) - 2 * boost::math::pow<4>(x) + x + 0.5, 6 * boost::math::pow<5>(x) - 8 * boost::math::pow<3>(x) + 1); }, 0.75, -20., 20., 52), boost::math::evaluation_error); +#endif // There is no root: BOOST_CHECK_THROW(boost::math::tools::halley_iterate([](double x) { return std::make_tuple(x * x + 1, 2 * x, 2); }, 10.0, -12.0, 12.0, 52), boost::math::evaluation_error); From 4e6c6f6be1816ea6bc53f6e9f4dd0151ca71e4e1 Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 11 Jul 2023 21:47:41 -0400 Subject: [PATCH 05/18] update newton bounds --- include/boost/math/distributions/kolmogorov_smirnov.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/boost/math/distributions/kolmogorov_smirnov.hpp b/include/boost/math/distributions/kolmogorov_smirnov.hpp index 0365b69838..6ae80632d3 100644 --- a/include/boost/math/distributions/kolmogorov_smirnov.hpp +++ b/include/boost/math/distributions/kolmogorov_smirnov.hpp @@ -382,7 +382,7 @@ inline RealType quantile(const kolmogorov_smirnov_distribution std::uintmax_t m = policies::get_max_root_iterations(); // and max iterations. return tools::newton_raphson_iterate(detail::kolmogorov_smirnov_quantile_functor(dist, p), - k, RealType(0), boost::math::tools::max_value(), get_digits, m); + k, RealType(0), RealType(1), get_digits, m); } // quantile template @@ -407,7 +407,7 @@ inline RealType quantile(const complemented2_type(dist, p), - k, RealType(0), boost::math::tools::max_value(), get_digits, m); + k, RealType(0), RealType(1), get_digits, m); } // quantile (complemented) template From dd5693fbeb8637e2fe530c959b1e80f7581b5fea Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 11 Jul 2023 21:48:39 -0400 Subject: [PATCH 06/18] added special cases --- .../detail/ibeta_inverse.hpp | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/include/boost/math/special_functions/detail/ibeta_inverse.hpp b/include/boost/math/special_functions/detail/ibeta_inverse.hpp index 79eb95c45f..031a7b4952 100644 --- a/include/boost/math/special_functions/detail/ibeta_inverse.hpp +++ b/include/boost/math/special_functions/detail/ibeta_inverse.hpp @@ -25,23 +25,16 @@ namespace boost{ namespace math{ namespace detail{ template struct temme_root_finder { - temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {} + temme_root_finder(const T t_, const T a_) : t(t_), a(a_) { + const T x_extrema = 1 / (1 + a); + assert(0 < x_extrema && x_extrema < 1); + } boost::math::tuple operator()(T x) { BOOST_MATH_STD_USING // ADL of std names T y = 1 - x; - if(y == 0) - { - T big = tools::max_value() / 4; - return boost::math::make_tuple(static_cast(-big), static_cast(-big)); - } - if(x == 0) - { - T big = tools::max_value() / 4; - return boost::math::make_tuple(static_cast(-big), big); - } T f = log(x) + a * log(y) + t; T f1 = (1 / x) - (a / (y)); return boost::math::make_tuple(f, f1); @@ -410,6 +403,14 @@ T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const Policy& pol) T lower = eta < mu ? cross : 0; T upper = eta < mu ? 1 : cross; T x = (lower + upper) / 2; + + // Special cases due to lack of numerical precision + if (cross == 1) { + return 1; + } else if (cross == 0) { + return 0; + } + x = tools::newton_raphson_iterate( temme_root_finder(u, mu), x, lower, upper, policies::digits() / 2); #ifdef BOOST_INSTRUMENT From acd0241f432c697ffbccf4de7fe3de0332162acf Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 11 Jul 2023 21:53:15 -0400 Subject: [PATCH 07/18] added more tests --- test/test_roots.cpp | 108 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 101 insertions(+), 7 deletions(-) diff --git a/test/test_roots.cpp b/test/test_roots.cpp index 351f40d1a2..32714cd33a 100644 --- a/test/test_roots.cpp +++ b/test/test_roots.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -437,7 +438,7 @@ void test_beta(T, const char* /* name */) test_inverses(ibeta_large_data); } -void test_newton_bracket() { +void test_3_attempt() { int newton_limits = static_cast(std::numeric_limits::digits * 0.6); const auto tol = ldexp(1, 1 - newton_limits); @@ -455,13 +456,20 @@ void test_newton_bracket() { }; const auto fn_test_bracket = [&](const auto fn_quartic) { - // Initial search direction has no root - BOOST_CHECK_CLOSE_FRACTION(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.5, 0.1, 1.1, newton_limits), tol); - // Initial search direction has root + // Local minimization finds root BOOST_CHECK_CLOSE_FRACTION(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.8, 0.1, 1.1, newton_limits), tol); - // Initial search direction has root + + // Local minimization fails to find root. But, the initial search direction points to the left. + // And f(x_initial) and f(x_left_bound) have opposite signs. So, we can left bracket and find + // a root. BOOST_CHECK_CLOSE_FRACTION(-1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.5, -1.1, 1.1, newton_limits), tol); - // Initial search direction has root + + // Local minimization fails to find root. Left bracketing is not attempted because f(x_initial) + // and f(x_left_bound) have the same sign. But, f(x_initial) and f(x_right_bound) have opposite + // signs. So, we can right bracket and find a root. + BOOST_CHECK_CLOSE_FRACTION(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.5, 0.1, 1.1, newton_limits), tol); + + // Two roots exist (-1 and 1). The return root is the one found with local minimization. BOOST_CHECK_CLOSE_FRACTION(1, boost::math::tools::newton_raphson_iterate(fn_quartic, 0.8, -1.1, 1.1, newton_limits), tol); }; @@ -469,6 +477,90 @@ void test_newton_bracket() { fn_test_bracket(quartic_neg); } + +void test_cusp() { + std::cout << std::endl; + std::cout << std::endl; + + // Returns a lambda with an extrema at x_center and a y offset ye + const auto fn_make_cusp = [](const double x_center, const double ye, const double scale) { + return [=](const double x) -> std::pair { + BOOST_MATH_STD_USING + const auto x_shift = x - x_center; + const auto y = sqrt(fabs(x_shift)) + ye; + const auto y_dot = boost::math::sign(x_shift) / (2 * sqrt(fabs(x_shift))); + return std::make_pair(y * scale, y_dot * scale); + }; + }; + + const auto fn_test_case = [&](const double x_center, const double ye, const double scale) { + const auto fn_case = fn_make_cusp(x_center, ye, scale); + + int newton_limits = static_cast(std::numeric_limits::digits * 0.6); + const auto tol = ldexp(1, 1 - newton_limits); + + // Bounds + const double ye_sqr = ye * ye; + const double xl = x_center - ye_sqr * 2.0; + const double xh = x_center + ye_sqr * 2.0; + const double x = xh - 1.0; + + using std::fabs; + + // No roots exist + if (0 < ye) { + BOOST_CHECK_THROW(boost::math::tools::newton_raphson_iterate(fn_case, x, xl, xh, newton_limits), boost::math::evaluation_error); + } else { + const auto x_root = boost::math::tools::newton_raphson_iterate(fn_case, x, xl, xh, newton_limits); + + // One root exists + if (ye == 0.0) { + BOOST_CHECK_CLOSE_FRACTION(x_center, x_root, tol); + + // Two roots exist (One is picked arbitrarily) + } else { + // x_root = x_center ± ye^2 + // x_root - x_center = ± ye^2 + // abs(x_root - x_center) = ye^2 + BOOST_CHECK_CLOSE_FRACTION(fabs(x_root - x_center), ye * ye, tol); + } + } + }; + + for (const double scale : {-1.0, +1.0}) { + for (double x_center = -4.0; x_center < 4.0; x_center += static_cast(1) / 7) { + for (double ye = -5.0; ye < 5.0; ye += 0.25) { + fn_test_case(x_center, ye, scale); + } + } + } +} + + +void test_solved_at_bound() { + const auto fn_unsolved = [](const double x) -> std::pair { + return std::make_pair(x * x + 1, 2 * x); + }; + const auto fn_solved_at_neg_1 = [&](const double x) -> std::pair { + if (x == -1) { return std::make_pair(0.0, 0.0); } + return fn_unsolved(x); + }; + const auto fn_solved_at_pos_1 = [&](const double x) -> std::pair { + if (x == +1) { return std::make_pair(0.0, 0.0); } + return fn_unsolved(x); + }; + + const double x_start = 0.5; + + BOOST_CHECK_NE(x_start, 0.0); + BOOST_CHECK_NE(x_start, -1.0); + BOOST_CHECK_NE(x_start, +1.0); + + BOOST_CHECK_EQUAL(-1, boost::math::tools::newton_raphson_iterate(fn_solved_at_neg_1, x_start, -1.0, 1.0, 100)); + BOOST_CHECK_EQUAL(+1, boost::math::tools::newton_raphson_iterate(fn_solved_at_pos_1, x_start, -1.0, 1.0, 100)); +} + + #if !defined(BOOST_NO_CXX11_AUTO_DECLARATIONS) && !defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) && !defined(BOOST_NO_CXX11_LAMBDAS) template void test_complex_newton() @@ -688,7 +780,9 @@ BOOST_AUTO_TEST_CASE( test_main ) test_beta(0.1, "double"); - test_newton_bracket(); + test_3_attempt(); + test_cusp(); + test_solved_at_bound(); #if !defined(BOOST_NO_CXX11_AUTO_DECLARATIONS) && !defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) && !defined(BOOST_NO_CXX11_LAMBDAS) test_complex_newton>(); From 89ac04c7b0d1aaa85b0f6c5ea05fae103d01e2c5 Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 11 Jul 2023 22:09:10 -0400 Subject: [PATCH 08/18] Refactored into state and solvers --- include/boost/math/tools/roots.hpp | 2856 +++++----------------------- 1 file changed, 448 insertions(+), 2408 deletions(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index b70befbb4f..47152212ee 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -216,228 +216,119 @@ inline std::pair bisect(F f, T min, T max, Tol tol) noexcept(policies::is_ } -#if 0 - - // - template - class Optional { - public: - Optional() : has_data_(false) {} - Optional(C c) : has_data_(true), c_(c) {} - - bool has_data() const { return has_data_; } - const C& c() const { - assert(has_data_); - return c_; - } - - void replace(C c) { c_ = c; } - - private: - bool has_data_; - C c_; - }; - - // Enum for bracketing states: case_1, case_2U, case_2B - enum class BracketState { case_1, case_2U, case_2B, case_2I }; - - // Enum for solve outputs - enum class OutputResult { success, failure }; - - // Solve output type - template - using SolveOutput = std::pair; - - - - - - -#endif - - - - - - - - - - - - - namespace detail { + ////// 1D Newton root finder explanation. ////// + // + // Consider the following three root finding problems: + // + // Case 1: Find a root of the function f(x) given bounds [x_l_orig, x_h_orig] and + // an initial evaluation (x_0, f(x_0), f'(x_0)). + // Case 2U: Find a root of the function f(x) given evaluations at a low point + // (x_l, f(x_l), f'(x_l)) and a high point (x_h, f(x_h), f'(x_h)), where + // and f(x_l) and f(x_h) have the same sign, but have slopes that indicate + // that the function is sloping towards zero between x_l and x_h, forming + // a U shape (hence the name). The U may or may not cross zero. + // Case 2B: Find a root of the function f(x) given evaluations at a low point + // (x_l, f(x_l), f'(x_l)) and a high point (x_h, f(x_h), f'(x_h)), where + // and f(x_l) and f(x_h) have opposite signs that BRACKET a root. + // + // All three cases operate on similar data, but require different algorithms to + // solve. These cases can also transitions between each other. Specifically... + // + // Case 1: The problem will remain in Case 1 as long as dx does not change sign. + // If dx changes sign, the problem will transition to Case 2U or Case 2B. + // Case 2U: The problem will transition to Case 2B if the criteria for Case 2B is + // met, i.e., f(x_l) and f(x_h) have opposite signs. + // Case 2B: Does not transition to other cases. + // + // + ////// Code organization. ////// + // + // The code for the root finding problem is divided into two parts: a state that + // contains data and solvers that are dataless. State is stored in the Root1D_State + // class. The following solver classes operate on the state: + // - SolverCase1 + // - SolverCase2U + // - SolverCase2B + // The process of finding a root begins by iterating with the solver for Case 1. + // Iterations continue until one of four following things happen: + // - Returns success (f(x) is zero, ... etc) + // - Transition to Case 2U + // - Transition to Case 2B + // - Evaluated limit without being able to transition to Case 2U or Case 2B + // Similarly, when iterating in Case 2U, iterations will continue until one of the + // following happens: + // - Transition to Case 2B + // - Returns success (f(x) is zero, ... etc) + // - Returns failure (converging to non-zero extrema) + // When iterating in Case 2B, iterations will continue until: + // - Returns success (f(x) is zero, ... etc) + // + // Enum for the possible cases for the root finding state machine + enum class CaseFlag { case_1, case_2U, case_2B, success, failure }; - -#if 1 - - // constexpr bool DEBUG_PRINT = false; - constexpr bool DEBUG_PRINT = true; - - - // Stores the last two Newton step sizes - template - class StepSizeHistory { - public: - StepSizeHistory() { reset(); } - - void reset() { - // NOTE: for T = boost::math::concepts::real_concept, both std::numeric_limits::infinity() - // and std::numeric_limits::max() are 0 - dx_p_ = static_cast(1) / 0; - dx_pp_ = static_cast(1) / 0; - } - - void push(T input) { - dx_pp_ = dx_p_; - dx_p_ = input; - } - - T dx_p() const { return dx_p_; } - T dx_pp() const { return dx_pp_; } - - private: - T dx_p_; - T dx_pp_; - }; - - // Stores x, f(x), and f'(x) + // Stores x, f(x), and f'(x) from a function evaluation template - class CacheOrder1 { + class Data_X01 { public: template - CacheOrder1(F f, T x) : x_(x), has_f0_f1_(true) { detail::unpack_tuple(f(x_), f0_, f1_); } - CacheOrder1(T x) : x_(x) , has_f0_f1_(false) {} - - T calc_dx() const { return Step().calc_dx(*this); } - - bool is_l() const { return 0 <= calc_dx(); } - bool is_h() const { return !is_l(); } - - bool is_B(const CacheOrder1& z) const { return !is_same_sign(z); } - bool is_I(const CacheOrder1& z) const { return !is_B(z) && !is_U(z); } - bool is_U(const CacheOrder1& z) const { - if (is_ordered(z)) { - return is_U_ordered(*this, z); - } else { - return is_U_ordered(z, *this); - } - } - std::pair, CacheOrder1> sort(const CacheOrder1& z) const { - if (is_ordered(z)) { - return std::make_pair(*this, z); - } else { - return std::make_pair(z, *this); - } - } - bool is_ordered(const CacheOrder1& z) const { return x() < z.x(); } - bool is_U_ordered(const CacheOrder1& l, const CacheOrder1& h) const { - const bool is_pos = 0 < l.f0(); + Data_X01(F f, T x) : x_(x) { detail::unpack_tuple(f(x_), f0_, f1_); } - T l_slope = l.f1(); - T h_slope = h.f1(); + // Store just a constant x value + Data_X01(T x) : x_(x), f0_(static_cast(1) / 0), f1_(static_cast(1) / 0) {} - if (!is_pos) { - l_slope *= -1; - h_slope *= -1; - } + T calc_dx() const { return Step().calc_dx(*this); } + T calc_dx_newton() const { return -f0_ / f1_; } - return (l_slope < 0) && (0 < h_slope); - } - bool is_smaller(const CacheOrder1& z) const { - using std::fabs; - return fabs(f0_) < fabs(z.f0()); - } + // This evaluation is a high bound if the Newton step is negative. This function checks the following + // without division: -f(x) / f'(x) < 0 + bool is_bound_h() const { return sign(f0_) == sign(f1_); } - bool is_same_sign(const CacheOrder1& z) const { return sign(f0_) == sign(z.f0()); } - bool is_same_slope(const CacheOrder1& z) const { return sign(f1_) == sign(z.f1()); } + // Returns true if *this and z BRACKET a root. For example between: f0(x_this) = 2 and f0(x_z) = -1 OR + // f0(x_this) = 0 and f0(x_z) = -1 + bool is_B(const Data_X01& z) const { return !is_same_sign_f0(z); } - void print() const { - if (!DEBUG_PRINT) return; + bool is_same_sign_f0(const Data_X01& z) const { return sign(f0_) == sign(z.f0()); } - std::cout << "x: " << x_; - if (has_f0_f1_) { - assert(has_f0_f1_); - std::cout << ", f0: " << f0_ << ", f1: " << f1_; - } - std::cout << std::endl; + CaseFlag identify_status_flag(const Data_X01& h) const { + if (is_B(h)) { return CaseFlag::case_2B; } + if (is_U(h)) { return CaseFlag::case_2U; } + return CaseFlag::failure; } - bool is_x_only() const { return !has_f0_f1_; } - bool is_x_f0_f1() const { return has_f0_f1_; } - - bool is_dx_small(const T& factor) const { - // |dx| ≤ |x| * factor - // |f0| / |f1| ≤ |x| * factor - // |f0| ≤ |x| * |f1| * factor - using std::fabs; - - const bool bool_eval = fabs(f0_) <= fabs(x_ * f1_) * factor; - const bool bool_eval_2 = fabs(f0_) <= fabs(x_ * f1_ * factor); - - if (DEBUG_PRINT) { - std::cout.precision(200); - std::cout << "fabs(f0_): " << fabs(f0_) << std::endl; - std::cout << "fabs(x_ * f1_): " << fabs(x_ * f1_) << std::endl; - std::cout << "factor: " << factor << std::endl; - std::cout << "bool_eval; " << bool_eval << std::endl; - std::cout << "bool_eval_2; " << bool_eval_2 << std::endl; - std::cout << "fabs(f0_) == 0 " << (fabs(f0_) == 0) << std::endl; - std::cout << "fabs(x_ * f1_) * factor == 0 " << (fabs(x_ * f1_) * factor == 0) << std::endl; - - // Print data type T to std::cout - std::cout << "T: " << typeid(T).name() << std::endl; - std::cout << std::endl; - } - - - return bool_eval; - } + // Makes debugging easier if required + void print() const { std::cout << "x: " << x_ << ", f0: " << f0_ << ", f1: " << f1_ << std::endl; } T x() const { return x_; } - T f0() const { assert(has_f0_f1_); return f0_; } // f(x) - T f1() const { assert(has_f0_f1_); return f1_; } // f'(x) + T f0() const { return f0_; } // f(x) + T f1() const { return f1_; } // f'(x) private: + bool is_sorted_x(const Data_X01& h) const { return x() <= h.x(); } + + // Returns true if *this and h form a U shape. + bool is_U(const Data_X01& h) const { + const auto& l = *this; + if (!l.is_sorted_x(h)) { return h.is_U(l); } + return (sign(l.f0()) * sign(l.f1()) == -1) && // Above zero and sloping down OR below zero and sloping up + (sign(h.f0()) * sign(h.f1()) == +1); // Above zero and sloping up OR below zero and sloping down + } + T x_; T f0_; // f(x) T f1_; // f'(x) - bool has_f0_f1_; }; - // - // + // Stores x and dx for a potential next root finding evaluation template - class NEXT { + class X_DX { public: - NEXT static FromX(T x) { - return NEXT(x); - } - - NEXT static TakeStep(const CacheOrder1& cx) { - const T dx = cx.calc_dx(); - return NEXT(cx.x() + dx, dx); - } - - NEXT static From_x_dx(T x, T dx) { - return NEXT(x, dx); - } + X_DX(T x, T dx) : x_(x), dx_(dx) {} T x() const { return x_; } T dx() const { return dx_; } - void print() const { - if (!DEBUG_PRINT) return; - std::cout << "X_DX_NEXT -- x: " << x_; - std::cout << ", dx: " << dx_; - std::cout << std::endl; - } - private: - NEXT(T x) : x_(x), dx_(static_cast(1) / 0) {} - NEXT(T x, T dx) : x_(x), dx_(dx) {} - T x_; T dx_; }; @@ -446,2328 +337,480 @@ namespace detail { template class StepNewton { public: - T calc_dx(const CacheOrder1& z) const { return -z.f0() / z.f1(); } - }; - - // Enum for status flags - enum class StatusFlag { normal, dx_small, out_of_iterations, - case_1, case_2U, case_2B, case_2I, - success, failure, stall, f0_is_zero, case_0 + T calc_dx(const Data_X01& z) const { return -z.f0() / z.f1(); } }; - - //////////////////////////////////////////////////////////////////////////////////////// - - // - // + // The state of the 1D root finding problem. This state is acted upon by one of + // several solvers template - class ContextF { - public: - ContextF(F f, T xl, T xh, int digits, std::uintmax_t max_iter) - : f_(f) - , factor_(static_cast(ldexp(1.0, 1 - digits))) - , max_iter_(max_iter) - , count_(max_iter) - , l_orig_(xl) - , h_orig_(xh) - { - max_iter_ = std::min(max_iter_, static_cast(400)); - count_ = std::min(count_, static_cast(400)); - } - - CacheOrder1 operator()(const NEXT& next) { - const T x = next.x(); - - if (count_ <= 0) { - static const char* function = "boost::math::tools::detail::ContextF<%1%>"; - policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::ContextF, last closest guess was %1%", x, boost::math::policies::policy<>()); - } - - if (DEBUG_PRINT) std::cout << "EVAL IS: " << count_ << ", at: " << x << std::endl; - - --count_; - return CacheOrder1(f_, x); - } - - T dx_p() const { return dx_history_.dx_p(); } - T dx_pp() const { return dx_history_.dx_pp(); } - std::uintmax_t iters_used() const { return max_iter_ - count_; } - - T factor() const { return factor_; } - std::uintmax_t max_iter() const { return max_iter_; } - std::uintmax_t count() const { return count_; } - StepSizeHistory& dx_history() { return dx_history_; } - + class Root1D_State { private: - F f_; - T factor_; - std::uintmax_t max_iter_; - std::uintmax_t count_; - StepSizeHistory dx_history_; - T l_orig_; - T h_orig_; - }; + enum class SideSelect { sign_f0, // Updates based on the sign of f(x) + sign_dx, // Updates based on the sign of dx + not_last, // Updates entry not updated last + l, // Updates the low bound + h }; // Updates the high bound + // Allows template dispatch of update_cache + template + class Val {}; - -//////////////////////////////////////////////////////////////////////////////////////// - - // - // - template - class BoundState { public: - BoundState(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) - : cache_l_(xl) - , cache_h_(xh) - , flag_(StatusFlag::case_1) - , context_(f, xl, xh, digits, max_iter) - { - const auto next = NEXT::FromX(x); - const auto cache_x_temp = context_(next); - update_cache_sign_dx(cache_x_temp); - } - - T calc_dx_true(const CacheOrder1& cx) const { - return cx.x() - cache_x().x(); - } - - // CacheOrder1 operator()(T x) { return context_(x); } - CacheOrder1 operator()(const NEXT& next) { return context_(next); } - - void print_dx_history() const { - if (!DEBUG_PRINT) return; - // std::cout << "dx_p: " << dx_p() << std::endl; - // std::cout << "dx_pp: " << dx_pp() << std::endl; + Root1D_State(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) + : data_(xl, xh) + , flag_(CaseFlag::case_1) + , evaluator_(f, digits, max_iter) + { + update_cache_sign_dx(x); } - void update_cache_l(const CacheOrder1& cx) { - if (DEBUG_PRINT) std::cout << "set dx_true to: " << calc_dx_true(cx) << std::endl; - if (DEBUG_PRINT) std::cout << "dx_newton: " << cx.calc_dx() << std::endl; - dx_history().push(calc_dx_true(cx)); - print_dx_history(); + // Syntactic sugar for type dispatch of update_cache + template + void update_cache_sign_f0(Args... args) { update_cache(Val(), args...); } + template + void update_cache_sign_dx(Args... args) { update_cache(Val(), args...); } + template + void update_cache_not_last(Args... args) { update_cache(Val(), args...); } + template + void update_cache_l(Args... args) { update_cache(Val(), args...); } + template + void update_cache_h(Args... args) { update_cache(Val(), args...); } - if (DEBUG_PRINT) std::cout << "UPDATE_CACHE_L to: " << cx.x() << std::endl; + // Filter arguments into preferred format + template + void update_cache(Val vs, const T x) { update_cache(vs, calc_x_dx(x)); } + template + void update_cache(Val vs, const X_DX& x_dx) { update_cache(vs, eval_next(x_dx), x_dx.dx()); } + template + void update_cache(Val vs, const Data_X01& cx) { update_cache(vs, cx, calc_dx_true(cx)); } - cache_l_ = cx; - is_last_eval_l_ = true; + // Resolve SideSelect + void update_cache(Val, const Data_X01& cx, const T dx) { + const bool is_bound_h = cx.is_same_sign_f0(cache_h()); // Updates high bound if f0 is same sign as f0_h + // TODO: change to update_cache + update_cache_dispatch(is_bound_h, cx, dx); } - void update_cache_h(const CacheOrder1& cx) { - if (DEBUG_PRINT) std::cout << "set dx_true to: " << calc_dx_true(cx) << std::endl; - if (DEBUG_PRINT) std::cout << "dx_newton: " << cx.calc_dx() << std::endl; - dx_history().push(calc_dx_true(cx)); - print_dx_history(); - - if (DEBUG_PRINT) std::cout << "UPDATE_CACHE_H to: " << cx.x() << std::endl; - - cache_h_ = cx; - is_last_eval_l_ = false; - } - - void update_cache_sign_f0(const CacheOrder1& cx) { - if (cx.is_same_sign(cache_l_)) { - update_cache_l(cx); - } else { - update_cache_h(cx); - } + void update_cache(Val, const Data_X01& cx, const T dx) { + const bool is_bound_h = cx.is_bound_h(); // Updates high bound if dx is negative + update_cache_dispatch(is_bound_h, cx, dx); } - - void update_cache_sign_dx(const CacheOrder1& cx) { - if (cx.is_l()) { - update_cache_l(cx); - } else { - update_cache_h(cx); - } + void update_cache(Val, const Data_X01& cx, const T dx) { + const bool is_bound_h = !data_.ind_last(); // Updates high bound if last eval was not high + update_cache_dispatch(is_bound_h, cx, dx); } - - bool is_converge_too_slow(const T& dx) const { - using std::fabs; - return fabs(dx_pp()) < fabs(dx * 4); + void update_cache(Val, const Data_X01& cx, const T dx) { + update_cache_dispatch(false, cx, dx); // Updates low bound } - - T midpoint() const { return (cache_l_.x() + cache_h_.x()) / 2; } - bool is_inbounds(T x) const { return (l_smart() <= x) && (x <= h_smart()); } - - T l_smart() const { return cache_l_.x(); } - T h_smart() const { return cache_h_.x(); } - - bool is_dx_small(const NEXT& next) const { - using std::fabs; - - const T side_1 = fabs(next.dx()); - const T side_2 = fabs(next.x()) * factor(); - - // if (DEBUG_PRINT) std::cout << "side_1: " << side_1 << std::endl; - // if (DEBUG_PRINT) std::cout << "side_2: " << side_2 << std::endl; - - return side_1 <= side_2; - - // return fabs(dx_p()) <= fabs(cache_x().x()) * factor(); + void update_cache(Val, const Data_X01& cx, const T dx) { + update_cache_dispatch(true, cx, dx); // Updates high bound } - void set_flag_if_dx_small(const NEXT& next) { - if (cache_x().f0() == 0) { - set_flag(StatusFlag::f0_is_zero); - } else if (is_dx_small(next)) { - set_flag(StatusFlag::dx_small); - } - } - void set_flag_if_count_zero() { - if (count() <= 0) { - set_flag(StatusFlag::out_of_iterations); - if (DEBUG_PRINT) std::cout << "OUT OF ITERATIONS: count: " << count() << std::endl; - } - } - void set_flag_if_stall(const NEXT& next) { - if (next.dx() == 0) { set_flag(StatusFlag::stall); } + // Update the cache + void update_cache_dispatch(const bool is_bound_h, const Data_X01& cx, const T dx) { + dx_history_.push(dx); + data_.set(cx, is_bound_h); + if (cx.f0() == 0) { flag_ = CaseFlag::success; } } - void materialize_l() { - if (is_last_eval_l_) { return; } - cache_l_ = (*this)(l_smart()); - is_last_eval_l_ = true; - } - void materialize_h() { - if (!is_last_eval_l_) { return; } - cache_h_ = (*this)(h_smart()); - is_last_eval_l_ = false; - } + void update_case_flag() { + if (flag() == CaseFlag::success) { return; } // Don't update if already successful - const CacheOrder1& cache_x() const { - if (is_last_eval_l_) { return cache_l_; } else { return cache_h_; } + // Find status flag if both caches are valid + if (num_valid_caches() == 2) { + flag_ = cache_l().identify_status_flag(cache_h()); + } } - int num_x_only() const { return cache_l_.is_x_only() + cache_h_.is_x_only(); } - - T get_limit() const { - assert(num_x_only() == 1); - return is_last_eval_l_ ? cache_h_.x() : cache_l_.x(); - } + T midpoint() const { return (cache_l().x() + cache_h().x()) / 2; } + bool is_inbounds(T x) const { return (x_l() <= x) && (x <= x_h()); } - void set_state_2() { - assert(num_x_only() == 0); + // NOTE: caches always contain valid x data, as both caches are initalized with the + // bounds of the root finding problem. + T x_l() const { return cache_l().x(); } + T x_h() const { return cache_h().x(); } - const auto& cache_l = cache_l_; - const auto& cache_h = cache_h_; + // Returns the last cache that was evaluated + const Data_X01& cache_last() const { return data_.cache_last(); } + + int num_x_only() const { return data_.num_x_only(); } + int num_valid_caches() const { return 2 - num_x_only(); } - if (cache_l.is_B(cache_h)) { - flag_ = StatusFlag::case_2B; - } else if (cache_l.is_U(cache_h)) { - flag_ = StatusFlag::case_2U; - } else { - assert(cache_l.is_I(cache_h)); - flag_ = StatusFlag::case_2I; - } + // Extrapolates the last evaluation to calculate a next x value, and the distance dx + // between the last x and this next x. + X_DX extrapolate_last_to_calc_x_dx() const { + const auto& cx = cache_last(); + const T x_next = cx.x() + cx.calc_dx(); + return calc_x_dx(x_next); + // NOTE: the code below is incorrect because (x + dx) - x != dx + // const T x = cx.x(); + // const T dx = cx.calc_dx(); + // return X_DX(x + dx, dx); } + X_DX calc_x_dx(T x) const { return X_DX(x, calc_dx_true(x)); } - void print() const { - if (!DEBUG_PRINT) return; - std::cout << "is_last_eval_l_: " << is_last_eval_l_ << std::endl; - std::cout << "l: "; cache_l_.print(); - std::cout << "h: "; cache_h_.print(); - std::cout << "flag: "; - print_root_status_flag(); - } + // Calculates the x and dx for the next step by subtracting the current x from the next x. + // This approach to calculating dx is aware of floating point precision issues for small dx. + T calc_dx_true(const Data_X01& cx) const { return calc_dx_true(cx.x()); } + T calc_dx_true(const T& x) const { return x - cache_last().x(); } - void print_root_status_flag() const { - if (!DEBUG_PRINT) return; + // Makes debugging easier if required + void print_status_flag() const { switch (flag()) { - case StatusFlag::normal: std::cout << "normal" << std::endl; break; - case StatusFlag::dx_small: std::cout << "dx_small" << std::endl; break; - case StatusFlag::out_of_iterations: std::cout << "out_of_iterations" << std::endl; break; - case StatusFlag::case_1: std::cout << "case_1" << std::endl; break; - case StatusFlag::case_2U: std::cout << "case_2U" << std::endl; break; - case StatusFlag::case_2B: std::cout << "case_2B" << std::endl; break; - case StatusFlag::case_2I: std::cout << "case_2I" << std::endl; break; - case StatusFlag::success: std::cout << "success" << std::endl; break; - case StatusFlag::failure: std::cout << "failure" << std::endl; break; - case StatusFlag::stall: std::cout << "stall" << std::endl; break; - case StatusFlag::f0_is_zero: std::cout << "f0_is_zero" << std::endl; break; - case StatusFlag::case_0: std::cout << "case_0" << std::endl; break; + case CaseFlag::case_1: std::cout << "case_1" << std::endl; break; + case CaseFlag::case_2U: std::cout << "case_2U" << std::endl; break; + case CaseFlag::case_2B: std::cout << "case_2B" << std::endl; break; + case CaseFlag::success: std::cout << "success" << std::endl; break; + case CaseFlag::failure: std::cout << "failure" << std::endl; break; } } - - void set_flag(StatusFlag flag) { flag_ = flag; } - - bool is_B() const { return cache_l_.is_B(cache_h_); } - bool is_U() const { return cache_l_.is_U(cache_h_); } - bool is_I() const { return cache_l_.is_I(cache_h_); } - - StepSizeHistory& dx_history() { return context_.dx_history(); } - std::uintmax_t count() const { return context_.count(); } - T factor() const { return context_.factor(); } - std::uintmax_t max_iter() const { return context_.max_iter(); } - T dx_p() const { return context_.dx_p(); } - T dx_pp() const { return context_.dx_pp(); } - std::uintmax_t iters_used() const { return context_.iters_used(); } - - NEXT create_next_for_x(T x) { - const T dx = x - cache_x().x(); - return NEXT::From_x_dx(x, dx); - } - const CacheOrder1& cache_l() const { return cache_l_; } - const CacheOrder1& cache_h() const { return cache_h_; } - StatusFlag flag() const { return flag_; } + // Returns true if the two caches bracket a root where the sign of f(x) changes + bool is_B() const { return cache_l().is_B(cache_h()); } - void set_x_next_value(T x) { x_next_value_ = x; } - T x_next_value() const { return x_next_value_; } + void set_flag(CaseFlag flag) { flag_ = flag; } + Data_X01 eval_next(const X_DX& x_dx) { return evaluator_(x_dx); } + void reset_dx_history() { dx_history_.reset(); } - private: - bool is_last_eval_l_; - CacheOrder1 cache_l_; - CacheOrder1 cache_h_; - StatusFlag flag_; - ContextF context_; - T x_next_value_; - }; + std::uintmax_t count() const { return evaluator_.count(); } + T factor() const { return evaluator_.factor(); } + T dx_pp() const { return dx_history_.dx_pp(); } + std::uintmax_t num_fn_evals() const { return evaluator_.num_fn_evals(); } -//////////////////////////////////////////////////////////////////////////////////////// + T get_active_bound() const { return data_.get_bound_from_unused_cache(); } - template - class Bound2B; - - template - class Bound2U; + const Data_X01& cache_l() const { return data_.cache_l(); } + const Data_X01& cache_h() const { return data_.cache_h(); } + CaseFlag flag() const { return flag_; } - template - class Bound1; + private: + // Holds Evaluation Data + class Root1D_Data { + public: + // Creates two caches whose x values are initialized with the bounds of the root finding + // problem, but whose f0 and f1 values are invalid. + Root1D_Data(T xl, T xh) + : v_cache_({Data_X01(xl), Data_X01(xh)}) + , v_is_x_only_({true, true}) + , ind_last_(true) {} -//////////////////////////////////////////////////////////////////////////////////////// + // Updates the cache using the bool is_h as an index. + void set(const Data_X01& cx, bool is_h) { + v_cache_[is_h] = cx; + v_is_x_only_[is_h] = false; + ind_last_ = is_h; + } - // - // - template - class BoundFn { - public: - void do_debug_printout(BoundState& b) { - b.print(); - if (DEBUG_PRINT) std::cout << std::endl; - } + // Returns the x value of the only cache that contains an original bound + T get_bound_from_unused_cache() const { + assert(num_x_only() == 1); + return v_is_x_only_[0] ? cache_l().x() : cache_h().x(); + } - void solve_next(BoundState& b) { - do_debug_printout(b); + const Data_X01& cache_last() const { return v_cache_[ind_last_]; } + const Data_X01& cache_l() const { return v_cache_[0]; } + const Data_X01& cache_h() const { return v_cache_[1]; } + + int num_x_only() const { return v_is_x_only_[0] + v_is_x_only_[1]; } - while (is_bound_still_valid(b)) { - const auto next = calc_next(b); - if (DEBUG_PRINT) next.print(); + bool ind_last() const { return ind_last_; } // Last evaluated index - b.set_x_next_value(next.x()); + private: + std::array, 2> v_cache_; + std::array v_is_x_only_; + bool ind_last_; + }; - if (is_stop(b, next)) { - if (DEBUG_PRINT) std::cout << "GOT STOP SIGNAL" << std::endl; - break; + // Stores the size of the last two steps. Calling reset will set all step sizes to infinity. + class StepSizeHistory { + public: + StepSizeHistory() { reset(); } + + void reset() { + // NOTE: for T = boost::math::concepts::real_concept, the following evaluate to 0: + // std::numeric_limits::infinity() and std::numeric_limits::infinity() + dx_p_ = static_cast(1) / 0; + dx_pp_ = static_cast(1) / 0; } - const auto cx = b(next); - this->update_cache(b, cx); - do_debug_printout(b); - } - if (DEBUG_PRINT) std::cout << "DONE WITH LOOP" << std::endl; - } + void push(T input) { + dx_pp_ = dx_p_; + dx_p_ = input; + } - NEXT calc_next(BoundState& b) { - // Attempt Newton step - const auto& cx = b.cache_x(); - const auto next = NEXT::TakeStep(cx); + T dx_pp() const { return dx_pp_; } - const bool c1 = !(b.is_inbounds(next.x())); - const bool c2 = b.is_converge_too_slow(next.dx()); + private: + T dx_p_; + T dx_pp_; + }; - if (c1) { - if (DEBUG_PRINT) std::cout << "NOT INBOUNDS x_next: " << next.x() << std::endl; - } - if (c2) { - if (DEBUG_PRINT) std::cout << "CONVERGE TOO SLOW dx_next: " << next.dx() << std::endl; + // Stores the problem setup for the 1D root finding algorithm + class Root1D_Evaluator { + public: + Root1D_Evaluator(F f, int digits, std::uintmax_t max_iter) + : f_(f) + , factor_(static_cast(ldexp(1.0, 1 - digits))) + , max_iter_(max_iter) + , count_(max_iter) {} + + Data_X01 operator()(const X_DX& x_dx) { + if (count_ <= 0) { + static const char* function = "boost::math::tools::detail::Root1D_Evaluator<%1%>"; + policies::raise_evaluation_error(function, "Ran out of iterations when finding root in boost::math::tools::detail::Root1D_Evaluator, attempted to evalutate %1%", x_dx.x(), boost::math::policies::policy<>()); + } + --count_; + return Data_X01(f_, x_dx.x()); } - // if (!(b.is_inbounds(next.x())) || b.is_converge_too_slow(next.dx())) { - if (c1 || c2) { - b.dx_history().reset(); - if (DEBUG_PRINT) std::cout << "TOOK BISECTION wanted: " << next.x() << ", dx: " << next.dx() << std::endl; - return this->calc_next_bisection(b); - } else { - if (DEBUG_PRINT) std::cout << "TOOK NEWTON -------------" << std::endl; - return next; - } - } + std::uintmax_t num_fn_evals() const { return max_iter_ - count_; } - NEXT calc_next_bisection_regular(BoundState& b) { - const auto x_mid = b.midpoint(); - if (x_mid == b.l_smart() || x_mid == b.h_smart()) { - b.set_flag(StatusFlag::dx_small); - } - // const auto next = NEXT::FromX(x_mid); - const auto next = b.create_next_for_x(x_mid); - return next; - } + T factor() const { return factor_; } + std::uintmax_t count() const { return count_; } + + private: + F f_; + T factor_; + std::uintmax_t max_iter_; + std::uintmax_t count_; + }; - virtual bool is_stop(BoundState& b, const NEXT next) = 0; - virtual NEXT calc_next_bisection(BoundState& b) = 0; - virtual void update_cache(BoundState& b, const CacheOrder1& cx) = 0; - virtual void print_class() const = 0; - // virtual void is_bound_still_valid() const = 0; - virtual bool is_bound_still_valid(const BoundState& b) const = 0; + // // // Data members // // // + Root1D_Data data_; + CaseFlag flag_; + StepSizeHistory dx_history_; + Root1D_Evaluator evaluator_; }; - // - // +//////////////////////////////////////////////////////////////////////////////////////// + + // The abstract base class for the 1D root finding algorithm template - class Bound1 : public BoundFn { + class SolverBase { public: - bool is_stop(BoundState& b, const NEXT next) override { - // if (b.num_x_only() != 2) { - b.set_flag_if_dx_small(next); - b.set_flag_if_count_zero(); - b.set_flag_if_stall(next); - // } - return b.flag() != StatusFlag::case_1; - } + // Solves the root finding problem for this status flag + T solve(Root1D_State& b) const { + X_DX x_dx = b.extrapolate_last_to_calc_x_dx(); + + while (is_solver_current(b)) { + // Calculate next x and dx + x_dx = b.extrapolate_last_to_calc_x_dx(); + if (is_need_bisection(b, x_dx)) { + b.reset_dx_history(); + x_dx = calc_next_bisection(b); + } - NEXT calc_next_bisection(BoundState& b) override { - const T x = b.cache_x().x(); - const T x_limit = b.get_limit(); - const T dx_limit = x_limit - x; - const T dx_8 = b.cache_x().calc_dx() * 1000; + if (is_exit_case_now(b, x_dx)) { break; } - using std::fabs; - if (fabs(dx_limit) < fabs(dx_8)) { - return b.create_next_for_x(x_limit); - // b.set_x_dx_next_for_x(x_limit); - // b.materialize_l(); - // b.materialize_h(); - // b.set_state_2(); - } else { - // update_cache(b, b(x + dx_8)); - // b.set_x_dx_next_for_x(x + dx_8); - return b.create_next_for_x(x + dx_8); + this->update_cache(b, x_dx); } + + return x_dx.x(); } - void update_cache(BoundState& b, const CacheOrder1& cx) override { - b.update_cache_sign_dx(cx); - if (b.num_x_only() == 0) { - b.set_state_2(); - } else { - if (b.get_limit() == cx.x()) { - b.set_flag(StatusFlag::case_2I); - } - } + bool is_need_bisection(Root1D_State& b, const X_DX& x_dx) const { + using std::fabs; + return !(b.is_inbounds(x_dx.x())) || + fabs(b.dx_pp()) < fabs(x_dx.dx() * 4); // Converges too slow + } + + bool is_dx_small(const Root1D_State& b, const X_DX x_dx) const { + using std::fabs; + return fabs(x_dx.dx()) <= fabs(x_dx.x()) * b.factor(); } - void print_class() const override { - if (!DEBUG_PRINT) return; - std::cout << "Bound1" << std::endl; } + bool is_solver_current(const Root1D_State& b) const { return b.flag() == flag_this(); } - bool is_bound_still_valid(const BoundState& b) const override { return b.flag() == StatusFlag::case_1; } + virtual bool is_exit_case_now(Root1D_State& b, const X_DX x_dx) const = 0; + virtual X_DX calc_next_bisection(Root1D_State& b) const = 0; + virtual void update_cache(Root1D_State& b, const X_DX& x_dx) const = 0; + virtual CaseFlag flag_this() const = 0; }; - // - // + // The solver for the case where only a single bracted (either low or high) has been + // evaluated. template - class Bound2B : public BoundFn { + class SolverCase1 : public SolverBase { public: - bool is_stop(BoundState& b, const NEXT next) override { - b.set_flag_if_dx_small(next); - b.set_flag_if_count_zero(); - b.set_flag_if_stall(next); - return b.flag() != StatusFlag::case_2B; + X_DX calc_next_bisection(Root1D_State& b) const override { + const T x_limit = b.get_active_bound(); + const auto next_limit = b.calc_x_dx(x_limit); + b.update_cache_not_last(next_limit); + b.update_case_flag(); + return next_limit; } - NEXT calc_next_bisection(BoundState& b) override { - return this->calc_next_bisection_regular(b); + void update_cache(Root1D_State&b, const X_DX& x_dx) const override { + b.update_cache_sign_dx(x_dx); + b.update_case_flag(); } - void update_cache(BoundState& b, const CacheOrder1& cx) override { - b.update_cache_sign_f0(cx); + bool is_exit_case_now(Root1D_State& b, const X_DX x_dx) const override { + if (this->is_dx_small(b, x_dx)) { + b.set_flag(CaseFlag::success); return true; } + return b.flag() != flag_this(); } - void print_class() const override { - if (!DEBUG_PRINT) return; - std::cout << "Bound2B" << std::endl; } - - bool is_bound_still_valid(const BoundState& b) const override { return b.flag() == StatusFlag::case_2B; } + CaseFlag flag_this() const override { return CaseFlag::case_1; } }; - // - // + // Base class for SolverCase2U and SolverCase2B template - class Bound2U : public BoundFn { + class SolverCase2Base : public SolverBase { public: - bool is_stop(BoundState& b, const NEXT next) override { - b.set_flag_if_dx_small(next); - b.set_flag_if_count_zero(); - return b.flag() != StatusFlag::case_2U; + X_DX calc_next_bisection(Root1D_State& b) const override { + const auto x_mid = b.midpoint(); + if (x_mid == b.x_l() || x_mid == b.x_h()) { + b.set_flag(this->flag_stall()); + } + return b.calc_x_dx(x_mid); } - NEXT calc_next_bisection(BoundState& b) override { - return this->calc_next_bisection_regular(b); - } + virtual CaseFlag flag_stall() const = 0; + }; + + // The solver for Case 2U + template + class SolverCase2U : public SolverCase2Base { + public: + void update_cache(Root1D_State& b, const X_DX& x_dx) const override { + const auto cx = b.eval_next(x_dx); - void update_cache(BoundState& b, const CacheOrder1& cx) override { - if (is_dx_bigger_than_prev(b, cx)) { - b.set_flag(StatusFlag::case_2I); + // Transition to Case 2B because will bracket a root with the new evaluation + if (b.cache_last().is_B(cx)) { + b.set_flag(CaseFlag::case_2B); + + // If the newton step is increasing this suggests that the solver is converging toward a + // local minimum that is not a root. This suggests failure. + } else if (is_newton_dx_bigger_than_prev(b, cx)) { + b.set_flag(CaseFlag::failure); } - b.update_cache_sign_dx(cx); + + b.update_cache_sign_dx(cx, x_dx.dx()); } - bool is_dx_bigger_than_prev(BoundState& b, const CacheOrder1& cx) const { + bool is_newton_dx_bigger_than_prev(Root1D_State& b, const Data_X01& cx) const { using std::fabs; - const auto& cache_same_side = cx.is_l() ? b.cache_l() : b.cache_h(); - const T dx_cx = cx.calc_dx(); - const T dx_same_side = cache_same_side.calc_dx(); - const bool is_dx_bigger = fabs(dx_same_side) < fabs(dx_cx); - return is_dx_bigger; + // NOTE: b.cache_last() isn't generally on the same side of the root as cx + const auto& cx_prev_side_eval = cx.is_bound_h() ? b.cache_h() : b.cache_l(); + return fabs(cx_prev_side_eval.calc_dx_newton()) < fabs(cx.calc_dx_newton()); // Uses Newton step + } + + bool is_exit_case_now(Root1D_State& b, const X_DX x_dx) const override { + if (x_dx.dx() == 0) { + b.set_flag(CaseFlag::failure); return true; + } + return b.flag() != flag_this(); } - void print_class() const override { - if (!DEBUG_PRINT) return; - std::cout << "Bound2U" << std::endl; } - - bool is_bound_still_valid(const BoundState& b) const override { return b.flag() == StatusFlag::case_2U; } + CaseFlag flag_this() const override { return CaseFlag::case_2U; } + CaseFlag flag_stall() const override { return CaseFlag::failure; } }; + // The solver for Case 2B + template + class SolverCase2B : public SolverCase2Base { + public: + void update_cache(Root1D_State& b, const X_DX& x_dx) const override { + b.update_cache_sign_f0(x_dx); + } -//////////////// -//////////////// -//////////////// + bool is_exit_case_now(Root1D_State& b, const X_DX x_dx) const override { + if (this->is_dx_small(b, x_dx)) { + std::cout << "dx is small --> success" << std::endl; + b.set_flag(CaseFlag::success); return true; } + return b.flag() != flag_this(); + } + CaseFlag flag_this() const override { return CaseFlag::case_2B; } + CaseFlag flag_stall() const override { return CaseFlag::success; } + }; - // template - // std::pair solve(BoundState& state) { +//////////////////////////////////////////////////////////////////////////////////////// + + // Iterates on the root finding problem using the solver associated with the current case flag. + // Returns the result. + template + std::pair solve_for_state(Root1D_State& state, T x_final) { + switch (state.flag()) { + case CaseFlag::case_1: return solve_for_state(state, SolverCase1().solve(state)); + case CaseFlag::case_2B: return solve_for_state(state, SolverCase2B().solve(state)); + case CaseFlag::case_2U: return solve_for_state(state, SolverCase2U().solve(state)); + case CaseFlag::failure: return std::make_pair(false, 0); + case CaseFlag::success: return std::make_pair(true, x_final); + } + return std::make_pair(false, 0); // Unreachable but suppresses compiler warning + } + template + std::pair solve_for_state(Root1D_State& state) { + return solve_for_state(state, state.cache_last().x()); + } - // state.is_finished(); - // } +//////////////////////////////////////////////////////////////////////////////////////// + // Makes three attempts to find a root: + // 1. Local minimization starting with the initial x + // 2. Bracketing the root in the direction of the initial dx if there is a sign flip + // between f(x) and f(x_bound in direction of dx) + // 3. Bracketing the root in the opposite direction of the initial dx if there is a + // sign flip between f(x) and f(x_bound in opposite direction of dx). + // + // If all three attempts fail, an error is thrown. + // + // template - std::pair solve(BoundState& state) { - - if (DEBUG_PRINT) { - std::cout << std::endl; - std::cout << "=========================================" << std::endl; - std::cout << "Count: " << state.count() << std::endl; - std::cout << "STATUS IS: "; - state.print_root_status_flag(); - std::cout << "-----------------------------------------" << std::endl; - std::cout << std::endl; - } + T solve_3_attempts(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { + // Create Root1D_State + Root1D_State state(f, x, xl, xh, digits, max_iter); - if (state.flag() == StatusFlag::case_1) { - auto solver_case_1 = Bound1(); - solver_case_1.solve_next(state); - return solve(state); + // Get CacheX + const auto cache_x_orig = state.cache_last(); - } else if (state.flag() == StatusFlag::case_2B) { - auto solver_case_2b = Bound2B(); - solver_case_2b.solve_next(state); - return solve(state); + auto fn_form_bracket_with_side = [&](const bool ind_bound){ + if (ind_bound) { + // NOTE: These two lines cannot be swapped in case xh is a solution. This is tested. + state.update_cache_l(cache_x_orig); // First to original cache + state.update_cache_h(xh); // Then set bound + } else { + // NOTE: These two lines cannot be swapped in case xl is a solution. This is tested. + state.update_cache_h(cache_x_orig); // First to original cache + state.update_cache_l(xl); // Then set bound + } - } else if (state.flag() == StatusFlag::case_2U) { - auto solver_case_2u = Bound2U(); - solver_case_2u.solve_next(state); - return solve(state); + // state.reset_dx_history(); // Not required for correctness + state.update_case_flag(); // Sets to case_2U, case_2B, or failure + }; - } else if (state.flag() == StatusFlag::case_2I) { - return std::make_pair(false, state.x_next_value()); + auto fn_return = [&](const Root1D_State& s, const std::pair p) { + max_iter = s.num_fn_evals(); // Update max_iter + return p.second; + }; - } else if (state.flag() == StatusFlag::stall) { - return std::make_pair(true, state.x_next_value()); + // Solve problem + const auto p = solve_for_state(state); - } else if (state.flag() == StatusFlag::dx_small) { - return std::make_pair(true, state.x_next_value()); + // If local minimization was successful, return root + if (p.first) { return fn_return(state, p); } - } else if (state.flag() == StatusFlag::f0_is_zero) { - return std::make_pair(true, state.x_next_value()); + // Attempt to bracket root in the direction of the initial guess + fn_form_bracket_with_side(!cache_x_orig.is_bound_h()); + if (state.is_B()) { return fn_return(state, solve_for_state(state)); } - } else if (state.flag() == StatusFlag::out_of_iterations) { - if (DEBUG_PRINT) std::cout << "ERROR OUT OF ITERATIONS" << std::endl; - assert(false); + // Attempt to bracket root in the opposite direction of the initial guess + fn_form_bracket_with_side(cache_x_orig.is_bound_h()); + if (state.is_B()) { return fn_return(state, solve_for_state(state)); } - } else { - - // std::cout << "NOT HANDLED: " << static_cast(state.flag()) << std::endl; - if (DEBUG_PRINT) std::cout << "NOT HANDLED: " << std::endl; - state.print_root_status_flag(); - assert(false); - } - - assert(false); - return std::make_pair(false, state.x_next_value()); - } - - - -///////////////////////// - - template - T solve_2(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { - // Create BoundState - BoundState state(f, x, xl, xh, digits, max_iter); - - - } - -///////////////////////// - - - - template - T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { - if (DEBUG_PRINT) std::cout << "================= newton_raphson_iterate ==================== with max_iter: " << max_iter << std::endl; - - if (DEBUG_PRINT) { - std::cout << "l_orig: " << xl << std::endl; - std::cout << "x_orig: " << x << std::endl; - std::cout << "h_orig: " << xh << std::endl; - } - - - - // Create BoundState -#if 1 - BoundState state(f, x, xl, xh, digits, max_iter); - -#else - BoundState state(f, xl, xh, digits, max_iter); - - // Calc NEXT - const auto next = NEXT(x); - - // Calc CacheX - const auto cache_x_orig = state(next); - - - - - // // Set up problem - // state.update_cache_sign_dx(cache_x_orig); - - // // auto solver_case_1 = Bound1(); - // Bound1().solve(state); -#endif - - - - - - - - - - - -////////////////////////////////////////////////////////////////////////// - - - // Get CacheX - const auto cache_x_orig = state.cache_x(); - - // Solve problem - auto p = solve(state); - - -////////////////////////////////////////////////////////////////////////// - // Define lambdas - auto fn_form_bracket_l = [&](){ - const auto next = NEXT::FromX(xl); - const auto cache_l_orig = state(next); - state.update_cache_l(cache_l_orig); - state.update_cache_h(cache_x_orig); - }; - auto fn_form_bracket_h = [&](){ - const auto next = NEXT::FromX(xh); - const auto cache_h_orig = state(next); - state.update_cache_l(cache_x_orig); - state.update_cache_h(cache_h_orig); - }; - auto fn_return = [&](const BoundState& s, const std::pair p) { - s.print_root_status_flag(); - assert(p.first); - max_iter = s.iters_used(); - return p.second; - }; -////////////////////////////////////////////////////////////////////////// - - // Do the natural - if (p.first) { - if (DEBUG_PRINT) std::cout << "succuss first time" << std::endl; - max_iter = state.iters_used(); - return fn_return(state, p); - } - if (DEBUG_PRINT) std::cout << "FAILURE first time" << std::endl << std::endl; - - - cache_x_orig.is_l() ? fn_form_bracket_h() : fn_form_bracket_l(); - state.set_state_2(); - if (state.is_B()) { - if (DEBUG_PRINT) std::cout << "BRACKETED PRIORITY 1" << std::endl; - p = solve(state); - return fn_return(state, p); - } else { - if (DEBUG_PRINT) std::cout << "COULD NOT BRACKET PRIORITY 1" << std::endl << std::endl; - } - - cache_x_orig.is_l() ? fn_form_bracket_l() : fn_form_bracket_h(); - state.set_state_2(); - if (state.is_B()) { - if (DEBUG_PRINT) std::cout << "BRACKETED PRIORITY 2" << std::endl; - p = solve(state); - return fn_return(state, p); - } else { - if (DEBUG_PRINT) std::cout << "COULD NOT BRACKET PRIORITY 2" << std::endl << std::endl; - } - - static const char* function = "boost::math::tools::detail::solve<%3%>"; - const T x_last = state.x_next_value(); - return policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::solve, last closest guess was %1%", x_last, boost::math::policies::policy<>()); - }; - -// from compile_test/gauss_kronrod_concept_test.cpp:11: -// ../../../libs/math/test/test_ibeta_inv.hpp(64): last checkpoint -// ../../../bin.v2/libs/math/test/test_bessel_airy_zeros.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/test_bessel_airy_zeros.o... -// test_inverse_gaussian -// cat "../../../bin.v2/libs/math/test/test_ibeta_inv_double.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/test_ibeta_inv_double.output" - -// test_hypergeometric_laplace_transform - - - -// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/test_vector_barycentric_rational.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/test_vector_barycentric_rational.o... -// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/cstdfloat_concept_check_1.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/compile_test/cstdfloat_concept_check_1.o... -// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/std_real_concept_check.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/compile_test/std_real_concept_check.o... -// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/std_real_concept_check_80.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/std_real_concept_check_80.o... -// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/std_real_concept_check_32.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/std_real_concept_check_32.o... -// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/std_real_concept_check_128.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/std_real_concept_check_128.o... -// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/cstdfloat_concept_check_2.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/compile_test/cstdfloat_concept_check_2.o... -// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/cstdfloat_concept_check_3.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/compile_test/cstdfloat_concept_check_3.o... -// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/cstdfloat_concept_check_4.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/compile_test/cstdfloat_concept_check_4.o... -// ...failed gcc.compile.c++ ../../../bin.v2/libs/math/test/std_real_concept_check_64.test/gcc-9/debug/link-static/threading-multi/visibility-hidden/std_real_concept_check_64.o... - - - -// ../../../b2 --report_level=detailed --show_progress=false | grep -E "test case|failed" - - - - -// ### ??? -// -// ### CONFIRMED TO BE REAL TESTS ### -// test_legendre -// test_ibeta_inv_double -// test_laguerre -// test_bessel_airy_zeros -// test_hermite -// -// testing.capture-output ../../../bin.v2/libs/math/example/legendre_stieltjes_example.test/gcc-9/debug/threading-multi/visibility-hidden/legendre_stieltjes_example.run - - - - - - -// Start here -#if 0 -//////////////////////////////////////////////////////////////////////////////////////// - - // Forward declare ContextF - template - class ContextF; - - - // - // - // - // - // - template - class BoundsState { - public: - virtual ~BoundsState() = default; - - void set_b(BoundState& b) { b_ = &b; } - - void do_step_eval_then_x_next(ContextF& helper_f) { - // Attempt Newton step - const auto cx = cache_x(); - const T dx_newton = cx.calc_dx(); - const T x_next = cx.x() + dx_newton; - - // Need to do bisection - if (!is_inbounds(x_next) || is_converge_too_slow(dx_newton)) { - calc_next_bisection(); - return; - } - - // Newton's step - update_cache(helper_f(x_next)); - return; - } - - bool is_converge_too_slow(const T& dx) const { - return fabs(b_->dx_pp()) < fabs(dx * 4); - } - - virtual const CacheOrder1& cache_x() const = 0; - virtual NEXT calc_next_bisection() = 0; - virtual void update_cache(const CacheOrder1& cache_x) = 0; - - bool is_inbounds(T x) const { return (l_smart() <= x) && (x <= h_smart()); } - virtual T l_smart() const = 0; - virtual T h_smart() const = 0; - - - protected: - ContextF* b_; - }; - - - - // - template - class ContextF { - public: - ContextF(F f, int digits, std::uintmax_t max_iter) - : f_(f) - , factor_(static_cast(ldexp(1.0, 1 - digits))) - , count_(max_iter) - , max_iter_(max_iter) {} - - ~ContextF() { - delete bounds_; - } - - CacheOrder1 operator()(T x) { - if (0 == count_) { - static const char* function = "boost::math::tools::detail::ContextF<%1%>"; - policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::ContextF, last closest guess was %1%", x, boost::math::policies::policy<>()); - } - --count_; - return CacheOrder1(f_, x); - } - - void transition_to_state(BoundsState* bounds) { - if (this->bounds_ != nullptr) - delete this->bounds_; - this->bounds_ = bounds; - this->bounds_->set_b(*this); - } - - SolveOutput solve(BoundsState* bounds) { - transition_to_state(bounds); - return solve(); - } - - T dx_pp() const { return dx_history_.dx_pp(); } - int iters_used() const { return max_iter_ - count_; } - - private: - SolveOutput solve() { - BOOST_MATH_STD_USING - std::cout << std::scientific << std::setprecision(6) << std::endl; - - while (count_) { - bounds_->do_step_eval_then_x_next(*this); - } - assert(false); - return std::make_pair(OutputResult::failure, 0); - } - -#if 0 - std::pair solve(Bounds& bounds, CacheOrder1 cache_x) { - BOOST_MATH_STD_USING - - std::cout << std::scientific << std::setprecision(6) << std::endl; - - T dx_true; - bool is_prev_newton = false; - - T x_next; - while (count_) { - // Calculate Newton step - const T dx_newton = cache_x.calc_dx(); - x_next = cache_x.x() + dx_newton; // Checked against bracket bounds below - - - const bool c1 = !bounds.is_inbounds(x_next); - if (c1 || fabs(dx_pp()) < fabs(dx_newton * 4)) { // Bisection halves step size each iteration - const auto p_bisect = bounds.handle_bisection(*this, cache_x); - if (bounds.is_case_2I()) { std::cout << "is_case_2I" << std::endl; return std::make_pair(false, 0); } - - const auto cache_x_next = p_bisect.second; - std::cout << "STEP: BISECTION: " << " -- wanted: " << x_next << " -- got: " << cache_x_next.x() << " -- bounds: " << c1 << " -- converge: " << (fabs(dx_pp()) < fabs(dx_newton * 4)) << std::endl; - x_next = cache_x_next.x(); - dx_history_.reset(); // Invalidate step size history - - dx_true = x_next - cache_x.x(); - - cache_x = cache_x_next; - is_prev_newton = false; - } else { - std::cout << "STEP: NEWTON" << std::endl; - dx_true = x_next - cache_x.x(); - cache_x = (*this)(x_next); // Update cache_x_ with x_next - is_prev_newton = true; - } - - std::cout << "AFTER STEP-- x_next: " << x_next << ", dx_true: " << dx_true << std::endl; - bounds.print(); - - assert(bounds.is_inbounds(x_next)); - - // std::cout << " UPDATED cache_x: " << std::endl; - // std::cout << " "; - // cache_x.print(); - // std::cout << " count: " << count_ << std::endl; - - // Determine if successful - const bool is_dx_small = fabs(cache_x.calc_dx()) < fabs(cache_x.x() * factor_); - if (is_dx_small) { - if (bounds.is_case_2B() || bounds.is_case_1()) { - std::cout << "is_case_2B && is_dx_small" << std::endl; - std::cout << "dx_true: " << fabs(dx_true) << ", cache_x.x(): " << fabs(cache_x.x()) << ", factor_: " << factor_ << std::endl; - break; // Tolerance met - } else { - bounds.print_case(); - std::cout << "dx_true: " << fabs(dx_true) << ", cache_x.x(): " << fabs(cache_x.x()) << ", factor_: " << factor_ << std::endl; - std::cout << "FAILURE dx_small" << std::endl; - return std::make_pair(false, 0); // Did not bracket root - } - } - - if (cache_x.f0() == 0) { - std::cout << "SUCCESS f0 = 0" << std::endl; - break; // Function is zero - } - - const bool is_numerical = (x_next == bounds.l_smart() || x_next == bounds.h_smart() || dx_true == 0); - if (is_numerical) { - if (bounds.is_case_2B()) { - std::cout << "is_case_2B && is_numerical" << std::endl; - break; // Tolerance met - } else { - std::cout << "x_next == bounds.l_smart(): " << (x_next == bounds.l_smart()) << std::endl; - std::cout << "x_next == bounds.h_smart(): " << (x_next == bounds.h_smart()) << std::endl; - std::cout << "dx_true == 0: " << (dx_true == 0) << std::endl; - std::cout << "FAILURE numerical" << std::endl; - return std::make_pair(false, 0); // Did not bracket root - } - } - - - if (bounds.is_case_2I()) { - std::cout << "FAILURE is_case_2I" << std::endl; - return std::make_pair(false, 0); - } - - bounds.update_state(cache_x); // Update bracket state - - if (bounds.is_case_2I()) { - std::cout << "FAILURE is_case_2I" << std::endl; - return std::make_pair(false, 0); - } - - std::cout << std::endl; - dx_history_.push(dx_newton); // Store step size - } - - return std::make_pair(true, cache_x.x()); - } -#endif - - private: - F f_; - T factor_; - std::uintmax_t count_; - std::uintmax_t max_iter_; - StepSizeHistory dx_history_; - BoundsState* bounds_; - }; - - - - template - class Bounds2; // : public BoundsState; - - template - class Bounds2B; // : public BoundsState; - - template - class Bounds2U; // : public BoundsState; - - template - class Bounds2I; // : public BoundsState; - - // - // - // - // - // - template - class Bounds2 : public BoundsState { - public: - static Bounds2* Identify(CacheOrder1 cache_a, CacheOrder1 cache_b) { - const auto cache_sorted = cache_a.sort(cache_b); - const auto cache_l = cache_sorted.first; - const auto cache_h = cache_sorted.second; - - if (cache_l.is_B(cache_h)) { - return new Bounds2B(cache_l, cache_h); - } else if (cache_l.is_U(cache_h)) { - return new Bounds2U(cache_l, cache_h); - } else { - return new Bounds2I(cache_l, cache_h); - } - } - - NEXT calc_next_bisection() override { - const auto midpoint = (cache_l_.x() + cache_h_.x()) / 2; - const auto cache_midpoint = helper_f(midpoint); - update_cache(cache_midpoint); - } - - void update_cache_sign_f0(const CacheOrder1& cache_x) { - if (cache_x.is_same_sign(cache_l_)) { - cache_l_ = cache_x; - } else { - cache_h_ = cache_x; - } - } - - void update_cache_sign_dx(const CacheOrder1& cache_x) { - if (cache_x.is_same_slope(cache_l_)) { // TODO: make better - cache_l_ = cache_x; - } else { - cache_h_ = cache_x; - } - } - - T l_smart() const override { return cache_l_.x(); } - T h_smart() const override { return cache_h_.x(); } - - const detail::CacheOrder1& cache_x() const override { - return (is_last_eval_l_) ? cache_l_ : cache_h_; - } - - bool is_B() const { return cache_l_.is_B(cache_h_); } - bool is_U() const { return cache_l_.is_U(cache_h_); } - bool is_I() const { return cache_l_.is_I(cache_h_); } - - const CacheOrder1& cache_l() { return cache_l_; } - const CacheOrder1& cache_h() { return cache_h_; } - - private: - bool is_last_eval_l_; - CacheOrder1 cache_l_; - CacheOrder1 cache_h_; - }; - - // - // - // - // - // - template - class Bounds2B : public Bounds2 { - public: - Bounds2B(const CacheOrder1& cache_l, const CacheOrder1& cache_h) - : Bounds2(cache_l, cache_h) - { - assert(this->is_B()); - } - - void update_cache(const CacheOrder1& cache_x) override { - update_cache_sign_f0(cache_x); - } - }; - - // - // - // - // - // - template - class Bounds2U : public Bounds2 { - public: - Bounds2U(const CacheOrder1& cache_l, const CacheOrder1& cache_h) - : Bounds2(cache_l, cache_h) - { - assert(this->is_U()); - } - - void update_cache(const CacheOrder1& cache_x) override { - update_cache_sign_dx(cache_x); - if (this->is_B()) { - this->b_->transition_to_state( - new Bounds2B(this->cache_l_, this->cache_h_) - ); - } - } - }; - - - // - // - // - // - // - template - class Bounds2I : public Bounds2 { - public: - Bounds2I(const CacheOrder1& cache_l, const CacheOrder1& cache_h) - : Bounds2(cache_l, cache_h) - { - assert(this->is_U()); - } - - void update_cache(const CacheOrder1& cache_x) override { - assert(false); - } - }; - - - // - // - // - // - // - template - class Bounds1 : public BoundsState { - public: - Bounds1(T l, const detail::CacheOrder1& cache_x, T h) - : cache_x_(cache_x) - { - limit_ = cache_x_.is_l() ? h : l; - } - - NEXT calc_next_bisection() override { - const auto cache_limit = (*this->b_)(limit_); - auto* bounds_new = Bounds2::Identify(cache_limit, cache_x_); - this->b_->transition_to_state(bounds_new); - } - - T l_smart() const override { return cache_x_.is_l() ? cache_x_.x() : limit_; } - T h_smart() const override { return cache_x_.is_h() ? cache_x_.x() : limit_; } - - void update_cache(const CacheOrder1& cache_x) override { cache_x_ = cache_x; } - - const CacheOrder1& cache_x() const override { return cache_x_; } - - private: - CacheOrder1 cache_x_; - T limit_; - }; - - - -/////////////////////////////////////////////////////////////////////// - - - template - T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { - // Create ContextF - ContextF b(f, digits, max_iter); - - // Create cache_x - const auto cache_x = b(x); - if (0 == cache_x.f0()) return x; // Check if root is already found - - // Create Bounds - auto bounds = Bounds1(xl, cache_x, xh); - - // Solve - const auto p = b.solve(&bounds); - - return p.second; - -#if 0 - if (p.first) { - std::cout << "succuss first time" << std::endl; - max_iter = b.iters_used(); - return p.second; - } else { - std::cout << "FAILURE first time" << std::endl << std::endl; - - bounds = Bounds(xl, cache_x, xh); - cache_x.is_l() ? bounds.eval_h_orig(b) : bounds.eval_l_orig(b); - if (bounds.is_case_2B()) { - std::cout << "BRACKETED PRIORITY 1" << std::endl; - return b.solve(bounds, cache_x).second; - } else { - std::cout << "COULD NOT BRACKET PRIORITY 1" << std::endl << std::endl; - } - - bounds = Bounds(xl, cache_x, xh); - cache_x.is_l() ? bounds.eval_l_orig(b) : bounds.eval_h_orig(b); - if (bounds.is_case_2B()) { - std::cout << "BRACKETED PRIORITY 2" << std::endl; - return b.solve(bounds, cache_x).second; - } else { - std::cout << "l: " << bounds.l_smart() << ", h: " << bounds.h_smart() << std::endl; - std::cout << "f0_l: " << bounds.cache_l().f0() << ", f0_h: " << bounds.cache_h().f0() << std::endl; - std::cout << "COULD NOT BRACKET PRIORITY 2" << std::endl << std::endl; - } - } -#endif - - static const char* function = "boost::math::tools::detail::ContextF<%1%>"; - policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::ContextF, last closest guess was %1%", x, boost::math::policies::policy<>()); - return 0; - }; - - - - - - -#if 0 - - - // - // - // - // - // - template - class Bounds2 { - public: - Bounds2(T l, const detail::CacheOrder1& cache_x, T h) - : var_l_(LimitOrig(l)), var_h_(LimitOrig(h)) - { - if (cache_x.is_l()) { - is_last_eval_l_ = true; - var_l_ = cache_x; - } else { - is_last_eval_l_ = false; - var_h_ = cache_x; - } - } - - template - void do_step_eval_then_x_next(H& helper_f) { - const auto cx = cache_x(); - const T x_next = cx.x() + cx.calc_dx(); - - if (!is_inbounds(x_next) || is_converge_too_slow()) { - return calc_next_bisection(helper_f); - } - - // Newton's step - update_cache(helper_f(x_next)); - } - - template - NEXT calc_next_bisection(H& helper_f) { - if (is_case_1()) { - materialize(helper_f); - return; - } - const T x_mid = (x_l() + x_h()) / 2; - update_cache_2(helper_f(x_mid)); - return; - } - - void update_cache(const detail::CacheOrder1& cache_x) { - if (is_case_1()) { - update_cache_1(cache_x); - } else { - update_cache_2(cache_x); - } - } - - void update_cache_1(const detail::CacheOrder1& cache_x) { - update_cache_using_dx(cache_x); - if (size() == 2) { - determine_state_2(); - } - } - - void update_cache_2(const detail::CacheOrder1& cache_x) { - if (is_case_2B()) { - update_cache_using_sign(cache_x); - } else { - assert(is_case_2U()); - update_cache_2U(cache_x); - } - } - - void update_cache_2U(const detail::CacheOrder1& cache_x) { - const T dx_l_orig = cache_l().calc_dx(); - const T dx_h_orig = cache_h().calc_dx(); - - update_cache_using_dx(cache_x); - - if (cache_l().is_B(cache_h())) { // Now bracket root! - state_ = BracketState::case_2B; - } else { - const T dx_l = cache_l().calc_dx(); - const T dx_h = cache_h().calc_dx(); - - if (dx_l_orig < dx_l || dx_h_orig < dx_h) { - state_ = BracketState::case_2I; - } - } - } - - template - void materialize(H& helper_f) { - // If var_l_ holds a LimitOrig, then we need to evaluate it - if (var_l_.is_limit_orig()) { - set_l(helper_f(var_l_.x())); - } - // If var_h_ holds a LimitOrig, then we need to evaluate it - if (var_h_.is_limit_orig()) { - set_h(helper_f(var_h_.x())); - } - determine_state_2(); - } - - void determine_state_2() { - assert(size() == 2); - if (cache_l().is_B(cache_h())) { - state_ = BracketState::case_2B; - } else if (cache_l().is_U(cache_h())) { - state_ = BracketState::case_2U; - } else { - state_ = BracketState::case_2I; - } - } - - void update_state_using_dx(const CacheOrder1& cache_x) { - if (cache_x.is_l()) { - var_l_ = cache_x; - } else { - var_h_ = cache_x; - } - } - - void update_state_using_sign(const CacheOrder1& cache_x) { - if (has_l()) { - if (cache_x.is_same_sign(cache_l())) { - var_l_ = cache_x; - } else { - var_h_ = cache_x; - } - } else { - if (cache_x.is_same_sign(cache_h())) { - var_h_ = cache_x; - } else { - var_l_ = cache_x; - } - } - } - -//////////////////////////////////////////////////////////////////////////////// - - detail::CacheOrder1 cache_x() const { return is_last_eval_l_ ? val_l_ : var_h_; } - - int size() const { return var_l_.is_cache() + var_h_.is_cache(); } - - bool is_inbounds(T x) const { return (l() <= x) && (x <= h()); } - - void set_l(const detail::CacheOrder1& cache_x) { - is_last_eval_l_ = true; - var_l_ = cache_x; - } - void set_h(const detail::CacheOrder1& cache_x) { - is_last_eval_l_ = false; - var_h_ = cache_x; - } - - T l() const { return var_l_.x(); } - T h() const { return var_h_.x(); } - - private: - BracketState state_; - bool is_last_eval_l_; - booth::variant, detail::CacheOrder1> var_l_; - booth::variant, detail::CacheOrder1> var_h_; - StepSizeHistory dx_history_; - }; - - - - - - - - - - - - - - - - - - - // - template - class Bounds { - public: - Bounds(T l, detail::CacheOrder1 cache_x, T h) - : state_(BracketState::case_1) - , l_orig_(l) - , h_orig_(h) - { - update_state_using_dx(cache_x); - } - - void update_state_using_dx(const CacheOrder1& cache_x) { - if (cache_x.is_l()) { - opt_cache_l_ = cache_x; - } else { - opt_cache_h_ = cache_x; - } - } - - void update_state_using_sign(const CacheOrder1& cache_x) { - if (has_l()) { - if (cache_x.is_same_sign(cache_l())) { - opt_cache_l_ = cache_x; - } else { - opt_cache_h_ = cache_x; - } - } else { - if (cache_x.is_same_sign(cache_h())) { - opt_cache_h_ = cache_x; - } else { - opt_cache_l_ = cache_x; - } - } - } - - void determine_state_exhaustive() { - if (cache_l().is_B(cache_h())) { - state_ = BracketState::case_2B; - } else if (cache_l().is_U(cache_h())) { - state_ = BracketState::case_2U; - } else { - state_ = BracketState::case_2I; - } - } - - bool update_state(const CacheOrder1& cache_x) { - if (is_case_1()) { - const bool c1 = cache_x.is_l() && has_h(); - const bool c2 = cache_x.is_h() && has_l(); - update_state_using_dx(cache_x); - if (c1 || c2) { - determine_state_exhaustive(); - } - } else if (is_case_2B()) { - update_state_using_sign(cache_x); - if (!cache_l().is_B(cache_h())) { - std::cout << std::endl; - std::cout << "SOMETHING IS WRONG WRONG WRONG" << std::endl; - print(); - cache_x.print(); - assert(false); - } - } else { - assert(is_case_2U()); - update_state_2U(cache_x); - } - return true; - } - - void update_state_2U(const CacheOrder1& cache_x) { - const T dx_l_orig = cache_l().calc_dx(); - const T dx_h_orig = cache_h().calc_dx(); - - update_state_using_dx(cache_x); - - if (cache_l().is_B(cache_h())) { // Now bracket root! - state_ = BracketState::case_2B; - } else { - const T dx_l = cache_l().calc_dx(); - const T dx_h = cache_h().calc_dx(); - - if (dx_l_orig < dx_l || dx_h_orig < dx_h) { - state_ = BracketState::case_2I; - } - } - } - - static T midpoint_smart(T l, T h) { - BOOST_MATH_STD_USING - - const T t_min = std::numeric_limits::min(); - const T t_max = std::numeric_limits::max(); - const T t_lowest = std::numeric_limits::lowest(); - const T t_inf = std::numeric_limits::infinity(); - - const T l_abs = fabs(l); - const T h_abs = fabs(h); - if (std::max(l_abs, h_abs) < 100 || // Numbers are small - (h - l) < std::min(l_abs, h_abs)) { // Numbers are close - return (l + h) / 2; - } - - // Want both numbers to have the same sign - if (sign(l) * sign(h) == -1) { return 0; } - - // If l is -Inf make it the lowest float - if (l == -t_inf) { - l = t_lowest; - } else if (h == t_inf) { - h = t_max; - } - - // If l is 0 make it the smallest float - if (l == 0) { - l = t_min; - } else if (h == 0) { - h = -t_min, l; - } - - return copysign(sqrt(l * h), l); - } - - template - std::pair> handle_bisection(H& helper_f, const CacheOrder1& cache_x) { - if (is_case_1()) { - const T x_next = cache_x.is_l() ? h_orig_ : l_orig_; - - CacheOrder1 cache_x_next; - - if (cache_x.is_l()) { - cache_x_next = helper_f(x_next); - opt_cache_h_ = cache_x_next; - } else { - cache_x_next = helper_f(x_next); - opt_cache_l_ = cache_x_next; - } - - determine_state_exhaustive(); - print_case(); - } - - // const T x_mid = (l_smart() + h_smart()) / 2; - const T x_mid = midpoint_smart(l_smart(), h_smart()); - - std::cout << "x_mid: " << x_mid << std::endl; - - const auto cache_x_next = helper_f(x_mid); - - std::cout << "post eval" << std::endl; - - return std::make_pair(true, cache_x_next); - } - - template - void eval_l_orig(H& helper_f) { - const auto cache_x = helper_f(l_orig_); - push_front(cache_x); - determine_state_exhaustive(); - } - template - void eval_h_orig(H& helper_f) { - const auto cache_x = helper_f(h_orig_); - push_back(cache_x); - determine_state_exhaustive(); - } - - void push_front(const CacheOrder1& cache_x) { - assert(size() == 1); - if (has_l()) { opt_cache_h_ = opt_cache_l_; } - opt_cache_l_ = cache_x; - } - void push_back(const CacheOrder1& cache_x) { - assert(size() == 1); - if (has_h()) { opt_cache_l_ = opt_cache_h_; } - opt_cache_h_ = cache_x; - } - - int size() const { return has_l() + has_h(); } - bool has_l() const { return opt_cache_l_.has_data(); } - bool has_h() const { return opt_cache_h_.has_data(); } - - const CacheOrder1& cache_l() const { return opt_cache_l_.c(); } - const CacheOrder1& cache_h() const { return opt_cache_h_.c(); } - - T l_smart() const { return has_l() ? cache_l().x() : l_orig_; } - T h_smart() const { return has_h() ? cache_h().x() : h_orig_; } - - bool is_inbounds(T x) const { return (l_smart() <= x) && (x <= h_smart()); } - - bool is_case_1() const { return state_ == BracketState::case_1; } - bool is_case_2U() const { return state_ == BracketState::case_2U; } - bool is_case_2B() const { return state_ == BracketState::case_2B; } - bool is_case_2I() const { return state_ == BracketState::case_2I; } - - std::string get_case_string() const { - if (is_case_1()) { return "case_1"; } - if (is_case_2U()) { return "case_2U"; } - if (is_case_2B()) { return "case_2B"; } - assert(is_case_2I()); - return "case_2I"; - } - - void print_case() const { std::cout << get_case_string() << std::endl; } - - void print() const { - std::cout << "BOUNDS: " << get_case_string() << std::endl; - std::cout << " l_orig: " << l_orig_ << std::endl; - if (has_l()) { cache_l().print(); } - if (has_h()) { cache_h().print(); } - std::cout << " h_orig: " << h_orig_ << std::endl; - } - - private: - BracketState state_; - T l_orig_; - detail::Optional> opt_cache_l_; - detail::Optional> opt_cache_h_; - T h_orig_; - }; - - template - T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { - // Create ContextF - ContextF helper_f(f, digits, max_iter); - - // Create cache_x - const auto cache_x = helper_f(x); - if (0 == cache_x.f0()) return x; // Check if root is already found - - // Create Bounds2 - auto bounds = Bounds1(xl, cache_x, xh); - - // Do Newton loop - const auto p = helper_f.solve(bounds, cache_x); - - if (p.first) { - std::cout << "succuss first time" << std::endl; - max_iter = helper_f.iters_used(); - return p.second; - } else { - std::cout << "FAILURE first time" << std::endl << std::endl; - - bounds = Bounds2(xl, cache_x, xh); - cache_x.is_l() ? bounds.eval_h_orig(helper_f) : bounds.eval_l_orig(helper_f); - if (bounds.is_case_2B()) { - std::cout << "BRACKETED PRIORITY 1" << std::endl; - return helper_f.solve(bounds, cache_x).second; - } else { - std::cout << "COULD NOT BRACKET PRIORITY 1" << std::endl << std::endl; - } - - bounds = Bounds2(xl, cache_x, xh); - cache_x.is_l() ? bounds.eval_l_orig(helper_f) : bounds.eval_h_orig(helper_f); - if (bounds.is_case_2B()) { - std::cout << "BRACKETED PRIORITY 2" << std::endl; - return helper_f.solve(bounds, cache_x).second; - } else { - std::cout << "l: " << bounds.l_smart() << ", h: " << bounds.h_smart() << std::endl; - std::cout << "f0_l: " << bounds.cache_l().f0() << ", f0_h: " << bounds.cache_h().f0() << std::endl; - std::cout << "COULD NOT BRACKET PRIORITY 2" << std::endl << std::endl; - } - } - - static const char* function = "boost::math::tools::detail::ContextF<%1%>"; - policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::ContextF, last closest guess was %1%", x, boost::math::policies::policy<>()); - return 0; - }; - -#if 0 - template - T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { - // Create ContextF - ContextF helper_f(f, digits, max_iter); - - // Create cache_x - const auto cache_x = helper_f(x); - if (0 == cache_x.f0()) return x; // Check if root is already found - - // Create Bounds - auto bounds = Bounds(xl, cache_x, xh); - - // Do Newton loop - const auto p = helper_f.solve(bounds, cache_x); - - if (p.first) { - std::cout << "succuss first time" << std::endl; - max_iter = helper_f.iters_used(); - return p.second; - } else { - std::cout << "FAILURE first time" << std::endl << std::endl; - - bounds = Bounds(xl, cache_x, xh); - cache_x.is_l() ? bounds.eval_h_orig(helper_f) : bounds.eval_l_orig(helper_f); - if (bounds.is_case_2B()) { - std::cout << "BRACKETED PRIORITY 1" << std::endl; - return helper_f.solve(bounds, cache_x).second; - } else { - std::cout << "COULD NOT BRACKET PRIORITY 1" << std::endl << std::endl; - } - - bounds = Bounds(xl, cache_x, xh); - cache_x.is_l() ? bounds.eval_l_orig(helper_f) : bounds.eval_h_orig(helper_f); - if (bounds.is_case_2B()) { - std::cout << "BRACKETED PRIORITY 2" << std::endl; - return helper_f.solve(bounds, cache_x).second; - } else { - std::cout << "l: " << bounds.l_smart() << ", h: " << bounds.h_smart() << std::endl; - std::cout << "f0_l: " << bounds.cache_l().f0() << ", f0_h: " << bounds.cache_h().f0() << std::endl; - std::cout << "COULD NOT BRACKET PRIORITY 2" << std::endl << std::endl; - } - } - - static const char* function = "boost::math::tools::detail::ContextF<%1%>"; - policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::ContextF, last closest guess was %1%", x, boost::math::policies::policy<>()); - return 0; - }; -#endif -#endif -#endif - -} // namespace detail - - - - -template -T newton_raphson_iterate(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) -{ - return detail::solve>(f, x, xl, xh, digits, max_iter); -} - -template -inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) -{ - std::uintmax_t m = (std::numeric_limits::max)(); - return newton_raphson_iterate(f, guess, min, max, digits, m); -} - - - -#if 1 - -template -T newton_raphson_iterate_classic(F f, T guess, T min, T max, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) -{ - BOOST_MATH_STD_USING - - static const char* function = "boost::math::tools::newton_raphson_iterate_classic<%1%>"; - if (min > max) - { - return policies::raise_evaluation_error(function, "Range arguments in wrong order in boost::math::tools::newton_raphson_iterate_classic(first arg=%1%)", min, boost::math::policies::policy<>()); - } - - T f0(0), f1, last_f0(0); - T result = guess; - - T factor = static_cast(ldexp(1.0, 1 - digits)); - T delta = tools::max_value(); - T delta1 = tools::max_value(); - T delta2 = tools::max_value(); - - // - // 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); - -#ifdef BOOST_MATH_INSTRUMENT - std::cout << "newton_raphson_iterate_classic, guess = " << guess << ", min = " << min << ", max = " << max - << ", digits = " << digits << ", max_iter = " << max_iter << "\n"; -#endif - - do { - last_f0 = f0; - delta2 = delta1; - delta1 = delta; - std::cout << std::endl; - std::cout << "EVAL IS: " << count << ", at: " << result << std::endl; - detail::unpack_tuple(f(result), f0, f1); - --count; - - std::cout << "result: " << result << ", f0: " << f0 << ", f1: " << f1; - - 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; - std::cout << ", delta: " << delta << std::endl; - } -#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) * 1.1f; // Protect against huge jumps! - //delta = sign(delta) * result; // Protect against huge jumps! Failed for negative result. https://github.com/boostorg/math/issues/216 - } - 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_classic, perhaps we have a local minima near current best guess of %1%", guess, boost::math::policies::policy<>()); - } - }while(count && (fabs(result * factor) < fabs(delta))); - - std::cout << "delta: " << delta << std::endl; - std::cout << "fabs(result * factor)" << fabs(result * factor) << std::endl; - std::cout << "fabs(delta) " << fabs(delta) << std::endl; - - max_iter -= count; - -#ifdef BOOST_MATH_INSTRUMENT - std::cout << "Newton Raphson required " << max_iter << " iterations\n"; -#endif - - return result; -} - -template -inline T newton_raphson_iterate_classic(F f, T guess, T min, T max, int digits) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) -{ - std::uintmax_t m = (std::numeric_limits::max)(); - return newton_raphson_iterate_classic(f, guess, min, max, digits, m); -} - - -#endif - - - - - - - - - - - - - -#else - - - -///////////////////////////////// - - - // Enforces bracketing when root finding with Newton's method - template - class BracketHelper { - // Stores the last two Newton step sizes - class StepSizeHistory { - public: - StepSizeHistory() { reset(); } - - void reset() { - dx_p_ = (std::numeric_limits::max)(); - dx_pp_ = (std::numeric_limits::max)(); - } - - void push(T input) { - dx_pp_ = dx_p_; - dx_p_ = input; - } - - T dx_p() const { return dx_p_; } - T dx_pp() const { return dx_pp_; } - - private: - T dx_p_; - T dx_pp_; - }; - - // - - - - - static Bounds create_bracket_candidate(ContextF& helper_f, const detail::CacheOrder1& cache_x, T bound) { - if (is_can_eval_bound(bound)) { - return Bounds(cache_x, helper_f(bound)); - } else { - return rollout_for_bracket_candidate(helper_f, cache_x, bound); - } - } - - static - - public: - static T solve(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { - // // Create ContextF - // ContextF helper_f(f, digits, max_iter); - - // // Create cache_x - // const auto cache_x = helper_f(x); - - // // Check if root is already found - // if (0 == cache_x.f0()) return x; - - // // Get preferred bound - // const T dx = Step().calc_dx(cache_x); - // const T bound_pref = (0 <= dx) ? xl : xh; - - // // Check if bound is valid - // const auto bc_pref = create_bracket_candidate(helper_f, cache_x, bound_pref); - - - - - const bool is_can_eval_bound = is_can_eval_bound(bound); - if (is_can_eval_bound) { - const auto cache_b = detail::CacheOrder1(f, bound); - const auto bc = Bounds(cache_x, cache_b); - if (bc.may_contain_root()) { - return solve_bc(f, bc, digits, max_iter, cache_x); - } else { - return std::make_pair(false, 0); - } - } - - - if (0 <= dx) { // Want to go lower - const auto p_l = solve_l(f, x, xl, xh, digits, max_iter, cache_x); - } else { - - } - } - - static std::pain solve_bound_prefered(F f, BranchCandidate bc, int digits, std::uintmax_t& max_iter) { - if (bc.can_bracket()) { - return solve_bracket(f, bc, digits, max_iter, cache_x); - } else if (bc.can_brent()) { - return solve_brent(f, bc, digits, max_iter, cache_x); - } else { - return std::make_pair(false, 0); - } - } - - static std::pair solve_l(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter, const detail::CacheOrder1& cache_x) { - const bool is_can_eval_bound(xl); - if (is_can_eval_bound) { - const auto bc = BranchCandidate(cache_x, detail::CacheOrder1(f, xl)); - if (bc.may_contain_root()) { - return solve_bc(f, x, xl, xh, digits, max_iter, cache_x); - } else { - return std::make_pair(false, 0); - } - } - - if (!is_can_eval_bound) { - // Find lower bound - } - - - } - - - - - // std::cout << "================= solve ====================" << std::endl; - - // const auto cache_debug = detail::CacheOrder1(f, x); - // std::cout << "x: " << cache_debug.x() << ", f0: " << cache_debug.f0() << ", f1: " << cache_debug.f1() << std::endl; - // // x: 0.56763, f0: -0.00218081, f1: 0.0721978 - - // const auto cache_debug_l = detail::CacheOrder1(f, xl); - // std::cout << "xl: " << cache_debug_l.x() << ", f0: " << cache_debug_l.f0() << ", f1: " << cache_debug_l.f1() << std::endl; - // // xl: 0, f0: -0.999001, f1: 0 - - // const auto cache_debug_y = detail::CacheOrder1(f, 1.0); - // std::cout << "y: " << cache_debug_y.x() << ", f0: " << cache_debug_y.f0() << ", f1: " << cache_debug_y.f1() << std::endl; - // // y: 1, f0: 0.000998974, f1: 1.64892e-07 - - // const auto cache_debug_z = detail::CacheOrder1(f, 2.0); - // std::cout << "z: " << cache_debug_z.x() << ", f0: " << cache_debug_z.f0() << ", f1: " << cache_debug_z.f1() << std::endl; - // // z: 2, f0: 0.000998974,flip f1: 2.88776e-33 - - - // std::cout << "xh: " << xh << std::endl; - // // xh: 3.40282e+38 - - - // Get maximum value of T - // const T x_limit = std::numeric_limits::max(); - - // std::cout << "x_limit: " << x_limit << std::endl; - // const bool bb = x_limit == xh; - // std::cout << bb << std::endl; - - // assert(false); - - const auto cache_debug_h = detail::CacheOrder1(f, xh); - std::cout << "xh: " << cache_debug_h.x() << ", f0: " << cache_debug_h.f0() << ", f1: " << cache_debug_h.f1() << std::endl; - // Did not finish - - - BracketHelper helper(f, x, xl, xh, digits, max_iter); - - T x_sol = helper.solve_impl(); - max_iter -= helper.count(); // Make max_iter the number of iterations used. - std::cout << "max_iter: " << max_iter << std::endl; - return x_sol; - } - - private: - // Same as the inputs to newton_raphson_iterate - BracketHelper(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) - : f_(f) - , cache_x_(f, x) - , cache_xl_(f, xl) - , cache_xh_(f, xh) - , count_(max_iter) - , factor_(static_cast(ldexp(1.0, 1 - digits))) // factor_ = 1.0 * 2^(1-digits) - { - if (xh < xl) { - static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; - policies::raise_evaluation_error(function, "Range arguments in wrong order in boost::math::tools::detail::BracketHelper(first arg=%1%)", xl, boost::math::policies::policy<>()); - } - - std::cout << "================= factor_: " << factor_ << "====================" << std::endl; - } - - - T calc_dx_next_default() const { return Step().calc_dx(cache_x_); } - T calc_x_next_default() const { return x() - calc_dx_next_default(); } - - bool can_bracket_l() const { return f0() * f0_l() <= 0; } - bool can_bracket_h() const { return f0() * f0_h() <= 0; } - - T solve_bracket_l() { - set_cache_H_to_cache_x(); // Evaluated point x becomes new high bound - return solve_bracket(); - } - T solve_bracket_h() { - set_cache_L_to_cache_x(); // Evaluated point x becomes new low bound - return solve_bracket(); - } - - bool want_to_go_l() const { return 0 <= calc_dx_next_default(); } // Because of negative sign in update formula x_next = x() - dx - - // - // - // - // - // - T solve_impl() { - if (0 == f0()) return x(); // Solved at initial x - - - } - - -#if 0 - T solve_impl() { - if (0 == f0()) return x(); // Solved at initial x - - // Calculate direction towards root - const T dx = Step().calc_dx(cache_x_); - const T x_next = x() - dx; // x_next calculated with step dx - - const bool can_bracket_l = f0() * f0_l() <= 0; // Sign flip between x and xl - const bool can_bracket_h = f0() * f0_h() <= 0; // Sign flip between x and xh - const bool want_to_go_l = 0 <= dx; // Because of negative sign in update formula x_next = x() - dx - - if (can_bracket_l && (want_to_go_l || !can_bracket_h)) { - set_cache_H_to_cache_x(); // Evaluated point x becomes new high bound - return solve_bracket(x_next); - } else if (can_bracket_h) { - set_cache_L_to_cache_x(); // Evaluated point x becomes new low bound - return solve_bracket(x_next); - } else { - static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; - return policies::raise_evaluation_error(function, "There appears to be no root to be found in boost::math::tools::detail::BracketHelper, perhaps we have a local minima near current best guess of %1%", x(), boost::math::policies::policy<>()); - } - } -#else - -//////////////////////////////////// - - T solve_impl() { - if (0 == f0()) return x(); // Solved at initial x - - std::cout << "SOLVE_IMPL" << std::endl; - - // Bracket root in direction of dx if possible - if (want_to_go_l() && can_bracket_l()) { - std::cout << "DIRECT --> SOLVE_BRACKET_L" << std::endl; - return solve_bracket_l(); - } else if (can_bracket_h()) { - std::cout << "DIRECT --> SOLVE_BRACKET_H" << std::endl; - return solve_bracket_h(); - } - - return solve_classic_with_fallback(); - } - - T solve_classic_with_fallback() { - std::cout << "Doing classic with fallback" << std::endl; - - // Copy caches - const auto cache_x_copy = cache_x_; - const auto cache_xh_copy = cache_xh_; - const auto cache_xl_copy = cache_xl_; - - bool is_success_classic = false; - const T root = solve_classic(&is_success_classic); - - if (is_success_classic) { - return root; - } else { - // Restore caches - cache_x_ = cache_x_copy; - cache_xh_ = cache_xh_copy; - cache_xl_ = cache_xl_copy; - - // Try bracketing - if (can_bracket_l()) { - return solve_bracket_l(); - } else if (can_bracket_h()) { - return solve_bracket_h(); - } - } - - static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; - return policies::raise_evaluation_error(function, "There appears to be no root to be found in boost::math::tools::detail::BracketHelper, perhaps we have a local minima near current best guess of %1%", x(), boost::math::policies::policy<>()); - } - - // - // The classic solver is used when the root is not bracketed. If the classic solver detects - // a sign flip it passes control to the bracket solver. - // - T solve_classic(bool* is_success) { - BOOST_MATH_STD_USING - static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; - - const bool same_sign_lh = 1 == sign(f0_l()) * sign(f0_h()); - const bool same_sign_lx = 1 == sign(f0_l()) * sign(f0()); - - if (!same_sign_lh || !same_sign_lx) { - *is_success = false; - // 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<>()); - return policies::raise_evaluation_error(function, "Not same sign boost::math::tools::detail::BracketHelper %1%", xl(), boost::math::policies::policy<>()); - } - - const bool is_negate = (1 == sign(f0_l())) ? false : true; - - const auto fn_unoriented = [&](T x) { - T brent_0; - T brent_1; - detail::unpack_tuple(f_(x), brent_0, brent_1); - std::cout << " BRENT eval at: " << x << ", f(x): " << brent_0 << std::endl; - return brent_0; - }; - - const auto fn_brent = [&](T x) { - T brent_0 = fn_unoriented(x); - return is_negate ? -brent_0 : brent_0; - }; - - - for (double i = 0.0; i <= 1.0; i += 0.0001) { - // std::cout << "fn_brent(" << i << "): "; // << fn_brent(i) << std::endl; - fn_brent(i); - } - - std::cout << std::endl; - - - const auto pair_x_fx = brent_find_minima(fn_brent, xl(), xh(), 1024); - const auto x_min = pair_x_fx.first; - const auto fx_min = fn_unoriented(x_min); - - std::cout << "x_min: " << x_min << std::endl; - std::cout << "fx_min: " << fx_min << std::endl; - - if (1 == sign(f0_l()) * sign(fx_min)) { - std::cout << "CLASSIC could not find root bracket" << std::endl; - std::cout << "x: " << x() << ", f0: " << f0() << ", f1: " << cache_x_.f1() << std::endl; - std::cout << "xl: " << xl() << ", f0_l: " << f0_l() << ", f1_l: " << f1_l() << std::endl; - std::cout << "xh: " << xh() << ", f0_h: " << f0_h() << ", f1_h: " << f1_h() << std::endl; - *is_success = false; - return x_min; // Return doesn't matter - } else { - *is_success = true; - std::cout << "CLASSIC found root bracket" << std::endl; - cache_x_ = detail::CacheOrder1(f_, x_min); // Update cache_x_ with x_next - - - - return (sign(fx_min) == 1) ? solve_bracket_h() : solve_bracket_l(); - } - } - -#endif - - T solve_bracket() { - BOOST_MATH_STD_USING - static const char* function = "boost::math::tools::detail::BracketHelper<%1%>"; - - std::cout << "SOLVE_BRACKET" - << "RealType: " << typeid(T).name() - << " -- count: " << count_ << "====================" << std::endl; - - - if (f0_l() == 0.0) { // Low bound is root - std::cout << "Low bound is root" << std::endl; - return xl(); - } else if (f0_h() == 0.0) { // High bound is root - std::cout << "High bound is root" << std::endl; - return xh(); - } else if (0 < f0_l() * f0_h()) { // Bounds do not bracket root - std::cout << "Bounds do not bracket root" << std::endl; - return policies::raise_evaluation_error(function, "Does not bracket boost::math::tools::detail::BracketHelper %1%", xl(), boost::math::policies::policy<>()); - } - - // Orient cache_l and cache_h so that f0_l is negative and f0_h is positive - if (0 < f0_l()) std::swap(cache_xl_, cache_xh_); - - T x_next; - do { - count_--; - - // Calculate Newton step - const T dx_newton = Step().calc_dx(cache_x_); - x_next = x() - dx_newton; // Checked against bracket bounds below - - if (!is_x_between_l_and_h(x_next) || - fabs(dx_pp()) < fabs(dx_newton * 4)) { // Bisection halves step size each iteration - - // if xl and xh have different signs set x_next to 0 - if (-1 == sign(xl()) * sign(xh())) { - x_next = 0; - } else { - const T x_S = (fabs(xl()) < fabs(xh())) ? xl() : xh(); - const T x_B = (fabs(xl()) < fabs(xh())) ? xh() : xl(); - const T SB_mag = fabs(xl() - xh()); - - if (x_S < SB_mag && 1 < fabs(x_S)) { // Spans orders of magnitude - x_next = copysign(sqrt(fabs(x_S)) * sqrt(fabs(x_B)), x_S); - } else { - x_next = (xl() + xh()) / 2; // Midpoint - } - } - - dx_history_.reset(); // Invalidate step size history - } - - const T dx = x_next - x(); - cache_x_ = detail::CacheOrder1(f_, x_next); // Update cache_x_ with x_next - - std::cout << "x: " << x() << ", dx: " << dx << ", dx_newton: " << dx_newton - << ", f1: " << cache_x_.f1() << ", f0: " << cache_x_.f0() - << ", xl: " << xl() << ", xh: " << xh() - << std::endl; - - if (fabs(dx) < fabs(x() * factor_)) { - std::cout << "EXIT: tolerance met" << std::endl; - break; } // Tolerance met - if (x_next == xl() || x_next == xh() || dx == 0) { - std::cout << "EXIT: numerical reason" << std::endl; - break; } - - - // Update bracket - (f0() < 0.0) ? set_cache_L_to_cache_x() : set_cache_H_to_cache_x(); - - dx_history_.push(dx); // Store step size - - if (count_ == 0) { - std::cout << "EXIT: max_iter reached" << std::endl; - break; - } - } while (count_); - - return x(); - } - - bool is_x_between_l_and_h(T x) const { return (std::min(xl(), xh()) <= x && x <= std::max(xl(), xh())); } - - void set_cache_L_to_cache_x() { cache_xl_ = cache_x_; } - void set_cache_H_to_cache_x() { cache_xh_ = cache_x_; } - - T x() const { return cache_x_.x(); } - T f0() const { return cache_x_.f0(); } - T xl() const { return cache_xl_.x(); } - T f0_l() const { return cache_xl_.f0(); } - T f1_l() const { return cache_xl_.f1(); } - T xh() const { return cache_xh_.x(); } - T f0_h() const { return cache_xh_.f0(); } - T f1_h() const { return cache_xh_.f1(); } - T dx_pp() const { return dx_history_.dx_pp(); } - std::uintmax_t count() const { return count_; } - - F f_; - detail::CacheOrder1 cache_x_; - detail::CacheOrder1 cache_xl_; - detail::CacheOrder1 cache_xh_; - std::uintmax_t count_; - T factor_; - StepSizeHistory dx_history_; - }; -} // namespace detail + static const char* function = "boost::math::tools::detail::solve<%3%>"; + const T x_last = state.cache_last().x(); + return policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::solve, last closest guess was %1%", x_last, boost::math::policies::policy<>()); + }; +} // namespace detail template T newton_raphson_iterate(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { - std::cout << "================= newton_raphson_iterate ====================" << std::endl; - return detail::BracketHelper>::solve(f, x, xl, xh, digits, max_iter); + return detail::solve_3_attempts>(f, x, xl, xh, digits, max_iter); } template @@ -2777,9 +820,6 @@ inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) noexcept return newton_raphson_iterate(f, guess, min, max, digits, m); } -#endif - - namespace detail { From 849e9d58509115ad421d774aaed1fa8973477ab3 Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 11 Jul 2023 23:08:21 -0400 Subject: [PATCH 09/18] clarity improvements in comments and code --- include/boost/math/tools/roots.hpp | 196 ++++++++++++++--------------- test/test_root_iterations.cpp | 7 +- test/test_roots.cpp | 4 - 3 files changed, 96 insertions(+), 111 deletions(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index 47152212ee..87a5a58d85 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -24,9 +24,6 @@ #include #include -#include -#include - namespace boost { namespace math { namespace tools { @@ -296,8 +293,10 @@ namespace detail { return CaseFlag::failure; } +#ifdef BOOST_MATH_INSTRUMENT // Makes debugging easier if required void print() const { std::cout << "x: " << x_ << ", f0: " << f0_ << ", f1: " << f1_ << std::endl; } +#endif T x() const { return x_; } T f0() const { return f0_; } // f(x) @@ -305,7 +304,7 @@ namespace detail { private: bool is_sorted_x(const Data_X01& h) const { return x() <= h.x(); } - + // Returns true if *this and h form a U shape. bool is_U(const Data_X01& h) const { const auto& l = *this; @@ -313,7 +312,7 @@ namespace detail { return (sign(l.f0()) * sign(l.f1()) == -1) && // Above zero and sloping down OR below zero and sloping up (sign(h.f0()) * sign(h.f1()) == +1); // Above zero and sloping up OR below zero and sloping down } - + T x_; T f0_; // f(x) T f1_; // f'(x) @@ -340,10 +339,8 @@ namespace detail { T calc_dx(const Data_X01& z) const { return -z.f0() / z.f1(); } }; -//////////////////////////////////////////////////////////////////////////////////////// - // The state of the 1D root finding problem. This state is acted upon by one of - // several solvers + // several solvers template class Root1D_State { private: @@ -353,62 +350,61 @@ namespace detail { l, // Updates the low bound h }; // Updates the high bound - // Allows template dispatch of update_cache + // Allows template dispatch of update_eval template class Val {}; public: - Root1D_State(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) + Root1D_State(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) : data_(xl, xh) , flag_(CaseFlag::case_1) - , evaluator_(f, digits, max_iter) - { - update_cache_sign_dx(x); + , evaluator_(f, digits, max_iter) + { + update_eval_sign_dx(x); } - // Syntactic sugar for type dispatch of update_cache + // Syntactic sugar for type dispatch of update_eval template - void update_cache_sign_f0(Args... args) { update_cache(Val(), args...); } + void update_eval_sign_f0(Args... args) { update_eval(Val(), args...); } template - void update_cache_sign_dx(Args... args) { update_cache(Val(), args...); } + void update_eval_sign_dx(Args... args) { update_eval(Val(), args...); } template - void update_cache_not_last(Args... args) { update_cache(Val(), args...); } + void update_eval_not_last(Args... args) { update_eval(Val(), args...); } template - void update_cache_l(Args... args) { update_cache(Val(), args...); } + void update_eval_l(Args... args) { update_eval(Val(), args...); } template - void update_cache_h(Args... args) { update_cache(Val(), args...); } + void update_eval_h(Args... args) { update_eval(Val(), args...); } // Filter arguments into preferred format template - void update_cache(Val vs, const T x) { update_cache(vs, calc_x_dx(x)); } + void update_eval(Val vs, const T x) { update_eval(vs, calc_x_dx(x)); } template - void update_cache(Val vs, const X_DX& x_dx) { update_cache(vs, eval_next(x_dx), x_dx.dx()); } + void update_eval(Val vs, const X_DX& x_dx) { update_eval(vs, eval_next(x_dx), x_dx.dx()); } template - void update_cache(Val vs, const Data_X01& cx) { update_cache(vs, cx, calc_dx_true(cx)); } + void update_eval(Val vs, const Data_X01& cx) { update_eval(vs, cx, calc_dx_true(cx)); } // Resolve SideSelect - void update_cache(Val, const Data_X01& cx, const T dx) { - const bool is_bound_h = cx.is_same_sign_f0(cache_h()); // Updates high bound if f0 is same sign as f0_h - // TODO: change to update_cache - update_cache_dispatch(is_bound_h, cx, dx); + void update_eval(Val, const Data_X01& cx, const T dx) { + const bool is_bound_h = cx.is_same_sign_f0(eval_h()); // Updates high bound if f0 is same sign as f0_h + update_eval(is_bound_h, cx, dx); } - void update_cache(Val, const Data_X01& cx, const T dx) { + void update_eval(Val, const Data_X01& cx, const T dx) { const bool is_bound_h = cx.is_bound_h(); // Updates high bound if dx is negative - update_cache_dispatch(is_bound_h, cx, dx); + update_eval(is_bound_h, cx, dx); } - void update_cache(Val, const Data_X01& cx, const T dx) { + void update_eval(Val, const Data_X01& cx, const T dx) { const bool is_bound_h = !data_.ind_last(); // Updates high bound if last eval was not high - update_cache_dispatch(is_bound_h, cx, dx); + update_eval(is_bound_h, cx, dx); } - void update_cache(Val, const Data_X01& cx, const T dx) { - update_cache_dispatch(false, cx, dx); // Updates low bound + void update_eval(Val, const Data_X01& cx, const T dx) { + update_eval(false, cx, dx); // Updates low bound } - void update_cache(Val, const Data_X01& cx, const T dx) { - update_cache_dispatch(true, cx, dx); // Updates high bound + void update_eval(Val, const Data_X01& cx, const T dx) { + update_eval(true, cx, dx); // Updates high bound } - // Update the cache - void update_cache_dispatch(const bool is_bound_h, const Data_X01& cx, const T dx) { + // Update the eval + void update_eval(const bool is_bound_h, const Data_X01& cx, const T dx) { dx_history_.push(dx); data_.set(cx, is_bound_h); if (cx.f0() == 0) { flag_ = CaseFlag::success; } @@ -417,30 +413,30 @@ namespace detail { void update_case_flag() { if (flag() == CaseFlag::success) { return; } // Don't update if already successful - // Find status flag if both caches are valid - if (num_valid_caches() == 2) { - flag_ = cache_l().identify_status_flag(cache_h()); + // Find status flag if both evals are valid + if (num_valid_evals() == 2) { + flag_ = eval_l().identify_status_flag(eval_h()); } } - T midpoint() const { return (cache_l().x() + cache_h().x()) / 2; } + T midpoint() const { return (eval_l().x() + eval_h().x()) / 2; } bool is_inbounds(T x) const { return (x_l() <= x) && (x <= x_h()); } - // NOTE: caches always contain valid x data, as both caches are initalized with the + // NOTE: evals always contain valid x data, as both evals are initalized with the // bounds of the root finding problem. - T x_l() const { return cache_l().x(); } - T x_h() const { return cache_h().x(); } + T x_l() const { return eval_l().x(); } + T x_h() const { return eval_h().x(); } + + // Returns the last eval that was evaluated + const Data_X01& eval_last() const { return data_.eval_last(); } - // Returns the last cache that was evaluated - const Data_X01& cache_last() const { return data_.cache_last(); } - int num_x_only() const { return data_.num_x_only(); } - int num_valid_caches() const { return 2 - num_x_only(); } + int num_valid_evals() const { return 2 - num_x_only(); } // Extrapolates the last evaluation to calculate a next x value, and the distance dx // between the last x and this next x. - X_DX extrapolate_last_to_calc_x_dx() const { - const auto& cx = cache_last(); + X_DX extrapolate_last_to_calc_x_dx() const { + const auto& cx = eval_last(); const T x_next = cx.x() + cx.calc_dx(); return calc_x_dx(x_next); // NOTE: the code below is incorrect because (x + dx) - x != dx @@ -450,11 +446,11 @@ namespace detail { } X_DX calc_x_dx(T x) const { return X_DX(x, calc_dx_true(x)); } - // Calculates the x and dx for the next step by subtracting the current x from the next x. - // This approach to calculating dx is aware of floating point precision issues for small dx. + // Calculates dx by evaluating x_next - x_prev. This is the actual floating point dx. T calc_dx_true(const Data_X01& cx) const { return calc_dx_true(cx.x()); } - T calc_dx_true(const T& x) const { return x - cache_last().x(); } + T calc_dx_true(const T& x) const { return x - eval_last().x(); } +#ifdef BOOST_MATH_INSTRUMENT // Makes debugging easier if required void print_status_flag() const { switch (flag()) { @@ -465,9 +461,10 @@ namespace detail { case CaseFlag::failure: std::cout << "failure" << std::endl; break; } } +#endif - // Returns true if the two caches bracket a root where the sign of f(x) changes - bool is_B() const { return cache_l().is_B(cache_h()); } + // Returns true if the two evals bracket a root where the sign of f(x) changes + bool is_B() const { return eval_l().is_B(eval_h()); } void set_flag(CaseFlag flag) { flag_ = flag; } Data_X01 eval_next(const X_DX& x_dx) { return evaluator_(x_dx); } @@ -478,46 +475,46 @@ namespace detail { T dx_pp() const { return dx_history_.dx_pp(); } std::uintmax_t num_fn_evals() const { return evaluator_.num_fn_evals(); } - T get_active_bound() const { return data_.get_bound_from_unused_cache(); } + T get_active_bound() const { return data_.get_bound_from_unused_eval(); } - const Data_X01& cache_l() const { return data_.cache_l(); } - const Data_X01& cache_h() const { return data_.cache_h(); } + const Data_X01& eval_l() const { return data_.eval_l(); } + const Data_X01& eval_h() const { return data_.eval_h(); } CaseFlag flag() const { return flag_; } private: // Holds Evaluation Data class Root1D_Data { public: - // Creates two caches whose x values are initialized with the bounds of the root finding + // Creates two evals whose x values are initialized with the bounds of the root finding // problem, but whose f0 and f1 values are invalid. Root1D_Data(T xl, T xh) - : v_cache_({Data_X01(xl), Data_X01(xh)}) + : v_eval_({Data_X01(xl), Data_X01(xh)}) , v_is_x_only_({true, true}) , ind_last_(true) {} - // Updates the cache using the bool is_h as an index. + // Updates the eval using the bool is_h as an index. void set(const Data_X01& cx, bool is_h) { - v_cache_[is_h] = cx; + v_eval_[is_h] = cx; v_is_x_only_[is_h] = false; ind_last_ = is_h; } - // Returns the x value of the only cache that contains an original bound - T get_bound_from_unused_cache() const { + // Returns the x value of the only eval that contains an original bound + T get_bound_from_unused_eval() const { assert(num_x_only() == 1); - return v_is_x_only_[0] ? cache_l().x() : cache_h().x(); + return v_is_x_only_[0] ? eval_l().x() : eval_h().x(); } - const Data_X01& cache_last() const { return v_cache_[ind_last_]; } - const Data_X01& cache_l() const { return v_cache_[0]; } - const Data_X01& cache_h() const { return v_cache_[1]; } + const Data_X01& eval_last() const { return v_eval_[ind_last_]; } + const Data_X01& eval_l() const { return v_eval_[0]; } + const Data_X01& eval_h() const { return v_eval_[1]; } int num_x_only() const { return v_is_x_only_[0] + v_is_x_only_[1]; } bool ind_last() const { return ind_last_; } // Last evaluated index private: - std::array, 2> v_cache_; + std::array, 2> v_eval_; std::array v_is_x_only_; bool ind_last_; }; @@ -583,8 +580,6 @@ namespace detail { Root1D_Evaluator evaluator_; }; -//////////////////////////////////////////////////////////////////////////////////////// - // The abstract base class for the 1D root finding algorithm template class SolverBase { @@ -601,9 +596,11 @@ namespace detail { x_dx = calc_next_bisection(b); } + // Maybe exit if (is_exit_case_now(b, x_dx)) { break; } - this->update_cache(b, x_dx); + // Update eval + this->update_eval(b, x_dx); } return x_dx.x(); @@ -624,25 +621,24 @@ namespace detail { virtual bool is_exit_case_now(Root1D_State& b, const X_DX x_dx) const = 0; virtual X_DX calc_next_bisection(Root1D_State& b) const = 0; - virtual void update_cache(Root1D_State& b, const X_DX& x_dx) const = 0; + virtual void update_eval(Root1D_State& b, const X_DX& x_dx) const = 0; virtual CaseFlag flag_this() const = 0; }; - // The solver for the case where only a single bracted (either low or high) has been - // evaluated. + // The solver for the case where only a single evulation (either low or high) exists. template class SolverCase1 : public SolverBase { public: X_DX calc_next_bisection(Root1D_State& b) const override { const T x_limit = b.get_active_bound(); const auto next_limit = b.calc_x_dx(x_limit); - b.update_cache_not_last(next_limit); + b.update_eval_not_last(next_limit); b.update_case_flag(); return next_limit; } - void update_cache(Root1D_State&b, const X_DX& x_dx) const override { - b.update_cache_sign_dx(x_dx); + void update_eval(Root1D_State&b, const X_DX& x_dx) const override { + b.update_eval_sign_dx(x_dx); b.update_case_flag(); } @@ -674,26 +670,26 @@ namespace detail { template class SolverCase2U : public SolverCase2Base { public: - void update_cache(Root1D_State& b, const X_DX& x_dx) const override { + void update_eval(Root1D_State& b, const X_DX& x_dx) const override { const auto cx = b.eval_next(x_dx); // Transition to Case 2B because will bracket a root with the new evaluation - if (b.cache_last().is_B(cx)) { + if (b.eval_last().is_B(cx)) { b.set_flag(CaseFlag::case_2B); - + // If the newton step is increasing this suggests that the solver is converging toward a // local minimum that is not a root. This suggests failure. } else if (is_newton_dx_bigger_than_prev(b, cx)) { b.set_flag(CaseFlag::failure); } - b.update_cache_sign_dx(cx, x_dx.dx()); + b.update_eval_sign_dx(cx, x_dx.dx()); } bool is_newton_dx_bigger_than_prev(Root1D_State& b, const Data_X01& cx) const { using std::fabs; - // NOTE: b.cache_last() isn't generally on the same side of the root as cx - const auto& cx_prev_side_eval = cx.is_bound_h() ? b.cache_h() : b.cache_l(); + // NOTE: b.eval_last() isn't generally on the same side of the root as cx + const auto& cx_prev_side_eval = cx.is_bound_h() ? b.eval_h() : b.eval_l(); return fabs(cx_prev_side_eval.calc_dx_newton()) < fabs(cx.calc_dx_newton()); // Uses Newton step } @@ -712,13 +708,12 @@ namespace detail { template class SolverCase2B : public SolverCase2Base { public: - void update_cache(Root1D_State& b, const X_DX& x_dx) const override { - b.update_cache_sign_f0(x_dx); + void update_eval(Root1D_State& b, const X_DX& x_dx) const override { + b.update_eval_sign_f0(x_dx); } bool is_exit_case_now(Root1D_State& b, const X_DX x_dx) const override { if (this->is_dx_small(b, x_dx)) { - std::cout << "dx is small --> success" << std::endl; b.set_flag(CaseFlag::success); return true; } return b.flag() != flag_this(); } @@ -727,14 +722,12 @@ namespace detail { CaseFlag flag_stall() const override { return CaseFlag::success; } }; -//////////////////////////////////////////////////////////////////////////////////////// - // Iterates on the root finding problem using the solver associated with the current case flag. // Returns the result. template std::pair solve_for_state(Root1D_State& state, T x_final) { switch (state.flag()) { - case CaseFlag::case_1: return solve_for_state(state, SolverCase1().solve(state)); + case CaseFlag::case_1: return solve_for_state(state, SolverCase1().solve(state)); case CaseFlag::case_2B: return solve_for_state(state, SolverCase2B().solve(state)); case CaseFlag::case_2U: return solve_for_state(state, SolverCase2U().solve(state)); case CaseFlag::failure: return std::make_pair(false, 0); @@ -744,11 +737,9 @@ namespace detail { } template std::pair solve_for_state(Root1D_State& state) { - return solve_for_state(state, state.cache_last().x()); + return solve_for_state(state, state.eval_last().x()); } -//////////////////////////////////////////////////////////////////////////////////////// - // Makes three attempts to find a root: // 1. Local minimization starting with the initial x // 2. Bracketing the root in the direction of the initial dx if there is a sign flip @@ -765,17 +756,17 @@ namespace detail { Root1D_State state(f, x, xl, xh, digits, max_iter); // Get CacheX - const auto cache_x_orig = state.cache_last(); + const auto eval_x_orig = state.eval_last(); auto fn_form_bracket_with_side = [&](const bool ind_bound){ if (ind_bound) { // NOTE: These two lines cannot be swapped in case xh is a solution. This is tested. - state.update_cache_l(cache_x_orig); // First to original cache - state.update_cache_h(xh); // Then set bound + state.update_eval_l(eval_x_orig); // First to original eval + state.update_eval_h(xh); // Then set bound } else { // NOTE: These two lines cannot be swapped in case xl is a solution. This is tested. - state.update_cache_h(cache_x_orig); // First to original cache - state.update_cache_l(xl); // Then set bound + state.update_eval_h(eval_x_orig); // First to original eval + state.update_eval_l(xl); // Then set bound } // state.reset_dx_history(); // Not required for correctness @@ -791,20 +782,20 @@ namespace detail { const auto p = solve_for_state(state); // If local minimization was successful, return root - if (p.first) { return fn_return(state, p); } + if (p.first) { return fn_return(state, p); } // Attempt to bracket root in the direction of the initial guess - fn_form_bracket_with_side(!cache_x_orig.is_bound_h()); + fn_form_bracket_with_side(!eval_x_orig.is_bound_h()); if (state.is_B()) { return fn_return(state, solve_for_state(state)); } // Attempt to bracket root in the opposite direction of the initial guess - fn_form_bracket_with_side(cache_x_orig.is_bound_h()); + fn_form_bracket_with_side(eval_x_orig.is_bound_h()); if (state.is_B()) { return fn_return(state, solve_for_state(state)); } static const char* function = "boost::math::tools::detail::solve<%3%>"; - const T x_last = state.cache_last().x(); + const T x_last = state.eval_last().x(); return policies::raise_evaluation_error(function, "Unable to bracket root in boost::math::tools::detail::solve, last closest guess was %1%", x_last, boost::math::policies::policy<>()); - }; + } } // namespace detail template @@ -820,7 +811,6 @@ inline T newton_raphson_iterate(F f, T guess, T min, T max, int digits) noexcept return newton_raphson_iterate(f, guess, min, max, digits, m); } - namespace detail { struct halley_step diff --git a/test/test_root_iterations.cpp b/test/test_root_iterations.cpp index 9f2b4d213d..fada5d9bb0 100644 --- a/test/test_root_iterations.cpp +++ b/test/test_root_iterations.cpp @@ -169,7 +169,7 @@ BOOST_AUTO_TEST_CASE( test_main ) // Newton next: iters = 1000; dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, guess / 2, result * 10, newton_limits, iters); - BOOST_CHECK_CLOSE_FRACTION(dr, result, ldexp(1.0, 1 - newton_limits)); + BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits::epsilon() * 2); BOOST_CHECK_LE(iters, 12); // Halley next: iters = 1000; @@ -192,7 +192,7 @@ BOOST_AUTO_TEST_CASE( test_main ) // Newton next: iters = 1000; dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters); - BOOST_CHECK_CLOSE_FRACTION(dr, result, ldexp(1.0, 1 - newton_limits)); + BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits::epsilon() * 2); BOOST_CHECK_LE(iters, 12); // Halley next: iters = 1000; @@ -238,7 +238,7 @@ BOOST_AUTO_TEST_CASE( test_main ) // Newton next: iters = 1000; dr = boost::math::tools::newton_raphson_iterate(cbrt_functor_deriv(arg), guess, result / 10, result * 10, newton_limits, iters); - BOOST_CHECK_CLOSE_FRACTION(dr, result, ldexp(1.0, 1 - newton_limits)); + BOOST_CHECK_CLOSE_FRACTION(dr, result, std::numeric_limits::epsilon() * 2); BOOST_CHECK_LE(iters, 5); // Halley next: iters = 1000; @@ -262,7 +262,6 @@ BOOST_AUTO_TEST_CASE( test_main ) #endif #define T double - # include "ibeta_small_data.ipp" for (unsigned i = 0; i < ibeta_small_data.size(); ++i) diff --git a/test/test_roots.cpp b/test/test_roots.cpp index 32714cd33a..53d1131a77 100644 --- a/test/test_roots.cpp +++ b/test/test_roots.cpp @@ -762,10 +762,6 @@ void test_failures() // There is no root: BOOST_CHECK_THROW(boost::math::tools::newton_raphson_iterate([](double x) { return std::make_pair(x * x + 1, 2 * x); }, 10.0, -12.0, 12.0, 52), boost::math::evaluation_error); BOOST_CHECK_THROW(boost::math::tools::newton_raphson_iterate([](double x) { return std::make_pair(x * x + 1, 2 * x); }, -10.0, -12.0, 12.0, 52), boost::math::evaluation_error); - // There is a root, but a bad guess takes us into a local minima: -#if 0 - BOOST_CHECK_THROW(boost::math::tools::newton_raphson_iterate([](double x) { return std::make_pair(boost::math::pow<6>(x) - 2 * boost::math::pow<4>(x) + x + 0.5, 6 * boost::math::pow<5>(x) - 8 * boost::math::pow<3>(x) + 1); }, 0.75, -20., 20., 52), boost::math::evaluation_error); -#endif // There is no root: BOOST_CHECK_THROW(boost::math::tools::halley_iterate([](double x) { return std::make_tuple(x * x + 1, 2 * x, 2); }, 10.0, -12.0, 12.0, 52), boost::math::evaluation_error); From 7a4c5617e8d93e0bfe74bc49a9dcc174df872f20 Mon Sep 17 00:00:00 2001 From: Ryan Date: Wed, 12 Jul 2023 00:24:32 -0400 Subject: [PATCH 10/18] fix inspect issues --- include/boost/math/special_functions/detail/ibeta_inverse.hpp | 2 +- include/boost/math/tools/roots.hpp | 2 +- test/test_roots.cpp | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/boost/math/special_functions/detail/ibeta_inverse.hpp b/include/boost/math/special_functions/detail/ibeta_inverse.hpp index 43ecbf4666..90471bfc73 100644 --- a/include/boost/math/special_functions/detail/ibeta_inverse.hpp +++ b/include/boost/math/special_functions/detail/ibeta_inverse.hpp @@ -27,7 +27,7 @@ struct temme_root_finder { temme_root_finder(const T t_, const T a_) : t(t_), a(a_) { const T x_extrema = 1 / (1 + a); - assert(0 < x_extrema && x_extrema < 1); + BOOST_MATH_ASSERT(0 < x_extrema && x_extrema < 1); } boost::math::tuple operator()(T x) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index 87a5a58d85..ba834f4ed3 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -501,7 +501,7 @@ namespace detail { // Returns the x value of the only eval that contains an original bound T get_bound_from_unused_eval() const { - assert(num_x_only() == 1); + BOOST_MATH_ASSERT(num_x_only() == 1); return v_is_x_only_[0] ? eval_l().x() : eval_h().x(); } diff --git a/test/test_roots.cpp b/test/test_roots.cpp index 53d1131a77..65d60449eb 100644 --- a/test/test_roots.cpp +++ b/test/test_roots.cpp @@ -519,8 +519,8 @@ void test_cusp() { // Two roots exist (One is picked arbitrarily) } else { - // x_root = x_center ± ye^2 - // x_root - x_center = ± ye^2 + // x_root = x_center +- ye^2 + // x_root - x_center = +- ye^2 // abs(x_root - x_center) = ye^2 BOOST_CHECK_CLOSE_FRACTION(fabs(x_root - x_center), ye * ye, tol); } From 90582fc9c5f5d71cd06cc3a96441c1994c431592 Mon Sep 17 00:00:00 2001 From: Ryan Date: Wed, 12 Jul 2023 14:19:40 -0400 Subject: [PATCH 11/18] fix CI errors: division by zero and array include --- include/boost/math/tools/roots.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index ba834f4ed3..beaffaeac4 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -11,6 +11,7 @@ #endif #include // test for multiprecision types in complex Newton +#include #include #include #include @@ -525,10 +526,9 @@ namespace detail { StepSizeHistory() { reset(); } void reset() { - // NOTE: for T = boost::math::concepts::real_concept, the following evaluate to 0: - // std::numeric_limits::infinity() and std::numeric_limits::infinity() - dx_p_ = static_cast(1) / 0; - dx_pp_ = static_cast(1) / 0; + // Need to static cast double infinity because for T = boost::math::concepts::real_concept, + // std::numeric_limits::infinity() evaluates to 0. + dx_p_ = dx_pp_ = static_cast(std::numeric_limits::infinity()); } void push(T input) { From 42b0bbce3dc6dea53333b2ac1acf6891c34442ef Mon Sep 17 00:00:00 2001 From: Ryan Date: Wed, 12 Jul 2023 15:03:12 -0400 Subject: [PATCH 12/18] direct init array to fix clang 6 error --- include/boost/math/tools/roots.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index beaffaeac4..fe2e6ab396 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -490,7 +490,7 @@ namespace detail { // problem, but whose f0 and f1 values are invalid. Root1D_Data(T xl, T xh) : v_eval_({Data_X01(xl), Data_X01(xh)}) - , v_is_x_only_({true, true}) + , v_is_x_only_({{true, true}}) , ind_last_(true) {} // Updates the eval using the bool is_h as an index. From 2cea88830d3e3de03e25f0302462afb14803c61f Mon Sep 17 00:00:00 2001 From: Ryan Date: Sat, 22 Jul 2023 11:49:36 -0400 Subject: [PATCH 13/18] fix another windows 1/0 compile error --- include/boost/math/tools/roots.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index fe2e6ab396..46ee489157 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -273,7 +273,10 @@ namespace detail { Data_X01(F f, T x) : x_(x) { detail::unpack_tuple(f(x_), f0_, f1_); } // Store just a constant x value - Data_X01(T x) : x_(x), f0_(static_cast(1) / 0), f1_(static_cast(1) / 0) {} + Data_X01(T x) + : x_(x) + , f0_(static_cast(std::numeric_limits::infinity())) + , f1_(static_cast(std::numeric_limits::infinity())) {} T calc_dx() const { return Step().calc_dx(*this); } T calc_dx_newton() const { return -f0_ / f1_; } From e94fa65794492e48b6b78463d0422069ad14d3a5 Mon Sep 17 00:00:00 2001 From: Ryan Date: Mon, 31 Jul 2023 11:01:54 -0400 Subject: [PATCH 14/18] Merge remote-tracking branch 'origin/develop' --- include/boost/math/special_functions/fibonacci.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/boost/math/special_functions/fibonacci.hpp b/include/boost/math/special_functions/fibonacci.hpp index 7aaf6dd854..2eb4cabaff 100644 --- a/include/boost/math/special_functions/fibonacci.hpp +++ b/include/boost/math/special_functions/fibonacci.hpp @@ -32,8 +32,8 @@ inline BOOST_CXX14_CONSTEXPR T unchecked_fibonacci(unsigned long long n) noexcep if (n <= 2) return n == 0 ? 0 : 1; /* * This is based on the following identities by Dijkstra: - * F(2*n) = F(n)^2 + F(n+1)^2 - * F(2*n+1) = (2 * F(n) + F(n+1)) * F(n+1) + * F(2*n-1) = F(n-1)^2 + F(n)^2 + * F(2*n) = (2*F(n-1) + F(n)) * F(n) * The implementation is iterative and is unrolled version of trivial recursive implementation. */ unsigned long long mask = 1; From 632a2590c424ab28cea792c7fe6c4d6c760942de Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 1 Aug 2023 22:59:22 -0400 Subject: [PATCH 15/18] add BitSpace for smarter bisection --- include/boost/math/tools/roots.hpp | 189 ++++++++++++++++++++++++++++- 1 file changed, 185 insertions(+), 4 deletions(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index 46ee489157..e2d720e57d 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -25,6 +25,8 @@ #include #include +#include + namespace boost { namespace math { namespace tools { @@ -343,6 +345,155 @@ namespace detail { T calc_dx(const Data_X01& z) const { return -z.f0() / z.f1(); } }; + ////// Motivation for the BitSpace class. ////// + // + // What's the best way to bisect between a lower bound (lb) and an upper + // bound (ub) during root finding? Let's consider options... + // + // Arithmetic bisection: + // - The natural choice, but it doesn't always work well. For example, if + // lb = 1.0 and ub = std::numeric_limits::max(), many bisections + // may be needed to converge if the root is near 1. + // + // Geometric bisection: + // - This approach performs much better for the example above, but it + // too has issues. For example, if lb = 0.0, it gets stuck at 0.0. + // It also fails if lb and ub have different signs. + // + // In addition to the limitations outlined above, neither of these approaches + // works if ub is infinity. We want a more robust way to handle bisection + // for general root finding problems. The BitSpace class does this. + // + // + ////// The BitSpace class. ////// + // + // On a conceptual level, the BitSpace class is designed to solve the + // following root finding problem. + // - A function f(x) has a single root x_solution somewhere in the interval + // [-infinity, +infinity]. For all values below x_solution f(x) is -1. + // For all values above x_solution f(x) is +1. The best way to root find + // this problem is to bisect in bit space. + // + // Efficient bit space bisection is possible because of the IEEE 754 standard. + // According to the standard, the bits in floating point numbers are partitioned + // into three parts: sign, exponent, and mantissa. Assuming constant sign, + // increasing numbers in bit space have increasing floating point values + // starting at zero, and ending at infinity! + // + // 0 | 0 00000000 00000000000000000000000 | positive zero + // 1.4013e-45 | 0 00000000 00000000000000000000001 | std::numeric_limits::denorm_min() + // 1.17549e-38 | 0 00000001 00000000000000000000000 | std::numeric_limits::min() + // 1.19209e-07 | 0 01101000 00000000000000000000000 | std::numeric_limits::epsilon() + // 1 | 0 01111111 00000000000000000000000 | positive one + // 3.40282e+38 | 0 11111110 11111111111111111111111 | std::numeric_limits::max() + // inf | 0 11111111 00000000000000000000000 | std::numeric_limits::infinity() + // nan | 0 11111111 10000000000000000000000 | std::numeric_limits::quiet_NaN() + // + // This class interprets the bits as either a float or a uint, depending on + // the context. + // + class BitSpace { + public: + BitSpace() = delete; + + template + static T calc_midpoint(T x, T X) { + // Sort x and X by magnitude + using std::fabs; + if (fabs(X) < fabs(x)) { + std::swap(x, X); + } + + // Return the arithmetic mean if x and X have the same sign and are within a factor of 2 of each other + if (sign(x) == sign(X) && fabs(X) < 2 * fabs(x)) { + return (x + X) / 2; + } + + return do_bit_space_mean_754(x, X); + } + + private: + static float do_bit_space_mean_754(float x, float X) { + return BitSpaceMidpoint754(x, X)(); + } + + template + static T do_bit_space_mean_754(T x, T X) { + double x_double = cast_to_double(x); + double X_double = cast_to_double(X); + return BitSpaceMidpoint754(x_double, X_double)(); + } + + // Static casting the float type T to double is not possible for T = + // boost::concepts::real_concept types. This helper template detects + // if T is a real_concept type indirectly through the presence of a + // value() member function. + template + struct has_value_function : std::false_type {}; + template + struct has_value_function().value())>> : std::true_type {}; + + // Cast type T to double if T is a boost::concepts::real_concept type. + template + static typename std::enable_if::value, double>::type cast_to_double(T x) { + return static_cast(x.value()); + } + template + static typename std::enable_if::value, double>::type cast_to_double(T x) { + return static_cast(x); + } + + // Does the bisection in bit space for IEEE 754 floating point numbers. + template + class BitSpaceMidpoint754 { + public: + // NOTE: x and X must be sorted by magnitude + BitSpaceMidpoint754(T x, T X) : x_(x), X_(X) { + static_assert(std::numeric_limits::is_iec559, "Type must be IEEE 754 floating point."); + static_assert(std::is_unsigned::value, "U must be an unsigned integer type."); + static_assert(sizeof(T) == sizeof(U), "Type and uint size must be the same."); + } + + T operator()() const { + // Convert x_ and X_ to magnitude and sign + using std::fabs; + const T x_mag = fabs(x_); + const T X_mag = fabs(X_); + const T sign_x = sign(x_); + const T sign_X = sign(X_); + + // Convert the magnitudes to bits + U bits_mag_x = float_to_uint(x_mag); + U bits_mag_X = float_to_uint(X_mag); + + // Calculate the average magnitude in bits + U bits_mag = (sign_x == sign_X) ? (bits_mag_X + bits_mag_x) : (bits_mag_X - bits_mag_x); + bits_mag = bits_mag >> 1; // Divide by 2 + + // Reconstruct upl_mean from average magnitude and sign of X + return uint_to_float(bits_mag) * sign_X; + } + + private: + // Convert a float to bits + static U float_to_uint(T x) { + U bits; + std::memcpy(&bits, &x, sizeof(U)); + return bits; + } + + // Convert bits to a float + static T uint_to_float(U bits) { + T x; + std::memcpy(&x, &bits, sizeof(T)); + return x; + } + + T x_; + T X_; // Has a larger magnitude than x_ + }; + }; + // The state of the 1D root finding problem. This state is acted upon by one of // several solvers template @@ -423,7 +574,7 @@ namespace detail { } } - T midpoint() const { return (eval_l().x() + eval_h().x()) / 2; } + // T midpoint() const { return (eval_l().x() + eval_h().x()) / 2; } bool is_inbounds(T x) const { return (x_l() <= x) && (x <= x_h()); } // NOTE: evals always contain valid x data, as both evals are initalized with the @@ -485,6 +636,8 @@ namespace detail { const Data_X01& eval_h() const { return data_.eval_h(); } CaseFlag flag() const { return flag_; } + void print_step_size_history() const { dx_history_.print(); } + private: // Holds Evaluation Data class Root1D_Data { @@ -539,6 +692,10 @@ namespace detail { dx_p_ = input; } + void print() const { + std::cout << "dx_p: " << dx_p_ << ", dx_pp: " << dx_pp_ << std::endl; + } + T dx_pp() const { return dx_pp_; } private: @@ -561,7 +718,16 @@ namespace detail { policies::raise_evaluation_error(function, "Ran out of iterations when finding root in boost::math::tools::detail::Root1D_Evaluator, attempted to evalutate %1%", x_dx.x(), boost::math::policies::policy<>()); } --count_; - return Data_X01(f_, x_dx.x()); + + // Set precision of cout to 20 digits + std::cout.precision(20); + + std::cout << "evaluated at x: " << x_dx.x() << "-- with dx: " << x_dx.dx() << " -- count: " << count_ << std::endl;\ + const auto data_out = Data_X01(f_, x_dx.x()); + + std::cout << "f0: " << data_out.f0() << ", f1: " << data_out.f1() << std::endl; + + return data_out; } std::uintmax_t num_fn_evals() const { return max_iter_ - count_; } @@ -599,6 +765,8 @@ namespace detail { x_dx = calc_next_bisection(b); } + b.print_step_size_history(); + // Maybe exit if (is_exit_case_now(b, x_dx)) { break; } @@ -628,7 +796,7 @@ namespace detail { virtual CaseFlag flag_this() const = 0; }; - // The solver for the case where only a single evulation (either low or high) exists. + // The solver for the case where only a single evaluation (either low or high) exists. template class SolverCase1 : public SolverBase { public: @@ -659,10 +827,21 @@ namespace detail { class SolverCase2Base : public SolverBase { public: X_DX calc_next_bisection(Root1D_State& b) const override { - const auto x_mid = b.midpoint(); + const auto x_mid = BitSpace::calc_midpoint(b.x_l(), b.x_h()); + + // const auto x_mid = b.midpoint(); if (x_mid == b.x_l() || x_mid == b.x_h()) { b.set_flag(this->flag_stall()); } + + // Assert that x_mid is in bounds + if (!b.is_inbounds(x_mid)) { + std::cout << "x_mid: " << x_mid << std::endl; + std::cout << "x_l: " << b.x_l() << std::endl; + std::cout << "x_h: " << b.x_h() << std::endl; + assert(false); + } + return b.calc_x_dx(x_mid); } @@ -755,6 +934,8 @@ namespace detail { // template T solve_3_attempts(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { + std::cout << "CONSTRUCTED" << std::endl; + // Create Root1D_State Root1D_State state(f, x, xl, xh, digits, max_iter); From 33319e7e7eb8af3c018b268fd08003f2e8c9bd9a Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 1 Aug 2023 23:29:45 -0400 Subject: [PATCH 16/18] Remove void_t because it requires C++17 --- include/boost/math/tools/roots.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index e2d720e57d..e83ea66e83 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -431,7 +431,7 @@ namespace detail { template struct has_value_function : std::false_type {}; template - struct has_value_function().value())>> : std::true_type {}; + struct has_value_function().value()))> : std::true_type {}; // Cast type T to double if T is a boost::concepts::real_concept type. template From 2756324185fc9fbe86382d456bc6ebe30914714c Mon Sep 17 00:00:00 2001 From: Ryan Date: Fri, 4 Aug 2023 16:52:53 -0400 Subject: [PATCH 17/18] refactor approach --- include/boost/math/tools/roots.hpp | 297 ++++++++++++++++++----------- test/test_roots.cpp | 84 ++++++++ 2 files changed, 266 insertions(+), 115 deletions(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index e83ea66e83..4991333926 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -25,8 +25,6 @@ #include #include -#include - namespace boost { namespace math { namespace tools { @@ -345,7 +343,7 @@ namespace detail { T calc_dx(const Data_X01& z) const { return -z.f0() / z.f1(); } }; - ////// Motivation for the BitSpace class. ////// + ////// Motivation for the Bisection namespace. ////// // // What's the best way to bisect between a lower bound (lb) and an upper // bound (ub) during root finding? Let's consider options... @@ -362,13 +360,14 @@ namespace detail { // // In addition to the limitations outlined above, neither of these approaches // works if ub is infinity. We want a more robust way to handle bisection - // for general root finding problems. The BitSpace class does this. - // + // for general root finding problems. That's what this namespace is for. // - ////// The BitSpace class. ////// + namespace Bisection { + + ////// The Midpoint754 class ////// // - // On a conceptual level, the BitSpace class is designed to solve the - // following root finding problem. + // On a conceptual level, this class is designed to solve the following root + // finding problem. // - A function f(x) has a single root x_solution somewhere in the interval // [-infinity, +infinity]. For all values below x_solution f(x) is -1. // For all values above x_solution f(x) is +1. The best way to root find @@ -376,9 +375,10 @@ namespace detail { // // Efficient bit space bisection is possible because of the IEEE 754 standard. // According to the standard, the bits in floating point numbers are partitioned - // into three parts: sign, exponent, and mantissa. Assuming constant sign, - // increasing numbers in bit space have increasing floating point values - // starting at zero, and ending at infinity! + // into three parts: sign, exponent, and mantissa. As long as the sign of the + // of the number stays the same, increasing numbers in bit space have increasing + // floating point values starting at zero, and ending at infinity! The table + // below shows select numbers for float (single precision). // // 0 | 0 00000000 00000000000000000000000 | positive zero // 1.4013e-45 | 0 00000000 00000000000000000000001 | std::numeric_limits::denorm_min() @@ -389,110 +389,196 @@ namespace detail { // inf | 0 11111111 00000000000000000000000 | std::numeric_limits::infinity() // nan | 0 11111111 10000000000000000000000 | std::numeric_limits::quiet_NaN() // - // This class interprets the bits as either a float or a uint, depending on - // the context. + // Negative values are similar, but the sign bit is set to 1. My keeping track of the possible + // sign flip, it can bisect numbers with different signs. // - class BitSpace { - public: - BitSpace() = delete; + template + class Midpoint754 { + private: + // Does the bisection in bit space for IEEE 754 floating point numbers. + // Infinities are allowed. It's assumed that neither x nor X is NaN. + static_assert(std::numeric_limits::is_iec559, "Type must be IEEE 754 floating point."); + static_assert(std::is_unsigned::value, "U must be an unsigned integer type."); + static_assert(sizeof(T) == sizeof(U), "Type and uint size must be the same."); + + // Convert float to bits + static U float_to_uint(T x) { + U bits; + std::memcpy(&bits, &x, sizeof(U)); + return bits; + } - template - static T calc_midpoint(T x, T X) { - // Sort x and X by magnitude + // Convert bits to float + static T uint_to_float(U bits) { + T x; + std::memcpy(&x, &bits, sizeof(T)); + return x; + } + + public: + static T solve(T x, T X) { using std::fabs; + + // Sort so that X has the larger magnitude if (fabs(X) < fabs(x)) { std::swap(x, X); } + + const T x_mag = std::fabs(x); + const T X_mag = std::fabs(X); + const T sign_x = sign(x); + const T sign_X = sign(X); - // Return the arithmetic mean if x and X have the same sign and are within a factor of 2 of each other - if (sign(x) == sign(X) && fabs(X) < 2 * fabs(x)) { - return (x + X) / 2; - } + // Convert the magnitudes to bits + U bits_mag_x = float_to_uint(x_mag); + U bits_mag_X = float_to_uint(X_mag); - return do_bit_space_mean_754(x, X); + // Calculate the average magnitude in bits + U bits_mag = (sign_x == sign_X) ? (bits_mag_X + bits_mag_x) : (bits_mag_X - bits_mag_x); + bits_mag = bits_mag >> 1; // Divide by 2 + + // Reconstruct upl_mean from average magnitude and sign of X + return uint_to_float(bits_mag) * sign_X; } + }; // class Midpoint754 + + template + class MidpointNon754 { private: - static float do_bit_space_mean_754(float x, float X) { - return BitSpaceMidpoint754(x, X)(); - } - - template - static T do_bit_space_mean_754(T x, T X) { - double x_double = cast_to_double(x); - double X_double = cast_to_double(X); - return BitSpaceMidpoint754(x_double, X_double)(); - } + static_assert(!std::is_same::value, "Need to use Midpoint754 solver when T is float"); + static_assert(!std::is_same::value, "Need to use Midpoint754 solver when T is double"); - // Static casting the float type T to double is not possible for T = - // boost::concepts::real_concept types. This helper template detects - // if T is a real_concept type indirectly through the presence of a - // value() member function. - template - struct has_value_function : std::false_type {}; - template - struct has_value_function().value()))> : std::true_type {}; + public: + static T solve(T x, T X) { + const T sx = sign(x); + const T sX = sign(X); + + // Sign flip return zero + if (sx * sX == -1) { return T(0.0); } - // Cast type T to double if T is a boost::concepts::real_concept type. - template - static typename std::enable_if::value, double>::type cast_to_double(T x) { - return static_cast(x.value()); - } - template - static typename std::enable_if::value, double>::type cast_to_double(T x) { - return static_cast(x); + // At least one is positive + if (0 < sx + sX) { return do_solve(x, X); } + + // At least one is negative + return -do_solve(-x, -X); } - // Does the bisection in bit space for IEEE 754 floating point numbers. - template - class BitSpaceMidpoint754 { + private: + struct EqZero { + EqZero(T x) { BOOST_MATH_ASSERT(x == 0 && "x must be zero."); } + }; + + struct EqInf { + EqInf(T x) { BOOST_MATH_ASSERT(x == static_cast(std::numeric_limits::infinity()) && "x must be infinity."); } + }; + + class PosFinite { public: - // NOTE: x and X must be sorted by magnitude - BitSpaceMidpoint754(T x, T X) : x_(x), X_(X) { - static_assert(std::numeric_limits::is_iec559, "Type must be IEEE 754 floating point."); - static_assert(std::is_unsigned::value, "U must be an unsigned integer type."); - static_assert(sizeof(T) == sizeof(U), "Type and uint size must be the same."); + PosFinite(T x) : x_(x) { + BOOST_MATH_ASSERT(0 < x && "x must be positive."); + BOOST_MATH_ASSERT(x < std::numeric_limits::infinity() && "x must be less than infinity."); } - T operator()() const { - // Convert x_ and X_ to magnitude and sign - using std::fabs; - const T x_mag = fabs(x_); - const T X_mag = fabs(X_); - const T sign_x = sign(x_); - const T sign_X = sign(X_); - - // Convert the magnitudes to bits - U bits_mag_x = float_to_uint(x_mag); - U bits_mag_X = float_to_uint(X_mag); + T value() const { return x_; } - // Calculate the average magnitude in bits - U bits_mag = (sign_x == sign_X) ? (bits_mag_X + bits_mag_x) : (bits_mag_X - bits_mag_x); - bits_mag = bits_mag >> 1; // Divide by 2 + private: + T x_; + }; - // Reconstruct upl_mean from average magnitude and sign of X - return uint_to_float(bits_mag) * sign_X; + // Two unknowns + static T do_solve(T x, T X) { + if (X < x) { + return do_solve(X, x); } - - private: - // Convert a float to bits - static U float_to_uint(T x) { - U bits; - std::memcpy(&bits, &x, sizeof(U)); - return bits; + + if (x == 0) { + return do_solve(EqZero(x), X); + } else if (x == static_cast(std::numeric_limits::infinity())) { + return static_cast(std::numeric_limits::infinity()); + } else { + return do_solve(PosFinite(x), X); } + } - // Convert bits to a float - static T uint_to_float(U bits) { - T x; - std::memcpy(&x, &bits, sizeof(T)); - return x; + // One unknowns + static T do_solve(EqZero x, T X) { + if (X == 0) { + return T(0.0); + } else if (X == static_cast(std::numeric_limits::infinity())) { + return T(1.0); + } else { + return do_solve(x, PosFinite(X)); } + } + static T do_solve(PosFinite x, T X) { + if (X == static_cast(std::numeric_limits::infinity())) { + return do_solve(x, EqInf(X)); + } else { + return do_solve(x, PosFinite(X)); + } + } - T x_; - T X_; // Has a larger magnitude than x_ - }; - }; + // Zero unknowns + template + static typename std::enable_if::is_specialized, T>::type + do_solve(PosFinite x, EqInf X) { + return do_solve(x, PosFinite(std::numeric_limits::max())); + } + template + static typename std::enable_if::is_specialized, T>::type + do_solve(PosFinite x, EqInf X) { + BOOST_MATH_ASSERT(false && "infinite bounds support requires specialization."); + return static_cast(std::numeric_limits::signaling_NaN()); + } + + template + static typename std::enable_if::is_specialized, U>::type + do_solve(EqZero x, PosFinite X) { + const auto get_smallest_value = []() { + const U denorm_min = std::numeric_limits::denorm_min(); + if (denorm_min != 0) { return denorm_min; } + + const U min = std::numeric_limits::min(); + if (min != 0) { return min; } + + BOOST_MATH_ASSERT(false && "denorm_min and min are both zero."); + return static_cast(std::numeric_limits::signaling_NaN()); + }; + + return do_solve(PosFinite(get_smallest_value()), X); + } + template + static typename std::enable_if::is_specialized, U>::type + do_solve(EqZero x, PosFinite X) { return X.value() / U(2); } + + static T do_solve(PosFinite x, PosFinite X) { + BOOST_MATH_ASSERT(x.value() <= X.value() && "x must be less than or equal to X."); + + const T xv = x.value(); + const T Xv = X.value(); + + // Take arithmetic mean if they are close enough + if (Xv < xv * 8) { return (Xv - xv) / 2 + xv; } // NOTE: avoids overflow + + // Take geometric mean if they are far apart + using std::sqrt; + return sqrt(xv) * sqrt(Xv); // NOTE: avoids overflow + } + }; // class MidpointNon754 + + template + static T calc_midpoint(T x, T X) { + return MidpointNon754::solve(x, X); + } + static float calc_midpoint(float x, float X) { + return Midpoint754::solve(x, X); + } + static double calc_midpoint(double x, double X) { + return Midpoint754::solve(x, X); + } + + } // namespace Bisection // The state of the 1D root finding problem. This state is acted upon by one of // several solvers @@ -616,6 +702,8 @@ namespace detail { case CaseFlag::failure: std::cout << "failure" << std::endl; break; } } + + void print_step_size_history() const { dx_history_.print(); } #endif // Returns true if the two evals bracket a root where the sign of f(x) changes @@ -636,8 +724,6 @@ namespace detail { const Data_X01& eval_h() const { return data_.eval_h(); } CaseFlag flag() const { return flag_; } - void print_step_size_history() const { dx_history_.print(); } - private: // Holds Evaluation Data class Root1D_Data { @@ -692,9 +778,11 @@ namespace detail { dx_p_ = input; } +#if defined(BOOST_MATH_INSTRUMENT) void print() const { std::cout << "dx_p: " << dx_p_ << ", dx_pp: " << dx_pp_ << std::endl; } +#endif T dx_pp() const { return dx_pp_; } @@ -719,15 +807,7 @@ namespace detail { } --count_; - // Set precision of cout to 20 digits - std::cout.precision(20); - - std::cout << "evaluated at x: " << x_dx.x() << "-- with dx: " << x_dx.dx() << " -- count: " << count_ << std::endl;\ - const auto data_out = Data_X01(f_, x_dx.x()); - - std::cout << "f0: " << data_out.f0() << ", f1: " << data_out.f1() << std::endl; - - return data_out; + return Data_X01(f_, x_dx.x()); } std::uintmax_t num_fn_evals() const { return max_iter_ - count_; } @@ -765,8 +845,6 @@ namespace detail { x_dx = calc_next_bisection(b); } - b.print_step_size_history(); - // Maybe exit if (is_exit_case_now(b, x_dx)) { break; } @@ -827,21 +905,12 @@ namespace detail { class SolverCase2Base : public SolverBase { public: X_DX calc_next_bisection(Root1D_State& b) const override { - const auto x_mid = BitSpace::calc_midpoint(b.x_l(), b.x_h()); + const auto x_mid = Bisection::calc_midpoint(b.x_l(), b.x_h()); - // const auto x_mid = b.midpoint(); if (x_mid == b.x_l() || x_mid == b.x_h()) { b.set_flag(this->flag_stall()); } - // Assert that x_mid is in bounds - if (!b.is_inbounds(x_mid)) { - std::cout << "x_mid: " << x_mid << std::endl; - std::cout << "x_l: " << b.x_l() << std::endl; - std::cout << "x_h: " << b.x_h() << std::endl; - assert(false); - } - return b.calc_x_dx(x_mid); } @@ -934,8 +1003,6 @@ namespace detail { // template T solve_3_attempts(F f, T x, T xl, T xh, int digits, std::uintmax_t& max_iter) noexcept(policies::is_noexcept_error_policy >::value&& BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { - std::cout << "CONSTRUCTED" << std::endl; - // Create Root1D_State Root1D_State state(f, x, xl, xh, digits, max_iter); diff --git a/test/test_roots.cpp b/test/test_roots.cpp index baf9d41849..0c8e17e0cb 100644 --- a/test/test_roots.cpp +++ b/test/test_roots.cpp @@ -24,7 +24,9 @@ #include #include +#include #include +#include #define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \ {\ @@ -771,6 +773,86 @@ void test_failures() #endif } +template +void test_bisect() { + // Creates a vector of test values that include: zero, denorm_min, min, epsilon, max, + // and infinity. Also includes powers of 2 ranging from 2^-15 to 2^+15 by powers of 3. + const auto fn_create_test_ladder = [](){ + std::vector v; + + const auto fn_push_back_x_scaled = [&v](const T& x) { + const auto fn_add_if_non_zero_and_finite = [&v](const T& x) { + if (x != 0 && boost::math::isfinite(x)) { // Avoids most duplications + v.push_back(x); + } + }; + + const T two_thirds = T(2) / T(3); + const T four_thirds = T(4) / T(3); + + v.push_back(x); + fn_add_if_non_zero_and_finite(x * T(0.5)); + fn_add_if_non_zero_and_finite(x * two_thirds); + fn_add_if_non_zero_and_finite(x * four_thirds); + fn_add_if_non_zero_and_finite(x * T(2.0)); + }; + + const auto fn_push_back_pm = [&](std::vector& v, const T& x) { + fn_push_back_x_scaled( x); + fn_push_back_x_scaled(-x); + }; + + fn_push_back_pm(v, T(0.0)); + fn_push_back_pm(v, std::numeric_limits::denorm_min()); + fn_push_back_pm(v, std::numeric_limits::min()); + fn_push_back_pm(v, std::numeric_limits::epsilon()); + const int test_exp_range = 5; + for (int i = -test_exp_range; i < (test_exp_range + 1); ++i) { + const int exponent = i * 3; + const T x = std::ldexp(1.0, exponent); + fn_push_back_pm(v, x); + } + fn_push_back_pm(v, std::numeric_limits::max()); + + // Real_concept doesn't have infinity + if (std::numeric_limits::has_infinity) { + fn_push_back_pm(v, std::numeric_limits::infinity()); + } + + std::sort(v.begin(), v.end()); + + const int num_zeros = std::count(v.begin(), v.end(), T(0.0)); + BOOST_CHECK_LE(2, num_zeros); // Positive and negative zero + + return v; + }; + + // Test for all combinations from the ladder + const auto v = fn_create_test_ladder(); + for (const T& x_i: v) { + for (const T& x_j: v) { + const T x_mid_ij = boost::math::tools::detail::Bisection::calc_midpoint(x_i, x_j); + + const T x_lo = std::min(x_i, x_j); + const T x_hi = std::max(x_i, x_j); + + BOOST_CHECK_LE(x_lo, x_mid_ij); + BOOST_CHECK_LE(x_mid_ij, x_hi); + } + } +} + +void test_bisect_all_cases() { + test_bisect(); + test_bisect(); + test_bisect(); + test_bisect(); + test_bisect(); + test_bisect(); + test_bisect(); + test_bisect(); +} + BOOST_AUTO_TEST_CASE( test_main ) { @@ -780,6 +862,8 @@ BOOST_AUTO_TEST_CASE( test_main ) test_cusp(); test_solved_at_bound(); + test_bisect_all_cases(); + // bug reports: boost::math::skew_normal_distribution<> dist(2.0, 1.0, -2.5); BOOST_CHECK(boost::math::isfinite(quantile(dist, 0.075))); From 07271f30d676c5f64e5b4e7d11ab05fc9cf602b4 Mon Sep 17 00:00:00 2001 From: Ryan Date: Fri, 4 Aug 2023 17:15:03 -0400 Subject: [PATCH 18/18] added () around min and max --- include/boost/math/tools/roots.hpp | 4 ++-- test/test_roots.cpp | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index 4991333926..6f096c5705 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -523,7 +523,7 @@ namespace detail { template static typename std::enable_if::is_specialized, T>::type do_solve(PosFinite x, EqInf X) { - return do_solve(x, PosFinite(std::numeric_limits::max())); + return do_solve(x, PosFinite((std::numeric_limits::max)())); } template static typename std::enable_if::is_specialized, T>::type @@ -539,7 +539,7 @@ namespace detail { const U denorm_min = std::numeric_limits::denorm_min(); if (denorm_min != 0) { return denorm_min; } - const U min = std::numeric_limits::min(); + const U min = (std::numeric_limits::min)(); if (min != 0) { return min; } BOOST_MATH_ASSERT(false && "denorm_min and min are both zero."); diff --git a/test/test_roots.cpp b/test/test_roots.cpp index 0c8e17e0cb..ac195e3985 100644 --- a/test/test_roots.cpp +++ b/test/test_roots.cpp @@ -804,7 +804,7 @@ void test_bisect() { fn_push_back_pm(v, T(0.0)); fn_push_back_pm(v, std::numeric_limits::denorm_min()); - fn_push_back_pm(v, std::numeric_limits::min()); + fn_push_back_pm(v, (std::numeric_limits::min)()); fn_push_back_pm(v, std::numeric_limits::epsilon()); const int test_exp_range = 5; for (int i = -test_exp_range; i < (test_exp_range + 1); ++i) { @@ -812,7 +812,7 @@ void test_bisect() { const T x = std::ldexp(1.0, exponent); fn_push_back_pm(v, x); } - fn_push_back_pm(v, std::numeric_limits::max()); + fn_push_back_pm(v, (std::numeric_limits::max)()); // Real_concept doesn't have infinity if (std::numeric_limits::has_infinity) { @@ -833,8 +833,8 @@ void test_bisect() { for (const T& x_j: v) { const T x_mid_ij = boost::math::tools::detail::Bisection::calc_midpoint(x_i, x_j); - const T x_lo = std::min(x_i, x_j); - const T x_hi = std::max(x_i, x_j); + const T x_lo = (std::min)(x_i, x_j); + const T x_hi = (std::max)(x_i, x_j); BOOST_CHECK_LE(x_lo, x_mid_ij); BOOST_CHECK_LE(x_mid_ij, x_hi);