Replace Boost Newton-Raphson with internal version

The Boost version's bisection fallback does not work properly.

boostorg/math#808

src/Domain/CoordinateMaps/TimeDependent/CubicScale.cpp

@@ -187,12 +187,9 @@ std::optional<std::array<double, Dim>> CubicScale<Dim>::inverse(
   // the CubicEquation solver: ~ 480 ns
   // boost implemented Newton-Raphson: ~ 280 ns
   // minimal Newton-Raphson from Numerical Recipes: ~ 255 ns
-  // Despite the minimal Newton-Raphson being more efficient than the boost
-  // version, here we utilize the boost implementation, as it includes
-  // additional checks for zero derivative, checks on bounds, and can implement
-  // bisection if necessary.
   const double scale_factor =
-      RootFinder::newton_raphson(cubic_and_deriv, initial_guess, 0.0, 1.0, 14) /
+      RootFinder::newton_raphson(cubic_and_deriv, initial_guess, 0.0, 1.0,
+                                 1.0e-14, 1.0e-14, 0.0) /
       target_dimensionless_radius;
 
   std::array<double, Dim> result{};

src/Domain/CoordinateMaps/UniformCylindricalEndcap.cpp

@@ -566,11 +566,11 @@ std::optional<std::array<double, 3>> UniformCylindricalEndcap::inverse(
           2.0 * function_at_rhobar_min / deriv_function_at_rhobar_min;
       try {
         // We have a good guess, so use newton_raphson
-        const size_t digits = 14;
         // use rhobar_min as maximum value, since that brackets the root.
         rhobar =
             RootFinder::newton_raphson(function_and_deriv_to_zero, rhobar_min,
-                                       new_rhobar_min, rhobar_min, digits);
+                                       new_rhobar_min, rhobar_min, 0.0, 1.0e-14,
+                                       0.0);
       } catch (std::exception&) {
         const double function_at_new_rhobar_min =
             function_to_zero(new_rhobar_min);
@@ -595,11 +595,11 @@ std::optional<std::array<double, 3>> UniformCylindricalEndcap::inverse(
           2.0 * function_at_rhobar_max / deriv_function_at_rhobar_max;
       try {
         // We have a good guess, so use newton_raphson
-        const size_t digits = 14;
         // use rhobar_max as minimum value, since that brackets the root.
         rhobar =
             RootFinder::newton_raphson(function_and_deriv_to_zero, rhobar_max,
-                                       rhobar_max, new_rhobar_max, digits);
+                                       rhobar_max, new_rhobar_max, 0.0, 1.0e-14,
+                                       0.0);
       } catch (std::exception&) {
         const double function_at_new_rhobar_max =
             function_to_zero(new_rhobar_max);

src/Evolution/Systems/RadiationTransport/M1Grey/M1Closure.cpp

@@ -64,9 +64,8 @@ void compute_closure_impl(
   static constexpr double small_velocity = 1.e-15;
   // Dimension of spatial tensors
   constexpr size_t spatial_dim = 3;
-  // Number of significant digits used in the rootfinding rooting
-  // used to find the closure factor
-  constexpr size_t root_find_number_of_digits = 6;
+  // Tolerance used in the root finding used to find the closure
+  // factor
   constexpr double root_find_tolerance = 1.e-6;
   Variables<
       tmpl::list<hydro::Tags::LorentzFactorSquared<DataVector>, MomentumSquared,
@@ -244,8 +243,8 @@ void compute_closure_impl(
         get(*closure_factor)[s] = 1.;
       } else {
         get(*closure_factor)[s] = RootFinder::newton_raphson(
-            zeta_j_sqr_minus_h_sqr, zeta_guess, 1.e-15, 1.,
-            root_find_number_of_digits);
+            zeta_j_sqr_minus_h_sqr, zeta_guess, 1.e-15, 1., 0.0,
+            root_find_tolerance, 0.0);
       }
 
       // Assemble output quantities:

src/NumericalAlgorithms/RootFinding/NewtonRaphson.hpp

@@ -6,9 +6,10 @@
 
 #pragma once
 
-#include <boost/math/tools/roots.hpp>
-#include <functional>
-#include <limits>
+#include <cmath>
+#include <cstddef>
+#include <optional>
+#include <utility>
 
 #include "DataStructures/DataVector.hpp"
 #include "Utilities/ErrorHandling/Exceptions.hpp"
@@ -17,102 +18,150 @@
 namespace RootFinder {
 /*!
  * \ingroup NumericalAlgorithmsGroup
- * \brief Finds the root of the function `f` with the Newton-Raphson method.
+ * \brief Finds the root of the function `f` with the Newton-Raphson
+ * method, falling back to bisection on poor convergence.
  *
- * `f` is a unary invokable that takes a `double` which is the current value at
- * which to evaluate `f`. `f` must return a `std::pair<double, double>` where
- * the first element is the function value and the second element is the
- * derivative of the function.  An example is below.
+ * `f` is a unary invokable that takes a `double` which is the current
+ * value at which to evaluate `f`. `f` must return a
+ * `std::pair<double, double>` where the first element is the function
+ * value and the second element is the derivative of the function.
+ * The method converges when the residual is smaller than or equal to
+ * \p residual_tolerance or when the step size (or bracket size when
+ * bisecting) is smaller than \p step_absolute_tolerance or the
+ * proposed result times \p step_relative_tolerance.
  *
  * \snippet Test_NewtonRaphson.cpp double_newton_raphson_root_find
  *
- * See the [Boost](http://www.boost.org/) documentation for more details.
- *
  * \requires Function `f` is invokable with a `double`
- * \note The parameter `digits` specifies the precision of the result in its
- * desired number of base-10 digits.
  *
  * \throws `convergence_error` if the requested precision is not met after
  * `max_iterations` iterations.
  */
 template <typename Function>
 double newton_raphson(const Function& f, const double initial_guess,
                       const double lower_bound, const double upper_bound,
-                      const size_t digits, const size_t max_iterations = 50) {
-  ASSERT(digits < std::numeric_limits<double>::digits10,
-         "The desired accuracy of " << digits
-                                    << " base-10 digits must be smaller than "
-                                       "the machine numeric limit of "
-                                    << std::numeric_limits<double>::digits10
-                                    << " base-10 digits.");
-
-  boost::uintmax_t max_iters = max_iterations;
-  // clang-tidy: internal boost warning, can't fix it.
-  const auto result = boost::math::tools::newton_raphson_iterate(  // NOLINT
-      f, initial_guess, lower_bound, upper_bound,
-      std::round(std::log2(std::pow(10, digits))), max_iters);
-  if (max_iters >= max_iterations) {
-    throw convergence_error(MakeString{}
-                            << "newton_raphson reached max iterations of "
-                            << max_iterations
-                            << " without converging. Best result is: " << result
-                            << " with residual " << f(result).first);
+                      const double residual_tolerance,
+                      const double step_absolute_tolerance,
+                      const double step_relative_tolerance,
+                      const size_t max_iterations = 50) {
+  // Check if a and b have the same sign.  Zero does not have the same
+  // sign as anything.
+  const auto same_sign = [](const double a, const double b) {
+    return a * b > 0.0;
+  };
+
+  ASSERT(residual_tolerance >= 0.0, "residual_tolerance must be non-negative.");
+  ASSERT(step_absolute_tolerance >= 0.0,
+         "step_absolute_tolerance must be non-negative.");
+  ASSERT(step_relative_tolerance >= 0.0,
+         "step_relative_tolerance must be non-negative.");
+
+  ASSERT(not same_sign(f(lower_bound).first, f(upper_bound).first) or
+         std::abs(f(upper_bound).first) < residual_tolerance or
+         std::abs(f(lower_bound).first) < residual_tolerance,
+         "Root not bracketed: "
+         "f(" << lower_bound << ") = " << f(lower_bound).first << "  "
+         "f(" << upper_bound << ") = " << f(upper_bound).first);
+
+  // Avoid evaluating at the bounds if we do not need to.  This will
+  // result in a non-optimal bracket until the region containing the
+  // root has been determined, but will avoid extra function
+  // evaluations if bisection is never needed.
+  std::optional<double> bracket_positive{};
+  std::optional<double> bracket_negative{};
+
+  double x = initial_guess;
+  double step = std::abs(upper_bound - lower_bound);
+
+  for (size_t i = 0; i < max_iterations; ++i) {
+    const auto [value, deriv] = f(x);
+    if (std::abs(value) <= residual_tolerance) {
+      return x;
+    }
+    if (value > 0.0) {
+      bracket_positive.emplace(x);
+    } else {
+      bracket_negative.emplace(x);
+    }
+
+    const auto current_bracket =
+        bracket_positive.has_value() and bracket_negative.has_value()
+            ? std::pair(*bracket_positive, *bracket_negative)
+            : std::pair(upper_bound, lower_bound);
+    if (same_sign((x - current_bracket.first) * deriv - value,
+                  (x - current_bracket.second) * deriv - value) or
+        std::abs(2.0 * value) > std::abs(step * deriv)) {
+      // Next guess not bracketed or converging slowly.  Perform a
+      // bisection.
+
+      // We need a bracket to bisect, so find any unknown bounds.
+      if (not(bracket_positive.has_value() and bracket_negative.has_value())) {
+        const double low_value = f(lower_bound).first;
+        if (std::abs(low_value) <= residual_tolerance) {
+          return lower_bound;
+        }
+
+        // Only one can be unset, because we set one of them right
+        // after we evaluated f(x).
+        if (not bracket_positive.has_value()) {
+          bracket_positive.emplace(low_value > 0.0 ? lower_bound : upper_bound);
+        } else {
+          bracket_negative.emplace(low_value > 0.0 ? upper_bound : lower_bound);
+        }
+      }
+      x = 0.5 * (*bracket_positive + *bracket_negative);
+      // Convenient value to allow combining the convergence checks
+      // for the two branches, and large enough that it won't trigger
+      // the slow-convergence check in the next iteration.
+      step = *bracket_positive - *bracket_negative;
+    } else {
+      // Convergence is fine.  Keep going with plain Newton-Raphson.
+      step = value / deriv;
+      x -= step;
+    }
+    if (std::abs(step) <= step_absolute_tolerance or
+        std::abs(step) <= step_relative_tolerance * std::abs(x)) {
+      return x;
+    }
   }
-  return result;
+
+  throw convergence_error(MakeString{}
+                          << "newton_raphson reached max iterations of "
+                          << max_iterations
+                          << " without converging. Best result is: " << x
+                          << " with residual " << f(x).first);
 }
 
 /*!
  * \ingroup NumericalAlgorithmsGroup
  * \brief Finds the root of the function `f` with the Newton-Raphson method on
  * each element in a `DataVector`.
  *
- * `f` is a binary invokable that takes a `double` as its first argument and a
- * `size_t` as its second. The `double` is the current value at which to
- * evaluate `f`, and the `size_t` is the current index into the `DataVector`s.
- *  `f` must return a `std::pair<double, double>` where the first element is
- * the function value and the second element is the derivative of the function.
- * Below is an example of how to root find different functions by indexing into
- * a lambda-captured `DataVector` using the `size_t` passed to `f`.
+ * Calls the `double` version of `newton_raphson` for each element of
+ * the input `DataVector`s, passing an additional `size_t` to \p f
+ * indicating the current index into the `DataVector`s.
  *
  * \snippet Test_NewtonRaphson.cpp datavector_newton_raphson_root_find
  *
- * See the [Boost](http://www.boost.org/) documentation for more details.
- *
  * \requires Function `f` be callable with a `double` and a `size_t`
- * \note The parameter `digits` specifies the precision of the result in its
- * desired number of base-10 digits.
  *
  * \throws `convergence_error` if, for any index, the requested precision is not
  * met after `max_iterations` iterations.
  */
 template <typename Function>
 DataVector newton_raphson(const Function& f, const DataVector& initial_guess,
                           const DataVector& lower_bound,
-                          const DataVector& upper_bound, const size_t digits,
+                          const DataVector& upper_bound,
+                          const double residual_tolerance,
+                          const double step_absolute_tolerance,
+                          const double step_relative_tolerance,
                           const size_t max_iterations = 50) {
-  ASSERT(digits < std::numeric_limits<double>::digits10,
-         "The desired accuracy of " << digits
-                                    << " base-10 digits must be smaller than "
-                                       "the machine numeric limit of "
-                                    << std::numeric_limits<double>::digits10
-                                    << " base-10 digits.");
-  const auto digits_binary = std::round(std::log2(std::pow(10, digits)));
-
   DataVector result_vector{lower_bound.size()};
   for (size_t i = 0; i < result_vector.size(); ++i) {
-    boost::uintmax_t max_iters = max_iterations;
-    // clang-tidy: internal boost warning, can't fix it.
-    result_vector[i] = boost::math::tools::newton_raphson_iterate(  // NOLINT
-        [&f, i](double x) { return f(x, i); }, initial_guess[i], lower_bound[i],
-        upper_bound[i], digits_binary, max_iters);
-    if (max_iters >= max_iterations) {
-      throw convergence_error(MakeString{}
-                              << "newton_raphson reached max iterations of "
-                              << max_iterations
-                              << " without converging. Best result is: "
-                              << result_vector[i] << " with residual "
-                              << f(result_vector[i], i).first);
-    }
+    result_vector[i] = newton_raphson(
+        [&f, i](const double x) { return f(x, i); }, initial_guess[i],
+        lower_bound[i], upper_bound[i], residual_tolerance,
+        step_absolute_tolerance, step_relative_tolerance, max_iterations);
   }
   return result_vector;
 }

src/NumericalAlgorithms/Spectral/Legendre.cpp

@@ -134,8 +134,10 @@ compute_collocation_points_and_weights<Basis::Legendre, Quadrature::Gauss>(
             newton_raphson_step,
             // Initial guess
             -cos((2. * j + 1.) * M_PI / (2. * poly_degree + 2.)),
-            // Lower and upper bound, and number of desired base-10 digits
-            -1., 1., 14);
+            -cos((2. * j) * M_PI / (2. * poly_degree + 2.)),
+            -cos((2. * j + 2.) * M_PI / (2. * poly_degree + 2.)),
+            // Tolerances
+            0.0, 1.0e-14, 0.0);
         const LegendrePolynomialAndDerivative L_and_dL(poly_degree + 1,
                                                        logical_coord);
         x[j] = logical_coord;
@@ -223,8 +225,8 @@ std::pair<DataVector, DataVector> compute_collocation_points_and_weights<
             // Initial guess
             -cos((j + 0.25) * M_PI / poly_degree -
                  0.375 / (poly_degree * M_PI * (j + 0.25))),
-            // Lower and upper bound, and number of desired base-10 digits
-            -1., 1., 14);
+            // Lower and upper bound, and tolerances.
+            -1., 1., 0.0, 1.0e-14, 0.0);
         const EvaluateQandL q_and_L(poly_degree, logical_coord);
         x[j] = logical_coord;
         x[poly_degree - j] = -logical_coord;

src/PointwiseFunctions/AnalyticSolutions/NewtonianEuler/RiemannProblem.cpp

@@ -180,7 +180,7 @@ RiemannProblem<Dim>::RiemannProblem(
     try {
       pressure_star_ = RootFinder::newton_raphson(
           f_of_p_and_deriv, guess_for_pressure_star, pressure_lower,
-          pressure_upper, -log10(pressure_star_tol_));
+          pressure_upper, 0.0, pressure_star_tol_, 0.0);
     } catch (std::exception& exception) {
       ERROR(
           "Failed to find p_* with Newton-Raphson root finder. Got "

src/PointwiseFunctions/AnalyticSolutions/Xcts/ConstantDensityStar.cpp

@@ -40,9 +40,9 @@ double compute_alpha(const double density, const double radius) {
                 pow<2>(pow_2_one_plus_a_square)};
       },
       // Choose initial guess for no particular reason
-      2. * sqrt(5.), sqrt(5.), std::numeric_limits<double>::max(),
+      1.0 / alpha_source, sqrt(5.), 1.0 / alpha_source,
       // Choose a precision of 14 base-10 digits for no particular reason
-      14);
+      0.0, 1.0e-14, 0.0);
 }
 
 template <typename DataType>

