prospicio_math/
roots.rs

1//! Root finding in one variable.
2//!
3//! Every solver here takes a bracket and returns a point inside it; none
4//! extrapolates. Callers own the bracket, because only they know the
5//! domain (an alpha in `[1e-6, 1e6]`, a probability in `(0, 1)`).
6
7/// The boundary in `(lo, hi)` of a predicate that holds below a point and
8/// fails above it, by bisection to full precision: the loop stops when
9/// the midpoint no longer lies strictly between the ends.
10///
11/// For `f(x) = target` with `f` increasing, pass `|x| f(x) < target`.
12/// If `below` is not monotone the result is some point where it changes.
13///
14/// ```
15/// use prospicio_math::roots::bisect;
16///
17/// let r = bisect(0.0, 2.0, |x| x * x < 2.0);
18/// assert!((r - 2f64.sqrt()).abs() < 1e-15);
19/// ```
20pub fn bisect(lo: f64, hi: f64, mut below: impl FnMut(f64) -> bool) -> f64 {
21    let (mut lo, mut hi) = (lo, hi);
22    // 2100 halvings span the whole f64 range; the exact stop comes first.
23    for _ in 0..2100 {
24        let mid = 0.5 * (lo + hi);
25        if mid <= lo || mid >= hi {
26            break;
27        }
28        if below(mid) {
29            lo = mid;
30        } else {
31            hi = mid;
32        }
33    }
34    0.5 * (lo + hi)
35}
36
37/// [`bisect`] on a log scale, for a bracket `0 < lo < hi` that spans orders
38/// of magnitude: each step splits at the geometric mean, so relative
39/// precision improves at the same rate everywhere in the bracket.
40///
41/// ```
42/// use prospicio_math::roots::bisect_log;
43///
44/// let r = bisect_log(1e-6, 1e6, |x| x.ln() < 3.0);
45/// assert!((r / 3f64.exp() - 1.0).abs() < 1e-14);
46/// ```
47pub fn bisect_log(lo: f64, hi: f64, mut below: impl FnMut(f64) -> bool) -> f64 {
48    let (mut lo, mut hi) = (lo, hi);
49    for _ in 0..2100 {
50        let mid = lo.sqrt() * hi.sqrt();
51        if mid <= lo || mid >= hi {
52            break;
53        }
54        if below(mid) {
55            lo = mid;
56        } else {
57            hi = mid;
58        }
59    }
60    lo.sqrt() * hi.sqrt()
61}
62
63/// A root of `f` between `x0` and `x1`, where `f0 = f(x0)` and `f1 = f(x1)`
64/// have opposite signs, by the Illinois variant of regula falsi.
65///
66/// Regula falsi keeps a bracket like bisection but steps to where the
67/// chord crosses zero; Illinois halves the stale end's value whenever the
68/// same end is kept twice, which restores superlinear convergence. On a
69/// smooth score it needs a handful of evaluations where bisection needs
70/// fifty. Returns an exact zero as soon as one is evaluated.
71///
72/// ```
73/// use prospicio_math::roots::illinois;
74///
75/// let f = |x: f64| x.exp() - 2.0;
76/// let r = illinois(0.0, 1.0, f(0.0), f(1.0), f);
77/// assert!((r - 2f64.ln()).abs() < 1e-15);
78/// ```
79pub fn illinois(x0: f64, x1: f64, f0: f64, f1: f64, mut f: impl FnMut(f64) -> f64) -> f64 {
80    let (mut x0, mut x1, mut f0, mut f1) = (x0, x1, f0, f1);
81    let positive_left = f0 > 0.0;
82    let mut side = 0;
83    for _ in 0..200 {
84        let x = (x0 * f1 - x1 * f0) / (f1 - f0);
85        if !(x > x0.min(x1) && x < x0.max(x1)) || (x1 - x0).abs() <= 1e-15 * x1.abs().max(1.0) {
86            break;
87        }
88        let fx = f(x);
89        if fx == 0.0 {
90            return x;
91        }
92        if (fx > 0.0) == positive_left {
93            x0 = x;
94            f0 = fx;
95            if side == -1 {
96                f1 *= 0.5;
97            }
98            side = -1;
99        } else {
100            x1 = x;
101            f1 = fx;
102            if side == 1 {
103                f0 *= 0.5;
104            }
105            side = 1;
106        }
107    }
108    0.5 * (x0 + x1)
109}
110
111#[cfg(test)]
112mod tests {
113    use super::*;
114
115    #[test]
116    fn bisect_reaches_full_precision() {
117        let r = bisect(1.0, 2.0, |x| x * x * x < 3.0);
118        let next = f64::from_bits(r.to_bits() + 1);
119        let prev = f64::from_bits(r.to_bits() - 1);
120        assert!(prev * prev * prev < 3.0 || next * next * next >= 3.0);
121        assert!((r - 3f64.cbrt()).abs() <= 2.0 * f64::EPSILON);
122    }
123
124    #[test]
125    fn bisect_handles_tiny_roots_in_wide_brackets() {
126        let r = bisect(1e-12, 1e3, |x| x < 3e-10);
127        assert!((r / 3e-10 - 1.0).abs() < 1e-14);
128    }
129
130    #[test]
131    fn bisect_log_handles_huge_brackets() {
132        let r = bisect_log(1e-300, 1e300, |x| x < 7e200);
133        assert!((r / 7e200 - 1.0).abs() < 1e-14);
134    }
135
136    #[test]
137    fn illinois_accepts_either_sign_order() {
138        let f = |x: f64| 2.0 - x * x;
139        let r = illinois(0.0, 3.0, f(0.0), f(3.0), f);
140        assert!((r - 2f64.sqrt()).abs() < 1e-14);
141        let g = |x: f64| x * x - 2.0;
142        let r = illinois(0.0, 3.0, g(0.0), g(3.0), g);
143        assert!((r - 2f64.sqrt()).abs() < 1e-14);
144    }
145
146    #[test]
147    fn illinois_is_fast_on_smooth_functions() {
148        let mut calls = 0;
149        let f = |x: f64| x.powi(5) - 7.0;
150        let r = illinois(0.0, 3.0, f(0.0), f(3.0), |x| {
151            calls += 1;
152            f(x)
153        });
154        assert!((r - 7f64.powf(0.2)).abs() < 1e-14);
155        assert!(calls < 30, "{calls} calls");
156    }
157}
prospicio_math/roots.rs

prospicio_math/
roots.rs
