polycool/
quadratic.rs

1// Copyright 2025 the Kurbo Authors
2// SPDX-License-Identifier: Apache-2.0 OR MIT
3
4use arrayvec::ArrayVec;
5
6use crate::Quadratic;
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
12impl Quadratic {
13    /// This is like [`Quadratic::eval`] but faster.
14    ///
15    /// It would be nice if we could just make `eval` like this, but I couldn't
16    /// figure out how, given the lack of specialization.
17    #[doc(hidden)]
18    pub fn eval_opt(&self, x: f64) -> f64 {
19        let [c0, c1, c2] = self.coeffs;
20        c0 + c1 * x + c2 * x * x
21    }
22
23    /// Returns the roots of this quadratic, in increasing order.
24    ///
25    /// Double-roots are only counted once.
26    pub fn roots(&self) -> ArrayVec<f64, 2> {
27        let &[c, b, a] = self.coeffs();
28        let disc = b * b - 4.0 * a * c;
29        if disc.is_finite() {
30            let mut ret = ArrayVec::new();
31            let mut push = |r: f64| {
32                if r.is_finite() {
33                    ret.push(r)
34                }
35            };
36            if disc > 0.0 {
37                let q = -0.5 * (b + disc.sqrt().copysign(b));
38                let r0 = q / a;
39                let r1 = c / q;
40                push(r0.min(r1));
41                push(r0.max(r1));
42            } else if disc == 0.0 {
43                let root = -0.5 * b / a;
44                if root.is_finite() {
45                    push(root);
46                } else if c == 0.0 {
47                    // This is kurbo's behavior: the intention is that if the
48                    // whole thing is zero, return zero as a single root. I'm
49                    // not sure I love it.
50                    //
51                    // Bear in mind that this branch is not *only* for the
52                    // identically zero case: if a == c == 0.0 and b * b
53                    // underflows then we will end up here. In that case,
54                    // zero is the only root.
55                    push(0.0);
56                }
57            } else {
58                // No roots.
59            }
60            ret
61        } else {
62            // At least one of the coefficients was too large and triggered
63            // overflow.
64            //
65            // The exponent of f64 maxes out at 1023, so scaling down by
66            // 2^{-512} is enough to ensure that squaring doesn't overflow. We
67            // do an extra factor of 2^{-3} for some wiggle room. This can't
68            // completely destroy all the coefficients: because of the overflow,
69            // we know that at least one of them was big.
70            let scale = 2.0f64.powi(-515);
71            // If we're infinite, just give up. (Otherwise, we'd stack overflow
72            // by repeatedly trying to rescale.)
73            if self.is_finite() {
74                (*self * scale).roots()
75            } else {
76                ArrayVec::new()
77            }
78        }
79    }
80
81    /// Returns the two distinct roots of this quadratic, but only if the
82    /// discriminant is positive.
83    pub fn positive_discriminant_roots(&self) -> Option<(f64, f64)> {
84        let &[c, b, a] = self.coeffs();
85        let disc = b * b - 4.0 * a * c;
86        if disc.is_finite() {
87            if disc > 0.0 {
88                let q = -0.5 * (b + disc.sqrt().copysign(b));
89                let r0 = q / a;
90                let r1 = c / q;
91                Some((r0.min(r1), r0.max(r1)))
92            } else {
93                None
94            }
95        } else {
96            self.positive_discriminant_roots_scaled()
97        }
98    }
99
100    #[cold]
101    fn positive_discriminant_roots_scaled(&self) -> Option<(f64, f64)> {
102        if self.is_finite() {
103            let scale = 2.0f64.powi(-515);
104            (*self * scale).positive_discriminant_roots()
105        } else {
106            None
107        }
108    }
109}
110
111#[cfg(test)]
112mod tests {
113    #[test]
114    fn root_evaluation() {
115        arbtest::arbtest(|u| {
116            let q = crate::arbitrary::quadratic(u)?;
117            // Arbitrary quadratics can have coefficients with wild magnitudes,
118            // so we need to adjust our error expectations accordingly.
119            let magnitude = q.max_abs_coefficient().max(1.0);
120
121            for r in q.roots() {
122                let y = q.eval(r);
123                // To evaluate the polynomial, we need to square r, so our error
124                // should be relative to the magnitude of r squared.
125                let r_magnitude = r.abs().max(1.0);
126                let threshold = r_magnitude * 1e-14 * r_magnitude * magnitude;
127                if y.is_finite() && threshold.is_finite() {
128                    assert!(y.abs() <= threshold);
129                }
130            }
131            Ok(())
132        })
133        .budget_ms(5_000);
134    }
135}
polycool/quadratic.rs

polycool/
quadratic.rs
