ogeom_math/
interval.rs

1//! Interval arithmetic for filtered predicates.
2//!
3//! An [`Interval`] encloses a real number the program cannot represent: every
4//! operation returns bounds that certainly contain the true result, obtained
5//! by computing with the IEEE operations (which are correctly rounded, so
6//! the true value lies within one ulp of the computed one) and then widening
7//! each bound one step outward. The enclosure is conservative, never wrong.
8//!
9//! The point of carrying bounds is [`Interval::certain_sign`]: a sign
10//! decision made through an interval is either *certain*, because zero lies
11//! outside the bounds, or honestly undecided, because it does not. That is
12//! the filter a predicate wants: answer fast when floating point can, and
13//! say so when it cannot, instead of reading rounding noise as a direction.
14//!
15//! The vocabulary is the arithmetic predicates need: add, subtract, multiply,
16//! negate, square, absolute value, square root, and division away from zero.
17//! Transcendentals are deliberately absent; the standard library does not
18//! state error bounds for them, and an enclosure that might not enclose is
19//! worse than none.
20
21use ogeom_core::predicates::Sign;
22
23/// A closed interval `[lo, hi]` certainly containing a true real value.
24#[derive(Debug, Clone, Copy, PartialEq)]
25pub struct Interval {
26    lo: f64,
27    hi: f64,
28}
29
30impl Interval {
31    /// The degenerate interval holding one exactly-represented value.
32    #[must_use]
33    pub const fn point(value: f64) -> Self {
34        Self {
35            lo: value,
36            hi: value,
37        }
38    }
39
40    /// An interval from stated bounds. Panics if `lo > hi` or either is NaN,
41    /// because such an "enclosure" encloses nothing.
42    #[must_use]
43    pub fn new(lo: f64, hi: f64) -> Self {
44        assert!(lo <= hi, "an interval needs ordered, comparable bounds");
45        Self { lo, hi }
46    }
47
48    /// A value with a symmetric absolute uncertainty.
49    #[must_use]
50    pub fn about(value: f64, radius: f64) -> Self {
51        assert!(radius >= 0.0, "an uncertainty is a magnitude");
52        Self::new(value - radius, value + radius)
53    }
54
55    /// The lower bound.
56    #[must_use]
57    pub const fn lo(&self) -> f64 {
58        self.lo
59    }
60
61    /// The upper bound.
62    #[must_use]
63    pub const fn hi(&self) -> f64 {
64        self.hi
65    }
66
67    /// The width of the enclosure.
68    #[must_use]
69    pub fn width(&self) -> f64 {
70        self.hi - self.lo
71    }
72
73    /// Whether the enclosure contains `value`.
74    #[must_use]
75    pub fn contains(&self, value: f64) -> bool {
76        self.lo <= value && value <= self.hi
77    }
78
79    /// The sign of the true value, where the bounds decide it: `None` means
80    /// zero lies inside the enclosure and floating point genuinely cannot
81    /// tell, which is an answer, not a failure.
82    #[must_use]
83    pub fn certain_sign(&self) -> Option<Sign> {
84        if self.lo > 0.0 {
85            Some(Sign::Positive)
86        } else if self.hi < 0.0 {
87            Some(Sign::Negative)
88        } else if self.lo == 0.0 && self.hi == 0.0 {
89            Some(Sign::Zero)
90        } else {
91            None
92        }
93    }
94
95    /// The negation, exact; negation never rounds.
96    #[must_use]
97    pub const fn neg(&self) -> Self {
98        Self {
99            lo: -self.hi,
100            hi: -self.lo,
101        }
102    }
103
104    /// The sum.
105    #[must_use]
106    pub fn add(&self, other: &Self) -> Self {
107        Self {
108            lo: (self.lo + other.lo).next_down(),
109            hi: (self.hi + other.hi).next_up(),
110        }
111    }
112
113    /// The difference.
114    #[must_use]
115    pub fn sub(&self, other: &Self) -> Self {
116        Self {
117            lo: (self.lo - other.hi).next_down(),
118            hi: (self.hi - other.lo).next_up(),
119        }
120    }
121
122    /// The product: the extremes over the four bound products, widened.
123    #[must_use]
124    pub fn mul(&self, other: &Self) -> Self {
125        let products = [
126            self.lo * other.lo,
127            self.lo * other.hi,
128            self.hi * other.lo,
129            self.hi * other.hi,
130        ];
131        let mut lo = products[0];
132        let mut hi = products[0];
133        for p in &products[1..] {
134            lo = lo.min(*p);
135            hi = hi.max(*p);
136        }
137        Self {
138            lo: lo.next_down(),
139            hi: hi.next_up(),
140        }
141    }
142
143    /// The square, tighter than `mul` with itself, because a square cannot
144    /// be negative even when the interval straddles zero.
145    #[must_use]
146    pub fn square(&self) -> Self {
147        let (a, b) = (self.lo * self.lo, self.hi * self.hi);
148        if self.lo <= 0.0 && self.hi >= 0.0 {
149            Self {
150                lo: 0.0,
151                hi: a.max(b).next_up(),
152            }
153        } else {
154            Self {
155                lo: a.min(b).next_down().max(0.0),
156                hi: a.max(b).next_up(),
157            }
158        }
159    }
160
161    /// The absolute value, exact.
162    #[must_use]
163    pub fn abs(&self) -> Self {
164        if self.lo >= 0.0 {
165            *self
166        } else if self.hi <= 0.0 {
167            self.neg()
168        } else {
169            Self {
170                lo: 0.0,
171                hi: self.hi.max(-self.lo),
172            }
173        }
174    }
175
176    /// The square root, for enclosures of non-negative values. A lower bound
177    /// pushed below zero by widening is clamped; the true value it encloses
178    /// was non-negative. An interval entirely below zero has no real root
179    /// and returns `None`.
180    #[must_use]
181    pub fn sqrt(&self) -> Option<Self> {
182        if self.hi < 0.0 {
183            return None;
184        }
185        let lo = if self.lo <= 0.0 {
186            0.0
187        } else {
188            self.lo.sqrt().next_down().max(0.0)
189        };
190        Some(Self {
191            lo,
192            hi: self.hi.sqrt().next_up(),
193        })
194    }
195
196    /// The quotient, defined only when the divisor certainly excludes zero.
197    #[must_use]
198    pub fn checked_div(&self, other: &Self) -> Option<Self> {
199        if other.lo <= 0.0 && other.hi >= 0.0 {
200            return None;
201        }
202        let quotients = [
203            self.lo / other.lo,
204            self.lo / other.hi,
205            self.hi / other.lo,
206            self.hi / other.hi,
207        ];
208        let mut lo = quotients[0];
209        let mut hi = quotients[0];
210        for q in &quotients[1..] {
211            lo = lo.min(*q);
212            hi = hi.max(*q);
213        }
214        Some(Self {
215            lo: lo.next_down(),
216            hi: hi.next_up(),
217        })
218    }
219}
220
221#[cfg(test)]
222#[allow(clippy::unwrap_used)]
223mod tests {
224    use super::*;
225    use proptest::prelude::*;
226
227    fn finite() -> impl Strategy<Value = f64> {
228        // Magnitudes a kernel actually computes with, both signs, zero included.
229        prop_oneof![
230            Just(0.0),
231            -1e12..1e12f64,
232            (-1.0..1.0f64).prop_map(|x| x * 1e-12),
233        ]
234    }
235
236    proptest! {
237        /// The defining law: an operation on point intervals encloses the
238        /// floating-point result, which is within half an ulp of the truth.
239        #[test]
240        fn point_operations_enclose_their_own_result(a in finite(), b in finite()) {
241            let (x, y) = (Interval::point(a), Interval::point(b));
242            prop_assert!(x.add(&y).contains(a + b));
243            prop_assert!(x.sub(&y).contains(a - b));
244            prop_assert!(x.mul(&y).contains(a * b));
245            prop_assert!(x.square().contains(a * a));
246            prop_assert!(x.abs().contains(a.abs()));
247            if a >= 0.0 {
248                prop_assert!(x.sqrt().unwrap().contains(a.sqrt()));
249            }
250            if b != 0.0 {
251                prop_assert!(x.checked_div(&y).unwrap().contains(a / b));
252            }
253        }
254
255        /// Containment is monotone: whatever holds a value before an
256        /// operation holds the operated value after.
257        #[test]
258        fn enclosures_stay_enclosures(a in finite(), b in finite(), r in 0.0..1e-6f64) {
259            let x = Interval::about(a, r);
260            let y = Interval::about(b, r);
261            // The true values are a and b themselves; every combination of
262            // them must land inside.
263            prop_assert!(x.add(&y).contains(a + b));
264            prop_assert!(x.mul(&y).contains(a * b));
265            prop_assert!(x.sub(&y).contains(a - b));
266        }
267    }
268
269    #[test]
270    fn signs_are_certain_only_away_from_zero() {
271        assert_eq!(
272            Interval::new(1e-300, 2e-300).certain_sign(),
273            Some(Sign::Positive)
274        );
275        assert_eq!(
276            Interval::new(-2.0, -1e-300).certain_sign(),
277            Some(Sign::Negative)
278        );
279        assert_eq!(Interval::point(0.0).certain_sign(), Some(Sign::Zero));
280        assert_eq!(Interval::new(-1e-300, 1e-300).certain_sign(), None);
281    }
282
283    #[test]
284    fn the_classic_rounding_case_is_enclosed() {
285        // 0.1 + 0.2 in doubles is famously not 0.3; the enclosure holds the
286        // computed sum and stays within a couple of ulps.
287        let z = Interval::point(0.1).add(&Interval::point(0.2));
288        assert!(z.contains(0.1 + 0.2));
289        assert!(z.width() <= 4.0 * f64::EPSILON);
290    }
291
292    #[test]
293    fn squares_of_straddling_intervals_start_at_zero() {
294        let s = Interval::new(-2.0, 3.0).square();
295        assert_eq!(s.lo(), 0.0);
296        assert!(s.contains(9.0) && s.contains(0.25));
297    }
298
299    #[test]
300    fn division_through_zero_refuses() {
301        assert!(
302            Interval::point(1.0)
303                .checked_div(&Interval::new(-1.0, 1.0))
304                .is_none()
305        );
306    }
307
308    /// The use the type exists for: a cross-product magnitude computed from
309    /// normals carrying a stated residual either certainly clears a floor or
310    /// the interval says the arithmetic cannot decide.
311    #[test]
312    fn a_filtered_gate_decision() {
313        let residual = 1e-12;
314        let (ax, ay) = (
315            Interval::about(1.0, residual),
316            Interval::about(0.0, residual),
317        );
318        let (bx, by) = (
319            Interval::about(1.0, residual),
320            Interval::about(1e-6, residual),
321        );
322        // The 2D cross product ax*by - ay*bx encloses the true sine-scale.
323        let cross = ax.mul(&by).sub(&ay.mul(&bx));
324        assert_eq!(cross.certain_sign(), Some(Sign::Positive));
325        // Shrink the angle to the residual's own scale and the sign is
326        // honestly undecided.
327        let by = Interval::about(1e-12, residual);
328        let cross = ax.mul(&by).sub(&ay.mul(&bx));
329        assert_eq!(cross.certain_sign(), None);
330    }
331}
ogeom_math/interval.rs

ogeom_math/
interval.rs
