polycool/cubic.rs
1// Copyright 2025 the Kurbo Authors
2// SPDX-License-Identifier: Apache-2.0 OR MIT
3
4use arrayvec::ArrayVec;
5
6use crate::{Cubic, Quadratic, different_signs};
7
8#[cfg(feature = "libm")]
9#[allow(unused_imports, reason = "unused if libm and std are both around")]
10use crate::libm_polyfill::FloatFuncs as _;
11
12// We assume that there is at least one element, at most 3 elements, and buf[1]
13// and buf[2] are already in the right order (if present).
14//
15// This might be faster than the stdlib sort for our special case,
16// but the real reason it's here is for no_std support.
17fn partial_sort<const M: usize>(buf: &mut ArrayVec<f64, M>) {
18 if buf.len() > 1 && buf[0] > buf[1] {
19 buf.swap(0, 1);
20 if buf.len() > 2 && buf[1] > buf[2] {
21 buf.swap(1, 2);
22 }
23 }
24}
25
26impl Cubic {
27 /// This is like [`Cubic::eval`] but faster.
28 ///
29 /// It would be nice if we could just make `eval` like this, but I couldn't
30 /// figure out how, given the lack of specialization.
31 ///
32 /// Marked `pub` for benchmarking.
33 #[doc(hidden)]
34 pub fn eval_opt(&self, x: f64) -> f64 {
35 let [c0, c1, c2, c3] = self.coeffs;
36 let xx = x * x;
37 let xxx = xx * x;
38 c0 + c1 * x + c2 * xx + c3 * xxx
39 }
40
41 /// Evaluate this cubic and its gradient at the same time, reusing some
42 /// of the intermediate computations.
43 ///
44 /// In micro-benchmarks, this is faster than evaluating separately. But
45 /// it doesn't help with the performance of Yuksel's algorithm, presumably
46 /// because the compiler is inlining and optimizing reuse of the
47 /// intermediate computations already.
48 ///
49 /// Marked `pub` for benchmarking.
50 #[doc(hidden)]
51 pub fn eval_with_deriv_opt(&self, deriv: &Quadratic, x: f64) -> (f64, f64) {
52 let [c0, c1, c2, c3] = self.coeffs;
53 let [d0, d1, d2] = deriv.coeffs;
54 let xx = x * x;
55 let xxx = xx * x;
56 (c0 + c1 * x + c2 * xx + c3 * xxx, d0 + d1 * x + d2 * xx)
57 }
58
59 /// Computes the critical points of this cubic, as long
60 /// as the discriminant of the derivative is positive.
61 /// The return values are in increasing order.
62 ///
63 /// Some corner cases worth noting:
64 /// - If the discriminant is zero, returns nothing. That is, we don't find
65 /// double-roots of the derivative. These aren't needed anyway, as these
66 /// critical points are used to construct bracketing intervals for the
67 /// roots. Double-roots of the derivative are inflection points of the
68 /// cubic, and they aren't needed to bracket roots.
69 /// - If the derivative is linear or close to it, we might
70 /// return +/- infinity as one of the roots.
71 /// - Unless some input is NaN, we don't return NaN.
72 fn critical_points(&self) -> Option<(f64, f64)> {
73 let a = 3.0 * self.coeffs[3];
74 let b_2 = self.coeffs[2];
75 let c = self.coeffs[1];
76 let disc_4 = b_2 * b_2 - a * c;
77
78 if !disc_4.is_finite() {
79 return self.rescaled_critical_points();
80 }
81
82 if disc_4 > 0.0 {
83 let q = -(b_2 + disc_4.sqrt().copysign(b_2));
84 let r0 = q / a;
85 let r1 = c / q;
86 Some((r0.min(r1), r0.max(r1)))
87 } else {
88 None
89 }
90 }
91
92 #[cold]
93 fn rescaled_critical_points(&self) -> Option<(f64, f64)> {
94 let scale = 2.0f64.powi(-515);
95 (*self * scale).critical_points()
96 }
97
98 fn one_root(
99 &self,
100 lower: f64,
101 upper: f64,
102 lower_val: f64,
103 upper_val: f64,
104 x_error: f64,
105 ) -> f64 {
106 let deriv = self.deriv();
107 if !deriv.is_finite() {
108 return f64::NAN;
109 }
110 crate::yuksel::find_root(
111 |x| self.eval_opt(x),
112 |x| deriv.eval_opt(x),
113 lower,
114 upper,
115 lower_val,
116 upper_val,
117 x_error,
118 )
119 }
120
121 fn first_root(self, lower: f64, upper: f64, x_error: f64) -> Option<f64> {
122 if let Some((x0, x1)) = self.critical_points() {
123 let possible_endpoints: [f64; 3] = [x0, x1, upper];
124 let mut last = lower;
125 let mut last_val = self.eval(last);
126 for x in possible_endpoints {
127 if x > last && x <= upper {
128 let val = self.eval(x);
129 if different_signs(last_val, val) {
130 return Some(self.one_root(last, x, last_val, val, x_error));
131 }
132
133 last = x;
134 last_val = val;
135 }
136 }
137 None
138 } else {
139 let lower_val = self.eval(lower);
140 let upper_val = self.eval(upper);
141 if different_signs(lower_val, upper_val) {
142 Some(self.one_root(lower, upper, lower_val, upper_val, x_error))
143 } else {
144 None
145 }
146 }
147 }
148
149 /// Computes all roots between `lower` and `upper`, to the desired accuracy.
150 ///
151 /// We make no guarantees about multiplicity. In fact, if there's a
152 /// double-root that isn't a triple-root (and therefore has no sign change
153 /// nearby) then there's a good chance we miss it altogether. This is
154 /// fine if you're using this root-finding to optimize a quartic, because
155 /// double-roots of the derivative aren't local extrema.
156 pub fn roots_between(self, lower: f64, upper: f64, x_error: f64) -> ArrayVec<f64, 3> {
157 let mut ret = ArrayVec::new();
158 let mut scratch = ArrayVec::new();
159 self.roots_between_with_buffer(lower, upper, x_error, &mut ret, &mut scratch);
160 ret
161 }
162
163 pub(crate) fn roots_between_with_buffer<const M: usize>(
164 self,
165 lower: f64,
166 upper: f64,
167 x_error: f64,
168 out: &mut ArrayVec<f64, M>,
169 _scratch: &mut ArrayVec<f64, M>,
170 ) {
171 if let Some(r) = self.first_root(lower, upper, x_error) {
172 out.push(r);
173 let quad = self.deflate(r);
174 if let Some((x0, x1)) = quad.positive_discriminant_roots() {
175 if lower <= x0 && x0 <= upper {
176 out.push(x0);
177 }
178 if lower <= x1 && x1 <= upper {
179 out.push(x1);
180 }
181
182 // `self.first_root` is supposed to return the smallest root in
183 // our interval, but it's possible it doesn't because it misses
184 // a double-root (or near-double-root).
185 if lower <= x0 && x0 < r {
186 partial_sort(out);
187 }
188 }
189 }
190 }
191}
192
193#[cfg(test)]
194mod tests {
195 use super::*;
196
197 // Embarrassingly, an early version gave NaN on x^3, because
198 // the Newton iterations ran into 0.0/0.0.
199 #[test]
200 fn x_cubed() {
201 let p = Cubic::new([0.0, 0.0, 0.0, 1.0]);
202 let roots = p.roots_between(-1.0, 1.0, 1e-6);
203 assert_eq!(roots.len(), 1);
204 assert!(roots[0].abs() <= 1e-6);
205 }
206}