feanor_math/algorithms/
newton.rs

1use std::cmp::Ordering;
2
3use oorandom::Rand64;
4
5use crate::algorithms::poly_gcd::squarefree_part::poly_squarefree_part_local;
6use crate::computation::DontObserve;
7use crate::rings::poly::dense_poly::*;
8use crate::homomorphism::*;
9use crate::field::FieldStore;
10use crate::divisibility::{DivisibilityRing, DivisibilityRingStore};
11use crate::integer::*;
12use crate::ring::*;
13use crate::MAX_PROBABILISTIC_REPETITIONS;
14use crate::rings::float_complex::{Complex64, Complex64Base, Complex64El};
15use crate::rings::poly::*;
16
17const NEWTON_MAX_SCALE: u32 = 10;
18const NEWTON_ITERATIONS: usize = 16;
19const ASSUME_RADIUS_TO_APPROX_RADIUS_FACTOR: f64 = 2.;
20
21#[derive(Debug)]
22#[stability::unstable(feature = "enable")]
23pub struct PrecisionError;
24
25#[stability::unstable(feature = "enable")]
26pub fn absolute_error_of_poly_eval<P>(poly_ring: P, f: &El<P>, poly_deg: usize, point: Complex64El, relative_error_point: f64) -> f64
27    where P: RingStore,
28        P::Type: PolyRing + DivisibilityRing,
29        <P::Type as RingExtension>::BaseRing: RingStore<Type = Complex64Base>
30{
31    let CC = Complex64::RING;
32    let mut current = point;
33    let mut current_relative_error = relative_error_point;
34    let mut total_error = 0.;
35    for i in 1..=poly_deg {
36        total_error += CC.abs(*poly_ring.coefficient_at(f, i)) * CC.abs(current) * current_relative_error;
37        // technically, we would have `(1 + current_relative_error)(1 + f64::EPSILON) - 1`, but we ignore `O(f64::EPSILON^2)` terms here
38        current_relative_error += f64::EPSILON;
39        CC.mul_assign(&mut current, point);
40    }
41    return total_error;
42}
43
44fn bound_distance_to_root<P>(approx_root: Complex64El, CCX: P, poly: &El<P>, poly_deg: usize) -> Result<f64, PrecisionError>
45    where P: RingStore,
46        P::Type: PolyRing + DivisibilityRing,
47        <P::Type as RingExtension>::BaseRing: RingStore<Type = Complex64Base>
48{
49    let CC = Complex64::RING;
50    let f = poly;
51    let f_prime = derive_poly(&CCX, &f);
52    
53    let approx_radius = (CC.abs(CCX.evaluate(&f, &approx_root, CC.identity())) + absolute_error_of_poly_eval(&CCX, &f, poly_deg, approx_root, 0.)) / CC.abs(CCX.evaluate(&f_prime, &approx_root, CC.identity()));
54    if !approx_radius.is_finite() {
55        return Err(PrecisionError);
56    }
57
58    let assume_radius = approx_radius * ASSUME_RADIUS_TO_APPROX_RADIUS_FACTOR;
59    // we bound `|f'(x + t) - f'(x)| <= epsilon` for `t` at most `assume_radius` using the taylor series
60    let mut abs_taylor_series_coeffs = Vec::new();
61    let mut current = derive_poly(&CCX, &f_prime);
62    for i in 0..poly_deg.saturating_sub(1) {
63        CCX.inclusion().mul_assign_map(&mut current, CC.from_f64(1. / (i as f64 + 1.)));
64        abs_taylor_series_coeffs.push(CC.abs(CCX.evaluate(&current, &approx_root, CC.identity())));
65        current = derive_poly(&CCX, &current);
66    }
67    let f_prime_bound = abs_taylor_series_coeffs.iter().enumerate()
68        .map(|(i, c)| {
69            c * assume_radius.powi(i as i32 + 1)
70        })
71        .sum::<f64>();
72
73    // The idea is as follows: We have `f(x + t) = f(x) + f'(g(t)) t` where `g(t)` is a value between `x` and `x + t`;
74    // Using `|f'(g(t))| >= |f'(x)| - f_prime_bound` and rearranging it gives `t <= |f(x)| / (|f'(x)| - f_prime_bound)`;
75    // For this to be at most `assume_radius`, it suffices to assume `f_prime_bound <= |f'(x)| - |f(x)|/R = |f'(x)| (1 - approx_radius / assume_radius)`
76    if f_prime_bound > CC.abs(CCX.evaluate(&f_prime, &approx_root, CC.identity())) * (ASSUME_RADIUS_TO_APPROX_RADIUS_FACTOR - 1.) / ASSUME_RADIUS_TO_APPROX_RADIUS_FACTOR {
77        return Err(PrecisionError);
78    }
79
80    return Ok(assume_radius);
81}
82
83fn newton_with_initial<P>(poly_ring: P, f: &El<P>, poly_deg: usize, initial: El<Complex64>) -> Result<(El<Complex64>, f64), PrecisionError>
84    where P: RingStore,
85        P::Type: PolyRing + DivisibilityRing,
86        <P::Type as RingExtension>::BaseRing: RingStore<Type = Complex64Base>
87{
88    let CC = Complex64::RING;
89    let f_prime = derive_poly(&poly_ring, f);
90    let mut current = initial;
91    for _ in 0..NEWTON_ITERATIONS {
92        current = CC.sub(current, CC.div(&poly_ring.evaluate(f, &current, CC.identity()), &poly_ring.evaluate(&f_prime, &current, CC.identity())));
93    }
94    return Ok((current, bound_distance_to_root(current, poly_ring, f, poly_deg)?));
95}
96
97fn find_approximate_complex_root_squarefree<P>(poly_ring: P, f: &El<P>, poly_deg: usize) -> Result<(El<Complex64>, f64), PrecisionError>
98    where P: RingStore,
99        P::Type: PolyRing + DivisibilityRing,
100        <P::Type as RingExtension>::BaseRing: RingStore<Type = Complex64Base>
101{
102    let mut rng = Rand64::new(1);
103    let CC = Complex64::RING;
104    assert!(poly_ring.terms(f).all(|(c, _): (&Complex64El, _)| CC.re(*c).is_finite() && CC.im(*c).is_finite()));
105    let f_prime = derive_poly(&poly_ring, f);
106    
107    let (approx_root, approx_radius) = (0..MAX_PROBABILISTIC_REPETITIONS).map(|_| {
108        let starting_point_unscaled = CC.add(
109            // this cast might wrap around i64::MAX, so also produces negative values
110            CC.from_f64((rng.rand_u64() as i64) as f64 / i64::MAX as f64),
111            CC.mul(Complex64::I, CC.from_f64((rng.rand_u64() as i64) as f64 / i64::MAX as f64))
112        );
113        let scale = (rng.rand_u64() % (2 * NEWTON_MAX_SCALE as u64)) as i32 - NEWTON_MAX_SCALE as i32;
114        let starting_point = CC.mul(starting_point_unscaled, CC.from_f64(2.0f64.powi(scale)));
115
116        let mut current = starting_point;
117        for _ in 0..NEWTON_ITERATIONS {
118            current = CC.sub(current, CC.div(&poly_ring.evaluate(f, &current, CC.identity()), &poly_ring.evaluate(&f_prime, &current, CC.identity())));
119        }
120
121        // we expect the root to lie within this radius of current
122        let approx_radius = CC.abs(poly_ring.evaluate(f, &current, CC.identity())) / CC.abs(poly_ring.evaluate(&f_prime, &current, CC.identity()));
123
124        return (current, approx_radius);
125    }).min_by(|(_, r1), (_, r2)| match (r1.is_finite(), r2.is_finite()) {
126        (false, false) => Ordering::Equal,
127        (true, false) => Ordering::Less,
128        (false, true) => Ordering::Greater,
129        (true, true) => f64::total_cmp(r1, r2)
130    }).unwrap();
131    if !approx_radius.is_finite() || approx_radius < 0. {
132        return Err(PrecisionError);
133    }
134
135    return Ok((approx_root, bound_distance_to_root(approx_root, poly_ring, f, poly_deg)?));
136}
137
138///
139/// Finds an approximation to a complex root of the given integer polynomial.
140/// 
141/// This function does not try to be as efficient as possible, but instead tries
142/// to avoid (or at least detect) as many numerical problems as possible.
143/// 
144/// The first return value is an approximation to the root of the polynomial, and the
145/// second return value is an upper bound to the distance to the exact root.
146/// 
147#[stability::unstable(feature = "enable")]
148pub fn find_approximate_complex_root<P>(ZZX: P, f: &El<P>) -> Result<(El<Complex64>, f64), PrecisionError>
149    where P: RingStore,
150        P::Type: PolyRing + DivisibilityRing,
151        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: IntegerRing
152{
153    assert!(ZZX.degree(f).unwrap_or(0) > 0);
154    let CC = Complex64::RING;
155    assert!(ZZX.checked_div(&poly_squarefree_part_local(&ZZX, ZZX.clone_el(f), DontObserve), f).is_some(), "polynomial must be square-free");
156    let CCX = DensePolyRing::new(CC, "X");
157    return find_approximate_complex_root_squarefree(&CCX, &CCX.lifted_hom(&ZZX, CC.can_hom(ZZX.base_ring()).unwrap()).map_ref(f), ZZX.degree(f).unwrap());
158}
159
160///
161/// Finds an approximation to all complex roots of the given integer polynomial.
162/// 
163/// This function does not try to be as efficient as possible, but instead tries
164/// to avoid (or at least detect) as many numerical problems as possible.
165/// However, the task of finding all roots has quite bad numerical stability, especially
166/// if some of the roots are close together. Hence, this is likely to fail for polynomials
167/// with large degrees (say > 100) or very large coefficients.
168/// 
169/// The first component of each returned tuple is an approximation to a root of the
170/// polynomial, and the second component is an upper bound to the distance to exact root.
171/// 
172#[stability::unstable(feature = "enable")]
173pub fn find_all_approximate_complex_roots<P>(ZZX: P, poly: &El<P>) -> Result<Vec<(El<Complex64>, f64)>, PrecisionError>
174    where P: RingStore,
175        P::Type: PolyRing + DivisibilityRing,
176        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: IntegerRing
177{
178    assert!(ZZX.degree(poly).unwrap_or(0) > 0);
179    let CC = Complex64::RING;
180    let ZZ_to_CC = CC.can_hom(ZZX.base_ring()).unwrap();
181    assert!(ZZX.checked_div(&poly_squarefree_part_local(&ZZX, ZZX.clone_el(poly), DontObserve), poly).is_some(), "polynomial must be square-free");
182    let CCX = DensePolyRing::new(CC, "X");
183    let ZZX_to_CCX = CCX.lifted_hom(&ZZX, ZZ_to_CC);
184
185    let d = ZZX.degree(poly).unwrap();
186    let f = ZZX_to_CCX.map_ref(poly);
187    let mut remaining_poly = CCX.clone_el(&f);
188    let mut result = Vec::new();
189    for i in 0..ZZX.degree(&poly).unwrap() {
190        let (next_root_initial, _) = find_approximate_complex_root_squarefree(&CCX, &remaining_poly, d - i)?;
191        let (next_root, distance) = newton_with_initial(&CCX, &f, d, next_root_initial)?;
192        if result.iter().any(|(prev_root, prev_distance)| CC.abs(CC.sub(*prev_root, next_root)) <= distance + prev_distance) {
193            return Err(PrecisionError);
194        }
195        result.push((next_root, distance));
196        let mut new_remaining_poly = CCX.zero();
197        for j in (1..=(d - i)).rev() {
198            let lc = *CCX.coefficient_at(&remaining_poly, j);
199            CCX.get_ring().add_assign_from_terms(&mut new_remaining_poly, [(lc, j - 1)]);
200            CCX.get_ring().add_assign_from_terms(&mut remaining_poly, [(CC.mul(lc, next_root), j - 1)]);
201        }
202        remaining_poly = new_remaining_poly;
203    }
204    // just some canonical order to make tests and debugging easier
205    result.sort_unstable_by(|(l, _), (r, _)| f64::total_cmp(&(CC.re(*l) + CC.im(*l) * 0.000001), &(CC.re(*r) + CC.im(*r) * 0.000001)));
206    return Ok(result);
207}
208
209#[cfg(test)]
210use crate::algorithms::cyclotomic::cyclotomic_polynomial;
211#[cfg(test)]
212use std::f64::consts::PI;
213#[cfg(test)]
214use crate::algorithms::eea::signed_gcd;
215#[cfg(test)]
216use crate::primitive_int::StaticRing;
217
218#[test]
219fn test_find_approximate_complex_root() {
220    let ZZ = BigIntRing::RING;
221    let ZZX = DensePolyRing::new(&ZZ, "X");
222    let CC = Complex64::RING;
223
224    let [f] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(2) + 1]);
225    let (root, radius) = find_approximate_complex_root(&ZZX, &f).unwrap();
226    assert!(radius <= 0.000000001);
227    assert!(
228        CC.abs(CC.sub(Complex64::I, root)) <= radius ||
229        CC.abs(CC.add(Complex64::I, root)) <= radius
230    );
231
232    let [f] = ZZX.with_wrapped_indeterminate(|X| [100000000 * X.pow_ref(2) + 1]);
233    let (root, radius) = find_approximate_complex_root(&ZZX, &f).unwrap();
234    assert!(radius <= 0.000000001);
235    assert!(
236        CC.abs(CC.sub(CC.mul(Complex64::I, CC.from_f64(0.0001)), root)) <= radius ||
237        CC.abs(CC.add(CC.mul(Complex64::I, CC.from_f64(0.0001)), root)) <= radius
238    );
239    
240    let f = cyclotomic_polynomial(&ZZX, 105);
241    let (root, radius) = find_approximate_complex_root(&ZZX, &f).unwrap();
242    assert!(radius <= 0.000000001);
243    let root_of_unity = |k, n| CC.exp(CC.mul(CC.from_f64(2.0 * PI * k as f64 / n as f64), Complex64::I));
244    assert!((0..105).filter(|k| signed_gcd(*k, 105, StaticRing::<i64>::RING) == 1).any(|k| CC.abs(CC.sub(root_of_unity(k, 105), root)) <= radius));
245    
246    let [f] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(4) - 6 * X.pow_ref(2) - 11]);
247    let (root, radius) = find_approximate_complex_root(&ZZX, &f).unwrap();
248    assert!(radius <= 0.000000001);
249    let expected = [
250        CC.from_f64(-2.7335207983477241859817441249894389670263693372196541057936155022),
251        CC.from_f64(2.7335207983477241859817441249894389670263693372196541057936155022),
252        CC.mul(CC.from_f64(-1.2133160985495821701841076107103126272177604441592580488679327999), Complex64::I),
253        CC.mul(CC.from_f64(1.2133160985495821701841076107103126272177604441592580488679327999), Complex64::I)
254    ];
255    assert!(expected.iter().any(|x| CC.abs(CC.sub(*x, root)) <= radius));
256}
257
258#[test]
259fn test_find_all_approximate_complex_roots() {
260    let ZZ = BigIntRing::RING;
261    let ZZX = DensePolyRing::new(&ZZ, "X");
262    let CC = Complex64::RING;
263
264    let [f] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(3) - 7 * X + 100]);
265    let expected = [
266        CC.from_f64(-5.1425347350362689414102462079341122827700721701123021518405977272),
267        CC.add(
268            CC.from_f64(2.5712673675181344707051231039670561413850360850561510759202988636),
269            CC.mul(CC.from_f64(-3.5824918179656616775077885147170618458138064112538190024361795659), Complex64::I)
270        ),
271        CC.add(
272            CC.from_f64(2.5712673675181344707051231039670561413850360850561510759202988636),
273            CC.mul(CC.from_f64(3.5824918179656616775077885147170618458138064112538190024361795659), Complex64::I)
274        )
275    ];
276    let actual = find_all_approximate_complex_roots(&ZZX, &f).unwrap();
277    for (expected, (actual, dist)) in expected.iter().copied().zip(actual.iter().copied()) {
278        assert!(dist < 0.000000001);
279        assert!(CC.abs(CC.sub(actual, expected)) <= dist);
280    }
281    
282    let root_of_unity = |k, n| CC.exp(CC.mul(CC.from_f64(2.0 * PI * k as f64 / n as f64), Complex64::I));
283    let f = cyclotomic_polynomial(&ZZX, 105);
284    let expected = (1..=52).rev().filter(|i| signed_gcd(*i, 105, StaticRing::<i64>::RING) == 1).flat_map(|i| [root_of_unity(105 - i, 105), root_of_unity(i, 105)]).chain([CC.one()]);
285    let actual = find_all_approximate_complex_roots(&ZZX, &f).unwrap();
286    assert_eq!(48, actual.len());
287    for (expected, (actual, dist)) in expected.zip(actual.iter().copied()) {
288        assert!(dist < 0.000000001);
289        assert!(CC.abs(CC.sub(actual, expected)) <= dist);
290    }
291}
feanor_math/algorithms/newton.rs

feanor_math/algorithms/
newton.rs
