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
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
use arrayvec::ArrayVec;
#[derive(Copy, Clone, Debug, PartialEq)]
pub struct Quadratic {
pub c0: f64,
pub c1: f64,
pub c2: f64,
}
impl std::ops::Mul<f64> for Quadratic {
type Output = Quadratic;
fn mul(self, rhs: f64) -> Self::Output {
Self {
c0: self.c0 * rhs,
c1: self.c1 * rhs,
c2: self.c2 * rhs,
}
}
}
impl std::ops::Div<f64> for Quadratic {
type Output = Quadratic;
fn div(self, rhs: f64) -> Self::Output {
Self {
c0: self.c0 / rhs,
c1: self.c1 / rhs,
c2: self.c2 / rhs,
}
}
}
impl Quadratic {
pub fn eval(&self, x: f64) -> f64 {
self.c0 + self.c1 * x + self.c2 * x * x
}
pub fn is_finite(&self) -> bool {
self.c0.is_finite() && self.c1.is_finite() && self.c2.is_finite()
}
pub fn roots(&self) -> ArrayVec<f64, 2> {
let a = self.c2;
let b = self.c1;
let c = self.c0;
let disc = b * b - 4.0 * a * c;
if disc.is_finite() {
let mut ret = ArrayVec::new();
let mut push = |r: f64| {
if r.is_finite() {
ret.push(r)
}
};
if disc > 0.0 {
let q = -0.5 * (b + disc.sqrt().copysign(b));
let r0 = q / a;
let r1 = c / q;
push(r0.min(r1));
push(r0.max(r1));
} else if disc == 0.0 {
let root = -0.5 * b / a;
if root.is_finite() {
push(root);
} else if c == 0.0 {
// This is kurbo's behavior: the intention is that if the
// whole thing is zero, return zero as a single root. I'm
// not sure I love it.
//
// Bear in mind that this branch is not *only* for the
// identically zero case: if a == c == 0.0 and b * b
// underflows then we will end up here. In that case,
// zero is the only root.
push(0.0);
}
} else {
// No roots.
}
ret
} else {
// At least one of the coefficients was too large and triggered
// overflow.
//
// The exponent of f64 maxes out at 1023, so scaling down by
// 2^{-512} is enough to ensure that squaring doesn't overflow. We
// do an extra factor of 2^{-3} for some wiggle room. This can't
// completely destroy all the coefficients: because of the overflow,
// we know that at least one of them was big.
let scale = 2.0f64.powi(-515);
// TODO: this can stack overflow if we're infinite. How should
// we handle that?
(*self * scale).roots()
}
}
pub fn positive_discriminant_roots(&self) -> Option<(f64, f64)> {
let a = self.c2;
let b = self.c1;
let c = self.c0;
let disc = b * b - 4.0 * a * c;
if disc.is_finite() {
if disc > 0.0 {
let q = -0.5 * (b + disc.sqrt().copysign(b));
let r0 = q / a;
let r1 = c / q;
Some((r0.min(r1), r0.max(r1)))
} else {
None
}
} else {
self.positive_discriminant_roots_scaled()
}
}
#[cold]
fn positive_discriminant_roots_scaled(&self) -> Option<(f64, f64)> {
if self.is_finite() {
let scale = 2.0f64.powi(-515);
(*self * scale).positive_discriminant_roots()
} else {
None
}
}
pub fn positive_discriminant_roots_no_overflow_check(&self) -> Option<(f64, f64)> {
let a = self.c2;
let b = self.c1;
let c = self.c0;
let disc = b * b - 4.0 * a * c;
if disc > 0.0 {
let q = -0.5 * (b + disc.sqrt().copysign(b));
let r0 = q / a;
let r1 = c / q;
Some((r0.min(r1), r0.max(r1)))
} else {
None
}
}
pub fn positive_discriminant_roots_no_overflow_check_half_b(&self) -> Option<(f64, f64)> {
let a = self.c2;
let b = self.c1;
let c = self.c0;
let disc = b * b - a * c;
if disc > 0.0 {
let q = -(b + disc.sqrt().copysign(b));
let r0 = q / a;
let r1 = c / q;
Some((r0.min(r1), r0.max(r1)))
} else {
None
}
}
}
#[cfg(test)]
mod tests {
#[test]
fn root_evaluation() {
arbtest::arbtest(|u| {
let q = crate::arbitrary::quadratic(u)?;
// Arbitrary quadratics can have coefficients with wild magnitudes,
// so we need to adjust our error expectations accordingly.
let magnitude = q.c0.abs().max(q.c1.abs()).max(q.c2.abs()).max(1.0);
for r in q.roots() {
let y = q.eval(r);
// To evaluate the polynomial, we need to square r, so our error
// should be relative to the magnitude of r squared.
let r_magnitude = r.abs().max(1.0);
let threshold = r_magnitude * 1e-14 * r_magnitude * magnitude;
if y.is_finite() && threshold.is_finite() {
assert!(y.abs() <= threshold);
}
}
Ok(())
})
.budget_ms(5_000);
}
}