polycool/
yuksel.rs

1// Copyright 2025 the Kurbo Authors
2// SPDX-License-Identifier: Apache-2.0 OR MIT
3
4//! An implementation of Yuksel's robust version of Newton's algorithm.
5
6fn different_signs(x: f64, y: f64) -> bool {
7    (x < 0.0) != (y < 0.0)
8}
9
10pub(crate) fn find_root<F: Fn(f64) -> f64, DF: Fn(f64) -> f64>(
11    f: F,
12    deriv: DF,
13    mut lower: f64,
14    mut upper: f64,
15    val_lower: f64,
16    val_upper: f64,
17    x_error: f64,
18) -> f64 {
19    if !val_lower.is_finite() || !val_upper.is_finite() {
20        return f64::NAN;
21    }
22    debug_assert!(different_signs(val_lower, val_upper));
23
24    let mut x = lower + (upper - lower) / 2.0;
25    let mut step = (upper - lower) / 2.0;
26
27    if step.abs() <= x_error {
28        return x;
29    }
30
31    while step.abs() > x_error && x.is_finite() {
32        // There's potentially a performance gain to be had by evaluating the
33        // derivative and the polynomial "jointly", so that subexpressions can
34        // be shared. At least for the cubic case, the compiler is smart enough
35        // to do that for us.
36        let deriv_x = deriv(x);
37        let val_x = f(x);
38
39        // By returning early on zero, we make it impossible for `step` to
40        // be NaN. This seems to benchmark a little better than checking for
41        // NaN after the fact.
42        if val_x == 0.0 {
43            return x;
44        }
45        let root_in_first_half = different_signs(val_lower, val_x);
46        if root_in_first_half {
47            upper = x;
48        } else {
49            lower = x;
50        }
51
52        debug_assert!(deriv_x.is_finite());
53        debug_assert!(val_x.is_finite());
54
55        step = -val_x / deriv_x;
56        let mut new_x = x + step;
57
58        if new_x <= lower || new_x >= upper {
59            new_x = lower + (upper - lower) / 2.0;
60
61            if new_x == upper || new_x == lower {
62                // This should be rare, but it happens if they ask for more
63                // accuracy than is reasonable. For example, suppose (because
64                // of large coefficients) the output value jumps from -1.0
65                // to 1.0 between adjacent floats and they ask for an output
66                // error of smaller than 0.5. Then we'll eventually shrink
67                // the search interval to a pair of adjacent floats and hit
68                // this case.
69                return new_x;
70            }
71        }
72        step = new_x - x;
73        x = new_x;
74    }
75    x
76}
polycool/yuksel.rs

polycool/
yuksel.rs
