zenith_float_num/
roots.rs

1//! Scalar real root finders on [`ExactNum`].
2
3use crate::defs::WORD_BIT_SIZE;
4use crate::Consts;
5use crate::ExactNum;
6use crate::RoundingMode;
7
8/// Maximum iterations for every method in this module.
9pub const ROOT_MAX_ITER: usize = 256;
10
11/// Default absolute tolerance exponent: `tol = 2^{ROOT_DEFAULT_TOL}`.
12pub const ROOT_DEFAULT_TOL: i32 = -256;
13
14fn work_p(p: usize) -> usize {
15    p.saturating_add(WORD_BIT_SIZE)
16}
17
18fn finite(x: &ExactNum) -> bool {
19    !x.is_nan() && !x.is_inf()
20}
21
22fn opposite_signs(a: &ExactNum, b: &ExactNum) -> bool {
23    if a.is_zero() || b.is_zero() || !finite(a) || !finite(b) {
24        return false;
25    }
26    a.is_positive() != b.is_positive()
27}
28
29fn below_tol(x: &ExactNum, tol: &ExactNum) -> bool {
30    x.is_zero() || x.cmp(tol) == Some(-1) || x.cmp(tol) == Some(0)
31}
32
33/// Suggested absolute tolerance `2^{ROOT_DEFAULT_TOL}` at precision `p`.
34pub fn root_default_tol(p: usize, rm: RoundingMode) -> ExactNum {
35    ExactNum::from_u8(1, p).ldexp(ROOT_DEFAULT_TOL, p, rm)
36}
37
38fn bisect_counted<F>(
39    mut f: F,
40    mut a: ExactNum,
41    mut b: ExactNum,
42    tol: &ExactNum,
43    p: usize,
44    rm: RoundingMode,
45    cc: &mut Consts,
46) -> Option<(ExactNum, usize)>
47where
48    F: FnMut(&ExactNum, usize, RoundingMode, &mut Consts) -> ExactNum,
49{
50    if !finite(&a) || !finite(&b) || !finite(tol) || a.cmp(&b) != Some(-1) {
51        return None;
52    }
53    let wrk = work_p(p);
54    let mut fa = f(&a, wrk, RoundingMode::None, cc);
55    let mut fb = f(&b, wrk, RoundingMode::None, cc);
56    if !opposite_signs(&fa, &fb) {
57        return None;
58    }
59    let two = ExactNum::from_u8(2, wrk);
60    for k in 1..=ROOT_MAX_ITER {
61        let mid = a
62            .add(&b, wrk, RoundingMode::None)
63            .div(&two, wrk, RoundingMode::None);
64        let fm = f(&mid, wrk, RoundingMode::None, cc);
65        if !finite(&fm) {
66            return None;
67        }
68        if opposite_signs(&fa, &fm) {
69            b = mid.clone();
70            fb = fm.clone();
71        } else {
72            a = mid.clone();
73            fa = fm.clone();
74        }
75        if fm.is_zero() || below_tol(&b.sub(&a, wrk, RoundingMode::None).abs(), tol) {
76            let mut out = mid;
77            let _ = out.set_precision(p, rm);
78            return Some((out, k));
79        }
80        let _ = fb;
81    }
82    None
83}
84
85/// Bisection on `[a, b]`. `None` if `f(a)` and `f(b)` do not have opposite signs.
86pub fn bisect<F>(
87    f: F,
88    a: &ExactNum,
89    b: &ExactNum,
90    tol: &ExactNum,
91    p: usize,
92    rm: RoundingMode,
93    cc: &mut Consts,
94) -> Option<ExactNum>
95where
96    F: FnMut(&ExactNum, usize, RoundingMode, &mut Consts) -> ExactNum,
97{
98    bisect_counted(f, a.clone(), b.clone(), tol, p, rm, cc).map(|(x, _)| x)
99}
100
101fn newton_counted<F, D>(
102    mut f: F,
103    mut df: D,
104    x0: ExactNum,
105    tol: &ExactNum,
106    max_iter: usize,
107    p: usize,
108    rm: RoundingMode,
109    cc: &mut Consts,
110) -> Option<(ExactNum, usize)>
111where
112    F: FnMut(&ExactNum, usize, RoundingMode, &mut Consts) -> ExactNum,
113    D: FnMut(&ExactNum, usize, RoundingMode, &mut Consts) -> ExactNum,
114{
115    if !finite(&x0) || !finite(tol) || max_iter == 0 || max_iter > ROOT_MAX_ITER {
116        return None;
117    }
118    let wrk = work_p(p);
119    let mut x = x0;
120    for k in 1..=max_iter {
121        let y = f(&x, wrk, RoundingMode::None, cc);
122        let d = df(&x, wrk, RoundingMode::None, cc);
123        if !finite(&y) || !finite(&d) || d.is_zero() {
124            return None;
125        }
126        let step = y.div(&d, wrk, RoundingMode::None);
127        x = x.sub(&step, wrk, RoundingMode::None);
128        if !finite(&x) {
129            return None;
130        }
131        if y.is_zero() || below_tol(&step.abs(), tol) {
132            let _ = x.set_precision(p, rm);
133            return Some((x, k));
134        }
135    }
136    None
137}
138
139/// Newton–Raphson from `x0` with exact derivative `df`.
140pub fn newton<F, D>(
141    f: F,
142    df: D,
143    x0: &ExactNum,
144    tol: &ExactNum,
145    max_iter: usize,
146    p: usize,
147    rm: RoundingMode,
148    cc: &mut Consts,
149) -> Option<ExactNum>
150where
151    F: FnMut(&ExactNum, usize, RoundingMode, &mut Consts) -> ExactNum,
152    D: FnMut(&ExactNum, usize, RoundingMode, &mut Consts) -> ExactNum,
153{
154    newton_counted(f, df, x0.clone(), tol, max_iter, p, rm, cc).map(|(x, _)| x)
155}
156
157fn illinois_counted<F>(
158    mut f: F,
159    mut a: ExactNum,
160    mut b: ExactNum,
161    tol: &ExactNum,
162    p: usize,
163    rm: RoundingMode,
164    cc: &mut Consts,
165) -> Option<(ExactNum, usize)>
166where
167    F: FnMut(&ExactNum, usize, RoundingMode, &mut Consts) -> ExactNum,
168{
169    if !finite(&a) || !finite(&b) || !finite(tol) || a.cmp(&b) != Some(-1) {
170        return None;
171    }
172    let wrk = work_p(p);
173    let two = ExactNum::from_u8(2, wrk);
174    let mut fa = f(&a, wrk, RoundingMode::None, cc);
175    let mut fb = f(&b, wrk, RoundingMode::None, cc);
176    if !opposite_signs(&fa, &fb) {
177        return None;
178    }
179    for k in 1..=ROOT_MAX_ITER {
180        let den = fb.sub(&fa, wrk, RoundingMode::None);
181        if den.is_zero() {
182            return None;
183        }
184        let c = a
185            .mul(&fb, wrk, RoundingMode::None)
186            .sub(
187                &b.mul(&fa, wrk, RoundingMode::None),
188                wrk,
189                RoundingMode::None,
190            )
191            .div(&den, wrk, RoundingMode::None);
192        let fc = f(&c, wrk, RoundingMode::None, cc);
193        if !finite(&c) || !finite(&fc) {
194            return None;
195        }
196        if fc.is_zero() || below_tol(&b.sub(&a, wrk, RoundingMode::None).abs(), tol) {
197            let mut out = c;
198            let _ = out.set_precision(p, rm);
199            return Some((out, k));
200        }
201        if opposite_signs(&fa, &fc) {
202            b = c;
203            fb = fc;
204            fa = fa.div(&two, wrk, RoundingMode::None);
205        } else {
206            a = c;
207            fa = fc;
208            fb = fb.div(&two, wrk, RoundingMode::None);
209        }
210    }
211    None
212}
213
214/// Illinois (modified regula falsi) on `[a, b]`.
215pub fn illinois<F>(
216    f: F,
217    a: &ExactNum,
218    b: &ExactNum,
219    tol: &ExactNum,
220    p: usize,
221    rm: RoundingMode,
222    cc: &mut Consts,
223) -> Option<ExactNum>
224where
225    F: FnMut(&ExactNum, usize, RoundingMode, &mut Consts) -> ExactNum,
226{
227    illinois_counted(f, a.clone(), b.clone(), tol, p, rm, cc).map(|(x, _)| x)
228}
229
230fn strictly_inside(x: &ExactNum, lo: &ExactNum, hi: &ExactNum) -> bool {
231    x.cmp(lo) == Some(1) && x.cmp(hi) == Some(-1)
232}
233
234fn brent_counted<F>(
235    mut f: F,
236    mut a: ExactNum,
237    mut b: ExactNum,
238    tol: &ExactNum,
239    p: usize,
240    rm: RoundingMode,
241    cc: &mut Consts,
242) -> Option<(ExactNum, usize)>
243where
244    F: FnMut(&ExactNum, usize, RoundingMode, &mut Consts) -> ExactNum,
245{
246    if !finite(&a) || !finite(&b) || !finite(tol) || a.cmp(&b) != Some(-1) {
247        return None;
248    }
249    let wrk = work_p(p);
250    let two = ExactNum::from_u8(2, wrk);
251    let mut fa = f(&a, wrk, RoundingMode::None, cc);
252    let mut fb = f(&b, wrk, RoundingMode::None, cc);
253    if !opposite_signs(&fa, &fb) {
254        return None;
255    }
256    let mut c = a.clone();
257    let mut fc = fa.clone();
258    for k in 1..=ROOT_MAX_ITER {
259        if below_tol(&b.sub(&a, wrk, RoundingMode::None).abs(), tol) || fb.is_zero() {
260            let mut out = b;
261            let _ = out.set_precision(p, rm);
262            return Some((out, k));
263        }
264        let den_ba = fb.sub(&fa, wrk, RoundingMode::None);
265        let mut s = if !den_ba.is_zero() {
266            a.mul(&fb, wrk, RoundingMode::None)
267                .sub(
268                    &b.mul(&fa, wrk, RoundingMode::None),
269                    wrk,
270                    RoundingMode::None,
271                )
272                .div(&den_ba, wrk, RoundingMode::None)
273        } else {
274            a.add(&b, wrk, RoundingMode::None)
275                .div(&two, wrk, RoundingMode::None)
276        };
277        if a.cmp(&c) != Some(0) && b.cmp(&c) != Some(0) {
278            let d1 = fa.sub(&fb, wrk, RoundingMode::None);
279            let d2 = fa.sub(&fc, wrk, RoundingMode::None);
280            let d3 = fb.sub(&fc, wrk, RoundingMode::None);
281            if !d1.is_zero() && !d2.is_zero() && !d3.is_zero() {
282                let t1 = a
283                    .mul(&fb, wrk, RoundingMode::None)
284                    .mul(&fc, wrk, RoundingMode::None)
285                    .div(
286                        &d1.mul(&d2, wrk, RoundingMode::None),
287                        wrk,
288                        RoundingMode::None,
289                    );
290                let t2 = b
291                    .mul(&fa, wrk, RoundingMode::None)
292                    .mul(&fc, wrk, RoundingMode::None)
293                    .div(
294                        &fb.sub(&fa, wrk, RoundingMode::None)
295                            .mul(&d3, wrk, RoundingMode::None),
296                        wrk,
297                        RoundingMode::None,
298                    );
299                let t3 = c
300                    .mul(&fa, wrk, RoundingMode::None)
301                    .mul(&fb, wrk, RoundingMode::None)
302                    .div(
303                        &fc.sub(&fa, wrk, RoundingMode::None).mul(
304                            &fc.sub(&fb, wrk, RoundingMode::None),
305                            wrk,
306                            RoundingMode::None,
307                        ),
308                        wrk,
309                        RoundingMode::None,
310                    );
311                let iqi = t1
312                    .add(&t2, wrk, RoundingMode::None)
313                    .add(&t3, wrk, RoundingMode::None);
314                if strictly_inside(&iqi, &a, &b) {
315                    s = iqi;
316                }
317            }
318        }
319        if !strictly_inside(&s, &a, &b) {
320            s = a
321                .add(&b, wrk, RoundingMode::None)
322                .div(&two, wrk, RoundingMode::None);
323        }
324        let fs = f(&s, wrk, RoundingMode::None, cc);
325        if !finite(&fs) {
326            return None;
327        }
328        c = b.clone();
329        fc = fb.clone();
330        if opposite_signs(&fa, &fs) {
331            b = s;
332            fb = fs;
333        } else {
334            a = s;
335            fa = fs;
336        }
337        if a.cmp(&b) == Some(1) {
338            core::mem::swap(&mut a, &mut b);
339            core::mem::swap(&mut fa, &mut fb);
340        }
341    }
342    None
343}
344
345/// Brent's method (bisection + secant + inverse quadratic) on `[a, b]`.
346pub fn brent<F>(
347    f: F,
348    a: &ExactNum,
349    b: &ExactNum,
350    tol: &ExactNum,
351    p: usize,
352    rm: RoundingMode,
353    cc: &mut Consts,
354) -> Option<ExactNum>
355where
356    F: FnMut(&ExactNum, usize, RoundingMode, &mut Consts) -> ExactNum,
357{
358    brent_counted(f, a.clone(), b.clone(), tol, p, rm, cc).map(|(x, _)| x)
359}
360
361#[cfg(test)]
362mod tests {
363    use super::*;
364    use crate::Consts;
365
366    fn gold_p() -> (usize, RoundingMode) {
367        (256, RoundingMode::ToEven)
368    }
369
370    #[test]
371    fn roots_bisect_newton_brent_illinois() {
372        let (p, rm) = gold_p();
373        let mut cc = Consts::new().expect("consts");
374        let tol = root_default_tol(p, rm);
375        let three = ExactNum::from_u8(3, p);
376        let four = ExactNum::from_u8(4, p);
377        let zero = ExactNum::new(p);
378        let one = ExactNum::from_u8(1, p);
379        let two = ExactNum::from_u8(2, p);
380        let pi = cc.pi(p, rm);
381
382        let (broot, biters) = bisect_counted(
383            |x, p, rm, cc| x.sin(p, rm, cc),
384            three.clone(),
385            four.clone(),
386            &tol,
387            p,
388            rm,
389            &mut cc,
390        )
391        .expect("bisect sin");
392        assert_eq!(broot.cmp(&pi), Some(0));
393
394        let nroot = newton(
395            |x, p, rm, _cc| x.mul(x, p, rm).sub(&two, p, rm),
396            |x, p, rm, _cc| two.mul(x, p, rm),
397            &one,
398            &tol,
399            ROOT_MAX_ITER,
400            p,
401            rm,
402            &mut cc,
403        )
404        .expect("newton sqrt2");
405        let s2 = two.sqrt(p, rm);
406        assert_eq!(nroot.cmp(&s2), Some(0));
407
408        let (rbrent, briters) = brent_counted(
409            |x, p, rm, cc| x.sin(p, rm, cc),
410            three.clone(),
411            four.clone(),
412            &tol,
413            p,
414            rm,
415            &mut cc,
416        )
417        .expect("brent sin");
418        assert_eq!(rbrent.cmp(&pi), Some(0));
419        assert!(briters < biters);
420
421        let iroot = illinois(
422            |x, p, rm, cc| x.sin(p, rm, cc),
423            &three,
424            &four,
425            &tol,
426            p,
427            rm,
428            &mut cc,
429        )
430        .expect("illinois sin");
431        assert_eq!(iroot.cmp(&pi), Some(0));
432
433        assert!(bisect(
434            |x, p, rm, cc| x.sin(p, rm, cc),
435            &zero,
436            &one,
437            &tol,
438            p,
439            rm,
440            &mut cc
441        )
442        .is_none());
443    }
444}
zenith_float_num/roots.rs

zenith_float_num/
roots.rs
