lox_core/math/

roots.rs

// SPDX-FileCopyrightText: 2024 Helge Eichhorn <git@helgeeichhorn.de>
//
// SPDX-License-Identifier: MPL-2.0

//! Root-finding algorithms: Steffensen, Newton, safeguarded Newton, and Brent
//! methods.

use lox_approx::approx_eq;
use thiserror::Error;

use crate::error::LoxError;
use crate::math::callback::{Callback, CallbackWithDerivative};
use crate::math::float::{abs, powi, sqrt};

/// Finds a root of `f` starting from an initial guess.
pub trait FindRoot {
    /// Finds a root of `f` starting from `initial_guess`.
    fn find(&self, f: impl Callback, initial_guess: f64) -> Result<f64, RootFinderError>;
}

/// Finds a root of `f` using both the function and its derivative.
pub trait FindRootWithDerivative {
    /// Finds a root of `f` using `derivative`, starting from `initial_guess`.
    fn find_with_derivative(
        &self,
        f: impl Callback,
        derivative: impl Callback,
        initial_guess: f64,
    ) -> Result<f64, RootFinderError>;
}

/// Finds a root of `f` within a bracket `(a, b)`.
pub trait FindBracketedRoot {
    /// Finds a root of `f` within `bracket`, reusing the function values at the
    /// bracket endpoints instead of evaluating them again.
    ///
    /// `values` must equal `(f(bracket.0), f(bracket.1))`.
    fn find_in_bracket_with_values(
        &self,
        f: impl Callback,
        bracket: (f64, f64),
        values: (f64, f64),
    ) -> Result<f64, RootFinderError>;

    /// Finds a root of `f` within the given `bracket`.
    fn find_in_bracket(
        &self,
        f: impl Callback,
        bracket: (f64, f64),
    ) -> Result<f64, RootFinderError> {
        let fa = f.call(bracket.0)?;
        let fb = f.call(bracket.1)?;
        self.find_in_bracket_with_values(f, bracket, (fa, fb))
    }
}

/// Finds a root of `f` within a bracket `(a, b)`, using its derivative.
pub trait FindBracketedRootWithDerivative {
    /// Finds a root of `f` within `bracket`, reusing the function values at the
    /// bracket endpoints instead of evaluating them again.
    ///
    /// `values` must equal `(f(bracket.0), f(bracket.1))`.
    fn find_in_bracket_with_derivative_values(
        &self,
        f: impl CallbackWithDerivative,
        bracket: (f64, f64),
        values: (f64, f64),
    ) -> Result<f64, RootFinderError>;

    /// Finds a root of `f` within the given `bracket`.
    fn find_in_bracket_with_derivative(
        &self,
        f: impl CallbackWithDerivative,
        bracket: (f64, f64),
    ) -> Result<f64, RootFinderError> {
        let (fa, _) = f.call(bracket.0)?;
        let (fb, _) = f.call(bracket.1)?;
        self.find_in_bracket_with_derivative_values(f, bracket, (fa, fb))
    }
}

/// Error returned by root-finding algorithms.
#[derive(Debug, Error)]
pub enum RootFinderError {
    /// The algorithm did not converge within the maximum number of iterations.
    #[error("not converged after {iterations} iterations at x = {x}, residual {residual}")]
    NotConverged {
        /// Number of iterations performed before giving up.
        iterations: u32,
        /// The best root estimate reached.
        x: f64,
        /// The residual `f(x)` at the best estimate.
        residual: f64,
    },
    /// The root is not within the given bracket.
    #[error("root not in bracket")]
    NotInBracket,
    /// The objective function returned a non-finite value.
    #[error("function returned a non-finite value ({value}) at x = {x}")]
    NonFinite {
        /// The point at which the function was evaluated.
        x: f64,
        /// The non-finite value returned.
        value: f64,
    },
    /// The derivative returned a non-finite value.
    #[error("derivative returned a non-finite value ({value}) at x = {x}")]
    NonFiniteDerivative {
        /// The point at which the derivative was evaluated.
        x: f64,
        /// The non-finite value returned.
        value: f64,
    },
    /// The iteration update diverged, e.g. due to a zero derivative at a
    /// stationary point or a vanishing update denominator.
    #[error("iteration step diverged at x = {x}")]
    DivergedStep {
        /// The iterate at which the update diverged.
        x: f64,
    },
    /// The objective function returned an error.
    #[error(transparent)]
    Callback(#[from] LoxError),
}

/// Evaluates `f` at `x`, mapping a callback failure or non-finite result to the
/// corresponding [`RootFinderError`].
fn eval_finite<F: Callback>(f: &F, x: f64) -> Result<f64, RootFinderError> {
    let value = f.call(x)?;
    if !value.is_finite() {
        return Err(RootFinderError::NonFinite { x, value });
    }
    Ok(value)
}

/// Evaluates `f` and its derivative at `x`, mapping a callback failure or a
/// non-finite function value to the corresponding [`RootFinderError`].
///
/// A non-finite *derivative* is not an error here: safeguarded methods fall
/// back to bisection when the derivative is unusable.
fn eval_finite_with_derivative<F: CallbackWithDerivative>(
    f: &F,
    x: f64,
) -> Result<(f64, f64), RootFinderError> {
    let (value, derivative) = f.call(x)?;
    if !value.is_finite() {
        return Err(RootFinderError::NonFinite { x, value });
    }
    Ok((value, derivative))
}

/// Steffensen's method for root-finding (derivative-free).
///
/// The tolerances bound the iteration step, not the residual: near steep
/// features the iteration can stall within tolerance where `f(x)` is not
/// small and report a false root. Prefer [`Brent`] when a bracket is
/// available.
#[derive(Debug, Copy, Clone, PartialEq)]
pub struct Steffensen {
    max_iter: u32,
    /// Absolute tolerance on the iteration step.
    abs_tol: f64,
    /// Relative tolerance on the iteration step.
    rel_tol: f64,
}

impl Default for Steffensen {
    fn default() -> Self {
        Self {
            max_iter: 1000,
            abs_tol: sqrt(f64::EPSILON),
            rel_tol: sqrt(f64::EPSILON),
        }
    }
}

impl Steffensen {
    /// Sets the maximum number of iterations.
    pub fn with_max_iter(mut self, max_iter: u32) -> Self {
        self.max_iter = max_iter;
        self
    }

