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
//! An implementation of Yuksel's robust version of Newton's algorithm.
fn different_signs(x: f64, y: f64) -> bool {
(x < 0.0) != (y < 0.0)
}
pub(crate) fn find_root<F: Fn(f64) -> f64, DF: Fn(f64) -> f64>(
f: F,
deriv: DF,
mut lower: f64,
mut upper: f64,
val_lower: f64,
val_upper: f64,
x_error: f64,
) -> f64 {
if !val_lower.is_finite() || !val_upper.is_finite() {
return f64::NAN;
}
debug_assert!(different_signs(val_lower, val_upper));
let mut x = lower + (upper - lower) / 2.0;
let mut step = (upper - lower) / 2.0;
if step.abs() <= x_error {
return x;
}
while step.abs() > x_error && x.is_finite() {
// There's potentially a performance gain to be had by evaluating the
// derivative and the polynomial "jointly", so that subexpressions can
// be shared. At least for the cubic case, the compiler is smart enough
// to do that for us.
let deriv_x = deriv(x);
let val_x = f(x);
let root_in_first_half = different_signs(val_lower, val_x);
if root_in_first_half {
upper = x;
} else {
lower = x;
}
debug_assert!(deriv_x.is_finite());
debug_assert!(val_x.is_finite());
step = -val_x / deriv_x;
let mut new_x = x + step;
if new_x <= lower || new_x >= upper {
new_x = lower + (upper - lower) / 2.0;
if new_x == upper || new_x == lower {
// This should be rare, but it happens if they ask for more
// accuracy than is reasonable. For example, suppse (because
// of large coefficients) the output value jumps from -1.0
// to 1.0 between adjacent floats and they ask for an output
// error of smaller than 0.5. Then we'll eventually shrink
// the search interval to a pair of adjacent floats and hit
// this case.
return new_x;
}
}
step = new_x - x;
x = new_x;
}
x
}