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<R>

The unified return type used by every 1D, polynomial-iterative, and N-D root finder. Defined in math/roots/root_finding_base.hpp. The template parameter R dispatches automatically: scalar T for 1D, vector V for N-D.

template<typename R>
struct RootFindingResult {
    using error_type = /* R if scalar, R::value_type if vector */;

    std::optional<R>  root;           // root approximation (scalar or vector) on success, std::nullopt on failure
    bool              converged;      // true when convergence criteria were met
    size_t            iterations;     // number of iterations actually performed
    error_type        error_estimate; // |residual| / |step| / bracket-width estimate
};

1D solvers return RootFindingResult<double>, N-D solvers return RootFindingResult<Vector<double>>.

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';
}

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:

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;
// Run output: iter=42  root=1.41421356237311
//   matches the true sqrt(2)=1.41421356237310 to 12 digits

regula_falsi

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

Behaviour: The basic false-position method (Regula Falsi). The line through the bracket endpoints meets the x-axis at $c$, which becomes the next trial point; the endpoint sharing the sign of $f(c)$ is replaced by $c$, keeping a sign-changing bracket. Convergence is always guaranteed (the root stays bracketed), but when $f$ is convex (or concave) one endpoint stays fixed — stagnation — and convergence can degrade to linear. For the stagnation-free production version see illinois_method.

Parameters:

Use: the plain Regula Falsi, for teaching/reference. In practice prefer the stagnation-free variants illinois_method (Illinois), king_method (Pegasus family), or anderson_bjork_method.

illinois_method

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

Behaviour: The basic false position (see regula_falsi) with the Illinois correction. Plain Regula Falsi can stagnate on convex/concave functions, degrading to linear order; when one endpoint fails to update for two consecutive steps, scaling fa or fb by weight $1/2$ avoids stagnation and restores super-linear convergence. For even faster variants see anderson_bjork_method and king_method (Pegasus).

Parameters:

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).

Parameters:

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.

Parameters:

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;
// Run output: iter=26  root=0.739085133215161
//   matches the Dottie number 0.7390851332151607 to 14 digits

alefeld_potra_shi

