rust_finprim/utils/
halley.rs

1use crate::error::RootFindingError;
2use crate::FinPrimError;
3use crate::FloatLike;
4
5/// Returns the numerator and denominator of the Halley's step (allows for some checking)
6#[inline(always)]
7pub(crate) fn step_halley<T: FloatLike>(f_x: T, f_prime_x: T, f_prime2_x: T) -> (T, T) {
8    let _2 = T::two();
9    let numerator = _2 * f_x * f_prime_x;
10    let denominator = _2 * f_prime_x * f_prime_x - f_x * f_prime2_x;
11    (numerator, denominator)
12}
13
14pub fn halley<T: FloatLike, F, D, D2>(
15    guess: T,
16    f: F,
17    f_prime: D,
18    f_prime2: D2,
19    tolerance: T,
20    max_iter: u16,
21) -> Result<T, FinPrimError<T>>
22where
23    F: Fn(T) -> T,
24    D: Fn(T) -> T,
25    D2: Fn(T) -> T,
26{
27    let mut x = guess;
28    let mut fx = f(x);
29    for _ in 0..max_iter {
30        if fx.abs() < tolerance {
31            return Ok(x);
32        }
33        let f_prime_x = f_prime(x);
34        let f_prime2 = f_prime2(x);
35
36        let (numerator, denominator) = step_halley(fx, f_prime_x, f_prime2);
37        if denominator.is_zero() {
38            return Err(FinPrimError::RootFindingError(RootFindingError::DivideByZero {
39                last_x: x,
40                last_fx: fx,
41            }));
42        }
43        x -= numerator / denominator;
44        fx = f(x);
45    }
46    Err(FinPrimError::RootFindingError(RootFindingError::FailedToConverge {
47        last_x: x,
48        last_fx: fx,
49    }))
50}
rust_finprim/utils/halley.rs

rust_finprim/utils/
halley.rs
