rust_poly/poly/roots/single_root/
halley.rs

1use super::{
2    line_search_accelerate, line_search_decelerate, multiplicity_lagouanelle, LazyDerivatives,
3};
4use crate::{
5    num::{Complex, One, Zero},
6    poly::roots,
7    roots::initial_guess::initial_guess_smallest,
8    util::{
9        self,
10        complex::c_from_f64,
11        doc_macros::{errors_no_converge, panic_t_from_f64, panic_t_from_int},
12    },
13    Poly, RealScalar,
14};
15
16/// Find a single root
17///
18/// # Returns
19/// - vector of roots (usually 1)
20/// - number of evaluations
21///
22/// # Errors
23#[doc = errors_no_converge!()]
24///
25/// # Panics
26#[doc = panic_t_from_f64!()]
27#[doc = panic_t_from_int!(r"usize")]
28#[allow(clippy::items_after_statements)]
29// TODO: remove this when there's a better alternative
30#[allow(clippy::type_complexity)]
31#[allow(clippy::similar_names)]
32pub fn halley<T: RealScalar>(
33    poly: &Poly<T>,
34    epsilon: T,
35    max_iter: Option<usize>,
36    initial_guess: Option<Complex<T>>,
37) -> std::result::Result<(Vec<Complex<T>>, u128), roots::Error<Vec<Complex<T>>>> {
38    let mut eval_counter = 0;
39    let mut guess = initial_guess.unwrap_or_else(|| initial_guess_smallest(poly));
40    let mut guess_old = guess.clone();
41    let mut guess_old_old = guess.clone();
42    let mut guess_delta_old = Complex::one();
43
44    // number of iterations without improvements after which we assume we're in
45    // a cycle set it to a prime number to avoid meta-cycles forming
46    const CYCLE_COUNT_THRESHOLD: usize = 17;
47    let mut cycle_counter = 0;
48    let mut best_guess = guess.clone();
49    let mut best_px_norm = guess.clone().norm_sqr();
50
51    let mut diffs = LazyDerivatives::new(poly);
52
53    // until convergence
54    for i in util::iterator::saturating_counter() {
55        let px = poly.eval(guess.clone());
56        log::trace!("{{current_guess: {guess}, error: {}}}", px.norm_sqr());
57
58        // stopping criterion 1: converged
59        if px.norm_sqr() <= epsilon {
60            return Ok((vec![guess], eval_counter));
61        }
62
63        // stopping criterion 2: no improvement predicted due to numeric precision
64        if i > 3
65            && super::stopping_criterion_garwick(guess.clone(), guess_old.clone(), guess_old_old)
66        {
67            return Ok((vec![guess], eval_counter));
68        }
69
70        // max iter exceeded
71        if max_iter.is_some_and(|max| i >= max) {
72            return Err(roots::Error::NoConverge(vec![guess]));
73        }
74
75        // check for cycles
76        if px.norm_sqr() >= best_px_norm {
77            cycle_counter += 1;
78            if cycle_counter > CYCLE_COUNT_THRESHOLD {
79                // arbitrary constants
80                const ROTATION_RADIANS: f64 = 0.925_024_5;
81                const SCALE: f64 = 5.0;
82
83                cycle_counter = 0;
84                log::trace!("cycle detected, backing off {{current_guess: {guess}, best_guess: {best_guess}}}");
85
86                // TODO: when const trait methods are supported, this should be
87                //       made fully const.
88                let backoff = c_from_f64(Complex::from_polar(SCALE, ROTATION_RADIANS));
89                // reverting to older base guess, but offset
90                guess = best_guess.clone() - guess_delta_old.clone() * backoff;
91            }
92        } else {
93            cycle_counter = 0;
94            best_guess = guess.clone();
95            best_px_norm = px.norm_sqr();
96        }
97
98        let pdx = diffs.get_nth_derivative(1).eval(guess.clone());
99        let pddx = diffs.get_nth_derivative(2).eval(guess.clone());
100        eval_counter += 2;
101        let denom = (pdx.clone() * pdx.clone()).scale(T::from_u8(2).expect("overflow"))
102            - px.clone() * pddx.clone();
103
104        let guess_delta = if denom.is_zero() || pdx.is_zero() {
105            // these are arbitrary, originally chosen by Madsen.
106            const ROTATION_RADIANS: f64 = 0.925_024_5;
107            const SCALE: f64 = 5.0;
108
109            log::trace!("local minimum, backing off");
110
111            // TODO: when const trait methods are supported, this should be
112            //       made fully const.
113            let backoff = c_from_f64(Complex::from_polar(SCALE, ROTATION_RADIANS));
114            <Complex<T> as std::ops::Mul>::mul(guess_delta_old.clone(), backoff)
115        } else {
116            let m = multiplicity_lagouanelle(px.clone(), pdx.clone(), pddx);
117            (m + Complex::one()) * (px.clone() * pdx) / denom
118        };
119
120        const EXPLODE_THRESHOLD: f64 = 5.0;
121        let guess_delta = if guess_delta.norm_sqr()
122            > guess_delta_old.norm_sqr() * T::from_f64(EXPLODE_THRESHOLD).expect("overflow")
123        {
124            // these are arbitrary, originally chosen by Madsen.
125            const ROTATION_RADIANS: f64 = 0.925_024_5;
126            const SCALE: f64 = 5.0;
127
128            log::trace!("exploding gradient, backing off");
129
130            // TODO: when const trait methods are supported, this should be
131            //       made fully const.
132            let backoff = c_from_f64(Complex::from_polar(SCALE, ROTATION_RADIANS))
133                .scale(guess_delta_old.norm_sqr() / guess_delta.norm_sqr());
134            guess_delta * backoff
135        } else {
136            guess_delta
137        };
138
139        eval_counter += 1;
140        let guess_new =
141            if poly.eval(guess.clone() - guess_delta.clone()).norm_sqr() >= px.norm_sqr() {
142                log::trace!("overshooting, shortening step");
143                let res = line_search_decelerate(poly, guess.clone(), guess_delta.clone());
144                eval_counter += res.1;
145                res.0
146            } else {
147                log::trace!("undershooting, lengthening step");
148                let res = line_search_accelerate(poly, guess.clone(), guess_delta.clone());
149                eval_counter += res.1;
150                res.0
151            };
152
153        guess_delta_old = guess_delta;
154        guess_old_old = guess_old;
155        guess_old = guess;
156        guess = guess_new;
157    }
158    unreachable!()
159}
160
161#[cfg(test)]
162mod test {
163    use crate::{num::One, roots::halley_deflate, util::__testing::check_roots, Poly64};
164
165    #[test]
166    pub fn degree_0() {
167        let mut p = Poly64::one();
168        let roots = halley_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
169        assert!(roots.is_empty());
170        assert!(p.is_one());
171    }
172
173    #[test]
174    fn degree_1() {
175        let roots_expected = vec![complex!(1.0)];
176        let mut p = crate::Poly::from_roots(&roots_expected);
177        let roots = halley_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
178        assert!(check_roots(roots, roots_expected, 1E-12));
179    }
180
181    #[test]
182    fn degree_2() {
183        let roots_expected = vec![complex!(1.0), complex!(2.0)];
184        let mut p = crate::Poly::from_roots(&roots_expected);
185        let roots = halley_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
186        assert!(check_roots(roots, roots_expected, 1E-12));
187    }
188
189    #[test]
190    fn degree_3() {
191        let roots_expected = vec![complex!(1.0), complex!(2.0), complex!(3.0)];
192        let mut p = crate::Poly::from_roots(&roots_expected);
193        let roots = halley_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
194        assert!(check_roots(roots, roots_expected, 1E-7));
195    }
196
197    #[test]
198    fn degree_3_complex() {
199        let roots_expected = vec![complex!(1.0), complex!(0.0, 1.0), complex!(0.0, -1.0)];
200        let mut p = crate::Poly::from_roots(&roots_expected);
201        let roots = halley_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
202        assert!(check_roots(roots, roots_expected, 1E-7));
203    }
204
205    #[test]
206    fn degree_5_multiplicity_3() {
207        let roots_expected = vec![
208            complex!(1.0),
209            complex!(2.0),
210            complex!(2.0),
211            complex!(2.0),
212            complex!(3.0),
213        ];
214        let mut p = crate::Poly::from_roots(&roots_expected);
215        let roots = halley_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
216        assert!(
217            check_roots(roots.clone(), roots_expected, 1E-3),
218            "{roots:?}"
219        );
220    }
221
222    #[test]
223    fn degree_5_2_zeros() {
224        let roots_expected = vec![
225            complex!(0.0),
226            complex!(0.0),
227            complex!(1.0),
228            complex!(2.0),
229            complex!(3.0),
230        ];
231        let mut p = crate::Poly::from_roots(&roots_expected);
232        let roots = halley_deflate(&mut p, Some(1E-14), Some(100), &[]).unwrap();
233        assert!(check_roots(roots, roots_expected, 1E-7));
234    }
235}
rust_poly/poly/roots/single_root/halley.rs

rust_poly/poly/roots/single_root/
halley.rs