    /// Sets the absolute tolerance on the iteration step.
    pub fn with_abs_tol(mut self, abs_tol: f64) -> Self {
        self.abs_tol = abs_tol;
        self
    }

    /// Sets the relative tolerance on the iteration step.
    pub fn with_rel_tol(mut self, rel_tol: f64) -> Self {
        self.rel_tol = rel_tol;
        self
    }
}

impl FindRoot for Steffensen {
    fn find(&self, f: impl Callback, initial_guess: f64) -> Result<f64, RootFinderError> {
        let mut p0 = initial_guess;
        let mut last: Option<(f64, f64)> = None;
        for _ in 0..self.max_iter {
            let fp0 = eval_finite(&f, p0)?;
            // An initial guess that is already a root is returned directly,
            // avoiding a 0/0 update.
            if fp0 == 0.0 {
                return Ok(p0);
            }
            last = Some((p0, fp0));
            let f1 = p0 + fp0;
            let ff1 = eval_finite(&f, f1)?;
            let f2 = f1 + ff1;
            let p = p0 - powi(f1 - p0, 2) / (f2 - 2.0 * f1 + p0);
            if !p.is_finite() {
                return Err(RootFinderError::DivergedStep { x: p0 });
            }
            if approx_eq!(p, p0, rtol <= self.rel_tol, atol <= self.abs_tol) {
                return Ok(p);
            }
            p0 = p;
        }
        // Report the last point where `f` was actually evaluated rather than
        // re-evaluating a stepped final iterate, which may lie outside the
        // callback's valid domain and turn non-convergence into a callback error
        // or panic. When `max_iter == 0` the loop never runs, so fall back to the
        // still-in-domain initial guess.
        let (x, residual) = match last {
            Some(pair) => pair,
            None => (p0, eval_finite(&f, p0)?),
        };
        Err(RootFinderError::NotConverged {
            iterations: self.max_iter,
            x,
            residual,
        })
    }
}

/// Newton-Raphson method for root-finding (requires derivative).
///
/// The tolerances bound the iteration step, not the residual: near steep
/// features the iteration can stall within tolerance where `f(x)` is not
/// small and report a false root. Prefer [`Brent`] when a bracket is
/// available.
#[derive(Debug, Copy, Clone, PartialEq)]
pub struct Newton {
    max_iter: u32,
    /// Absolute tolerance on the iteration step.
    abs_tol: f64,
    /// Relative tolerance on the iteration step.
    rel_tol: f64,
}

impl Default for Newton {
    fn default() -> Self {
        Self {
            max_iter: 50,
            abs_tol: sqrt(f64::EPSILON),
            rel_tol: sqrt(f64::EPSILON),
        }
    }
}

impl Newton {
    /// Sets the maximum number of iterations.
    pub fn with_max_iter(mut self, max_iter: u32) -> Self {
        self.max_iter = max_iter;
        self
    }

    /// Sets the absolute tolerance on the iteration step.
    pub fn with_abs_tol(mut self, abs_tol: f64) -> Self {
        self.abs_tol = abs_tol;
        self
    }

    /// Sets the relative tolerance on the iteration step.
    pub fn with_rel_tol(mut self, rel_tol: f64) -> Self {
        self.rel_tol = rel_tol;
        self
    }
}

impl FindRootWithDerivative for Newton {
    fn find_with_derivative(
        &self,
        f: impl Callback,
        derivative: impl Callback,
        initial_guess: f64,
    ) -> Result<f64, RootFinderError> {
        let mut p0 = initial_guess;
        let mut last: Option<(f64, f64)> = None;
        for _ in 0..self.max_iter {
            let fx = eval_finite(&f, p0)?;
            // An initial guess that is already a root is returned directly,
            // avoiding a 0/0 update at a stationary point.
            if fx == 0.0 {
                return Ok(p0);
            }
            last = Some((p0, fx));
            let dfx = derivative.call(p0)?;
            if !dfx.is_finite() {
                return Err(RootFinderError::NonFiniteDerivative { x: p0, value: dfx });
            }
            let p = p0 - fx / dfx;
            if !p.is_finite() {
                return Err(RootFinderError::DivergedStep { x: p0 });
            }
            if approx_eq!(p, p0, rtol <= self.rel_tol, atol <= self.abs_tol) {
                return Ok(p);
            }
            p0 = p;
        }
        // Report the last point where `f` was actually evaluated rather than
        // re-evaluating a stepped final iterate, which may lie outside the
        // callback's valid domain and turn non-convergence into a callback error
        // or panic. When `max_iter == 0` the loop never runs, so fall back to the
        // still-in-domain initial guess.
        let (x, residual) = match last {
            Some(pair) => pair,
            None => (p0, eval_finite(&f, p0)?),
        };
        Err(RootFinderError::NotConverged {
            iterations: self.max_iter,
            x,
            residual,
        })
    }
}

/// Brent's method for bracketed root-finding.
///
/// The tolerances bound the root location `x`, not the residual `f(x)`: the
/// returned root is accurate to within `abs_tol + rel_tol * |x|`, independent of
/// how the objective is scaled.
#[derive(Debug, Copy, Clone, PartialEq)]
pub struct Brent {
    max_iter: u32,
    /// Absolute tolerance on the root location.
    abs_tol: f64,
    /// Relative tolerance on the root location.
    rel_tol: f64,
}

impl Default for Brent {
    fn default() -> Self {
        Self {
            max_iter: 100,
            abs_tol: 1e-6,
            rel_tol: sqrt(f64::EPSILON),
        }
    }
}

impl Brent {
    /// Sets the maximum number of iterations.
    pub fn with_max_iter(mut self, max_iter: u32) -> Self {
        self.max_iter = max_iter;
        self
    }

