rust_poly/poly/roots/
initial_guess.rs

1use std::cmp::Ordering;
2
3use itertools::Itertools;
4
5use crate::{
6    num::{Complex, Zero},
7    scalar::Rational,
8    util::{
9        casting::usize_to_f64,
10        complex::{c_arg, c_exp, c_from_f64, c_min, c_neg, c_powf},
11        doc_macros::{panic_t_from_f64, panic_t_from_int},
12    },
13    Poly, RealScalar,
14};
15
16// TODO: initial guesses Bini (see allow(unused) functions below)
17
18/// Guess close to the root with the smallest magnitude ([Madsen 1973](https://doi.org/10.1007/BF01933524))
19///
20/// # Panics
21#[doc = panic_t_from_int!(r"usize")]
22#[allow(clippy::module_name_repetitions)]
23pub fn initial_guess_smallest<T: RealScalar>(poly: &Poly<T>) -> Complex<T> {
24    debug_assert!(poly.is_normalized());
25    debug_assert!(poly.len_raw() >= 2);
26
27    let small = Rational::recip(T::from_u16(1_000).expect("overflow"));
28    let p_diff = poly.clone().diff();
29    let mut pz = poly.eval(Complex::zero());
30    let mut pdz = p_diff.eval(Complex::zero());
31
32    // avoid divide by zero
33    if pdz.norm_sqr() < small {
34        pz += small.clone();
35        pdz += small;
36    }
37
38    let theta = c_arg(c_neg(pz) / pdz);
39    let mut iter_coeffs = poly.0.iter();
40    let a0 = iter_coeffs.next().expect("infallible");
41
42    let mut guess = iter_coeffs
43        .zip(1..)
44        .map(|(ak, k)| {
45            c_powf(
46                c_exp(Complex::i().scale(theta.clone())).scale((a0 / ak).norm_sqr()),
47                T::one() / T::from_usize(k).expect("overflow"),
48            )
49        })
50        .reduce(c_min)
51        .expect("infallible")
52        .scale(Rational::recip(T::from_u8(2).expect("overflow")));
53
54    if guess.im.is_zero() {
55        // add a small constant because some methods can't converge to
56        // complex roots if the initial guess is real
57        guess += Complex::i().scale(Rational::recip(T::from_u16(1_000).expect("overflow")));
58    }
59    guess
60}
61
62/// Equidistant points around a circle.
63///
64/// The bias parameter controls the radius of the circle. A bias of 0 means the
65/// lower bound is used, whereas a bias of 1 means the upper bound is used. This
66/// parameter can be extrapolated.
67///
68/// The perturbation parameter controls the amount of randomization. At 0 no
69/// randomization is applied. At 1, the radius is picked uniformly between the
70/// lower bound and upper bound, and the angle has a uniform random offset of
71/// `pi / n_odd` where `n_odd` is the degree of the polynomial, rounded up to
72/// the next odd number. In essence, the annulus that contains all the roots is
73/// partitioned into equal slices and then a guess is picked at random in each
74/// of these slices. This parameter can be extrapolated above 1.
75///
76/// # Panics
77#[doc = panic_t_from_f64!()]
78pub fn initial_guesses_circle<T: RealScalar>(
79    poly: &Poly<T>,
80    bias: T,
81    seed: u64,
82    perturbation: T,
83    out: &mut [Complex<T>],
84) {
85    let n = out.len();
86
87    // ensuring n is odd makes the points always asymmetrical, even with even
88    // number of roots. This is important as some methods may get stuck if the
89    // guesses are symmetrical.
90    let n_odd = if n % 2 == 0 { n + 1 } else { n };
91
92    let mut rng = fastrand::Rng::with_seed(seed);
93    let angle_increment = std::f64::consts::TAU / (usize_to_f64(n_odd));
94    let low = lower_bound(poly);
95    let high = upper_bound(poly);
96    let span = high.clone() - low.clone();
97    let radius = high * bias.clone() + low.clone() * (T::one() - bias);
98    let mut angle_accumulator = 0.0;
99    for y in out {
100        let angle = T::from_f64(angle_accumulator).expect("overflow")
101            + T::from_f64(rng.f64().mul_add(angle_increment, -(angle_increment / 2.0)))
102                .expect("overflow")
103                * perturbation.clone();
104        let radius = radius.clone() * (T::one() - perturbation.clone())
105            + (T::from_f64(rng.f64()).expect("overflow") * span.clone() + low.clone())
106                * perturbation.clone();
107        *y = c_from_f64(Complex::from_polar(
108            radius.to_f64().expect("overflow"),
109            angle.to_f64().expect("overflow"),
110        ));
111        angle_accumulator += angle_increment;
112    }
113}
114
115/// The radius of a disk containing all the roots
116///
117/// Uses Deutsch's simple formula \[[McNamee 2005](https://www.researchgate.net/publication/228745231_A_comparison_of_a_priori_bounds_on_real_or_complex_roots_of_polynomials)\]
118fn upper_bound<T: RealScalar>(poly: &Poly<T>) -> T {
119    debug_assert!(
120        poly.degree_raw() >= 2,
121        "upper bound of small degree polynomials is not supported, use explicit solver"
122    );
123    debug_assert!(
124        poly.is_monic(),
125        "Deuthsch's formula requires the polynomial to be monic"
126    );
127
128    let n = poly.len_raw();
129
130    let next_last = poly.0[poly.len_raw() - 2].clone();
131    let coeffs_iter = poly.0.iter().take(n - 2);
132    let coeffs_iter_shifted = poly.0.iter().skip(1).take(n - 2);
133    let max_term = coeffs_iter
134        .zip(coeffs_iter_shifted)
135        .map(|(num, denom)| num / denom)
136        .map(|z| Complex::norm_sqr(&z))
137        .reduce(|acc, z| if z > acc { z } else { acc })
138        .expect("infallible");
139    next_last.norm_sqr() + max_term
140}
141
142/// The radius of a disk containing none of the roots
143fn lower_bound<T: RealScalar>(poly: &Poly<T>) -> T {
144    let mut this = Poly::from_complex_vec(poly.0.iter().cloned().rev().collect_vec());
145    this.make_monic();
146    upper_bound(&this).recip()
147}
148
149/// 2D cross product of OA and OB vectors, i.e. z-component of their 3D cross product.
150/// Returns a positive value, if OAB makes a counter-clockwise turn,
151/// negative for clockwise turn, and zero if the points are collinear.
152///
153/// [From Wiki Books](https://web.archive.org/web/20240617105108/https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain#Python)
154#[allow(unused)]
155fn cross_2d<T: RealScalar>(o: (T, T), a: (T, T), b: (T, T)) -> T {
156    (a.0 - o.0.clone()) * (b.1 - o.1.clone()) - (a.1 - o.1) * (b.0 - o.0)
157}
158
159/// Extract upper envelope of the convex hull of a set of points
160///
161/// [From Wiki Books](https://web.archive.org/web/20240617105108/https://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain#Python)
162#[allow(unused)]
163// TODO: should return an option instead of panic
164fn upper_convex_envelope<T: RealScalar>(points: &mut [(T, T)]) -> Vec<(T, T)> {
165    points.sort_by(
166        |a, b| match (a.0.partial_cmp(&b.0), a.1.partial_cmp(&b.1)) {
167            (None, _) | (Some(Ordering::Equal), None) => panic!("cannot order NaNs"),
168            (Some(Ordering::Equal), Some(ord)) | (Some(ord), _) => ord,
169        },
170    );
171
172    let mut upper: Vec<(T, T)> = vec![];
173    for p in points.iter_mut().rev() {
174        while upper.len() >= 2
175            && cross_2d(
176                upper[upper.len() - 2].clone(),
177                upper[upper.len() - 1].clone(),
178                p.clone(),
179            ) <= T::zero()
180        {
181            upper.pop();
182        }
183        upper.push(p.clone());
184    }
185    upper
186}
rust_poly/poly/roots/initial_guess.rs

rust_poly/poly/roots/
initial_guess.rs
