Root Finding

Why Root Finding Matters

A sensor's calibration model says that the measured response is $g(x)=x^2$; the desired response is $2$. Finding the input means solving $f(x)=g(x)-2=0$. A table can locate the answer roughly, but an iterative solver can refine it to a specified tolerance while reporting what evidence supports the answer. The same pattern appears when locating a simulation event, inverting a monotone physical model, or solving a nonlinear equation inside a larger algorithm.

This lesson concerns a real scalar equation $f(x)=0$. It develops bisection, Newton–Raphson, secant, regula falsi, and Brent–Dekker. You will distinguish methods that preserve a sign-changing interval from methods that extrapolate from local information, trace their different updates, implement all five, and check both residual and root error on a known example. A system $F(\mathbf{x})=\mathbf{0}$ needs a different solver; scalar Newton requires no matrix elimination or pivoting.

Core Concepts and Terminology

The input is a callable real function, starting information, an absolute position tolerance $\tau_x\gt0$, a residual tolerance $\tau_f\ge0$, and an iteration budget $K$. Newton additionally needs $f'$. The output should include the estimate, its evaluated residual, how many new trial points were evaluated, and a termination reason. A floating-point result without a reason conceals failures such as stagnation or an exhausted budget.

For a bracket, continuity on the entire closed interval and opposite nonzero endpoint signs imply at least one root by the intermediate value theorem. They do not imply uniqueness. A strictly monotone function has at most one root. The implementation compares signs directly, using (fa < 0) != (fb < 0) after endpoint checks, rather than multiplying residuals and risking overflow or underflow.

For $f(x)=x^2-2$ on $[1,2]$, the endpoint residuals are $-1$ and $2$. The midpoint is $3/2$, the line through the two endpoint samples crosses the axis at $4/3$, and the tangent at $2$ crosses at $3/2$. These are different geometric constructions even when two happen to produce the same first point.

A small residual alone does not bound position error without additional information. If $\lvert f'(x)\rvert\ge\mu\gt0$ throughout the segment joining $\hat{x}$ to a root, the mean value theorem gives $\lvert\hat{x}-r\rvert\le\lvert f(\hat{x})\rvert/\mu$. For a nearly flat function, the same residual can correspond to a much larger position error.

How Root Finding Works, Step by Step

1. Validate the numerical contract. Require finite starts and tolerances, $\tau_x\gt0$, $\tau_f\ge0$, and an integer $K\ge0$. The programs also require the initial difference $b-a$ to be representable and finite. Check every function and derivative evaluation for a finite result. Continuity and derivative correctness are mathematical assumptions the program cannot establish by sampling.

2. Check initial residuals before rejecting a bracket. If an initial value already satisfies $\lvert f(x)\rvert\le\tau_f$, return it with a residual status. Otherwise bisection, regula falsi, and Brent require opposite endpoint signs. Secant needs two initial guesses without that sign requirement. Newton uses only the second start, called $b$ in the shared interface.

3. Propose a new trial point. Each method supplies a different rule:

• Bisection: $s=a+(b-a)/2$. The midpoint halves the interval regardless of local curvature.
• Newton–Raphson: $s=x-f(x)/f'(x)$. Stop with a failure status if the derivative is zero or the proposed point is nonfinite; this open method can leave the initial interval.
• Secant: $s=b-f(b)(b-a)/(f(b)-f(a))$. Retain the last two iterates after evaluation, even when their signs agree. An equal or nonfinite residual difference prevents this update.
• Regula falsi: use the same line intersection as secant, but retain the two opposite-sign endpoints after evaluation. The code computes an algebraically equivalent normalized weight formula to reduce residual overflow.
• Brent–Dekker: use inverse quadratic interpolation when three stored function values are pairwise distinct, otherwise secant; reject an unsafe proposal and bisect. It retains a bracket throughout.

For the normalized false-position formula, let $M=\max(\lvert f(a)\rvert,\lvert f(b)\rvert)$, $p=\lvert f(a)\rvert/M$, and $q=\lvert f(b)\rvert/M$. Then $s=a+(b-a)p/(p+q)$. Opposite signs make this a convex combination of the endpoints in exact arithmetic. A rounded endpoint or a nonfinite proposal triggers midpoint fallback in the bracketing implementations.

For three samples $(a,f_a)$, $(b,f_b)$, $(c,f_c)$ with distinct residuals, inverse quadratic interpolation proposes

