axiolid_numeric/
root.rs

1//! Bracketed scalar root finding with guaranteed convergence.
2//!
3//! [`find_root`] is Brent's method (inverse quadratic interpolation, the
4//! secant step, and bisection), as in Brent, *Algorithms for Minimization
5//! without Derivatives* (1973), chapter 4, and OCCT `math_BracketedRoot`.
6//! An interpolation step is accepted only if it is less than half the step
7//! before last and stays inside the bracket; otherwise the step bisects.
8//! Convergence is therefore guaranteed for any continuous function with a
9//! sign change, in at most about the square of the bisection count and in
10//! practice far fewer.
11//!
12//! The answer is a bracket, not a point: the function changes sign across
13//! `[lower, upper]` (or vanishes exactly at `x`), so for a continuous
14//! function a root lies inside and `|x - root| <= upper - lower`.
15
16use crate::error::{finite, tolerance, NumericError, NumericResult, Status};
17
18/// A root located inside a sign-change bracket.
19#[derive(Debug, Clone, Copy, PartialEq)]
20#[must_use]
21pub struct Root {
22    /// Best estimate: the bracket end with the smaller `|f|`.
23    pub x: f64,
24    /// `f(x)`.
25    pub value: f64,
26    /// Lower end of the final bracket.
27    pub lower: f64,
28    /// Upper end of the final bracket. The function changes sign across
29    /// `[lower, upper]`, or `lower == upper == x` with `value == 0`.
30    pub upper: f64,
31    /// Iterations performed.
32    pub iterations: usize,
33    /// Function evaluations performed.
34    pub evaluations: usize,
35    /// [`Status::Converged`], or [`Status::BudgetExhausted`] if the
36    /// iteration cap was hit first; the bracket is valid either way.
37    pub status: Status,
38}
39
40impl Root {
41    /// A bound on `|x - root|`: the final bracket width.
42    pub fn error_bound(&self) -> f64 {
43        self.upper - self.lower
44    }
45}
46
47/// Iteration cap. Bisection from any finite bracket reaches the stopping
48/// width in under 2200 halvings, and Brent's method rarely needs more steps
49/// than bisection.
50const MAX_ITERATIONS: usize = 10_000;
51
52/// Find a root of `f` in `[lower, upper]` by Brent's method.
53///
54/// `f(lower)` and `f(upper)` must differ in sign (or one of them be zero).
55/// Iteration stops when the bracket is no wider than
56/// `4 * eps * |x| + x_tolerance`, or when its ends are adjacent floats;
57/// `x_tolerance = 0` therefore asks for full precision.
58///
59/// # Errors
60///
61/// Refuses non-finite bounds or tolerance, `lower > upper`, a negative
62/// tolerance, bracket ends of the same sign
63/// ([`NumericError::NoSignChange`]), and any non-finite value of `f`.
64pub fn find_root<F>(mut f: F, lower: f64, upper: f64, x_tolerance: f64) -> NumericResult<Root>
65where
66    F: FnMut(f64) -> f64,
67{
68    finite(lower, "lower bound")?;
69    finite(upper, "upper bound")?;
70    tolerance(x_tolerance, "x tolerance")?;
71    if lower > upper {
72        return Err(NumericError::InvalidArgument {
73            name: "bracket",
74            reason: "must have lower <= upper",
75        });
76    }
77    let mut evaluations = 0usize;
78    let mut eval = |x: f64| -> NumericResult<f64> {
79        evaluations += 1;
80        let v = f(x);
81        if v.is_finite() {
82            Ok(v)
83        } else {
84            Err(NumericError::NonFiniteEvaluation {
85                name: "function",
86                at: Some(x),
87            })
88        }
89    };
90
91    let (mut a, mut b) = (lower, upper);
92    let mut fa = eval(a)?;
93    let mut fb = eval(b)?;
94    if fa == 0.0 {
95        return Ok(exact(a, 0, evaluations));
96    }
97    if fb == 0.0 {
98        return Ok(exact(b, 0, evaluations));
99    }
100    if fa.signum() == fb.signum() {
101        return Err(NumericError::NoSignChange {
102            f_lower: fa,
103            f_upper: fb,
104        });
105    }
106
107    let (mut c, mut fc) = (a, fa);
108    let mut d = b - a;
109    let mut e = d;
110    let mut iterations = 0usize;
111    loop {
112        if fb.signum() == fc.signum() {
113            c = a;
114            fc = fa;
115            d = b - a;
116            e = d;
117        }
118        if fc.abs() < fb.abs() {
119            a = b;
120            b = c;
121            c = a;
122            fa = fb;
123            fb = fc;
124            fc = fa;
125        }
126        let tol1 = 2.0 * f64::EPSILON * b.abs() + 0.5 * x_tolerance;
127        let xm = 0.5 * (c - b);
128        let midpoint = b + xm;
129        if xm.abs() <= tol1 || fb == 0.0 || midpoint == b || midpoint == c {
130            break;
131        }
132        if iterations >= MAX_ITERATIONS {
133            return Ok(finish(
134                b,
135                fb,
136                c,
137                iterations,
138                evaluations,
139                Status::BudgetExhausted,
140            ));
141        }
142        iterations += 1;
143
144        if e.abs() >= tol1 && fa.abs() > fb.abs() {
145            let s = fb / fa;
146            let (mut p, mut q);
147            if a == c {
148                p = 2.0 * xm * s;
149                q = 1.0 - s;
150            } else {
151                let qq = fa / fc;
152                let r = fb / fc;
153                p = s * (2.0 * xm * qq * (qq - r) - (b - a) * (r - 1.0));
154                q = (qq - 1.0) * (r - 1.0) * (s - 1.0);
155            }
156            if p > 0.0 {
157                q = -q;
158            } else {
159                p = -p;
160            }
161            if 2.0 * p < (3.0 * xm * q - (tol1 * q).abs()).min((e * q).abs()) {
162                e = d;
163                d = p / q;
164            } else {
165                d = xm;
166                e = d;
167            }
168        } else {
169            d = xm;
170            e = d;
171        }
172        a = b;
173        fa = fb;
174        b += if d.abs() > tol1 { d } else { tol1.copysign(xm) };
175        fb = eval(b)?;
176    }
177    if fb == 0.0 {
178        return Ok(exact(b, iterations, evaluations));
179    }
180    Ok(finish(b, fb, c, iterations, evaluations, Status::Converged))
181}
182
183fn exact(x: f64, iterations: usize, evaluations: usize) -> Root {
184    Root {
185        x,
186        value: 0.0,
187        lower: x,
188        upper: x,
189        iterations,
190        evaluations,
191        status: Status::Converged,
192    }
193}
194
195fn finish(b: f64, fb: f64, c: f64, iterations: usize, evaluations: usize, status: Status) -> Root {
196    Root {
197        x: b,
198        value: fb,
199        lower: b.min(c),
200        upper: b.max(c),
201        iterations,
202        evaluations,
203        status,
204    }
205}
axiolid_numeric/root.rs

axiolid_numeric/
root.rs