    /// Sets the absolute tolerance on the root location.
    pub fn with_abs_tol(mut self, abs_tol: f64) -> Self {
        self.abs_tol = abs_tol;
        self
    }

    /// Sets the relative tolerance on the root location.
    pub fn with_rel_tol(mut self, rel_tol: f64) -> Self {
        self.rel_tol = rel_tol;
        self
    }
}

impl FindBracketedRoot for Brent {
    fn find_in_bracket_with_values(
        &self,
        f: impl Callback,
        bracket: (f64, f64),
        values: (f64, f64),
    ) -> Result<f64, RootFinderError> {
        let mut fblk = 0.0;
        let mut xblk = 0.0;
        let (mut xpre, mut xcur) = bracket;
        let mut spre = 0.0;
        let mut scur = 0.0;

        let (mut fpre, mut fcur) = values;

        if !fpre.is_finite() {
            return Err(RootFinderError::NonFinite {
                x: xpre,
                value: fpre,
            });
        }
        if !fcur.is_finite() {
            return Err(RootFinderError::NonFinite {
                x: xcur,
                value: fcur,
            });
        }

        // An endpoint that is exactly a root is returned directly.
        if fpre == 0.0 {
            return Ok(xpre);
        }
        if fcur == 0.0 {
            return Ok(xcur);
        }

        // The endpoints must straddle the root. Comparing the sign bits of the
        // two finite, non-zero values avoids the underflow and NaN hazards of
        // testing the sign of their product.
        if fpre.is_sign_negative() == fcur.is_sign_negative() {
            return Err(RootFinderError::NotInBracket);
        }

        for _ in 0..self.max_iter {
            // Compare sign bits rather than the product to avoid the underflow
            // that loses the sign of very small opposite-sign residuals.
            if fpre.is_sign_negative() != fcur.is_sign_negative() {
                xblk = xpre;
                fblk = fpre;
                spre = xcur - xpre;
                scur = xcur - xpre;
            }

            if abs(fblk) < abs(fcur) {
                xpre = xcur;
                xcur = xblk;
                xblk = xpre;
                fpre = fcur;
                fcur = fblk;
                fblk = fpre;
            }

            let delta = (self.abs_tol + self.rel_tol * abs(xcur)) / 2.0;
            let sbis = (xblk - xcur) / 2.0;

            if fcur == 0.0 || abs(sbis) < delta {
                return Ok(xcur);
            }

            if abs(spre) > delta && abs(fcur) < abs(fpre) {
                let stry = if approx_eq!(xpre, xblk, rtol <= self.rel_tol) {
                    // interpolate
                    -fcur * (xcur - xpre) / (fcur - fpre)
                } else {
                    // extrapolate
                    let dpre = (fpre - fcur) / (xpre - xcur);
                    let dblk = (fblk - fcur) / (xblk - xcur);
                    -fcur * (fblk * dblk - fpre * dpre) / (dblk * dpre * (fblk - fpre))
                };

                if 2.0 * abs(stry) < abs(spre).min(3.0 * abs(sbis) - delta) {
                    spre = scur;
                    scur = stry;
                } else {
                    // bisect
                    spre = sbis;
                    scur = sbis;
                }
            } else {
                // bisect
                spre = sbis;
                scur = sbis;
            }

            xpre = xcur;
            fpre = fcur;

            if abs(scur) > delta {
                xcur += scur
            } else {
                xcur += if sbis > 0.0 { delta } else { -delta };
            }

            fcur = eval_finite(&f, xcur)?;
        }

        Err(RootFinderError::NotConverged {
            iterations: self.max_iter,
            x: xcur,
            residual: fcur,
        })
    }
}

/// Safeguarded ("rtsafe") Newton method for bracketed root-finding.
///
/// Takes a Newton step when it lands inside the current bracket and shrinks the
/// interval quickly enough; otherwise bisects. This retains the guaranteed
/// convergence of bisection while gaining Newton's quadratic rate near the
/// root. A zero, non-finite, or otherwise unhelpful derivative simply yields a
/// bisection step, so the method never fails as long as a valid bracket is
/// maintained.
///
/// The tolerances bound the root location `x`, not the residual `f(x)`: the
/// returned root is accurate to within `abs_tol + rel_tol * |x|`, independent of
/// how the objective is scaled.
///
/// # References
///
/// - Press et al., *Numerical Recipes*, 3rd ed., §9.4 ("Newton-Raphson Method
///   Using Derivative", `rtsafe`).
#[derive(Debug, Copy, Clone, PartialEq)]
pub struct BracketedNewton {
    max_iter: u32,
    /// Absolute tolerance on the root location.
    abs_tol: f64,
    /// Relative tolerance on the root location.
    rel_tol: f64,
}

impl Default for BracketedNewton {
    fn default() -> Self {
        Self {
            max_iter: 100,
            abs_tol: 1e-6,
            rel_tol: sqrt(f64::EPSILON),
        }
    }
}

impl BracketedNewton {
    /// Sets the maximum number of iterations.
    pub fn with_max_iter(mut self, max_iter: u32) -> Self {
        self.max_iter = max_iter;
        self
    }

    /// Sets the absolute tolerance on the root location.
    pub fn with_abs_tol(mut self, abs_tol: f64) -> Self {
        self.abs_tol = abs_tol;
        self
    }