$$ s=\frac{a f_b f_c}{(f_a-f_b)(f_a-f_c)} +\frac{b f_a f_c}{(f_b-f_a)(f_b-f_c)} +\frac{c f_a f_b}{(f_c-f_a)(f_c-f_b)}. $$

The code scales all three ordinates by their maximum magnitude before using this formula. A zero denominator produced by rounding rejects the interpolation; scaling does not make every extreme input numerically safe.

4. Evaluate and update the appropriate state. After evaluating $f(s)$, return on the residual criterion if it passes. In a bracketing method, replace the endpoint whose sign agrees with $f(s)$; in secant, shift the two latest iterates; in Newton, replace the single iterate. The bracketed code swaps endpoint labels when necessary so $\lvert f(b)\rvert\le\lvert f(a)\rvert$. Thus $b$ is the better endpoint and can lie to the left of $a$.

5. Apply the method's stopping rule. A bracket width at most $\tau_x$ permits returning either endpoint with position error at most $\tau_x$ to some enclosed root. The code returns the better endpoint, rather than the midpoint, so its certificate uses the full width. An open-method step at most $\tau_x$ produces a step status; that status is a reason to inspect the residual, not proof of a root. A nonrepresentable update, an unchanged trial point, or a zero slope denominator yields stalled; exhausting the budget yields limit.

The Brent implementation uses explicit history guards. Let $c$ be the previous best endpoint and $d$ its predecessor. A proposal must lie strictly between $b$ and $a+(b-a)/4$. Its distance from $b$ must be less than half the previous comparison distance: $\lvert b-c\rvert$ after bisection, or $\lvert c-d\rvert$ after interpolation. If that comparison distance is already at most $\tau_x$, bisect instead.

This teaching implementation adds a progress watchdog to those Brent–Dekker guards. Record a checkpoint width; reset the checkpoint whenever the bracket halves or a midpoint is used. After two accepted interpolation steps without halving the checkpoint width, force a midpoint next. Consequently, every block of at most three bracket updates halves a checkpoint interval in exact arithmetic. This conservative variant makes the worst-case bound easy to verify; its exact iterations need not match a library's Brent routine.

flowchart TD
    A["Validate inputs and evaluate starts"] --> B{"Initial residual accepted?"}
    B -->|Yes| R["Return estimate and termination reason"]
    B -->|No| C{"Bracketing method?"}
    C -->|Yes| D["Require opposite signs; check width"]
    C -->|No| O{"Budget remains?"}
    D --> F{"Width accepted?"}
    F -->|Yes| R
    F -->|No| P{"Budget remains?"}
    O -->|No| R
    O -->|Yes| E["Use Newton or secant proposal"]
    P -->|No| R
    P -->|Yes| G["Choose trial; safeguard interpolation"]
    E --> H{"Finite new point and valid slope?"}
    G --> H
    H -->|No| R
    H -->|Yes| I["Evaluate trial residual"]
    I --> J{"Residual accepted?"}
    J -->|Yes| R
    J -->|No| L["Update bracket or open iterates"]
    L --> M{"Width or open step accepted?"}
    M -->|Yes| R
    M -->|No| N{"Budget remains?"}
    N -->|No| R
    N -->|Yes| C

The diagram assumes each evaluation is finite. A nonfinite function or derivative value raises an exception, rather than continuing with invalid signs. The budget counts newly evaluated trial points, excluding initial function values and derivative calls.

Worked Example with a Full Trace

Solve $x^2-2=0$ using the initial interval $[1,2]$. First trace bisection all the way to its configured stop with $\tau_x=1/8$, $\tau_f=0$, and $K=20$. The exact root is $r=\sqrt{2}\approx1.4142135623730951$. Endpoint labels below are written in increasing coordinate order for readability, independently of the code's residual-based swaps.

Now stop: the better endpoint is $11/8$, because $7/64\lt1/4$. The result is bracket after three trial evaluations. Its actual position error is about $0.039214$, below the certified $0.125$ width, even though the residual $7/64$ is not close to zero. Returning the midpoint $23/16$ would give a half-width bound, but would be a different return policy.

To trace the other update rules completely without dozens of rows, use $\tau_x=10^{-12}$ and the deliberately loose residual target $\tau_f=10^{-2}$. Newton starts at $2$; the other methods start at $1,2$. In the following table every new trial is listed, and each run stops on the final listed residual.

Secant and regula falsi agree for the first two trials. After that, secant uses $4/3,7/5$, which both have negative residuals; regula falsi keeps $7/5,2$, preserving opposite signs. Brent first uses secant, then inverse quadratic interpolation, then secant because the stored third ordinate duplicates an endpoint ordinate. In this trace, no midpoint fallback is needed.

