axiolid_predicates/
orient3.rs

1//! `orient3d`: which side of a plane a point lies on.
2//!
3//! Returns the sign of the 3x3 determinant
4//!
5//! ```text
6//! | ax-dx  ay-dy  az-dz |
7//! | bx-dx  by-dy  bz-dz |
8//! | cx-dx  cy-dy  cz-dz |
9//! ```
10//!
11//! Positive means `d` sees `a, b, c` counter-clockwise, i.e. `d` is below the
12//! plane under the right-hand rule. Zero means the four points are exactly
13//! coplanar -- the case that decides tetrahedralisation, convex hull facets,
14//! and whether a boolean surface passes through a vertex.
15//!
16//! Same filtered cascade as `orient2d`: cheap f64 with an error bound, then
17//! exact expansion arithmetic when the bound cannot exclude zero.
18
19use axiolid_core::Point3;
20use axiolid_guarantees::{Certified, Precision, Sign};
21
22use crate::arithmetic::{
23    expansion_product, expansion_sign, expansion_sum, grow_expansion, negate_expansion,
24    scale_expansion,
25};
26use crate::expansion::{two_diff, two_product};
27use crate::orient3_dyadic::orient3d_exact_dyadic;
28
29/// Machine epsilon for binary64.
30const EPSILON: f64 = f64::EPSILON / 2.0;
31
32/// Relative error bound for the 3x3 determinant filter.
33///
34/// The determinant is a sum of three 2x2 cofactor products. Propagating the
35/// `(1 + eps)` model over that expression gives `7*eps`; the second-order term
36/// absorbs the remainder so this is a true upper bound.
37///
38/// The constant is deliberately conservative. A mutation probe lowering it to
39/// `3*eps` does not fail the suite -- the filter stays sound at that value for
40/// the inputs generated -- while `0.05*eps` is caught immediately by
41/// `near_degenerate_cases_recover_a_definite_sign`. The margin between 3 and 7
42/// buys nothing measurable in throughput (the escalation-rate gates show the
43/// filter settles clean data either way) and costs nothing, so the derivation's
44/// value is kept rather than the empirically minimal one: correctness here is
45/// argued from the error model, not tuned against a test suite.
46const ORIENT3D_ERROR_FACTOR: f64 = (7.0 + 56.0 * EPSILON) * EPSILON;
47
48/// Orientation of `d` relative to the plane through `a`, `b`, `c`.
49///
50/// Returns an exact certified sign for finite binary64 coordinates. The usual
51/// expansion path handles inputs whose intermediates remain representable;
52/// other finite exponent ranges fall back to a fixed-size exact dyadic
53/// accumulator. Non-finite input returns [`Certified::Uncertain`] rather than
54/// guessing.
55#[must_use]
56pub fn orient3d(a: Point3, b: Point3, c: Point3, d: Point3) -> Certified {
57    match orient3d_filter(a, b, c, d) {
58        Certified::Certain { sign, .. } => Certified::exact_sign(sign),
59        // Non-exhaustive enum: anything we do not recognise must escalate.
60        _ => orient3d_exact(a, b, c, d)
61            .or_else(|| orient3d_exact_dyadic(a, b, c, d))
62            .map_or(
63                Certified::Uncertain {
64                    attempted: Precision::Exact,
65                },
66                Certified::exact_sign,
67            ),
68    }
69}
70
71/// The fast filter alone, exposed so escalation can be measured.
72#[must_use]
73pub fn orient3d_filter(a: Point3, b: Point3, c: Point3, d: Point3) -> Certified {
74    let (adx, ady, adz) = (a.x - d.x, a.y - d.y, a.z - d.z);
75    let (bdx, bdy, bdz) = (b.x - d.x, b.y - d.y, b.z - d.z);
76    let (cdx, cdy, cdz) = (c.x - d.x, c.y - d.y, c.z - d.z);
77
78    let differences = [adx, ady, adz, bdx, bdy, bdz, cdx, cdy, cdz];
79    if !differences.iter().all(|value| {
80        value.is_finite() && (*value == 0.0 || (value.abs() >= 1.0e-90 && value.abs() <= 1.0e90))
81    }) {
82        return Certified::Uncertain {
83            attempted: Precision::F64,
84        };
85    }
86
87    let bdxcdy = bdx * cdy;
88    let cdxbdy = cdx * bdy;
89    let cdxady = cdx * ady;
90    let adxcdy = adx * cdy;
91    let adxbdy = adx * bdy;
92    let bdxady = bdx * ady;
93
94    let determinant = adz * (bdxcdy - cdxbdy) + bdz * (cdxady - adxcdy) + cdz * (adxbdy - bdxady);
95
96    // The bound tracks operand magnitudes, so it scales with the model's units
97    // instead of assuming a coordinate range.
98    let permanent = (bdxcdy.abs() + cdxbdy.abs()) * adz.abs()
99        + (cdxady.abs() + adxcdy.abs()) * bdz.abs()
100        + (adxbdy.abs() + bdxady.abs()) * cdz.abs();
101
102    Certified::from_filter(
103        determinant,
104        ORIENT3D_ERROR_FACTOR * permanent,
105        Precision::F64,
106    )
107}
108
109/// Exact sign of the 3x3 determinant.
110///
111/// Every coordinate difference includes the tail recovered by [`two_diff`].
112/// Cofactors and their remaining products operate on complete expansions, so
113/// the result is the determinant of the input binary64 coordinates rather than
114/// of rounded coordinate differences.
115#[must_use]
116fn orient3d_exact(a: Point3, b: Point3, c: Point3, d: Point3) -> Option<Sign> {
117    let (adx, adx_error) = two_diff(a.x, d.x);
118    let (ady, ady_error) = two_diff(a.y, d.y);
119    let (adz, adz_error) = two_diff(a.z, d.z);
120    let (bdx, bdx_error) = two_diff(b.x, d.x);
121    let (bdy, bdy_error) = two_diff(b.y, d.y);
122    let (bdz, bdz_error) = two_diff(b.z, d.z);
123    let (cdx, cdx_error) = two_diff(c.x, d.x);
124    let (cdy, cdy_error) = two_diff(c.y, d.y);
125    let (cdz, cdz_error) = two_diff(c.z, d.z);
126
127    let differences = [[adx, ady, adz], [bdx, bdy, bdz], [cdx, cdy, cdz]];
128    let errors = [
129        [adx_error, ady_error, adz_error],
130        [bdx_error, bdy_error, bdz_error],
131        [cdx_error, cdy_error, cdz_error],
132    ];
133    if !exact_components_are_representable(&differences, &errors) {
134        return None;
135    }
136
137    if errors.iter().flatten().all(|error| *error == 0.0) {
138        return Some(orient3d_exact_differences(
139            differences[0],
140            differences[1],
141            differences[2],
142        ));
143    }
144
145    let [[adx, ady, adz], [bdx, bdy, bdz], [cdx, cdy, cdz]] = differences;
146    let [[adx_error, ady_error, adz_error], [bdx_error, bdy_error, bdz_error], [cdx_error, cdy_error, cdz_error]] =
147        errors;
148
149    let adx = difference_expansion(adx, adx_error);
150    let ady = difference_expansion(ady, ady_error);
151    let adz = difference_expansion(adz, adz_error);
152    let bdx = difference_expansion(bdx, bdx_error);
153    let bdy = difference_expansion(bdy, bdy_error);
154    let bdz = difference_expansion(bdz, bdz_error);
155    let cdx = difference_expansion(cdx, cdx_error);
156    let cdy = difference_expansion(cdy, cdy_error);
157    let cdz = difference_expansion(cdz, cdz_error);
158
159    let bc = orient3d_expansion_cofactor(&bdx, &cdy, &cdx, &bdy);
160    let ca = orient3d_expansion_cofactor(&cdx, &ady, &adx, &cdy);
161    let ab = orient3d_expansion_cofactor(&adx, &bdy, &bdx, &ady);
162
163    let total = expansion_sum(
164        &expansion_sum(&expansion_product(&bc, &adz), &expansion_product(&ca, &bdz)),
165        &expansion_product(&ab, &cdz),
166    );
167    Some(expansion_sign(&total))
168}
169
170/// Whether the ordinary expansion evaluator can represent every intermediate.
171///
172/// Components above exponent 300 could overflow a three-factor term. The
173/// monomial check separately proves that each two- and three-factor product
174/// stays above binary64's underflow floor.
175#[must_use]
176fn exact_components_are_representable(differences: &[[f64; 3]; 3], errors: &[[f64; 3]; 3]) -> bool {
177    differences
178        .iter()
179        .flatten()
180        .chain(errors.iter().flatten())
181        .all(|value| value.is_finite() && (*value == 0.0 || highest_bit_exponent(*value) <= 300))
182        && determinant_terms_are_representable(differences, errors)
183}
184
185#[must_use]
186fn determinant_terms_are_representable(
187    differences: &[[f64; 3]; 3],
188    errors: &[[f64; 3]; 3],
189) -> bool {
190    let mut least_bits = [[None; 3]; 3];
191    for row in 0..3 {
192        for column in 0..3 {
193            least_bits[row][column] = [differences[row][column], errors[row][column]]
194                .into_iter()
195                .filter(|value| *value != 0.0)
196                .map(least_significant_bit_exponent)
197                .min();
198        }
199    }
200
201    [
202        [(1, 0), (2, 1), (0, 2)],
203        [(2, 0), (1, 1), (0, 2)],
204        [(2, 0), (0, 1), (1, 2)],
205        [(0, 0), (2, 1), (1, 2)],
206        [(0, 0), (1, 1), (2, 2)],
207        [(1, 0), (0, 1), (2, 2)],
208    ]
209    .into_iter()
210    .all(|term| {
211        let exponents = term.map(|(row, column)| least_bits[row][column]);
212        let [Some(first), Some(second), Some(third)] = exponents else {
213            return true;
214        };
215        first + second >= -1074 && first + second + third >= -1074
216    })
217}
218
219#[must_use]
220fn highest_bit_exponent(value: f64) -> i32 {
221    let bits = value.abs().to_bits();
222    let encoded_exponent = ((bits >> 52) & 0x7ff) as i32;
223    if encoded_exponent == 0 {
224        let fraction = bits & ((1_u64 << 52) - 1);
225        -1074 + (63 - fraction.leading_zeros() as i32)
226    } else {
227        encoded_exponent - 1023
228    }
229}
230
231#[must_use]
232fn least_significant_bit_exponent(value: f64) -> i32 {
233    let bits = value.abs().to_bits();
234    let encoded_exponent = ((bits >> 52) & 0x7ff) as i32;
235    let fraction = bits & ((1_u64 << 52) - 1);
236    if encoded_exponent == 0 {
237        -1074 + fraction.trailing_zeros() as i32
238    } else {
239        let significand = (1_u64 << 52) | fraction;
240        encoded_exponent - 1023 - 52 + significand.trailing_zeros() as i32
241    }
242}
243
244/// Evaluate the determinant when all coordinate differences are already exact.
245#[must_use]
246fn orient3d_exact_differences(a: [f64; 3], b: [f64; 3], c: [f64; 3]) -> Sign {
247    let bc = orient3d_cofactor(b[0], c[1], c[0], b[1]);
248    let ca = orient3d_cofactor(c[0], a[1], a[0], c[1]);
249    let ab = orient3d_cofactor(a[0], b[1], b[0], a[1]);
250    let total = expansion_sum(
251        &expansion_sum(&scale_expansion(&bc, a[2]), &scale_expansion(&ca, b[2])),
252        &scale_expansion(&ab, c[2]),
253    );
254    expansion_sign(&total)
255}
256
257#[must_use]
258fn difference_expansion(difference: f64, error: f64) -> Vec<f64> {
259    let mut expansion = Vec::new();
260    if error != 0.0 {
261        expansion.push(error);
262    }
263    if difference != 0.0 || expansion.is_empty() {
264        expansion.push(difference);
265    }
266    expansion
267}
268
269#[must_use]
270fn orient3d_expansion_cofactor(p: &[f64], q: &[f64], r: &[f64], s: &[f64]) -> Vec<f64> {
271    expansion_sum(
272        &expansion_product(p, q),
273        &negate_expansion(&expansion_product(r, s)),
274    )
275}
276
277/// Exact `p*q - r*s` as an expansion.
278///
279/// Shared with the Delaunay predicates, which need the same 2x2 minors.
280#[must_use]
281pub(crate) fn orient3d_cofactor(p: f64, q: f64, r: f64, s: f64) -> Vec<f64> {
282    let (pq, pq_err) = two_product(p, q);
283    let (rs, rs_err) = two_product(r, s);
284    // Subtract by adding the negation; the four terms are combined by the
285    // carry-propagating grow so the result stays a valid expansion.
286    let e = grow_expansion(&[pq_err], -rs_err);
287    let e = grow_expansion(&e, pq);
288    grow_expansion(&e, -rs)
289}
axiolid_predicates/orient3.rs

axiolid_predicates/
orient3.rs