    /// Sets the relative tolerance on the root location.
    pub fn with_rel_tol(mut self, rel_tol: f64) -> Self {
        self.rel_tol = rel_tol;
        self
    }
}

impl FindBracketedRootWithDerivative for BracketedNewton {
    fn find_in_bracket_with_derivative_values(
        &self,
        f: impl CallbackWithDerivative,
        bracket: (f64, f64),
        values: (f64, f64),
    ) -> Result<f64, RootFinderError> {
        let (x1, x2) = bracket;
        let (f1, f2) = values;

        if !f1.is_finite() {
            return Err(RootFinderError::NonFinite { x: x1, value: f1 });
        }
        if !f2.is_finite() {
            return Err(RootFinderError::NonFinite { x: x2, value: f2 });
        }

        // An endpoint that is exactly a root is returned directly.
        if f1 == 0.0 {
            return Ok(x1);
        }
        if f2 == 0.0 {
            return Ok(x2);
        }

        // The endpoints must straddle the root. Comparing sign bits avoids the
        // underflow that loses the sign of very small opposite-sign residuals.
        if f1.is_sign_negative() == f2.is_sign_negative() {
            return Err(RootFinderError::NotInBracket);
        }

        // Orient the bracket so that `f(xl) < 0 < f(xh)`.
        let (mut xl, mut xh) = if f1.is_sign_negative() {
            (x1, x2)
        } else {
            (x2, x1)
        };

        let mut rts = 0.5 * (x1 + x2);
        let mut dx_old = abs(x2 - x1);
        let mut dx = dx_old;
        let (mut fx, mut dfx) = eval_finite_with_derivative(&f, rts)?;
        // The midpoint may already be the root, in which case the derivative is
        // never used (and could give a 0/0 step if it also vanishes there).
        if fx == 0.0 {
            return Ok(rts);
        }

        for _ in 0..self.max_iter {
            let delta = self.abs_tol + self.rel_tol * abs(rts);

            // Bisect when the derivative is unusable, when the Newton iterate
            // would leave the bracket, or when it is not shrinking the interval
            // fast enough; otherwise take the Newton step.
            let bisect = !dfx.is_finite()
                || ((rts - xh) * dfx - fx) * ((rts - xl) * dfx - fx) > 0.0
                || abs(2.0 * fx) > abs(dx_old * dfx);

            if bisect {
                dx_old = dx;
                dx = 0.5 * (xh - xl);
                rts = xl + dx;
                // The bisection step is below the representable resolution.
                if xl == rts {
                    return Ok(rts);
                }
            } else {
                dx_old = dx;
                dx = fx / dfx;
                let prev = rts;
                rts -= dx;
                if prev == rts {
                    return Ok(rts);
                }
            }

            if abs(dx) < delta {
                return Ok(rts);
            }

            (fx, dfx) = eval_finite_with_derivative(&f, rts)?;
            if fx == 0.0 {
                return Ok(rts);
            }
            // Maintain the sign-oriented bracket around the new iterate.
            if fx.is_sign_negative() {
                xl = rts;
            } else {
                xh = rts;
            }
        }

        Err(RootFinderError::NotConverged {
            iterations: self.max_iter,
            x: rts,
            residual: fx,
        })
    }
}

#[cfg(test)]
mod tests {
    use alloc::string::ToString;
    use core::f64::consts::PI;
    use lox_approx::assert_approx_eq;

    use super::*;
    use crate::error::BoxedError;
    use crate::math::callback::Fallible;
    use crate::math::float::{cos, sin};

    type Result = core::result::Result<f64, BoxedError>;

    #[test]
    fn test_newton_kepler() {
        fn mean_to_ecc(mean: f64, eccentricity: f64) -> core::result::Result<f64, RootFinderError> {
            let newton = Newton::default();
            newton.find_with_derivative(
                |e: f64| e - eccentricity * sin(e) - mean,
                |e: f64| 1.0 - eccentricity * cos(e),
                mean,
            )
        }
        let act = mean_to_ecc(PI / 2.0, 0.3).expect("should converge");
        assert_approx_eq!(act, 1.85846841205333, rtol <= 1e-8);
    }

    #[test]
    fn test_newton_cubic() {
        let newton = Newton::default();
        let act = newton
            .find_with_derivative(
                |x: f64| powi(x, 3) + 4.0 * powi(x, 2) - 10.0,
                |x: f64| 2.0 * powi(x, 2) + 8.0 * x,
                1.5,
            )
            .expect("should converge");
        assert_approx_eq!(act, 1.3652300134140969, rtol <= 1e-8);
    }

    #[test]
    fn test_newton_exact_root_initial_guess() {
        // f(x) = x^2, f'(x) = 2x. At x = 0 both are zero; the guess is already
        // the root and must be returned rather than producing 0/0 = NaN.
        let newton = Newton::default();
        let act = newton
            .find_with_derivative(|x: f64| powi(x, 2), |x: f64| 2.0 * x, 0.0)
            .expect("guess is already the root");
        assert_eq!(act, 0.0);
    }

    #[test]
    fn test_newton_zero_derivative_diverges() {
        // f(x) = x^2 + 1 has no real root; at x = 0 the derivative is zero, so
        // the update diverges and must be reported as a diverged step.
        let newton = Newton::default();
        let err = newton
            .find_with_derivative(|x: f64| powi(x, 2) + 1.0, |x: f64| 2.0 * x, 0.0)
            .unwrap_err();
        assert!(matches!(err, RootFinderError::DivergedStep { x } if x == 0.0));
    }

    #[test]
    fn test_newton_large_root_relative_tolerance() {
        // A root at 1e8 is unreachable by an absolute step tolerance of
        // sqrt(EPSILON); the relative tolerance lets it converge.
        let newton = Newton::default();
        let act = newton
            .find_with_derivative(|x: f64| powi(x, 2) - 1e16, |x: f64| 2.0 * x, 9e7)
            .expect("should converge");
        assert_approx_eq!(act, 1e8, rtol <= 1e-9);
    }