Implementation in Python, C++, and Java

Each program implements the same five methods and the same guards using ordinary double-precision arithmetic. The nested arrays in the driver contain starting pairs; root finding itself has no matrix rows to pivot. The explicit state swaps serve the bracket's best-endpoint convention. The programs use no numerical library and each panel runs independently.

The shared solver takes a method name, $f$, $f'$, two starts, tolerances, and a budget. Only Newton calls $f'$, and only Newton ignores the first start. Valid method strings are bisection, newton, secant, regula, and brent. A result contains x, fx, iterations, and status; the statuses are residual, bracket, step, stalled, and limit. Invalid inputs, a missing sign change after initial residual checks, or nonfinite function/derivative values raise an exception. Proposals that overflow are rejected or cause a stalled status, as appropriate to the method.

For all five sample runs, use $f(x)=x^2-2$, $f'(x)=2x$, starts $1,2$, $\tau_x=\tau_f=10^{-12}$, and $K=200$. This is a tighter configuration than the hand traces. The driver verifies both $\lvert f(\hat{x})\rvert\le10^{-10}$ and $\lvert\hat{x}-\sqrt{2}\rvert\le10^{-10}$. The second check is possible here because the reference root is known independently.

import math
from dataclasses import dataclass


@dataclass
class Result:
    x: float
    fx: float
    iterations: int
    status: str


def value(f, x):
    y = f(x)
    if not math.isfinite(y):
        raise ValueError("nonfinite function or derivative value")
    return y


def opposite(u, v):
    return (u < 0.0) != (v < 0.0)


def solve(method, f, df, a, b, xtol=1e-12, ftol=1e-12, limit=200):
    names = ("bisection", "newton", "secant", "regula", "brent")
    if (method not in names or not all(map(math.isfinite, (a, b, xtol, ftol)))
            or xtol <= 0.0 or ftol < 0.0 or not isinstance(limit, int) or limit < 0
            or not math.isfinite(b - a)):
        raise ValueError("invalid method, endpoints, tolerances, or budget")
    # Newton uses b as its sole initial guess; a is unused.
    if method == "newton":
        x, fx = b, value(f, b)
        if abs(fx) <= ftol:
            return Result(x, fx, 0, "residual")
        for k in range(1, limit + 1):
            slope = value(df, x)
            if slope == 0.0:
                return Result(x, fx, k - 1, "stalled")
            s = x - fx / slope
            if not math.isfinite(s) or s == x:
                return Result(x, fx, k - 1, "stalled")
            fs = value(f, s)
            if abs(fs) <= ftol:
                return Result(s, fs, k, "residual")
            if abs(s - x) <= xtol:
                return Result(s, fs, k, "step")
            x, fx = s, fs
        return Result(x, fx, limit, "limit")

    fa, fb = value(f, a), value(f, b)
    if abs(fa) <= ftol:
        return Result(a, fa, 0, "residual")
    if abs(fb) <= ftol:
        return Result(b, fb, 0, "residual")
    if method == "secant":
        for k in range(1, limit + 1):
            den = fb - fa
            if den == 0.0 or not math.isfinite(den):
                return Result(b, fb, k - 1, "stalled")
            s = b - fb * ((b - a) / den)
            if not math.isfinite(s) or s == b:
                return Result(b, fb, k - 1, "stalled")
            fs = value(f, s)
            if abs(fs) <= ftol:
                return Result(s, fs, k, "residual")
            if abs(s - b) <= xtol:
                return Result(s, fs, k, "step")
            a, fa, b, fb = b, fb, s, fs
        return Result(b, fb, limit, "limit")

    if a == b or not opposite(fa, fb):
        raise ValueError("bracketing methods require opposite endpoint signs")
    if abs(fa) < abs(fb):
        a, b, fa, fb = b, a, fb, fa
    if abs(b - a) <= xtol:
        return Result(b, fb, 0, "bracket")
    c, fc, d = a, fa, a
    bisected, interpolations = True, 0
    checkpoint = abs(b - a)
    for k in range(1, limit + 1):
        if abs(b - a) <= xtol:
            return Result(b, fb, k - 1, "bracket")
        mid = a + (b - a) / 2.0
        s, use_mid = mid, method == "bisection"
        if method == "regula":
            scale = max(abs(fa), abs(fb))
            p, q = abs(fa) / scale, abs(fb) / scale
            s = a + (b - a) * (p / (p + q))
        elif method == "brent":
            if fa != fc and fb != fc:
                # Normalize ordinates before inverse quadratic interpolation.
                scale = max(abs(fa), abs(fb), abs(fc))
                u, v, w = fa / scale, fb / scale, fc / scale
                den0, den1, den2 = (u-v)*(u-w), (v-u)*(v-w), (w-u)*(w-v)
                s = math.nan
                if den0 != 0.0 and den1 != 0.0 and den2 != 0.0:
                    s = a*(v*w/den0) + b*(u*w/den1) + c*(u*v/den2)
            else:
                scale = max(abs(fa), abs(fb))
                u, v = fa / scale, fb / scale
                s = b - (v / (v - u)) * (b - a)
            gate = a + 0.25 * (b - a)
            inside = min(gate, b) < s < max(gate, b)
            old = abs(b - c) if bisected else abs(c - d)
            use_mid = (not inside or abs(s-b) >= 0.5*old
                       or old <= xtol or interpolations >= 2)
        if (use_mid or not math.isfinite(s)
                or not min(a, b) < s < max(a, b)):
            s, use_mid = mid, True
        if s == a or s == b:
            return Result(b, fb, k - 1, "stalled")
        fs = value(f, s)
        if abs(fs) <= ftol:
            return Result(s, fs, k, "residual")
        d, c, fc = c, b, fb
        if opposite(fa, fs):
            b, fb = s, fs
        else:
            a, fa = s, fs
        if abs(fa) < abs(fb):
            a, b, fa, fb = b, a, fb, fa
        bisected = use_mid
        if use_mid or abs(b - a) <= 0.5 * checkpoint:
            interpolations, checkpoint = 0, abs(b - a)
        else:
            interpolations += 1
        if abs(b - a) <= xtol:
            return Result(b, fb, k, "bracket")
    return Result(b, fb, limit, "limit")


