Root Finding — Numerical Root-Finding Algorithms

Overview

The calx 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 Type

1D root-finding functions return RootFindingResult<T>.

Member	Description
`value`	Approximate root
`converged`	Whether convergence was achieved
`iterations`	Number of iterations
`error_estimate`	Estimated error

Convergence is controlled by ConvergenceCriteria<T> (tolerance, max iterations).

1D: Bracketing Methods

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

Function	Convergence	Description
`bisection(f, a, b)`	Linear	Bisection method. Simplest and most reliable
`false_position(f, a, b)`	Superlinear	False position (Regula Falsi)
`ridders_method(f, a, b)`	$\sqrt{2}$	Ridders' method. Faster than bisection
`brent_method(f, a, b)`	Superlinear	Brent's method. Best practical choice
`toms748(f, a, b)`	$\approx 1.84$	TOMS Algorithm 748. Improved Brent

1D: Non-bracketing Methods

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

Function	Convergence	Description
`newton_raphson(f, f', x0)`	Quadratic	Newton-Raphson. Requires derivative
`secant_method(f, x0, x1)`	$\approx 1.618$	Secant method. No derivative needed

1D: Higher-Order Methods

Function	Convergence	Description
`halley_method(f, f', f'', x0)`	Cubic	Halley's method. Requires 2nd derivative
`schroder_method(f, f', f'', x0)`	Cubic	Schröder's method. Robust for multiple roots
`king_method(f, f', x0)`	4th order	King's method. 4th order with only 1st derivative
`popovski_method(f, f', x0)`	4th order	Popovski's method

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 calx::Vector (which is reserved for linear algebra).

Low Degree (Closed-Form)

Function	Description
`solveLinear(p)`	Degree 1: $ax + b = 0$
`solveQuadratic(p)`	Degree 2: discriminant-based (cancellation-safe)
`solveCubic(p)`	Degree 3: Cardano's formula + trigonometric method
`solveQuartic(p)`	Degree 4: Ferrari's method

High Degree (Iterative)

Function	Description
`jenkinsTraub(p)`	Jenkins-Traub. Most numerically stable
`laguerre(p)`	Laguerre's method. Cubic convergence
`durandKernerAberth(p)`	Durand-Kerner-Aberth. Simultaneous iteration for all roots
`bairstow(p)`	Bairstow's method. Finds complex roots using only real arithmetic

Unified Dispatcher

auto roots = findAllRoots(p);  // Automatically selects optimal algorithm by degree

findAllRoots selects the algorithm based on degree:

Degree 1–4: closed-form (solveLinear / Quadratic / Cubic / Quartic)
Degree 5+: Jenkins-Traub

N-D: Multivariate Root Finding

Find zeros of $\mathbf{F}: \mathbb{R}^n \to \mathbb{R}^n$.

Function	Description
`newton_nd(F, J, x0)`	Multivariate Newton. Requires Jacobian
`broyden(F, x0)`	Broyden's method. Jacobian-free (quasi-Newton)

Example

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

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

    // Polynomial: x^3 - 6x^2 + 11x - 6 = (x-1)(x-2)(x-3)
    Polynomial<double> p = {-6.0, 11.0, -6.0, 1.0};
    auto roots = findAllRoots(p);
    for (auto& r : roots)
        std::cout << r << std::endl;
    // 1, 2, 3

    // Degree 5+: Jenkins-Traub
    Polynomial<double> q = {1.0, 0.0, 0.0, 0.0, 0.0, 1.0};  // x^5 + 1
    auto roots5 = findAllRoots(q);
    for (auto& r : roots5)
        std::cout << r << 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