    #[test]
    fn test_steffensen_exact_root_initial_guess() {
        // f(x) = x^2 - 4 has a root at 2; the guess is already the root.
        let steffensen = Steffensen::default();
        let act = steffensen
            .find(|x: f64| powi(x, 2) - 4.0, 2.0)
            .expect("guess is already the root");
        assert_eq!(act, 2.0);
    }

    #[test]
    fn test_steffensen_zero_denominator_diverges() {
        // A constant non-zero residual makes the Aitken denominator vanish; the
        // update diverges and must be reported as a diverged step.
        let steffensen = Steffensen::default();
        let err = steffensen.find(|_x: f64| 1.0, 0.0).unwrap_err();
        assert!(matches!(err, RootFinderError::DivergedStep { x } if x == 0.0));
    }

    #[test]
    fn test_steffensen_cubic() {
        let steffensen = Steffensen::default();
        let act = steffensen
            .find(|x: f64| powi(x, 3) + 4.0 * powi(x, 2) - 10.0, 1.5)
            .expect("should converge");
        assert_approx_eq!(act, 1.3652300134140969, rtol <= 1e-8);
    }

    #[test]
    fn test_brent_cubic() {
        let brent = Brent::default();
        let act = brent
            .find_in_bracket(|x: f64| powi(x, 3) + 4.0 * powi(x, 2) - 10.0, (1.0, 1.5))
            .expect("should converge");
        assert_approx_eq!(act, 1.3652300134140969, rtol <= 1e-8);
    }

    #[test]
    #[should_panic(expected = "derivative failed")]
    fn test_newton_kepler_callback_error() {
        let newton = Newton::default();
        newton
            .find_with_derivative(
                |e: f64| e,
                Fallible(|_e: f64| -> Result { Err("derivative failed".into()) }),
                1.0,
            )
            .unwrap();
    }

    #[test]
    #[should_panic(expected = "f failed")]
    fn test_steffensen_cubic_error() {
        let steffensen = Steffensen::default();
        // function errors immediately
        steffensen
            .find(Fallible(|_x| -> Result { Err("f failed".into()) }), 1.0)
            .unwrap();
    }

    #[test]
    #[should_panic(expected = "negative x")]
    fn test_brent_cubic_error() {
        let brent = Brent::default();
        // error at bracket endpoint, then during iteration
        brent
            .find_in_bracket(
                Fallible(|x: f64| -> Result {
                    if x.is_sign_negative() {
                        Err("negative x".into())
                    } else {
                        Ok(x * x - 2.0)
                    }
                }),
                (-1.0, 2.0),
            )
            .unwrap();
    }

    #[test]
    fn test_find_in_bracket_with_values_reuses_endpoints() {
        use core::cell::Cell;

        let f = |x: f64| powi(x, 3) + 4.0 * powi(x, 2) - 10.0;
        let (a, b) = (1.0, 1.5);
        let fa = powi(a, 3) + 4.0 * powi(a, 2) - 10.0;
        let fb = powi(b, 3) + 4.0 * powi(b, 2) - 10.0;

        // The value-reuse entry point agrees with the recomputing one.
        let brent = Brent::default();
        let via_values = brent
            .find_in_bracket_with_values(f, (a, b), (fa, fb))
            .expect("should converge");
        let via_recompute = brent.find_in_bracket(f, (a, b)).expect("should converge");
        assert_approx_eq!(via_values, via_recompute, rtol <= 1e-12);

        // The supplied endpoints are not re-evaluated.
        let count = Cell::new(0usize);
        let counting = |x: f64| {
            if x == a || x == b {
                count.set(count.get() + 1);
            }
            powi(x, 3) + 4.0 * powi(x, 2) - 10.0
        };
        Brent::default()
            .find_in_bracket_with_values(counting, (a, b), (fa, fb))
            .expect("should converge");
        assert_eq!(count.get(), 0, "endpoints must not be re-evaluated");
    }

    #[test]
    fn test_brent_rejects_non_finite_endpoint() {
        let brent = Brent::default();
        let err = brent
            .find_in_bracket_with_values(|_x: f64| 1.0, (0.0, 1.0), (f64::NAN, 1.0))
            .unwrap_err();
        assert!(matches!(err, RootFinderError::NonFinite { .. }));
    }

    #[test]
    fn test_brent_rejects_same_sign_underflowing_bracket() {
        // Same-sign endpoints whose product underflows to 0.0 must still be
        // rejected rather than accepted as a bracket (and returned as a root).
        let brent = Brent::default();
        let err = brent
            .find_in_bracket(|x: f64| 1e-200 * (x + 1.0), (0.0, 1.0))
            .unwrap_err();
        assert!(matches!(err, RootFinderError::NotInBracket));
    }

    #[test]
    fn test_brent_scale_independent() {
        // A heavily down-scaled objective: f(0) = -1e-6 must not be mistaken for
        // a root by a residual tolerance. The true root is at x = 1e6.
        let brent = Brent::default();
        let act = brent
            .find_in_bracket(|x: f64| 1e-12 * (x - 1e6), (0.0, 2e6))
            .expect("should converge");
        assert_approx_eq!(act, 1e6, rtol <= 1e-5);
    }

    #[test]
    fn test_brent_builder_tolerances() {
        // Custom tolerances set via the builder still locate the root.
        let brent = Brent::default()
            .with_abs_tol(1e-2)
            .with_rel_tol(1e-8)
            .with_max_iter(50);
        let act = brent
            .find_in_bracket(|x: f64| powi(x, 3) + 4.0 * powi(x, 2) - 10.0, (1.0, 1.5))
            .expect("should converge");
        assert_approx_eq!(act, 1.3652300134140969, rtol <= 1e-2);
    }