def main():
    def f(x):
        return x*x - 2.0

    def df(x):
        return 2.0*x

    # Each row is [left/start0, right/start1]; no matrix pivoting is needed.
    starts = [[1.0, 2.0] for _ in range(5)]
    methods = ["bisection", "newton", "secant", "regula", "brent"]
    ok = True
    for method, pair in zip(methods, starts):
        r = solve(method, f, df, pair[0], pair[1])
        good = abs(r.fx) <= 1e-10 and abs(r.x - math.sqrt(2.0)) <= 1e-10
        ok = ok and good and r.status in ("residual", "bracket")
        print(f"{method}: root = {r.x:.6f}, check = {'PASS' if good else 'FAIL'}")
    print(f"all methods = {'PASS' if ok else 'FAIL'}")


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>
using Fn = std::function<double(double)>;
struct Result { double x, fx; int iterations; std::string status; };

double value(const Fn& f, double x) {
    double y = f(x);
    if (!std::isfinite(y)) throw std::invalid_argument("nonfinite function or derivative value");
    return y;
}
bool opposite(double u, double v) { return (u < 0.0) != (v < 0.0); }

Result solve(const std::string& method, const Fn& f, const Fn& df,
             double a, double b, double xtol=1e-12, double ftol=1e-12, int limit=200) {
    bool known = method=="bisection" || method=="newton" || method=="secant"
                 || method=="regula" || method=="brent";
    if (!known || !std::isfinite(a) || !std::isfinite(b) || !std::isfinite(xtol)
        || !std::isfinite(ftol) || xtol<=0 || ftol<0 || limit<0 || !std::isfinite(b-a))
        throw std::invalid_argument("invalid method, endpoints, tolerances, or budget");
    // Newton uses b as its sole initial guess; a is unused.
    if (method=="newton") {
        double x=b, fx=value(f,b);
        if (std::abs(fx)<=ftol) return {x,fx,0,"residual"};
        for (int k=1; k<=limit; ++k) {
            double slope=value(df,x);
            if (slope==0) return {x,fx,k-1,"stalled"};
            double s=x-fx/slope;
            if (!std::isfinite(s) || s==x) return {x,fx,k-1,"stalled"};
            double fs=value(f,s);
            if (std::abs(fs)<=ftol) return {s,fs,k,"residual"};
            if (std::abs(s-x)<=xtol) return {s,fs,k,"step"};
            x=s; fx=fs;
        }
        return {x,fx,limit,"limit"};
    }
    double fa=value(f,a), fb=value(f,b);
    if (std::abs(fa)<=ftol) return {a,fa,0,"residual"};
    if (std::abs(fb)<=ftol) return {b,fb,0,"residual"};
    if (method=="secant") {
        for (int k=1; k<=limit; ++k) {
            double den=fb-fa;
            if (den==0 || !std::isfinite(den)) return {b,fb,k-1,"stalled"};
            double s=b-fb*((b-a)/den);
            if (!std::isfinite(s) || s==b) return {b,fb,k-1,"stalled"};
            double fs=value(f,s);
            if (std::abs(fs)<=ftol) return {s,fs,k,"residual"};
            if (std::abs(s-b)<=xtol) return {s,fs,k,"step"};
            a=b; fa=fb; b=s; fb=fs;
        }
        return {b,fb,limit,"limit"};
    }
    if (a==b || !opposite(fa,fb))
        throw std::invalid_argument("bracketing methods require opposite endpoint signs");
    if (std::abs(fa)<std::abs(fb)) { std::swap(a,b); std::swap(fa,fb); }
    if (std::abs(b-a)<=xtol) return {b,fb,0,"bracket"};
    double c=a, fc=fa, d=a;
    bool bisected=true;
    int interpolations=0;
    double checkpoint=std::abs(b-a);
    for (int k=1; k<=limit; ++k) {
        if (std::abs(b-a)<=xtol) return {b,fb,k-1,"bracket"};
        double mid=a+(b-a)/2.0, s=mid;
        bool use_mid=method=="bisection";
        if (method=="regula") {
            double scale=std::max(std::abs(fa),std::abs(fb));
            double p=std::abs(fa)/scale, q=std::abs(fb)/scale;
            s=a+(b-a)*(p/(p+q));
        } else if (method=="brent") {
            if (fa!=fc && fb!=fc) {
                double scale=std::max({std::abs(fa),std::abs(fb),std::abs(fc)});
                double u=fa/scale, v=fb/scale, w=fc/scale;
                double den0=(u-v)*(u-w), den1=(v-u)*(v-w), den2=(w-u)*(w-v);
                s=std::numeric_limits<double>::quiet_NaN();
                if (den0!=0 && den1!=0 && den2!=0)
                    s=a*(v*w/den0)+b*(u*w/den1)+c*(u*v/den2);
            } else {
                double scale=std::max(std::abs(fa),std::abs(fb));
                double u=fa/scale, v=fb/scale;
                s=b-(v/(v-u))*(b-a);
            }
            double gate=a+0.25*(b-a);
            bool inside=std::min(gate,b)<s && s<std::max(gate,b);
            double old=bisected ? std::abs(b-c) : std::abs(c-d);
            use_mid=!inside || std::abs(s-b)>=0.5*old || old<=xtol || interpolations>=2;
        }
        if (use_mid || !std::isfinite(s) || !(std::min(a,b)<s && s<std::max(a,b))) {
            s=mid; use_mid=true;
        }
        if (s==a || s==b) return {b,fb,k-1,"stalled"};
        double fs=value(f,s);
        if (std::abs(fs)<=ftol) return {s,fs,k,"residual"};
        d=c; c=b; fc=fb;
        if (opposite(fa,fs)) { b=s; fb=fs; }
        else { a=s; fa=fs; }
        if (std::abs(fa)<std::abs(fb)) { std::swap(a,b); std::swap(fa,fb); }
        bisected=use_mid;
        if (use_mid || std::abs(b-a)<=0.5*checkpoint) {
            interpolations=0; checkpoint=std::abs(b-a);
        } else ++interpolations;
        if (std::abs(b-a)<=xtol) return {b,fb,k,"bracket"};
    }
    return {b,fb,limit,"limit"};
}

