root1d 0.4.1

One dimensional root finding algorithms
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140

use std::f64::consts::{PI, E};

// Functions use in the "Algorithms 748" paper.

fn f1(x: f64) -> f64 { x.sin() - x / 2. }

fn f2(x: f64) -> f64 {
    let mut s = 0.;
    for i in 1..=20 {
        let i = i as f64;
        s += (2. * i - 5.).powi(2) / (x - i * i).powi(3);
    }
    -2. * s
}

fn f3<const A: i32, const B: i32>(x: f64) -> f64 {
    A as f64 * x * (B as f64 * x).exp()
}

fn f4_0<const N:i32>(x: f64) -> f64 { x.powi(N) - 0.2 }
fn f4_1<const N:i32>(x: f64) -> f64 { x.powi(N) - 1. }

fn f5(x: f64) -> f64 { x.sin() - 0.5 }

fn f6<const N:i32>(x: f64) -> f64 {
    let n = (-N) as f64;
    2. * x * n.exp() - 2. * (n*x).exp() + 1.
}

fn f7<const N:i32>(x: f64) -> f64 {
    (1. + (1. - N as f64).powi(2)) * x - (1. - N as f64 * x).powi(2)
}

fn f8<const N:i32>(x: f64) -> f64 { x * x - (1. - x).powi(N) }

fn f9<const N:i32>(x: f64) -> f64 {
    (1. + (1. - N as f64).powi(4)) * x - (1. - N as f64 * x).powi(4)
}

fn f10<const N:i32>(x: f64) -> f64 {
    ((-N) as f64 * x).exp() * (x - 1.) + x.powi(N)
}

fn f11<const N:i32>(x: f64) -> f64 {
    let n = N as f64;
    (n * x - 1.) / ((n - 1.) * x)
}

fn f12<const N:i32>(x: f64) -> f64 {
    let n1 = 1. / N as f64;
    x.powf(n1) - (N as f64).powf(n1)
}

fn f13(x: f64) -> f64 {
    if x == 0. { 0. }  else { x / (1. / (x * x)).exp() }
}

fn f14<const N:i32>(x: f64) -> f64 {
    let n20 = N as f64 / 20.;
    if x >= 0. { n20 * (x / 1.5 + x.sin() - 1.) } else { -n20 }
}

fn f15<const N:i32>(x: f64) -> f64 {
    if x > 2e-3 / (1 + N) as f64 {
        E - 1.859
    } else if x >= 0. {
        (0.5 * (N + 1) as f64 * x * 1e3).exp() - 1.859
    } else {
        -0.859
    }
}

/// Return the list of test problems in the “Algorithm 748” paper.
#[inline]
pub fn toms748_test_problems() -> [Vec<(fn(f64) -> f64, f64, f64)>; 15] {
    macro_rules! f { ($f: expr) => { $f as fn(f64) -> f64 } }
    macro_rules! vec_gen {
        ($f: ident $($n: literal)+, $a: expr, $b:expr) => {
            vec![$((f!($f::<$n>), $a, $b),)+]
        }
    }
    [vec![(f!(f1), 0.5 * PI, PI)],
     vec![(f!(f2),   1. + 1e-9,   4. - 1e-9), // n = 1
          (f!(f2),   4. + 1e-9,   9. - 1e-9),
          (f!(f2),   9. + 1e-9,  16. - 1e-9),
          (f!(f2),  16. + 1e-9,  25. - 1e-9),
          (f!(f2),  25. + 1e-9,  36. - 1e-9), // n = 5
          (f!(f2),  36. + 1e-9,  49. - 1e-9),
          (f!(f2),  49. + 1e-9,  64. - 1e-9),
          (f!(f2),  64. + 1e-9,  81. - 1e-9),
          (f!(f2),  81. + 1e-9, 100. - 1e-9),
          (f!(f2), 100. + 1e-9, 121. - 1e-9)], // n = 10
     vec![(f!(f3::< -40, -1>), -9., 31.),
          (f!(f3::<-100, -2>), -9., 31.),
          (f!(f3::<-200, -3>), -9., 31.)],
     { let mut t = vec_gen!(f4_0 4 6 8 10 12, 0., 5.);
       t.append(&mut vec_gen!(f4_1 4 6 8 10 12, 0., 5.));
       t.append(&mut vec_gen!(f4_1 8 10 12 14, -0.95, 4.05));
       t },
     vec![(f!(f5), 0., 1.5)],
     vec_gen!(f6 1 2 3 4 5 20 40 60 80 100, 0., 1.),
     vec_gen!(f7 5 10 20, 0., 1.),
     vec_gen!(f8 2 5 10 15 20, 0., 1.),
     vec_gen!(f9 1 2 4 5 8 15 20, 0., 1.),
     vec_gen!(f10 1 5 10 15 20, 0., 1.),
     vec_gen!(f11 2 5 15 20, 0.01, 1.),
     vec_gen!(f12 2 3 4 5 6 7 9 11 13 15 17 19 21 23 25 27 29 31 33,
              1., 100.),
     vec![(f!(f13), -1., 4.)],
     vec_gen!(f14 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
              22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40,
              -1e4, PI / 2.),
     vec_gen!(f15 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
              38 39 40 100 200 300 400 500 600 700 800 900 1000,
              -1e4, 1e-4)
    ]
}



/// Tolerances for the "roots" crate.
pub struct Tol {
    pub rtol: f64,
    pub atol: f64,
    pub maxiter: usize,
}

impl roots::Convergency<f64> for Tol {
    fn is_root_found(&mut self, y: f64) -> bool {
        y == 0.
    }

    fn is_converged(&mut self, a: f64, b: f64) -> bool {
        (a - b).abs() <= self.rtol * a.abs().max(b.abs()) + self.atol
    }

    fn is_iteration_limit_reached(&mut self, iter: usize) -> bool {
        iter > self.maxiter
    }
}