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. It never certifies a
11//! sign the dynamic filter would reject, because it is a strictly larger
12//! bound. A caller that exceeds the declared range gets `None` rather than a
13//! silently invalid answer.
14
15use axiolid_core::{Point2, Point3};
16use axiolid_guarantees::Sign;
17
18/// Machine epsilon for binary64.
19const EPSILON: f64 = f64::EPSILON / 2.0;
20
21/// Precomputed error bounds for a declared coordinate range.
22///
23/// Construct once per model, reuse across every predicate call.
24#[derive(Debug, Clone, Copy, PartialEq)]
25pub struct StaticFilter {
26    /// Maximum absolute coordinate value this filter is valid for.
27    bound: f64,
28    /// Absolute error bound for the `orient2d` determinant.
29    orient2d: f64,
30    /// Absolute error bound for the `orient3d` determinant.
31    orient3d: f64,
32}
33
34impl StaticFilter {
35    /// Build bounds valid for coordinates with `|x| <= bound`.
36    ///
37    /// Returns `None` for a non-finite or non-positive bound, and for a bound
38    /// so large that the derived error bound is not finite: in both cases no
39    /// sign could be certified, and returning a filter that always defers
40    /// would hide the configuration error.
41    #[must_use]
42    pub fn new(bound: f64) -> Option<Self> {
43        if !bound.is_finite() || bound <= 0.0 {
44            return None;
45        }
46        // A coordinate difference is at most 2*bound, so a 2x2 determinant
47        // term is at most (2*bound)^2 and a 3x3 term (2*bound)^3.
48        let span = 2.0 * bound;
49        let orient2d = (3.0 + 16.0 * EPSILON) * EPSILON * (2.0 * span * span);
50        let orient3d = (7.0 + 56.0 * EPSILON) * EPSILON * (6.0 * span * span * span);
51        if !orient2d.is_finite() || !orient3d.is_finite() {
52            return None;
53        }
54        Some(Self {
55            bound,
56            orient2d,
57            orient3d,
58        })
59    }
60
61    /// The coordinate range this filter was built for.
62    #[must_use]
63    pub const fn bound(self) -> f64 {
64        self.bound
65    }
66
67    /// Whether every coordinate of a 2D point is inside the declared range.
68    #[must_use]
69    fn covers2(self, p: Point2) -> bool {
70        p.x.abs() <= self.bound && p.y.abs() <= self.bound
71    }
72
73    /// Whether every coordinate of a 3D point is inside the declared range.
74    #[must_use]
75    fn covers3(self, p: Point3) -> bool {
76        p.x.abs() <= self.bound && p.y.abs() <= self.bound && p.z.abs() <= self.bound
77    }
78}
79
80impl StaticFilter {
81    /// Try to settle `orient2d` with the precomputed bound.
82    ///
83    /// `None` means "not settled": either a point lies outside the declared
84    /// range, or the determinant is too small for the static bound. The caller
85    /// must fall back to the full predicate, which is always correct.
86    #[must_use]
87    pub fn orient2d(self, a: Point2, b: Point2, c: Point2) -> Option<Sign> {
88        if !(self.covers2(a) && self.covers2(b) && self.covers2(c)) {
89            return None;
90        }
91        let determinant = (a.x - c.x) * (b.y - c.y) - (a.y - c.y) * (b.x - c.x);
92        decide(determinant, self.orient2d)
93    }
94
95    /// Try to settle `orient3d` with the precomputed bound.
96    #[must_use]
97    pub fn orient3d(self, a: Point3, b: Point3, c: Point3, d: Point3) -> Option<Sign> {
98        if !(self.covers3(a) && self.covers3(b) && self.covers3(c) && self.covers3(d)) {
99            return None;
100        }
101        let (adx, ady, adz) = (a.x - d.x, a.y - d.y, a.z - d.z);
102        let (bdx, bdy, bdz) = (b.x - d.x, b.y - d.y, b.z - d.z);
103        let (cdx, cdy, cdz) = (c.x - d.x, c.y - d.y, c.z - d.z);
104        let determinant = adz * (bdx * cdy - cdx * bdy)
105            + bdz * (cdx * ady - adx * cdy)
106            + cdz * (adx * bdy - bdx * ady);
107        decide(determinant, self.orient3d)
108    }
109}
110
111/// Certify a sign only when the magnitude strictly exceeds the bound.
112///
113/// Strict, not `>=`: a determinant exactly equal to its own error bound could
114/// have true value zero, so the sign is not proven.
115#[inline]
116#[must_use]
117fn decide(determinant: f64, bound: f64) -> Option<Sign> {
118    if determinant > bound {
119        Some(Sign::Positive)
120    } else if determinant < -bound {
121        Some(Sign::Negative)
122    } else {
123        None
124    }
125}
axiolid_predicates/static_filter.rs

axiolid_predicates/
static_filter.rs