template<concepts::OrderedField T>
RootFindingResult<T> alefeld_potra_shi(
    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$.

Parameters:

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

anderson_bjork_method

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

Behaviour: Improved Regula Falsi. When the same endpoint stagnates for two consecutive iterations, the stagnant function value is rescaled by the weight $m = 1 - f_c / f_b$ (falling back to $m = 1/2$ if $m \le 0$) before the next interpolation. This achieves faster super-linear convergence than the Illinois method (fixed weight $1/2$), comparable to King-Pegasus.

Parameters:

Use when: a fast bracketing super-linear method is desired over plain regula_falsi. Reference: Anderson & Björck, "A new high order method of regula falsi type for computing a root of an equation", BIT (1973).

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:

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;
// Run output: 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:

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;
// Run output: iter=6  root=0.739085133215161
//   matches the Dottie number 0.7390851332151607 to 14 digits

fixed_point_iteration

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

Behaviour: iterates $x_{n+1} = g(x_n)$ to find a fixed point $x^* = g(x^*)$ (i.e. a solution of $x = g(x)$). Note that $g$ is the iteration map itself, not the $f$ of $f(x)=0$. If $g$ is a contraction near $x^*$ ($|g'(x^*)| < 1$) it converges linearly; in particular it is at least quadratic when $g'(x^*) = 0$. It diverges when $|g'(x^*)| > 1$.

Parameters:

Note: when $g$ is not a contraction ($|g'(x^*)| \ge 1$) it diverges and returns converged = false after reaching the iteration limit (the best estimate is still returned). newton_raphson is the special case of fixed-point iteration with $g(x) = x - f(x)/f'(x)$.

// g(x) = cos(x), fixed point = Dottie number
auto r = fixed_point_iteration<double>(
    [](double x) { return std::cos(x); }, 0.0);
if (r.converged) std::cout << *r.root;
// g(x) = cos x is a contraction (|g'(x)| = |sin x| < 1), so it converges to x* = 0.7390851332151607

steffensen_method

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

How it works: update $x_{n+1} = x_n - \frac{f(x_n)^2}{f(x_n + f(x_n)) - f(x_n)}$. The denominator $f(x+f)-f(x)$ approximates $f(x)\,f'(x)$, giving Newton's quadratic convergence with no derivative (Aitken $\Delta^2$ acceleration applied to fixed-point iteration). Unlike the secant method it keeps only one point. Two $f$ evaluations per iteration.

Parameters:

Note: stalls when the denominator is near zero → converged=false. No global convergence guarantee; diverges if the start is far from the root.

// f(x) = x^2 - 2, root = √2, derivative-free
auto r = steffensen_method<double>(
    [](double x) { return x * x - 2.0; }, 1.5);
if (r.converged) std::cout << *r.root;
// Result: iter=4  root=1.414213562373095

sidi_method

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

How it works: uses the derivative at the newest point of the Newton interpolating polynomial $p$ through the last $K+1$ points: $x_{n+1} = x_n - f(x_n)/p'(x_n)$. With memory $=2$ it matches the secant method ($\varphi \approx 1.618$); larger values raise the order of convergence (toward $2$). Derivative-free.

Parameters:

// f(x) = x^2 - 2, root = √2, quadratic interpolation (memory=3)
auto r = sidi_method<double>(
    [](double x) { return x * x - 2.0; }, 1.0, 2.0,
    ConvergenceCriteria<double>(), 3);
if (r.converged) std::cout << *r.root;
// Result: iter=6  root=1.414213562373095

inverse_quadratic_interpolation

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

How it works: evaluates the quadratic interpolant of the inverse function $x = g(f)$ through 3 points at $f = 0$ ($x = \sum_i x_i \prod_{j \neq i} f_j/(f_i - f_j)$). This is Brent's method's internal step made available standalone. Derivative-free, super-linear ($\approx 1.84$), but unstable when the function values cluster, so it is normally paired with a safeguard (Brent).

Parameters:

Note: interpolation breaks down when the 3 $f$ values cluster → converged=false.

// f(x) = x^2 - 2, root = √2
auto r = inverse_quadratic_interpolation<double>(
    [](double x) { return x * x - 2.0; }, 1.0, 1.3, 1.7);
if (r.converged) std::cout << *r.root;
// Result: iter=5  root=1.414213562373095

muller_method

template<concepts::OrderedField T>
RootFindingResult<Complex<T>> muller_method(
    const std::function<Complex<T>(Complex<T>)>& f,
    Complex<T> z0, Complex<T> z1, Complex<T> z2,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

How it works: steps to the root (nearest the newest point) of the parabola through 3 points. Because it proceeds in complex arithmetic even when the discriminant is negative, it can reach complex roots from real-ish initial points (useful for polynomial complex roots). Order $\approx 1.84$, derivative-free. $f$ must be evaluable at complex arguments.

Parameters:

Note: the return type is RootFindingResult<Complex<T>> (the root is complex). Which complex root it converges to depends on the initial points.

// z^2 + 1 = 0 → ±i  (complex root from real-ish starts)
std::function<Complex<double>(Complex<double>)> g =
    [](Complex<double> z) { return z * z + Complex<double>(1.0); };
auto r = muller_method<double>(g, {1.0}, {0.5}, {0.0, 0.3});
if (r.converged) std::cout << r.root->re << " + " << r.root->im << "i";
// Result: iter=2  root=0 + 1i  (= i)

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:

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;
// Run output: 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.

Parameters:

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: After a Regula Falsi step, if the same endpoint stagnates the fixed endpoint's value is rescaled $f_1 \leftarrow f_1 \cdot f_2 / (f_2 + f)$ before re-interpolating. Lineage: Illinois (fixed weight $1/2$) → Pegasus (Dowell & Jarratt, 1972; weight $f_2/(f_2+f)$) → King (1973). King's method improves Pegasus by removing one or two of its slower substeps for higher asymptotic efficiency (order $\approx 1.839$ using only first-order divided differences), so the name king_method is not a mislabel — it denotes King's improvement of Pegasus. Faster super-linear convergence than the Illinois method.

References: R. F. King, "An improved Pegasus method for root finding", BIT 13, 423–427 (1973); M. Dowell & P. Jarratt, "The Pegasus method for computing the root of an equation", BIT 12, 503–508 (1972).

Parameters:

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.

Parameters:

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

householder_method

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

How it works: the order-$d{=}3$ member of the Householder family. It extends Newton ($d{=}1$) and Halley ($d{=}2$) to 4th-order convergence: $x_{n+1} = x_n - \frac{6 f f'^2 - 3 f^2 f''}{6 f'^3 - 6 f f' f'' + f^2 f'''}$. Requires the 1st–3rd derivatives.

Parameters:

Note: diverges when the denominator $6 f'^3 - 6 f f' f'' + f^2 f'''$ is near zero. Pays off over Newton/Halley when the 3rd derivative is cheap (symbolic / AD).

// f(x) = x^3 - 2, root = 2^(1/3)
auto r = householder_method<double>(
    [](double x) { return x*x*x - 2.0; },   // f
    [](double x) { return 3.0*x*x; },        // f'
    [](double x) { return 6.0*x; },          // f''
    [](double)   { return 6.0; },            // f'''
    1.5);
if (r.converged) std::cout << *r.root;
// Result: iter=3  root=1.259921049894873  (= 2^(1/3))

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.

Parameters:

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 << ' ';
// Run output: "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.

Parameters:

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);
// Run output: 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:

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.

Parameters:

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.

Parameters:

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

weierstrass

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

How it works: the classic Durand-Kerner (Weierstrass) method that updates all roots simultaneously. Unlike durandKernerAberth it carries no Aberth correction (plain simultaneous iteration, quadratic convergence): $z_i \leftarrow z_i - \frac{P(z_i)}{\mathrm{lc}\cdot\prod_{j \neq i}(z_i - z_j)}$. Starting from complex seeds $(0.4+0.9i)^i$, it converges to all roots (including complex conjugate pairs) when there are no repeated roots.

Parameters:

Note: near repeated/clustered roots the denominator $\prod_{j \neq i}(z_i - z_j)$ approaches zero and convergence slows. For more numerical robustness use durandKernerAberth (the Aberth-corrected variant).

// (x-2)(x^2+1) = x^3 - 2x^2 + x - 2 → 2, ±i
Polynomial<double> p({-2, 1, -2, 1});
auto roots = weierstrass(p);
// Result: roots = { 2, i, -i }

graeffe

template<typename T>
std::vector<Complex<T>> graeffe(
    const Polynomial<T>& poly,
    int maxSquarings = 6,
    T eps = std::numeric_limits<T>::epsilon() * T(100));

How it works: Graeffe root-squaring. It repeatedly forms $p(x)\,p(-x) = s(x^2)$, squaring the roots $r_i \to r_i^{2^m}$, and reads the moduli of separated roots from the monic coefficient ratios: $|r_i| \approx \left|\frac{a_{n-i}}{a_{n-i+1}}\right|^{1/2^m}$. Signs of real roots are recovered from $P(\pm|r_i|)$, then a few Newton steps polish each root. Because the coefficients grow doubly-exponentially, each squaring is renormalized to monic and squaring stops just before overflow.

Parameters:

Scope: for real roots with separated moduli. Equal-modulus roots, complex conjugate pairs and repeated roots cannot be separated by basic Graeffe (a historical method). For general complex roots use jenkinsTraub / durandKernerAberth / weierstrass.

// (x-1)(x-2)(x-4) = x^3 - 7x^2 + 14x - 8 → 4, 2, 1
Polynomial<double> p({-8, 14, -7, 1});
auto roots = graeffe(p);
// Result (descending modulus): roots = { 4, 2, 1 }

lehmerSchur

template<typename T>
std::vector<Complex<T>> lehmerSchur(
    const Polynomial<T>& poly,
    T tol = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxDescend = 200);

How it works: covers the disk containing all roots (Cauchy bound) with nine half-radius sub-disks (center + 8 around) and applies the Schur-Cohn test (Lehmer's unit-disk test) to each, descending into a sub-disk that contains a root. Once the radius is small enough it polishes the center with Newton's method and deflates, repeating degree-many times for all roots. The test iterates the exact Schur transform $T q = \overline{a_0}\,q - a_n\,q^*$ until the degree drops (a root lies inside whenever $|a_0|^2 - |a_n|^2 < 0$). Subdivision supplies the seed and Newton supplies the accuracy, so complex roots (including on the unit circle) are resolved.

Parameters:

Use when: you want geometric localization of roots, or to reliably separate roots on the unit circle / clustered complex roots. For the most numerically stable general solver use jenkinsTraub.

// x^4 - 1 = 0 → 1, -1, i, -i  (complex roots on the unit circle)
Polynomial<double> p({-1, 0, 0, 0, 1});
auto roots = lehmerSchur(p);
// Result: roots = { 1, -1, i, -i }

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 RootFindingResult<V> (a vector-typed specialisation of the same unified template used by 1D solvers).

Function reference

newton_raphson_nd

template<concepts::OrderedField T, typename V, typename M>
    requires concepts::VectorOf<V, T> && concepts::MatrixOf<M, T>
RootFindingResult<V> newton_raphson_nd(
    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:

Returns: RootFindingResult<V> 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_nd<double, Vector<double>, Matrix<double>>(F, J, x0);
if (r.converged && r.root) {
    const auto& v = *r.root;
    std::cout << v[0] << ", " << v[1];
}
// Run output: 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>
RootFindingResult<V> broyden_method(
    const std::function<V(const V&)>& F,
    const V& x0,
    std::optional<M> J0 = std::nullopt,
    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}$).

Parameters:

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

newton_krylov

template<concepts::OrderedField T>
RootFindingResult<Vector<T>> newton_krylov(
    const std::function<Vector<T>(const Vector<T>&)>& F,
    const Vector<T>& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>(),
    size_t krylov_restart = 30,
    T inner_tol = T(1e-3));

How it works: Jacobian-free Newton-Krylov (JFNK). The user supplies only $F$ (no analytic Jacobian and no explicit numerical Jacobian). Each Newton step solves $\mathbf{J}\,\boldsymbol\delta = -\mathbf{F}$ with GMRES (Krylov), but never forms the Jacobian: it uses only the directional finite difference $\mathbf{J}\,\mathbf{v} \approx (\mathbf{F}(\mathbf{x}+\varepsilon\mathbf{v}) - \mathbf{F}(\mathbf{x}))/\varepsilon$ (matrix-free). Backtracking damping adds global convergence.

Parameters:

Use when: large nonlinear systems where the Jacobian cannot be formed explicitly / is too large. broyden_method also needs only $F$, but this solves each step's linear system in a Krylov subspace, which is more memory-efficient in high dimensions.

// x² + y² = 1, x - y = 0 → (1/√2, 1/√2)   (only F is needed)
using V = Vector<double>;
auto F = [](const V& v) -> V {
    return V({v[0]*v[0] + v[1]*v[1] - 1.0, v[0] - v[1]});
};
auto r = newton_krylov<double>(F, V({0.5, 0.5}));
// Result: x = 0.7071067811865476, y = 0.7071067811865476

newton_fd

template<concepts::OrderedField T>
RootFindingResult<Vector<T>> newton_fd(
    const std::function<Vector<T>(const Vector<T>&)>& F,
    const Vector<T>& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>(),
    T fd_step = std::sqrt(std::numeric_limits<T>::epsilon()));

How it works: finite-difference Newton. The user supplies only $F$. Each Newton step builds the Jacobian by forward differences $J[:,j] = (F(x + h_j e_j) - F(x))/h_j$ as an explicit n×n matrix (sangi's Matrix) and solves $\mathbf{J}\,\boldsymbol\delta = -\mathbf{F}$ directly (LU). It is a thin wrapper that feeds a numerical Jacobian to newton_raphson_nd.

Parameters:

vs. newton_krylov: newton_fd holds $J$ explicitly (memory $O(n^2)$) and solves the linear system directly, so it is robust at small/medium scale. For large systems where $J$ cannot be formed, use newton_krylov (matrix-free).

// x² + y² = 1, x - y = 0 → (1/√2, 1/√2)   (only F is needed; J is built internally)
using V = Vector<double>;
auto F = [](const V& v) -> V {
    return V({v[0]*v[0] + v[1]*v[1] - 1.0, v[0] - v[1]});
};
auto r = newton_fd<double>(F, V({0.5, 0.5}));
// Result: x = 0.7071067811865476, y = 0.7071067811865476

levenberg_marquardt

template<concepts::OrderedField T, typename V, typename M>
RootFindingResult<V> 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.

Parameters:

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>
    requires concepts::VectorOf<V, T> && concepts::MatrixOf<M, T>
RootFindingResult<V> powell_hybrid(
    const std::function<V(const V&)>& F,
    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. The Jacobian is approximated internally by finite differences, so the user supplies only F. Same family as MINPACK's hybrd routine (SciPy's default fsolve).

Parameters:

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