axiolid_predicates/
static_filter.rs

1//! Static filters: bounds computed once from a coordinate range.
2//!
3//! The dynamic filter in each predicate computes a *permanent* -- a sum of
4//! absolute products -- on every call, which costs roughly as much as the
5//! determinant itself. When a caller can state an upper bound on coordinate
6//! magnitude up front (a model's bounding box, a quantisation grid), the error
7//! bound can be precomputed and the per-call work drops to one comparison.
8//!
9//! The trade is coverage, not correctness: a static bound is necessarily
10//! looser than a per-input one, so it defers more often. Its bound covers
11//! every rounding error of the f64 evaluation, so a sign it certifies is the
12//! exact sign. It does not defer wherever the dynamic filter does: the dynamic
13//! `orient2d` filter declines differences outside [1e-90, 1e90], which the
14//! static one may still settle. A caller that exceeds the declared range gets
15//! `None` rather than a silently invalid answer, and a declared range so small
16//! that a bound would sit below the normal range, where products lose their
17//! relative accuracy, is refused by [`StaticFilter::new`] (#429).
18
19use axiolid_core::{Point2, Point3};
20use axiolid_guarantees::Sign;
21
22/// Machine epsilon for binary64.
23const EPSILON: f64 = f64::EPSILON / 2.0;
24
25/// Precomputed error bounds for a declared coordinate range.
26///
27/// Construct once per model, reuse across every predicate call.
28#[derive(Debug, Clone, Copy, PartialEq)]
29pub struct StaticFilter {
30    /// Maximum absolute coordinate value this filter is valid for.
31    bound: f64,
32    /// Absolute error bound for the `orient2d` determinant.
33    orient2d: f64,
34    /// Absolute error bound for the `orient3d` determinant.
35    orient3d: f64,
36}
37
38impl StaticFilter {
39    /// Build bounds valid for coordinates with `|x| <= bound`.
40    ///
41    /// Returns `None` for a non-finite or non-positive bound, and for a bound
42    /// so large that the derived error bound is not finite: in both cases no
43    /// sign could be certified, and returning a filter that always defers
44    /// would hide the configuration error. A bound so small that a derived
45    /// error bound falls below the normal range is refused as well: the
46    /// determinant's products then underflow, the bound no longer covers the
47    /// rounding error, and a bound of zero would certify every nonzero sign.
48    #[must_use]
49    pub fn new(bound: f64) -> Option<Self> {
50        if !bound.is_finite() || bound <= 0.0 {
51            return None;
52        }
53        // A coordinate difference is at most 2*bound, so a 2x2 determinant
54        // term is at most (2*bound)^2 and a 3x3 term (2*bound)^3.
55        let span = 2.0 * bound;
56        let orient2d = (3.0 + 16.0 * EPSILON) * EPSILON * (2.0 * span * span);
57        let orient3d = (7.0 + 56.0 * EPSILON) * EPSILON * (6.0 * span * span * span);
58        let usable = |bound: f64| bound.is_finite() && bound >= f64::MIN_POSITIVE;
59        if !usable(orient2d) || !usable(orient3d) {
60            return None;
61        }
62        Some(Self {
63            bound,
64            orient2d,
65            orient3d,
66        })
67    }
68
69    /// The coordinate range this filter was built for.
70    #[must_use]
71    pub const fn bound(self) -> f64 {
72        self.bound
73    }
74
75    /// Whether every coordinate of a 2D point is inside the declared range.
76    #[must_use]
77    fn covers2(self, p: Point2) -> bool {
78        p.x.abs() <= self.bound && p.y.abs() <= self.bound
79    }
80
81    /// Whether every coordinate of a 3D point is inside the declared range.
82    #[must_use]
83    fn covers3(self, p: Point3) -> bool {
84        p.x.abs() <= self.bound && p.y.abs() <= self.bound && p.z.abs() <= self.bound
85    }
86}
87
88impl StaticFilter {
89    /// Try to settle `orient2d` with the precomputed bound.
90    ///
91    /// `None` means "not settled": either a point lies outside the declared
92    /// range, or the determinant is too small for the static bound. The caller
93    /// must fall back to the full predicate, which is always correct.
94    #[must_use]
95    pub fn orient2d(self, a: Point2, b: Point2, c: Point2) -> Option<Sign> {
96        if !(self.covers2(a) && self.covers2(b) && self.covers2(c)) {
97            return None;
98        }
99        let determinant = (a.x - c.x) * (b.y - c.y) - (a.y - c.y) * (b.x - c.x);
100        decide(determinant, self.orient2d)
101    }
102
103    /// Try to settle `orient3d` with the precomputed bound.
104    #[must_use]
105    pub fn orient3d(self, a: Point3, b: Point3, c: Point3, d: Point3) -> Option<Sign> {
106        if !(self.covers3(a) && self.covers3(b) && self.covers3(c) && self.covers3(d)) {
107            return None;
108        }
109        let (adx, ady, adz) = (a.x - d.x, a.y - d.y, a.z - d.z);
110        let (bdx, bdy, bdz) = (b.x - d.x, b.y - d.y, b.z - d.z);
111        let (cdx, cdy, cdz) = (c.x - d.x, c.y - d.y, c.z - d.z);
112        let determinant = adz * (bdx * cdy - cdx * bdy)
113            + bdz * (cdx * ady - adx * cdy)
114            + cdz * (adx * bdy - bdx * ady);
115        decide(determinant, self.orient3d)
116    }
117}
118
119/// Certify a sign only when the magnitude strictly exceeds the bound.
120///
121/// Strict, not `>=`: a determinant exactly equal to its own error bound could
122/// have true value zero, so the sign is not proven.
123#[inline]
124#[must_use]
125fn decide(determinant: f64, bound: f64) -> Option<Sign> {
126    if determinant > bound {
127        Some(Sign::Positive)
128    } else if determinant < -bound {
129        Some(Sign::Negative)
130    } else {
131        None
132    }
133}
axiolid_predicates/static_filter.rs

axiolid_predicates/
static_filter.rs
