Root Finding — Numerical Root-Finding Algorithms

Overview

The sangi root-finding module is organized into three categories.

• 1D root finding — find real $x$ satisfying $f(x) = 0$. Over 10 methods from bisection to Brent
• Polynomial roots — find all roots (including complex) of $p(x) = 0$. Closed-form for low degree, iterative for higher
• N-D root finding — find vector $\mathbf{x}$ satisfying $\mathbf{F}(\mathbf{x}) = \mathbf{0}$

Build

// Include all modules at once
#include <math/roots/root_finding.hpp>

// Or include individually
#include <math/roots/root_finding_1d.hpp>   // 1D root finding
#include <math/roots/root_finding_nd.hpp>   // N-D root finding
#include <math/roots/polynomial_roots.hpp>  // Polynomial roots

Header-only. No library linkage required.

Result Types and Convergence Criteria

RootFindingResult<T>

Returned by all 1D root finders and polynomial iterative finders. Defined in math/roots/root_finding_base.hpp.

template<typename T>
struct RootFindingResult {
    std::optional<T> root;           // root approximation on success, std::nullopt on failure
    bool             converged;      // true when convergence criteria were met
    size_t           iterations;     // number of iterations actually performed
    T                error_estimate; // |residual| / |step| / interval-width estimate
};

Usage:

auto r = brent_method(f, a, b);
if (r.converged && r.root) {
    double x = *r.root;
    // ... use x
} else {
    // failed (no sign change / max-iter / numerical singularity / ...)
    std::cerr << "no root: iter=" << r.iterations
              << " err=" << r.error_estimate << '\n';
}

N-D root finders return SystemRootFindingResult<V, T> with the same shape (the root field is std::optional<V> instead of std::optional<T>).

ConvergenceCriteria<T>

Holds the convergence thresholds. Every root finder takes it as a trailing argument (defaulted).

template<concepts::OrderedField T>
struct ConvergenceCriteria {
    T      abs_ftol = std::numeric_limits<T>::epsilon() * 100;   // |f(x)|  absolute tol
    T      abs_xtol = std::numeric_limits<T>::epsilon() * 100;   // |Δx|    absolute tol
    T      rel_ftol = std::numeric_limits<T>::epsilon() * 1000;  // |f|     relative tol
    T      rel_xtol = std::numeric_limits<T>::epsilon() * 1000;  // Δx      relative tol
    size_t max_iterations = 100;                                  // iteration cap
};

Example:

// Loose criteria for fast cutoff
ConvergenceCriteria<double> loose;
loose.abs_ftol = 1e-6;
loose.max_iterations = 50;
auto r = newton_raphson(f, df, x0, loose);

// Tight criteria with arbitrary precision
ConvergenceCriteria<Float> tight;
tight.abs_ftol = Float("1e-50");
tight.max_iterations = 200;
auto r2 = brent_method<Float>(f, a, b, tight);

1D: Bracketing Methods

Require an interval $[a, b]$ where $f(a)$ and $f(b)$ have opposite signs. Guaranteed to converge.

Function reference

bisection