int main() {
    Fn f=[](double x){ return x*x-2.0; };
    Fn df=[](double x){ return 2.0*x; };
    // Rows are [left/start0, right/start1]; no matrix pivoting is needed.
    std::vector<std::vector<double>> starts(5,std::vector<double>{1.0,2.0});
    std::vector<std::string> methods={"bisection","newton","secant","regula","brent"};
    bool ok=true;
    std::cout<<std::fixed<<std::setprecision(6);
    for (int i=0; i<5; ++i) {
        Result r=solve(methods[i],f,df,starts[i][0],starts[i][1]);
        bool good=std::abs(r.fx)<=1e-10 && std::abs(r.x-std::sqrt(2.0))<=1e-10;
        ok=ok && good && (r.status=="residual" || r.status=="bracket");
        std::cout<<methods[i]<<": root = "<<r.x<<", check = "<<(good?"PASS":"FAIL")<<'\n';
    }
    std::cout<<"all methods = "<<(ok?"PASS":"FAIL")<<'\n';
    return 0;
}

import java.util.Locale;
import java.util.function.DoubleUnaryOperator;

public class Main {
    static class Result {
        double x, fx;
        int iterations;
        String status;
        Result(double x, double fx, int iterations, String status) {
            this.x=x; this.fx=fx; this.iterations=iterations; this.status=status;
        }
    }
    static double value(DoubleUnaryOperator f, double x) {
        double y=f.applyAsDouble(x);
        if (!Double.isFinite(y)) throw new IllegalArgumentException("nonfinite function or derivative value");
        return y;
    }
    static boolean opposite(double u, double v) { return (u<0.0)!=(v<0.0); }

