rustyqlib/core/solvers/
halley.rs

1//! Halley's method: cubic convergence using the first and second
2//! derivatives. Used to polish the inverse normal CDF
3//! ([`inv_norm_cdf`](crate::core::utils::inv_norm_cdf)).
4
5use super::solver_1d::{Root, Solver1d};
6
7/// Halley iteration `x - 2 f f' / (2 f'^2 - f f'')`.
8pub fn halley(
9    cfg: &Solver1d,
10    f: impl Fn(f64) -> f64,
11    df: impl Fn(f64) -> f64,
12    d2f: impl Fn(f64) -> f64,
13    x0: f64,
14) -> Root {
15    let mut x = x0;
16    for i in 0..cfg.max_iter {
17        let fx = f(x);
18        if fx.abs() <= cfg.tol {
19            return Root { x, iterations: i, converged: true };
20        }
21        let dfx = df(x);
22        let denom = 2.0 * dfx * dfx - fx * d2f(x);
23        if denom == 0.0 || !denom.is_finite() {
24            return Root { x, iterations: i, converged: false };
25        }
26        x -= 2.0 * fx * dfx / denom;
27    }
28    Root { x, iterations: cfg.max_iter, converged: f(x).abs() <= cfg.tol }
29}
30
31#[cfg(test)]
32mod tests {
33    use super::*;
34
35    fn f(x: f64) -> f64 {
36        x * x - 2.0
37    }
38    fn df(x: f64) -> f64 {
39        2.0 * x
40    }
41    fn d2f(_: f64) -> f64 {
42        2.0
43    }
44
45    #[test]
46    fn finds_sqrt_two() {
47        let r = halley(&Solver1d::default(), f, df, d2f, 1.0);
48        assert!(r.converged && (r.x - std::f64::consts::SQRT_2).abs() < 1e-12, "{r:?}");
49    }
50
51    #[test]
52    fn converges_in_fewer_iterations_than_newton() {
53        let cfg = Solver1d::new(1e-14, 200);
54        let n = super::super::newton_raphson::newton_raphson(&cfg, f, df, 100.0);
55        let h = halley(&cfg, f, df, d2f, 100.0);
56        assert!(n.converged && h.converged);
57        assert!(h.iterations <= n.iterations, "halley {} vs newton {}", h.iterations, n.iterations);
58    }
59}
rustyqlib/core/solvers/halley.rs

rustyqlib/core/solvers/
halley.rs