    #[test]
    fn test_newton_builder_max_iter() {
        // A single iteration is not enough to converge on the cubic root, and
        // the configured cap is reported back in the error.
        let newton = Newton::default().with_max_iter(1);
        let err = newton
            .find_with_derivative(
                |x: f64| powi(x, 3) + 4.0 * powi(x, 2) - 10.0,
                |x: f64| 2.0 * powi(x, 2) + 8.0 * x,
                1.5,
            )
            .unwrap_err();
        assert!(matches!(
            err,
            RootFinderError::NotConverged { iterations: 1, .. }
        ));
    }

    #[test]
    fn test_brent_tiny_opposite_sign_residuals() {
        // Residuals so small their product underflows to a signed zero; the
        // in-loop bracketing must recognise the sign change from the sign bits
        // rather than the product, which would otherwise return the endpoint.
        let brent = Brent::default();
        let act = brent
            .find_in_bracket(|x: f64| 1e-200 * (x - 0.5), (0.0, 1.0))
            .expect("should converge");
        assert_approx_eq!(act, 0.5, rtol <= 1e-8);
    }

    #[test]
    fn test_brent_interior_non_finite() {
        // Finite, opposite-sign endpoints but a non-finite value at an interior
        // point must be reported rather than silently accepted as a root.
        let calls = core::cell::Cell::new(0u32);
        let brent = Brent::default();
        let err = brent
            .find_in_bracket(
                |x: f64| {
                    let n = calls.get();
                    calls.set(n + 1);
                    // The first two calls evaluate the bracket endpoints.
                    if n >= 2 { f64::NAN } else { x - 0.5 }
                },
                (0.0, 1.0),
            )
            .unwrap_err();
        assert!(matches!(err, RootFinderError::NonFinite { .. }));
    }

    #[test]
    fn test_newton_not_converged_does_not_re_evaluate() {
        // The final iterate steps outside the callback's valid domain. Reaching
        // the iteration cap must still yield NotConverged, not a callback error
        // from re-evaluating the un-checked final iterate.
        let newton = Newton::default().with_max_iter(1);
        let err = newton
            .find_with_derivative(
                Fallible(|x: f64| -> Result {
                    if x > 5.0 {
                        Err("out of domain".into())
                    } else {
                        Ok(x * x - 2.0)
                    }
                }),
                |_x: f64| 0.1,
                1.0,
            )
            .unwrap_err();
        assert!(matches!(err, RootFinderError::NotConverged { .. }));
    }

    #[test]
    fn test_newton_zero_max_iter_reports_real_residual() {
        // With no iterations the solver still reports a finite residual measured
        // at the initial guess, not NaN.
        let newton = Newton::default().with_max_iter(0);
        let err = newton
            .find_with_derivative(|x: f64| x * x - 2.0, |x: f64| 2.0 * x, 1.0)
            .unwrap_err();
        match err {
            RootFinderError::NotConverged { x, residual, .. } => {
                assert_eq!(x, 1.0);
                assert_approx_eq!(residual, -1.0, atol <= 1e-12);
            }
            other => panic!("expected NotConverged, got {other:?}"),
        }
    }

    #[test]
    fn test_steffensen_zero_max_iter_reports_real_residual() {
        let steffensen = Steffensen::default().with_max_iter(0);
        let err = steffensen.find(|x: f64| x * x - 2.0, 1.0).unwrap_err();
        match err {
            RootFinderError::NotConverged { x, residual, .. } => {
                assert_eq!(x, 1.0);
                assert_approx_eq!(residual, -1.0, atol <= 1e-12);
            }
            other => panic!("expected NotConverged, got {other:?}"),
        }
    }

    #[test]
    fn test_root_finder_error_display() {
        let not_converged = RootFinderError::NotConverged {
            iterations: 7,
            x: 1.5,
            residual: -0.25,
        };
        assert_eq!(
            not_converged.to_string(),
            "not converged after 7 iterations at x = 1.5, residual -0.25"
        );
        assert_eq!(
            RootFinderError::NotInBracket.to_string(),
            "root not in bracket"
        );
        assert_eq!(
            RootFinderError::NonFinite {
                x: 1.5,
                value: f64::INFINITY
            }
            .to_string(),
            "function returned a non-finite value (inf) at x = 1.5"
        );
        assert_eq!(
            RootFinderError::NonFiniteDerivative {
                x: 1.5,
                value: f64::NAN
            }
            .to_string(),
            "derivative returned a non-finite value (NaN) at x = 1.5"
        );
        assert_eq!(
            RootFinderError::DivergedStep { x: 0.5 }.to_string(),
            "iteration step diverged at x = 0.5"
        );
    }

    #[test]
    fn test_steffensen_builder_tolerances() {
        // Custom tolerances set via the builder still locate the root.
        let steffensen = Steffensen::default()
            .with_abs_tol(1e-4)
            .with_rel_tol(1e-8)
            .with_max_iter(100);
        let act = steffensen
            .find(|x: f64| powi(x, 3) + 4.0 * powi(x, 2) - 10.0, 1.5)
            .expect("should converge");
        assert_approx_eq!(act, 1.3652300134140969, rtol <= 1e-3);
    }

    #[test]
    fn test_newton_builder_tolerances() {
        let newton = Newton::default().with_abs_tol(1e-4).with_rel_tol(1e-8);
        let act = newton
            .find_with_derivative(
                |x: f64| powi(x, 3) + 4.0 * powi(x, 2) - 10.0,
                |x: f64| 2.0 * powi(x, 2) + 8.0 * x,
                1.5,
            )
            .expect("should converge");
        assert_approx_eq!(act, 1.3652300134140969, rtol <= 1e-3);
    }