    static Result solve(String method, DoubleUnaryOperator f, DoubleUnaryOperator df,
                        double a, double b, double xtol, double ftol, int limit) {
        boolean known=method.equals("bisection") || method.equals("newton") || method.equals("secant")
                      || method.equals("regula") || method.equals("brent");
        if (!known || !Double.isFinite(a) || !Double.isFinite(b) || !Double.isFinite(xtol)
            || !Double.isFinite(ftol) || xtol<=0 || ftol<0 || limit<0 || !Double.isFinite(b-a))
            throw new IllegalArgumentException("invalid method, endpoints, tolerances, or budget");
        // Newton uses b as its sole initial guess; a is unused.
        if (method.equals("newton")) {
            double x=b, fx=value(f,b);
            if (Math.abs(fx)<=ftol) return new Result(x,fx,0,"residual");
            for (int k=1; k<=limit; ++k) {
                double slope=value(df,x);
                if (slope==0) return new Result(x,fx,k-1,"stalled");
                double s=x-fx/slope;
                if (!Double.isFinite(s) || s==x) return new Result(x,fx,k-1,"stalled");
                double fs=value(f,s);
                if (Math.abs(fs)<=ftol) return new Result(s,fs,k,"residual");
                if (Math.abs(s-x)<=xtol) return new Result(s,fs,k,"step");
                x=s; fx=fs;
            }
            return new Result(x,fx,limit,"limit");
        }
        double fa=value(f,a), fb=value(f,b);
        if (Math.abs(fa)<=ftol) return new Result(a,fa,0,"residual");
        if (Math.abs(fb)<=ftol) return new Result(b,fb,0,"residual");
        if (method.equals("secant")) {
            for (int k=1; k<=limit; ++k) {
                double den=fb-fa;
                if (den==0 || !Double.isFinite(den)) return new Result(b,fb,k-1,"stalled");
                double s=b-fb*((b-a)/den);
                if (!Double.isFinite(s) || s==b) return new Result(b,fb,k-1,"stalled");
                double fs=value(f,s);
                if (Math.abs(fs)<=ftol) return new Result(s,fs,k,"residual");
                if (Math.abs(s-b)<=xtol) return new Result(s,fs,k,"step");
                a=b; fa=fb; b=s; fb=fs;
            }
            return new Result(b,fb,limit,"limit");
        }
        if (a==b || !opposite(fa,fb))
            throw new IllegalArgumentException("bracketing methods require opposite endpoint signs");
        if (Math.abs(fa)<Math.abs(fb)) {
            double t=a; a=b; b=t; t=fa; fa=fb; fb=t;
        }
        if (Math.abs(b-a)<=xtol) return new Result(b,fb,0,"bracket");
        double c=a, fc=fa, d=a;
        boolean bisected=true;
        int interpolations=0;
        double checkpoint=Math.abs(b-a);
        for (int k=1; k<=limit; ++k) {
            if (Math.abs(b-a)<=xtol) return new Result(b,fb,k-1,"bracket");
            double mid=a+(b-a)/2.0, s=mid;
            boolean useMid=method.equals("bisection");
            if (method.equals("regula")) {
                double scale=Math.max(Math.abs(fa),Math.abs(fb));
                double p=Math.abs(fa)/scale, q=Math.abs(fb)/scale;
                s=a+(b-a)*(p/(p+q));
            } else if (method.equals("brent")) {
                if (fa!=fc && fb!=fc) {
                    double scale=Math.max(Math.max(Math.abs(fa),Math.abs(fb)),Math.abs(fc));
                    double u=fa/scale, v=fb/scale, w=fc/scale;
                    double den0=(u-v)*(u-w), den1=(v-u)*(v-w), den2=(w-u)*(w-v);
                    s=Double.NaN;
                    if (den0!=0 && den1!=0 && den2!=0)
                        s=a*(v*w/den0)+b*(u*w/den1)+c*(u*v/den2);
                } else {
                    double scale=Math.max(Math.abs(fa),Math.abs(fb));
                    double u=fa/scale, v=fb/scale;
                    s=b-(v/(v-u))*(b-a);
                }
                double gate=a+0.25*(b-a);
                boolean inside=Math.min(gate,b)<s && s<Math.max(gate,b);
                double old=bisected ? Math.abs(b-c) : Math.abs(c-d);
                useMid=!inside || Math.abs(s-b)>=0.5*old || old<=xtol || interpolations>=2;
            }
            if (useMid || !Double.isFinite(s) || !(Math.min(a,b)<s && s<Math.max(a,b))) {
                s=mid; useMid=true;
            }
            if (s==a || s==b) return new Result(b,fb,k-1,"stalled");
            double fs=value(f,s);
            if (Math.abs(fs)<=ftol) return new Result(s,fs,k,"residual");
            d=c; c=b; fc=fb;
            if (opposite(fa,fs)) { b=s; fb=fs; }
            else { a=s; fa=fs; }
            if (Math.abs(fa)<Math.abs(fb)) {
                double t=a; a=b; b=t; t=fa; fa=fb; fb=t;
            }
            bisected=useMid;
            if (useMid || Math.abs(b-a)<=0.5*checkpoint) {
                interpolations=0; checkpoint=Math.abs(b-a);
            } else ++interpolations;
            if (Math.abs(b-a)<=xtol) return new Result(b,fb,k,"bracket");
        }
        return new Result(b,fb,limit,"limit");
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);
        DoubleUnaryOperator f=x -> x*x-2.0, df=x -> 2.0*x;
        // Rows are [left/start0, right/start1]; no matrix pivoting is needed.
        double[][] starts={ {1.0,2.0}, {1.0,2.0}, {1.0,2.0}, {1.0,2.0}, {1.0,2.0} };
        String[] methods={"bisection","newton","secant","regula","brent"};
        boolean ok=true;
        for (int i=0; i<5; ++i) {
            Result r=solve(methods[i],f,df,starts[i][0],starts[i][1],1e-12,1e-12,200);
            boolean good=Math.abs(r.fx)<=1e-10 && Math.abs(r.x-Math.sqrt(2.0))<=1e-10;
            ok=ok && good && (r.status.equals("residual") || r.status.equals("bracket"));
            System.out.printf("%s: root = %.6f, check = %s%n",methods[i],r.x,good?"PASS":"FAIL");
        }
        System.out.println("all methods = "+(ok?"PASS":"FAIL"));
    }
}

