diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index 791ec1a7e5..6f096c5705 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 @@ -213,116 +214,845 @@ inline std::pair bisect(F f, T min, T max, Tol tol) noexcept(policies::is_ } -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()))) -{ - BOOST_MATH_STD_USING +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 }; + + // Stores x, f(x), and f'(x) from a function evaluation + template + class Data_X01 { + public: + template + 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(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_; } + + // 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_); } + + // 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); } + + bool is_same_sign_f0(const Data_X01& z) const { return sign(f0_) == sign(z.f0()); } + + 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; + } - 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<>()); - } +#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) + 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) + }; - T f0(0), f1, last_f0(0); - T result = guess; + // Stores x and dx for a potential next root finding evaluation + template + class X_DX { + public: + X_DX(T x, T dx) : x_(x), dx_(dx) {} - T factor = static_cast(ldexp(1.0, 1 - digits)); - T delta = tools::max_value(); - T delta1 = tools::max_value(); - T delta2 = tools::max_value(); + T x() const { return x_; } + T dx() const { return dx_; } + private: + T x_; + T dx_; + }; + + // Calculates the step with Newton's method + template + class StepNewton { + public: + T calc_dx(const Data_X01& z) const { return -z.f0() / z.f1(); } + }; + + ////// 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... + // + // 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. // - // 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. + // 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. // - T max_range_f = 0; - T min_range_f = 0; + // 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. That's what this namespace is for. + // + namespace Bisection { - std::uintmax_t count(max_iter); + ////// The Midpoint754 class ////// + // + // 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 + // 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. 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() + // 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() + // + // 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. + // + 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; + } -#ifdef BOOST_MATH_INSTRUMENT - std::cout << "Newton_raphson_iterate, guess = " << guess << ", min = " << min << ", max = " << max - << ", digits = " << digits << ", max_iter = " << max_iter << "\n"; -#endif + // Convert bits to float + static T uint_to_float(U bits) { + T x; + std::memcpy(&x, &bits, sizeof(T)); + return x; + } - 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); + 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); + + // 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; } - else - { - delta = f0 / f1; + }; // class Midpoint754 + + + template + class MidpointNon754 { + private: + 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"); + + 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); } + + // At least one is positive + if (0 < sx + sX) { return do_solve(x, X); } + + // At least one is negative + return -do_solve(-x, -X); } -#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. - delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; - // reset delta1/2 so we don't take this branch next time round: - delta1 = 3 * delta; - delta2 = 3 * delta; + + 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: + 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 value() const { return x_; } + + private: + T x_; + }; + + // Two unknowns + static T do_solve(T x, T X) { + if (X < x) { + return do_solve(X, x); + } + + 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); + } } - guess = result; - result -= delta; - if (result <= min) - { - delta = 0.5F * (guess - min); - result = guess - delta; - if ((result == min) || (result == max)) - break; + + // 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)); + } } - else if (result >= max) - { - delta = 0.5F * (guess - max); - result = guess - delta; - if ((result == min) || (result == max)) - break; + 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)); + } } - // Update brackets: - if (delta > 0) - { - max = guess; - max_range_f = f0; + + // 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)())); } - else - { - min = guess; - min_range_f = f0; + 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()); } - // - // Sanity check that we bracket the root: - // - if (max_range_f * min_range_f > 0) + + 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 + template + class Root1D_State { + private: + 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_eval + template + class Val {}; + + public: + 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) { - 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<>()); + update_eval_sign_dx(x); } - }while(count && (fabs(result * factor) < fabs(delta))); - max_iter -= count; + // Syntactic sugar for type dispatch of update_eval + template + void update_eval_sign_f0(Args... args) { update_eval(Val(), args...); } + template + void update_eval_sign_dx(Args... args) { update_eval(Val(), args...); } + template + void update_eval_not_last(Args... args) { update_eval(Val(), args...); } + template + void update_eval_l(Args... args) { update_eval(Val(), args...); } + template + void update_eval_h(Args... args) { update_eval(Val(), args...); } + + // Filter arguments into preferred format + template + void update_eval(Val vs, const T x) { update_eval(vs, calc_x_dx(x)); } + template + void update_eval(Val vs, const X_DX& x_dx) { update_eval(vs, eval_next(x_dx), x_dx.dx()); } + template + void update_eval(Val vs, const Data_X01& cx) { update_eval(vs, cx, calc_dx_true(cx)); } + + // Resolve SideSelect + 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_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_eval(is_bound_h, cx, 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_eval(is_bound_h, cx, dx); + } + void update_eval(Val, const Data_X01& cx, const T dx) { + update_eval(false, cx, dx); // Updates low bound + } + void update_eval(Val, const Data_X01& cx, const T dx) { + update_eval(true, cx, dx); // Updates high bound + } + + // 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; } + } + + void update_case_flag() { + if (flag() == CaseFlag::success) { return; } // Don't update if already successful + + // 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 (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 + // bounds of the root finding problem. + 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(); } + + int num_x_only() const { return data_.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 = 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 + // 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)); } + + // 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 - eval_last().x(); } #ifdef BOOST_MATH_INSTRUMENT - std::cout << "Newton Raphson required " << max_iter << " iterations\n"; + // Makes debugging easier if required + void print_status_flag() const { + switch (flag()) { + 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 print_step_size_history() const { dx_history_.print(); } #endif - return result; + // 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); } + void reset_dx_history() { dx_history_.reset(); } + + 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_eval(); } + + 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 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_eval_({Data_X01(xl), Data_X01(xh)}) + , v_is_x_only_({{true, true}}) + , ind_last_(true) {} + + // Updates the eval using the bool is_h as an index. + void set(const Data_X01& cx, bool is_h) { + v_eval_[is_h] = cx; + v_is_x_only_[is_h] = false; + ind_last_ = is_h; + } + + // Returns the x value of the only eval that contains an original bound + T get_bound_from_unused_eval() const { + BOOST_MATH_ASSERT(num_x_only() == 1); + return v_is_x_only_[0] ? eval_l().x() : eval_h().x(); + } + + 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_eval_; + std::array v_is_x_only_; + bool ind_last_; + }; + + // Stores the size of the last two steps. Calling reset will set all step sizes to infinity. + class StepSizeHistory { + public: + StepSizeHistory() { reset(); } + + void reset() { + // 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) { + dx_pp_ = dx_p_; + 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_; } + + private: + T dx_p_; + T dx_pp_; + }; + + // 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()); + } + + std::uintmax_t num_fn_evals() const { return max_iter_ - count_; } + + 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_; + }; + + // // // 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 SolverBase { + public: + // 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); + } + + // Maybe exit + if (is_exit_case_now(b, x_dx)) { break; } + + // Update eval + this->update_eval(b, x_dx); + } + + return x_dx.x(); + } + + 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(); + } + + bool is_solver_current(const Root1D_State& b) const { return b.flag() == flag_this(); } + + 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_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 evaluation (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_eval_not_last(next_limit); + b.update_case_flag(); + return next_limit; + } + + void update_eval(Root1D_State&b, const X_DX& x_dx) const override { + b.update_eval_sign_dx(x_dx); + b.update_case_flag(); + } + + 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(); + } + + CaseFlag flag_this() const override { return CaseFlag::case_1; } + }; + + // Base class for SolverCase2U and SolverCase2B + template + class SolverCase2Base : public SolverBase { + public: + X_DX calc_next_bisection(Root1D_State& b) const override { + const auto x_mid = Bisection::calc_midpoint(b.x_l(), b.x_h()); + + if (x_mid == b.x_l() || x_mid == b.x_h()) { + b.set_flag(this->flag_stall()); + } + + return b.calc_x_dx(x_mid); + } + + virtual CaseFlag flag_stall() const = 0; + }; + + // The solver for Case 2U + template + class SolverCase2U : public SolverCase2Base { + public: + 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.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_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.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 + } + + 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(); + } + + 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_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)) { + 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; } + }; + + // 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.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 + // 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 + 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); + + // Get CacheX + 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_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_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 + state.update_case_flag(); // Sets to case_2U, case_2B, or failure + }; + + auto fn_return = [&](const Root1D_State& s, const std::pair p) { + max_iter = s.num_fn_evals(); // Update max_iter + return p.second; + }; + + // Solve problem + const auto p = solve_for_state(state); + + // If local minimization was successful, return root + if (p.first) { return fn_return(state, p); } + + // Attempt to bracket root in the direction of the initial guess + 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(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.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 +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_3_attempts>(f, x, xl, xh, digits, max_iter); } template diff --git a/test/test_roots.cpp b/test/test_roots.cpp index e4e68e24b8..ac195e3985 100644 --- a/test/test_roots.cpp +++ b/test/test_roots.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -23,7 +24,9 @@ #include #include +#include #include +#include #define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \ {\ @@ -437,6 +440,129 @@ void test_beta(T, const char* /* name */) test_inverses(ibeta_large_data); } +void test_3_attempt() { + 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) { + // 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); + + // 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); + + // 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); + }; + + fn_test_bracket(quartic_pos); + 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() @@ -638,8 +764,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: - 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); // 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); @@ -649,11 +773,97 @@ 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 ) { test_beta(0.1, "double"); + test_3_attempt(); + 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)));