template<concepts::OrderedField T>
RootFindingResult<T> bisection(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: Halves the interval $[a, b]$ each step and keeps the side where $f$ changes sign. Linear convergence (error halves per iteration).

Parameters: f univariate real function / a, b bracket endpoints ($f(a) \cdot f(b) < 0$ required) / criteria convergence thresholds.

Returns: RootFindingResult<T>. If the interval has no sign change, returns converged = false, root = std::nullopt.

Use when: robustness is paramount and speed is secondary. The worst-case iteration count is predictable: log2((b-a)/tol).

auto r = bisection<double>([](double x){ return x*x - 2.0; }, 1.0, 2.0);
if (r.converged) std::cout << *r.root;
// Execution-verified (build & run): iter=42  root=1.41421356237311
//   matches the true sqrt(2)=1.41421356237310 to 12 digits

false_position

template<concepts::OrderedField T>
RootFindingResult<T> false_position(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: Linear interpolation through the bracket endpoints; the x-intercept is the next trial point. On concave functions one endpoint can stagnate (Illinois / Anderson-Björk style modifications are not currently implemented).

Parameters: same as bisection. $f(a) \cdot f(b) < 0$ required.

ridders_method

template<concepts::OrderedField T>
RootFindingResult<T> ridders_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: Combines the bracket midpoint $m$ with an exponentially weighted false-position step. Each iteration evaluates $f$ twice in exchange for $\sqrt{2}$ super-linear order (effective rate ~ 1.84).

Use when: faster than bisection, simpler than Brent. A good middle ground when function evaluations are cheap.

brent_method

template<concepts::OrderedField T>
RootFindingResult<T> brent_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: Brent (1973) hybrid that switches between inverse-quadratic interpolation, the secant method, and bisection based on internal state. Any non-shrinking step is forced to a bisection so both speed and guaranteed convergence are retained.

Use when: the recommended default for 1D root finding. This is the same algorithm as SciPy's brentq. As long as the endpoints have opposite signs, the function always returns a root.

auto r = brent_method<double>([](double x){ return std::cos(x) - x; }, 0.0, 1.0);
if (r.converged) std::cout << *r.root;
// Execution-verified (build & run): iter=26  root=0.739085133215161
//   matches the Dottie number 0.7390851332151607 to 14 digits

toms748

template<concepts::OrderedField T>
RootFindingResult<T> toms748(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: Alefeld-Potra-Shi (1995) TOMS Algorithm 748. Combines inverse cubic interpolation with double-length steps, achieving a higher effective rate (~ 1.84) than Brent on smooth $f$.

Use when: a guaranteed bracketing method that is competitive with Brent on smooth problems. Drop-in replacement for SciPy's brentq.

1D: Non-bracketing Methods

Start from initial value $x_0$ (and optionally derivatives). Faster convergence but initial-value dependent.

Function reference

newton_raphson

template<concepts::OrderedField T>
RootFindingResult<T> newton_raphson(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: $x_{n+1} = x_n - f(x_n) / f'(x_n)$. Quadratic convergence near a simple root ($|x_{n+1} - x^*| \le C |x_n - x^*|^2$); degrades to linear at multiple roots.

Parameters: f function / df derivative / x0 initial guess / criteria convergence settings.

Note: triggers a near-stationary detection when $|f'(x)| < \epsilon \cdot 10$ and returns converged = false. A poor initial guess can converge to a different root or diverge; prefer brent_method when a bracket is available.

// f(x) = x^2 - 2 by Newton's method
auto r = newton_raphson<double>(
    [](double x) { return x*x - 2.0; },
    [](double x) { return 2.0 * x; },
    1.0);
if (r.converged) std::cout << *r.root;
// Execution-verified (build & run): iter=5  root=1.4142135623731
//   matches sqrt(2)=1.41421356237310 to 15 digits — 5 iterations show quadratic convergence

secant_method

template<concepts::OrderedField T>
RootFindingResult<T> secant_method(
    const std::function<T(T)>& f,
    T x0, T x1,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: $x_{n+1} = x_n - f(x_n)\, (x_n - x_{n-1}) / (f(x_n) - f(x_{n-1}))$. Uses a finite-difference approximation of the derivative; slightly slower than Newton (golden-ratio order $\approx 1.618$).

Parameters: f function / x0, x1 two initial points (need not bracket the root, but should not be too close) / criteria convergence settings.

Note: stalls when $|f(x_n) - f(x_{n-1})| < \epsilon$ (slope ≈ 0); returns converged = false.

// f(x) = x - cos(x), root = Dottie number
auto r = secant_method<double>(
    [](double x) { return x - std::cos(x); }, 0.0, 1.0);
if (r.converged) std::cout << *r.root;
// Execution-verified (build & run): iter=6  root=0.739085133215161
//   matches the Dottie number 0.7390851332151607 to 14 digits

1D: Higher-Order Methods

Function reference

halley_method

template<concepts::OrderedField T>
RootFindingResult<T> halley_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: $x_{n+1} = x_n - 2 f f' / (2 f'^2 - f f'')$. Householder order 1, lifting Newton's quadratic convergence to cubic. Needs the 2nd derivative; roughly twice the cost per iteration of Newton but with fewer iterations.

Parameters: f, df, d2f function and its first two derivatives / x0 initial guess / criteria convergence settings.

Note: triggers near-singular detection when the denominator $2 f'^2 - f f''$ becomes too small.

// f(x) = x^3 - 2, root = 2^(1/3)
auto r = halley_method<double>(
    [](double x) { return x*x*x - 2.0; },
    [](double x) { return 3.0 * x * x; },
    [](double x) { return 6.0 * x; },
    1.0);
if (r.converged) std::cout << *r.root;
// Execution-verified (build & run): iter=3  root=1.25992104989487
//   matches 2^(1/3)=1.2599210498948732 to 14 digits — 3 iterations showcase cubic convergence

schroder_method

template<concepts::OrderedField T>
RootFindingResult<T> schroder_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: $x_{n+1} = x_n - f f' / (f'^2 - f f'')$. Schröder (1870) designed this to retain quadratic convergence even at multiple roots, whereas plain Newton degrades to linear.

Use when: dealing with multiple roots of polynomials; more stable than Halley in that case. For simple roots Halley (cubic order) is faster.

king_method

template<concepts::OrderedField T>
RootFindingResult<T> king_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: King-Pegasus (R.F.King, BIT 1973). After taking a Regula Falsi step, if the same endpoint stagnates the algorithm rescales $f_1 \leftarrow f_1 \cdot f_2 / (f_2 + f)$ before the next interpolation. Faster super-linear convergence than the Illinois method.

Parameters: f function / a, b endpoints ($f(a) \cdot f(b) < 0$ required) / criteria convergence settings.

popovski_method

template<concepts::OrderedField T>
RootFindingResult<T> popovski_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: Two-step update (Newton step + correction term) reaching 4th-order convergence using only the 1st derivative. King-Werner / Ostrowski family. Trades the 2nd derivative for 2 function evaluations per iteration.

Use when: the 2nd derivative is expensive (e.g. AD chains) but function evaluations are cheap.

Polynomial Roots

All functions return std::vector<Complex<T>> (real roots have imaginary part = 0). Roots are a collection of complex numbers, not elements of a vector space, so the result uses the STL std::vector rather than sangi::Vector (which is reserved for linear algebra). The input Polynomial<T> stores coefficients in ascending order (p[0] = constant term).

Low Degree (Closed-Form)

High Degree (Iterative)

Function reference

solveLinear / solveQuadratic / solveCubic / solveQuartic

template<typename T> std::vector<Complex<T>> solveLinear(const Polynomial<T>& p);
template<typename T> std::vector<Complex<T>> solveQuadratic(const Polynomial<T>& p);
template<typename T> std::vector<Complex<T>> solveCubic(const Polynomial<T>& p);
template<typename T> std::vector<Complex<T>> solveQuartic(const Polynomial<T>& p);

Behaviour: closed-form solutions for the respective degrees. Degree mismatch triggers an assertion.

Cancellation handling: solveQuadratic computes the larger-magnitude root directly and derives the other via Vieta's formula when the discriminant is small (Numerical Recipes recommendation; see also Okumura, Algorithm Encyclopedia, p.205).

// x^2 - 5x + 6 = (x-2)(x-3) = 0
Polynomial<double> p({6.0, -5.0, 1.0});  // ascending: c0=6, c1=-5, c2=1
auto roots = solveQuadratic(p);
for (auto& r : roots) std::cout << r.re << ' ';
// Execution-verified (build & run): "3 2" (matches true {2, 3} exactly)

jenkinsTraub

template<typename T>
std::vector<Complex<T>> jenkinsTraub(const Polynomial<T>& poly);

Behaviour: Jenkins-Traub (1970) RPOLY algorithm (three-stage shift-shift-shift). Computes roots one at a time and deflates the polynomial. The de-facto reference standard for the past several decades.

Use when: the recommended default for medium-degree polynomials (5-19). solvePolynomial dispatches here for that range.

// x^6 - 1 = 0 (the six 6th roots of unity)
Polynomial<double> p({-1.0, 0, 0, 0, 0, 0, 1.0});
auto roots = jenkinsTraub(p);
// Execution-verified (build & run): 6 roots = {±1, ±exp(±iπ/3)}
//   max |p(root)| = 4.97e-16 (16-digit precision)

laguerre

template<typename T>
std::vector<Complex<T>> laguerre(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

Behaviour: Laguerre + deflation. Each step computes $L(z) = -n / (G \pm \sqrt{(n-1)(n H - G^2)})$ where $G = p'/p$ and $H = G^2 - p''/p$. Globally convergent to all roots from any initial guess, with cubic order.

Parameters: eps convergence tolerance / maxIter iteration cap.

durandKernerAberth

template<typename T>
std::vector<Complex<T>> durandKernerAberth(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

Behaviour: Aberth-Ehrlich simultaneous iteration ($z_i \leftarrow z_i - (p(z_i)/p'(z_i)) / (1 - (p(z_i)/p'(z_i)) \cdot \sum_{j \ne i} 1/(z_i - z_j))$) followed by a Newton polish against the original polynomial. Because all roots are updated together, deflation error does not accumulate.

Use when: the recommended default for degree ≥ 20. solvePolynomial dispatches here for that range.

bairstow

template<typename T>
std::vector<Complex<T>> bairstow(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

Behaviour: divides the polynomial by a trial quadratic factor $x^2 + rx + s$ and uses Newton's method to refine $r, s$ until the remainder vanishes, then extracts the quadratic factor. Complex (conjugate-pair) roots emerge naturally from real-only arithmetic. Reference: Press et al., Numerical Recipes §9.5.

Use when: targeting embedded environments that prefer to avoid complex arithmetic, or when complex roots should be handled as paired units.

Unified Dispatcher

template<typename T>
std::vector<Complex<T>> solvePolynomial(
    const Polynomial<T>& p,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

solvePolynomial picks the best algorithm by degree:

• Degree 1: solveLinear (closed-form)
• Degree 2: solveQuadratic (closed-form, cancellation-safe)
• Degree 3: solveCubic (Cardano + trig)
• Degree 4: solveQuartic (Ferrari)
• Degree 5-19: jenkinsTraub (one root at a time + deflation, high accuracy)
• Degree ≥ 20: durandKernerAberth (simultaneous roots + Newton polish, avoids deflation drift)

N-D: Multivariate Root Finding

Find zeros of $\mathbf{F}: \mathbb{R}^n \to \mathbb{R}^n$. The return type is SystemRootFindingResult<V, T>, the vector-valued analogue of RootFindingResult.

Function reference

newton_raphson_system

template<concepts::OrderedField T, typename V, typename M>
    requires concepts::VectorOf<V, T> && concepts::MatrixOf<M, T>
SystemRootFindingResult<V, T> newton_raphson_system(
    const std::function<V(const V&)>& F,
    const std::function<M(const V&)>& J,
    const V& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: each iteration $\mathbf{x}_{n+1} = \mathbf{x}_n - J(\mathbf{x}_n)^{-1} F(\mathbf{x}_n)$. On a singular Jacobian the diagonal is regularised by abs_xtol * 10 and solve is retried.

Parameters: F system as a vector-valued function / J Jacobian / x0 initial vector / criteria convergence settings.

Returns: SystemRootFindingResult<V, T> with root of type std::optional<V>.

Convergence: $\|F(\mathbf{x})\|_2 < \mathrm{abs\_ftol}$ or $\|\Delta\mathbf{x}\|_2$ below the prescribed tolerance.

Example — find the intersection of the unit circle $x^2 + y^2 = 1$ with the line $y = x$. Pass sangi::Matrix<T>::solve() as the concept-satisfying solver by also including <math/linalg/solvers.hpp>.

#include <math/roots/root_finding.hpp>
#include <math/linalg/solvers.hpp>   // body of Matrix::solve() lives here

using namespace sangi;

auto F = [](const Vector<double>& v) {
    Vector<double> r(2);
    r[0] = v[0]*v[0] + v[1]*v[1] - 1.0;   // x^2 + y^2 - 1
    r[1] = v[0] - v[1];                   // x - y
    return r;
};
auto J = [](const Vector<double>& v) {
    Matrix<double> m(2, 2);
    m(0,0) = 2*v[0]; m(0,1) = 2*v[1];
    m(1,0) = 1.0;    m(1,1) = -1.0;
    return m;
};
Vector<double> x0(2); x0[0] = 0.5; x0[1] = 0.5;
auto r = newton_raphson_system<double, Vector<double>, Matrix<double>>(F, J, x0);
if (r.converged && r.root) {
    const auto& v = *r.root;
    std::cout << v[0] << ", " << v[1];
}
// Execution-verified (build & run): iter=5  root=(0.707106781186548, 0.707106781186548)
//   matches 1/sqrt(2)=0.7071067811865476 to 15 digits

broyden_method

template<concepts::OrderedField T, typename V, typename M>
SystemRootFindingResult<V, T> broyden_method(
    const std::function<V(const V&)>& F,
    const V& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: "good" Broyden quasi-Newton. Avoids re-computing the Jacobian by maintaining a rank-1 update ($B_{n+1} = B_n + \frac{(\Delta F - B_n \Delta \mathbf{x}) \Delta \mathbf{x}^\top}{\|\Delta \mathbf{x}\|^2}$).

Use when: the analytic Jacobian is unavailable or too expensive. $B_0$ defaults to the identity and is improved by the updates.

levenberg_marquardt

template<concepts::OrderedField T, typename V, typename M>
SystemRootFindingResult<V, T> levenberg_marquardt(
    const std::function<V(const V&)>& F,
    const std::function<M(const V&)>& J,
    const V& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: solves $(J^\top J + \lambda I) \Delta\mathbf{x} = -J^\top F$. The damping $\lambda$ is decreased on successful steps and increased on failed ones (Marquardt-Levenberg update rule). Interpolates between Newton and steepest descent.

Use when: the initial guess is far from a root, the Jacobian is near-singular, or the problem is best framed as residual minimisation rather than pure equation solving.

powell_hybrid

template<concepts::OrderedField T, typename V, typename M>
SystemRootFindingResult<V, T> powell_hybrid(
    const std::function<V(const V&)>& F,
    const std::function<M(const V&)>& J,
    const V& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

Behaviour: Powell (1970) dogleg trust-region. Follows a dogleg path that combines the Newton direction with the steepest-descent direction inside a trust region. Same family as MINPACK's hybrd routine (SciPy's default fsolve).

Use when: the recommended default for N-D root finding. Combines Newton-like local convergence with global convergence from the trust region.

Example

#include <math/roots/root_finding.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // 1D: find root of f(x) = x^2 - 2 using Brent's method
    auto result = brent_method<double>(
        [](double x) { return x * x - 2.0; },
        1.0, 2.0);
    if (result.converged && result.root)
        std::cout << "sqrt(2) = " << *result.root << std::endl;
    // Output: sqrt(2) = 1.41421...

    // Polynomial: x^3 - 6x^2 + 11x - 6 = (x-1)(x-2)(x-3)
    Polynomial<double> p({-6.0, 11.0, -6.0, 1.0});  // ascending
    auto roots = solvePolynomial(p);
    for (auto& r : roots)
        std::cout << r.re << std::endl;
    // Output: 1, 2, 3 (order is implementation-dependent)

    // Degree 5+: solvePolynomial dispatches to jenkinsTraub automatically
    Polynomial<double> q({1.0, 0.0, 0.0, 0.0, 0.0, 1.0});  // x^5 + 1
    auto roots5 = solvePolynomial(q);
    for (auto& r : roots5)
        std::cout << r.re << " + " << r.im << "i" << std::endl;
}

Related Mathematical Background

The following articles explain the mathematical concepts underlying the Roots module.

• Bisection Method — The most fundamental root-finding algorithm
• Newton's Method — Quadratically convergent iteration
• Secant Method — Derivative-free superlinear convergence
• Advanced Polynomial Root-Finding — Aberth-Ehrlich, companion matrix methods
• Nth Root via Precision Doubling — Multi-precision nth root computation
• Zimmermann Recursive Square Root — Fast multi-precision square root