The identical printed output is:

• bisection: root = 1.414214, check = PASS
• newton: root = 1.414214, check = PASS
• secant: root = 1.414214, check = PASS
• regula: root = 1.414214, check = PASS
• brent: root = 1.414214, check = PASS
• all methods = PASS

With this configuration, the trial counts are $38$, $5$, $6$, $16$, and $5$, respectively; all five stop on the residual criterion. Six printed decimal places are presentation rounding, not the internal accuracy target. The driver also requires a residual or bracket termination for the aggregate PASS, rather than accepting a small step as sufficient evidence.

Complexity, Correctness, and Invariants

Let $C_f$ be the cost of evaluating $f$, $C_d$ the cost of evaluating $f'$, $W$ the initial bracket width, and $\varepsilon=\tau_x$. Let $N=\max(0,\lceil\log_2(W/\varepsilon)\rceil)$ for $W\gt0$. Function values are cached, so an ordinary update uses one new function evaluation. Initial sampling costs one evaluation for Newton and up to two for the others. All implementations use $O(1)$ auxiliary state beyond the function's own workspace and the fixed sample driver.

There is no distribution-free “average iteration count”: it depends on the function, its scaling, and the starts. Under local quadratic convergence, Newton needs $O(\log\log(1/\varepsilon))$ further updates to reach a small position error; the secant method has the same asymptotic iteration form with a different convergence order. These are local statements, not global guarantees. A multiple root, inaccurate derivative, or a distant initial guess changes the conclusion. The table's Brent bound belongs to the implemented watchdog variant, rather than every algorithm called Brent.