src/PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.cpp

@@ -101,7 +101,8 @@ double kerr_schild_areal_radius_from_isotropic(const double isotropic_radius,
                 isotropic_radius,
             kerr_schild_isotropic_radius_from_areal_deriv(areal_radius, mass));
       },
-      isotropic_radius, 0., std::numeric_limits<double>::max(), 12);
+      isotropic_radius, isotropic_radius, isotropic_radius + mass, 0.0, 1.0e-12,
+      0.0);
 }
 
 DataVector kerr_schild_areal_radius_from_isotropic(
@@ -113,10 +114,8 @@ DataVector kerr_schild_areal_radius_from_isotropic(
                 isotropic_radius[i],
             kerr_schild_isotropic_radius_from_areal_deriv(areal_radius, mass));
       },
-      isotropic_radius, make_with_value<DataVector>(isotropic_radius, 0.),
-      make_with_value<DataVector>(isotropic_radius,
-                                  std::numeric_limits<double>::max()),
-      12);
+      isotropic_radius, isotropic_radius, isotropic_radius + mass, 0.0, 1.0e-12,
+      0.0);
 }
 }  // namespace
 

src/PointwiseFunctions/Hydro/EquationsOfState/Enthalpy.cpp

@@ -442,11 +442,11 @@ double Enthalpy<LowDensityEoS>::rest_mass_density_from_enthalpy(
           1.0 / density * evaluate_coefficients(derivative_coefficients_, x);
       return std::make_pair(f, df);
     };
-    const size_t digits = 14;
+    const double tolerance = 1.0e-14;
     const double intial_guess = 0.5 * (minimum_density_ + maximum_density_);
     const auto root_from_lambda = RootFinder::newton_raphson(
-        f_df_lambda, intial_guess, reference_density_, maximum_density_,
-        digits);
+        f_df_lambda, intial_guess, reference_density_, maximum_density_, 0.0,
+        tolerance, 0.0);
     return root_from_lambda;
   }
 }

src/PointwiseFunctions/Hydro/EquationsOfState/Spectral.cpp

@@ -281,10 +281,11 @@ double Spectral::rest_mass_density_from_enthalpy(
                         (density * density) * this->gamma(x);
       return std::make_pair(f, df);
     };
-    const size_t digits = 14;
+    const double tolerance = 1.0e-14;
     const double intial_guess = 0.5 * (reference_density_ + upper_density);
     const auto root_from_lambda = RootFinder::newton_raphson(
-        f_df_lambda, intial_guess, reference_density_, upper_density, digits);
+        f_df_lambda, intial_guess, reference_density_, upper_density, 0.0,
+        tolerance, 0.0);
     return root_from_lambda;
   }
 }