regit-svi 2.0.0

// Copyright 2026 Regit.io — Nicolas Koenig
// SPDX-License-Identifier: Apache-2.0

//! Bracket-preserving Brent root finder.

/// Successful bracket-preserving root solve with numerical evidence.
#[derive(Debug, Clone, Copy, PartialEq)]
pub(crate) struct BrentRoot {
    root: f64,
    lower: f64,
    upper: f64,
    residual: f64,
    evaluations: usize,
    exact: bool,
}

impl BrentRoot {
    pub(crate) const fn root(self) -> f64 {
        self.root
    }

    pub(crate) const fn lower(self) -> f64 {
        self.lower
    }

    pub(crate) const fn upper(self) -> f64 {
        self.upper
    }

    pub(crate) const fn residual(self) -> f64 {
        self.residual
    }

    pub(crate) const fn evaluations(self) -> usize {
        self.evaluations
    }

    pub(crate) const fn exact(self) -> bool {
        self.exact
    }
}

/// Finds an approximate root inside a sign-changing bracket by Brent's method.
#[must_use]
#[allow(clippy::many_single_char_names)] // Names follow Brent's published state variables.
#[cfg(test)]
pub(crate) fn brent_root<F>(
    f: F,
    a: f64,
    b: f64,
    tolerance: f64,
    max_iterations: usize,
) -> Option<f64>
where
    F: Fn(f64) -> f64,
{
    brent_root_with_evidence(f, a, b, tolerance, max_iterations).map(BrentRoot::root)
}

/// Finds a root and retains the final bracket, residual, and evaluation count.
#[must_use]
#[allow(clippy::many_single_char_names)] // Names follow Brent's published state variables.
pub(crate) fn brent_root_with_evidence<F>(
    f: F,
    a: f64,
    b: f64,
    tolerance: f64,
    max_iterations: usize,
) -> Option<BrentRoot>
where
    F: Fn(f64) -> f64,
{
    let mut a = a;
    let mut b = b;
    let mut fa = f(a);
    let mut fb = f(b);
    let mut evaluations = 2_usize;
    if !fa.is_finite() || !fb.is_finite() {
        return None;
    }
    if fa == 0.0 {
        return Some(BrentRoot {
            root: a,
            lower: a,
            upper: a,
            residual: 0.0,
            evaluations,
            exact: true,
        });
    }
    if fb == 0.0 {
        return Some(BrentRoot {
            root: b,
            lower: b,
            upper: b,
            residual: 0.0,
            evaluations,
            exact: true,
        });
    }
    if fa.is_sign_positive() == fb.is_sign_positive() {
        return None;
    }
    if fa.abs() < fb.abs() {
        core::mem::swap(&mut a, &mut b);
        core::mem::swap(&mut fa, &mut fb);
    }

    let mut c = a;
    let mut fc = fa;
    let mut d = a;
    let mut used_bisection = true;
    for _ in 0..max_iterations {
        if (b - a).abs() <= tolerance || fb == 0.0 {
            return Some(BrentRoot {
                root: b,
                lower: a.min(b),
                upper: a.max(b),
                residual: fb.abs(),
                evaluations,
                exact: fb == 0.0,
            });
        }
        let mut candidate = if (fa - fc).abs() > f64::EPSILON && (fb - fc).abs() > f64::EPSILON {
            a * fb * fc / ((fa - fb) * (fa - fc))
                + b * fa * fc / ((fb - fa) * (fb - fc))
                + c * fa * fb / ((fc - fa) * (fc - fb))
        } else {
            b - fb * (b - a) / (fb - fa)
        };
        let interpolation_bound = (3.0 * a + b) / 4.0;
        let lower = interpolation_bound.min(b);
        let upper = interpolation_bound.max(b);
        let reject_interpolation = !(lower..=upper).contains(&candidate)
            || (used_bisection && (candidate - b).abs() >= (b - c).abs() / 2.0)
            || (!used_bisection && (candidate - b).abs() >= (c - d).abs() / 2.0)
            || (used_bisection && (b - c).abs() < tolerance)
            || (!used_bisection && (c - d).abs() < tolerance);
        if reject_interpolation {
            candidate = f64::midpoint(a, b);
            used_bisection = true;
        } else {
            used_bisection = false;
        }

        let f_candidate = f(candidate);
        evaluations = evaluations.saturating_add(1);
        if !f_candidate.is_finite() {
            return None;
        }
        d = c;
        c = b;
        fc = fb;
        if fa.is_sign_positive() == f_candidate.is_sign_positive() {
            a = candidate;
            fa = f_candidate;
        } else {
            b = candidate;
            fb = f_candidate;
        }
        if fa.abs() < fb.abs() {
            core::mem::swap(&mut a, &mut b);
            core::mem::swap(&mut fa, &mut fb);
        }
    }
    // Exhausting the budget is not convergence evidence. Callers may retain
    // their original sign-changing bracket, but must not report an endpoint
    // as a refined root.
    None
}

#[cfg(test)]
#[allow(clippy::expect_used)]
mod tests {
    use super::*;

    #[test]
    fn finds_sqrt_two() {
        let root = brent_root(|x| x * x - 2.0, 0.0, 2.0, 1e-14, 200)
            .expect("sqrt(2) fixture has a finite sign-changing bracket");
        assert!((root - 2.0_f64.sqrt()).abs() < 1e-12);
    }

    #[test]
    fn finds_cubic_root() {
        let root = brent_root(|x| x * x * x - x - 2.0, 1.0, 2.0, 1e-14, 200)
            .expect("cubic fixture has a finite sign-changing bracket");
        assert!((root - 1.521_379_706_804_567_6).abs() < 1e-10);
    }

    #[test]
    fn accepts_endpoint_root() {
        let root =
            brent_root(|x| x - 3.0, 3.0, 5.0, 1e-12, 50).expect("endpoint is an exact finite root");
        assert!((root - 3.0).abs() < 1e-15);
    }

    #[test]
    fn rejects_invalid_brackets_and_values() {
        assert!(brent_root(|x| x * x + 1.0, -1.0, 1.0, 1e-12, 50).is_none());
        assert!(brent_root(|_| f64::NAN, -1.0, 1.0, 1e-12, 50).is_none());
        assert!(brent_root(|x| x - 0.25, 0.0, 1.0, 1e-12, 0).is_none());
        assert!(brent_root_with_evidence(|x| x - 0.25, 0.0, 1.0, 1e-12, 0).is_none());
    }

    #[test]
    fn finds_transcendental_root() {
        let root = brent_root(|x| x.cos() - x, 0.0, 1.0, 1e-14, 200)
            .expect("fixed-point fixture has a finite sign-changing bracket");
        assert!((root - 0.739_085_133_215_160_6).abs() < 1e-10);
    }
}