Bracket proof. Suppose $f(a)$ and $f(b)$ have opposite signs and $f$ is continuous. If $f(s)$ is nonzero, it matches exactly one endpoint sign. Replacing that endpoint preserves opposite signs and retains at least one root. Bisection halves the width. For the Brent variant, either interpolation halves the checkpoint interval within two updates or the following update bisects an interval no wider than that checkpoint. Inductively the checkpoint width halves within every three updates. The best-residual swap only renames endpoints and does not change the interval.

Open-method convergence. At a sufficiently smooth simple root $r$, Taylor expansion gives Newton's local error relation $e_{k+1}\approx[f''(r)/(2f'(r))]e_k^2$, where $e_k=x_k-r$. Secant uses a changing slope approximation and typically satisfies $e_{k+1}\approx C e_k e_{k-1}$, leading to the golden-ratio order. Neither relation supplies a global interval invariant, so a step criterion is reported separately from a residual criterion.

These certificates assume reliable function signs. Ordinary floating-point evaluation can suffer cancellation, underflow, or model error; the code's finite-value checks do not create a rigorous interval-arithmetic proof. If adjacent floating-point numbers prevent a new interior trial, the solver reports stalled instead of claiming to have reached an impossible tolerance.

Variants, Pitfalls, and Optimizations

Choose a method to match the available information. Use bisection for a transparent width guarantee when function evaluations are cheap. Use a bracketed hybrid such as Brent when evaluations are more expensive and interpolation is useful. Newton is attractive when a reliable derivative and a good starting point are available; secant avoids derivatives but gives up the bracket. For an equation known to be monotone, establish a bracket by bounded expansion within the model's domain before solving; terminate the expansion if no valid bracket is found.

A concrete Newton cycle. For $f(x)=x^3-2x+2$, $f'(x)=3x^2-2$. Starting at $0$ gives $1$; starting from that $1$ gives $0$ again. Every arithmetic operation is finite, yet Newton cycles forever without an iteration limit. The same function is continuous and has opposite signs at $-2$ and $-1$, so bisection or Brent can use that bracket. The two approaches can behave differently on the same equation.

Endpoint locking in false position. For $f(x)=x^{10}-1$ on $[0,2]$, the first false-position point is $2/1024=1/512$, very close to the lower endpoint. The high positive endpoint can remain fixed through many updates while the lower endpoint moves slowly. Illinois and Pegasus variants reweight the retained endpoint's interpolation ordinate; when implementing them, keep actual residuals separately for sign checks and stopping. Those variants are not implemented in the three programs.

Repeated roots. A bracket cannot detect the even-multiplicity root of $(x-1)^2$ on $[0,2]$ from endpoint signs, because both values are positive. Ordinary Newton typically converges only linearly to a multiple root. If the multiplicity $m$ is known, the modified update $x_{k+1}=x_k-m f(x_k)/f'(x_k)$ can restore quadratic local convergence under appropriate smoothness. Root isolation, derivative information, or polynomial-specific tools may be needed to find all roots; a single scalar run returns only one candidate.

Discontinuities and scaling. The function $1/x$ changes sign across $[-1,1]$ and has no zero; continuity fails at $0$. Bisection's first evaluation there is undefined. A discontinuous step can even shrink a sign interval while every residual stays far from zero. Conversely, $10^{-12}(x-100)$ at $x=0$ has residual $10^{-10}$ despite position error $100$. Choose residual tolerances in the model's output units, use derivative information to interpret them, and add relative position tolerances when the unknown spans many orders of magnitude. The supplied code uses absolute tolerances for clarity.

Numerical implementation details. Cache endpoint values, compare signs without multiplication, normalize interpolation ordinates, and reject invalid proposals. A computed zero derivative is tested exactly in these programs; there is no universal “small derivative” threshold because multiplying the entire equation by a constant changes derivative magnitude without changing its roots. Model-specific derivative checks or a safeguarded Newton variant can be useful. Fused multiply-add, extended precision, or interval arithmetic can improve difficult evaluations, but their use may change cross-language iteration paths.

Practice Problems and Real-World Use

Begin with equations whose roots are known independently, then use cases that break a method's assumptions. For each experiment, record function evaluations, derivative evaluations, termination reason, residual, and an error bound when available. Comparing only iteration counts hides the cost of derivatives and expensive function models.

For further reading, SciPy's Brent documentation describes the bracketed combination of bisection and inverse interpolation, and its Newton documentation documents derivative-based and secant interfaces. Brent's original research record points to the work underlying the hybrid method. Library stopping rules and safeguards vary; compare their documented contract with the explicit one used here.