    #[test]
    fn test_steffensen_non_finite_function_value() {
        // A non-finite value at the initial guess is reported immediately.
        let steffensen = Steffensen::default();
        let err = steffensen.find(|_x: f64| f64::INFINITY, 1.0).unwrap_err();
        assert!(matches!(err, RootFinderError::NonFinite { .. }));
    }

    #[test]
    fn test_steffensen_non_finite_aitken_probe() {
        // The first sample is finite, but the displaced Aitken probe f(p0 + f(p0))
        // is non-finite and must be reported.
        let steffensen = Steffensen::default();
        let err = steffensen
            .find(|x: f64| if x >= 2.0 { f64::INFINITY } else { x }, 1.5)
            .unwrap_err();
        assert!(matches!(err, RootFinderError::NonFinite { .. }));
    }

    #[test]
    fn test_steffensen_not_converged_reports_last_evaluation() {
        // One iteration cannot converge; the error reports the point that was
        // actually evaluated (the initial guess) and its residual.
        let steffensen = Steffensen::default().with_max_iter(1);
        let err = steffensen
            .find(|x: f64| powi(x, 3) + 4.0 * powi(x, 2) - 10.0, 1.5)
            .unwrap_err();
        match err {
            RootFinderError::NotConverged {
                iterations: 1,
                x,
                residual,
            } => {
                assert_eq!(x, 1.5);
                assert_approx_eq!(residual, 2.375, atol <= 1e-9);
            }
            other => panic!("expected NotConverged, got {other:?}"),
        }
    }

    #[test]
    fn test_newton_non_finite_function_value() {
        let newton = Newton::default();
        let err = newton
            .find_with_derivative(|_x: f64| f64::INFINITY, |_x: f64| 1.0, 1.0)
            .unwrap_err();
        assert!(matches!(err, RootFinderError::NonFinite { .. }));
    }

    #[test]
    fn test_newton_non_finite_derivative_value() {
        // A finite function value but a non-finite derivative is reported
        // against the derivative, not the objective.
        let newton = Newton::default();
        let err = newton
            .find_with_derivative(|x: f64| x, |_x: f64| f64::INFINITY, 1.0)
            .unwrap_err();
        assert!(matches!(
            err,
            RootFinderError::NonFiniteDerivative { x, .. } if x == 1.0
        ));
    }

    #[test]
    fn test_brent_rejects_non_finite_second_endpoint() {
        // The first endpoint is finite; the second is not.
        let brent = Brent::default();
        let err = brent
            .find_in_bracket_with_values(|_x: f64| 1.0, (0.0, 1.0), (1.0, f64::NAN))
            .unwrap_err();
        assert!(matches!(err, RootFinderError::NonFinite { .. }));
    }

    #[test]
    fn test_brent_endpoint_is_root() {
        let brent = Brent::default();
        // The lower endpoint is exactly the root.
        let lo = brent
            .find_in_bracket_with_values(|x: f64| x, (0.0, 1.0), (0.0, 1.0))
            .expect("lower endpoint is the root");
        assert_eq!(lo, 0.0);
        // The upper endpoint is exactly the root.
        let hi = brent
            .find_in_bracket_with_values(|x: f64| x - 1.0, (0.0, 1.0), (-1.0, 0.0))
            .expect("upper endpoint is the root");
        assert_eq!(hi, 1.0);
    }

    #[test]
    fn test_brent_not_converged() {
        // A single iteration cannot narrow a wide bracket below the tolerance.
        let brent = Brent::default().with_max_iter(1);
        let err = brent
            .find_in_bracket(|x: f64| powi(x, 3) - 0.5, (-1e6, 1e6))
            .unwrap_err();
        assert!(matches!(
            err,
            RootFinderError::NotConverged { iterations: 1, .. }
        ));
    }

    #[test]
    fn test_zero_max_iter_non_finite_initial_guess() {
        // With no iterations the initial guess is still evaluated; a non-finite
        // value there is reported as NonFinite, not NotConverged.
        let newton = Newton::default().with_max_iter(0);
        let err = newton
            .find_with_derivative(|_x: f64| f64::INFINITY, |_x: f64| 1.0, 1.0)
            .unwrap_err();
        assert!(matches!(err, RootFinderError::NonFinite { .. }));
    }

    #[test]
    fn test_brent_converges_on_stiff_function() {
        // A high-degree monomial has extreme curvature near its root, which
        // exercises Brent's switching between interpolation and bisection.
        let brent = Brent::default();
        let root = brent
            .find_in_bracket(|x: f64| powi(x, 15) - 0.5, (0.0, 1.0))
            .expect("should converge");
        assert!(abs(powi(root, 15) - 0.5) < 1e-3);
    }

    #[test]
    fn test_brent_in_loop_callback_error() {
        // Endpoints evaluate cleanly, but an interior evaluation fails; the
        // callback error must propagate.
        let brent = Brent::default();
        let err = brent
            .find_in_bracket_with_values(
                Fallible(|x: f64| -> Result {
                    if x == -1.0 || x == 2.0 {
                        Ok(x * x - 2.0)
                    } else {
                        Err("interior failure".into())
                    }
                }),
                (-1.0, 2.0),
                (-1.0, 2.0),
            )
            .unwrap_err();
        assert!(matches!(err, RootFinderError::Callback(_)));
    }

    // -----------------------------------------------------------------------
    // BracketedNewton (safeguarded Newton / rtsafe)
    // -----------------------------------------------------------------------

    #[test]
    fn test_bracketed_newton_cubic() {
        let solver = BracketedNewton::default();
        let act = solver
            .find_in_bracket_with_derivative(
                |x: f64| {
                    (
                        powi(x, 3) + 4.0 * powi(x, 2) - 10.0,
                        3.0 * powi(x, 2) + 8.0 * x,
                    )
                },
                (1.0, 1.5),
            )
            .expect("should converge");
        assert_approx_eq!(act, 1.3652300134140969, rtol <= 1e-8);
    }

    #[test]
    fn test_bracketed_newton_matches_brent() {
        // On a smooth function both solvers must land on the same root.
        let f = |x: f64| powi(x, 3) + 4.0 * powi(x, 2) - 10.0;
        let bracket = (1.0, 2.0);

        let newton = BracketedNewton::default()
            .find_in_bracket_with_derivative(
                |x: f64| {
                    (
                        powi(x, 3) + 4.0 * powi(x, 2) - 10.0,
                        3.0 * powi(x, 2) + 8.0 * x,
                    )
                },
                bracket,
            )
            .expect("should converge");
        let brent = Brent::default()
            .find_in_bracket(f, bracket)
            .expect("should converge");
        assert_approx_eq!(newton, brent, rtol <= 1e-9);
    }

    #[test]
    fn test_bracketed_newton_reuses_endpoints() {
        use core::cell::Cell;

        let (a, b) = (1.0, 1.5);
        let fa = powi(a, 3) + 4.0 * powi(a, 2) - 10.0;
        let fb = powi(b, 3) + 4.0 * powi(b, 2) - 10.0;

        let count = Cell::new(0usize);
        let counting = |x: f64| {
            if x == a || x == b {
                count.set(count.get() + 1);
            }
            powi(x, 3) + 4.0 * powi(x, 2) - 10.0
        };
        BracketedNewton::default()
            .find_in_bracket_with_derivative_values(
                |x: f64| (counting(x), 3.0 * powi(x, 2) + 8.0 * x),
                (a, b),
                (fa, fb),
            )
            .expect("should converge");
        assert_eq!(count.get(), 0, "endpoints must not be re-evaluated");
    }

    #[test]
    fn test_bracketed_newton_endpoint_is_root() {
        let solver = BracketedNewton::default();
        let lo = solver
            .find_in_bracket_with_derivative_values(|x: f64| (x, 1.0), (0.0, 1.0), (0.0, 1.0))
            .expect("lower endpoint is the root");
        assert_eq!(lo, 0.0);
    }

    #[test]
    fn test_bracketed_newton_rejects_non_bracket() {
        let solver = BracketedNewton::default();
        let err = solver
            .find_in_bracket_with_derivative_values(
                |x: f64| (x * x + 1.0, 2.0 * x),
                (1.0, 2.0),
                (2.0, 5.0),
            )
            .unwrap_err();
        assert!(matches!(err, RootFinderError::NotInBracket));
    }

    #[test]
    fn test_bracketed_newton_zero_derivative_bisects() {
        // A derivative that is zero at the midpoint would give a 0/0 Newton
        // step; the safeguard must fall back to bisection and still converge.
        // f(x) = x^3 has f'(0) = 0 exactly at the bracket midpoint (-1, 1).
        let solver = BracketedNewton::default();
        let act = solver
            .find_in_bracket_with_derivative(|x: f64| (powi(x, 3), 3.0 * powi(x, 2)), (-1.0, 1.0))
            .expect("should converge via bisection fallback");
        assert_approx_eq!(act, 0.0, atol <= 1e-6);
    }

    #[test]
    fn test_bracketed_newton_non_finite_derivative_bisects() {
        // A derivative that blows up must not fail the solve: the safeguard
        // bisects instead, keeping the bracket and converging on the root.
        let solver = BracketedNewton::default();
        let act = solver
            .find_in_bracket_with_derivative(|x: f64| (x - 0.5, f64::INFINITY), (0.0, 1.0))
            .expect("should converge despite non-finite derivative");
        assert_approx_eq!(act, 0.5, atol <= 1e-6);
    }

    #[test]
    fn test_bracketed_newton_scale_independent() {
        // A heavily down-scaled objective must not be mistaken for a root by a
        // residual tolerance; the true root is at x = 1e6.
        let solver = BracketedNewton::default();
        let act = solver
            .find_in_bracket_with_derivative(|x: f64| (1e-12 * (x - 1e6), 1e-12), (0.0, 2e6))
            .expect("should converge");
        assert_approx_eq!(act, 1e6, rtol <= 1e-5);
    }

    #[test]
    fn test_bracketed_newton_rejects_non_finite_endpoint() {
        let solver = BracketedNewton::default();
        let err = solver
            .find_in_bracket_with_derivative_values(
                |_x: f64| (1.0, 1.0),
                (0.0, 1.0),
                (f64::NAN, 1.0),
            )
            .unwrap_err();
        assert!(matches!(err, RootFinderError::NonFinite { .. }));
    }

    #[test]
    fn test_bracketed_newton_not_converged() {
        // A single iteration cannot narrow a wide bracket below the tolerance.
        let solver = BracketedNewton::default().with_max_iter(1);
        let err = solver
            .find_in_bracket_with_derivative(
                |x: f64| (powi(x, 3) - 0.5, 3.0 * powi(x, 2)),
                (-1e6, 1e6),
            )
            .unwrap_err();
        assert!(matches!(
            err,
            RootFinderError::NotConverged { iterations: 1, .. }
        ));
    }

    #[test]
    fn test_bracketed_newton_in_loop_callback_error() {
        // Endpoints evaluate cleanly, but an interior evaluation fails; the
        // callback error must propagate.
        struct FailsInterior;

        impl CallbackWithDerivative for FailsInterior {
            fn call(&self, x: f64) -> core::result::Result<(f64, f64), LoxError> {
                if x == -1.0 || x == 2.0 {
                    Ok((x * x - 2.0, 2.0 * x))
                } else {
                    Err("interior failure".into())
                }
            }
        }

        let solver = BracketedNewton::default();
        let err = solver
            .find_in_bracket_with_derivative_values(FailsInterior, (-1.0, 2.0), (-1.0, 2.0))
            .unwrap_err();
        assert!(matches!(err, RootFinderError::Callback(_)));
    }
}