mesh: exact predicates (Shewchuk adaptive, symbolic perturbation, indirect LPI/TPI) (#8)

* mesh: exact predicates (Shewchuk adaptive orient/incircle/insphere), symbolic perturbation, indirect LPI/TPI predicates

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* mesh: predicate benchmark (filter hit rate, CAD-degenerate inputs); record the predicate decision

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

* mesh: simplify predicates (destructuring over fallbacks, shared-minor exact insphere, cheaper expansion product)

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

---------

Co-authored-by: jyviko <[email protected]>
Co-authored-by: Claude Opus 5.5 <[email protected]>

ATTRIBUTION.md

@@ -23,3 +23,12 @@ is committed and listed here with its exports.
 | xtask/src/assets.rs, xtask/src/licences.rs, xtask/src/test.rs, xtask/src/main.rs | SolidCraft `xtask/src/` | MIT OR Apache-2.0 |
 | clippy.toml, rustfmt.toml | SolidCraft | MIT OR Apache-2.0 |
 | crates/geometry/src/tagged_mesh.rs (package and part reading, transforms, units, caps) | SolidCraft `crates/io/src/threemf.rs` | MIT OR Apache-2.0 |
+
+## Test oracles
+
+Used only by tests (dev-dependencies, not shipped); no code from them is copied or adapted.
+
+| Used in | Oracle | Licence |
+|---|---|---|
+| crates/mesh/src/predicates/tests.rs, crates/mesh/tests/predicates_bench.rs | `robust` (georust), a port of Shewchuk's public-domain `predicates.c` | MIT OR Apache-2.0 |
+| crates/mesh/src/predicates/tests.rs, crates/mesh/src/predicates/dyadic.rs (tests) | `num-bigint`, `num-rational`, `num-traits` (exact rationals) | MIT OR Apache-2.0 |

Cargo.lock

@@ -481,6 +481,10 @@ name = "modelcraft-mesh"
 version = "0.1.0"
 dependencies = [
  "modelcraft-model",
+ "num-bigint",
+ "num-rational",
+ "num-traits",
+ "robust",
  "thiserror 2.0.21",
 ]
 
@@ -493,6 +497,16 @@ dependencies = [
  "thiserror 2.0.21",
 ]
 
+[[package]]
+name = "num-bigint"
+version = "0.4.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c89e69e7e0f03bea5ef08013795c25018e101932225a656383bd384495ecc367"
+dependencies = [
+ "num-integer",
+ "num-traits",
+]
+
 [[package]]
 name = "num-complex"
 version = "0.4.6"
@@ -502,6 +516,26 @@ dependencies = [
  "num-traits",
 ]
 
+[[package]]
+name = "num-integer"
+version = "0.1.47"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7ce2d95d4b3734dc35aa2f45e1aa22cd416814592a4f9d9205e11affd5b8e10b"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-rational"
+version = "0.4.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824"
+dependencies = [
+ "num-bigint",
+ "num-integer",
+ "num-traits",
+]
+
 [[package]]
 name = "num-traits"
 version = "0.2.19"

ROADMAP.md

@@ -6,7 +6,7 @@ LS-PrePost, with hexahedral meshing as the centre of gravity, plus its own linea
 | Milestone | Scope | Status |
 |---|---|---|
 | M0 Foundation | Workspace, quality gates (`cargo xtask ci`), FE model, command engine with rollback and panic guard, CLI, geometry boundary on SolidCraft's kernel and tessellator | done |
-| M1 Surface preparation | Geometry in: open `.solidcraft` designs, rebuilt by SolidCraft, with persistent face names (done); face-tagged 3MF files (`geometry.import.mesh`, docs/formats/tagged-mesh.md, done); robust adaptive-precision predicates; weld, closed/manifold/orientation checks and feature edges on the tessellated boundary; mesh quality metrics | in progress |
+| M1 Surface preparation | Geometry in: open `.solidcraft` designs, rebuilt by SolidCraft, with persistent face names (done); face-tagged 3MF files (`geometry.import.mesh`, docs/formats/tagged-mesh.md, done); exact predicates (Shewchuk adaptive, symbolic perturbation, indirect LPI/TPI, done); weld, closed/manifold/orientation checks and feature edges on the tessellated boundary; mesh quality metrics | in progress |
 | M2 Tetrahedral meshing | Bowyer–Watson Delaunay, robust constrained Delaunay boundary recovery (Diazzi et al. 2023), Delaunay refinement to a size field, ODT smoothing, Klingner–Shewchuk sliver removal; linear and quadratic tets | planned |
 | M3 Hexahedral meshing (core) | Octree/grid dual all-hex (valid on any input), then hex-dominant by tet recombination, then polycube for blocky parts; see docs/meshing/literature.md | planned |
 | M4 Model definition | Materials, sections, sets and surfaces, boundary conditions, loads, steps, output requests | planned |

crates/mesh/Cargo.toml

@@ -11,3 +11,11 @@ workspace = true
 [dependencies]
 modelcraft-model = { workspace = true }
 thiserror = { workspace = true }
+
+[dev-dependencies]
+# Test oracles only (never linked into ModelCraft): exact rationals, and georust's port of
+# Shewchuk's public-domain predicates.c as a differential reference.
+num-bigint = "0.4"
+num-rational = { version = "0.4", default-features = false, features = ["num-bigint", "std"] }
+num-traits = "0.2"
+robust = "1.2"

crates/mesh/src/lib.rs

@@ -20,3 +20,5 @@ pub enum Error {
 }
 
 pub type Result<T> = std::result::Result<T, Error>;
+
+pub mod predicates;

crates/mesh/src/predicates/adaptive.rs

@@ -0,0 +1,449 @@
+//! Shewchuk's adaptive-precision `orient2d`, `orient3d`, `incircle` and `insphere`
+//! (Shewchuk 1997, §4 and §5), implemented from the paper.
+//!
+//! Each predicate evaluates its determinant in plain doubles and accepts the sign when the
+//! magnitude exceeds the paper's certified error bound (stage A, the filter). Otherwise it
+//! escalates: stage B computes the determinant of the rounded coordinate differences exactly as
+//! an expansion and accepts it under a tighter bound, or outright when the differences were exact
+//! (zero tails, the common case for grid-like CAD input); stage C (`orient2d`, `orient3d`) adds
+//! the first-order tail correction; the last stage evaluates the determinant exactly from the
+//! input coordinates with expansion arithmetic. Shewchuk's further hand-scheduled stages between
+//! C and exact are not reproduced: they only make rare inputs faster.
+//!
+//! The error bounds assume IEEE binary64 round-to-nearest-even arithmetic without overflow or
+//! underflow (§5); the coordinate domain of [`Point2`] / [`Point3`] rules both out.
+
+use super::eft::two_diff_tail;
+use super::expansion::{Exp, Expansion, cross_diff};
+use super::ring::{det2, det3, sq2, sq3};
+use super::{Point2, Point3, Sign};
+
+/// `ε = 2^-53`, half a unit in the last place of 1 (the machine epsilon of §2.1).
+const EPS: f64 = f64::EPSILON / 2.0;
+const RESULT_ERR: f64 = (3.0 + 8.0 * EPS) * EPS;
+const CCW_A: f64 = (3.0 + 16.0 * EPS) * EPS;
+const CCW_B: f64 = (2.0 + 12.0 * EPS) * EPS;
+const CCW_C: f64 = (9.0 + 64.0 * EPS) * EPS * EPS;
+const O3D_A: f64 = (7.0 + 56.0 * EPS) * EPS;
+const O3D_B: f64 = (3.0 + 28.0 * EPS) * EPS;
+const O3D_C: f64 = (26.0 + 288.0 * EPS) * EPS * EPS;
+const ICC_A: f64 = (10.0 + 96.0 * EPS) * EPS;
+const ICC_B: f64 = (4.0 + 48.0 * EPS) * EPS;
+const ISP_A: f64 = (16.0 + 224.0 * EPS) * EPS;
+const ISP_B: f64 = (5.0 + 72.0 * EPS) * EPS;
+
+#[inline(always)]
+fn certain(det: f64, errbound: f64) -> bool {
+    det > errbound || -det > errbound
+}
+
+// ---------------------------------------------------------------------------------------------
+// orient2d
+// ---------------------------------------------------------------------------------------------
+
+/// Stage A of [`orient2d`]: the sign when the double-precision filter certifies it.
+#[inline]
+pub fn orient2d_filter(a: &Point2, b: &Point2, c: &Point2) -> Option<Sign> {
+    orient2d_a(a, b, c).ok()
+}
+
+/// `Ok(sign)` when stage A decides, else `Err(detsum)` for the adaptive stages.
+#[inline(always)]
+fn orient2d_a(a: &Point2, b: &Point2, c: &Point2) -> Result<Sign, f64> {
+    let ([ax, ay], [bx, by], [cx, cy]) = (a.0, b.0, c.0);
+    let det_left = (ax - cx) * (by - cy);
+    let det_right = (ay - cy) * (bx - cx);
+    let det = det_left - det_right;
+    // When the two products have opposite signs (or one is zero) the subtraction cannot cancel
+    // and the rounded result has the exact sign.
+    let det_sum = if det_left > 0.0 {
+        if det_right <= 0.0 {
+            return Ok(Sign::of(det));
+        }
+        det_left + det_right
+    } else if det_left < 0.0 {
+        if det_right >= 0.0 {
+            return Ok(Sign::of(det));
+        }
+        -det_left - det_right
+    } else {
+        return Ok(Sign::of(det));
+    };
+    if certain(det, CCW_A * det_sum) { Ok(Sign::of(det)) } else { Err(det_sum) }
+}
+
+/// The orientation of `a, b, c`: [`Sign::Positive`] when they turn counterclockwise,
+/// [`Sign::Negative`] when clockwise, [`Sign::Zero`] when collinear. Exact. It is the sign of
+/// `det[a - c; b - c]`.
+#[inline]
+pub fn orient2d(a: &Point2, b: &Point2, c: &Point2) -> Sign {
+    match orient2d_a(a, b, c) {
+        Ok(sign) => sign,
+        Err(det_sum) => orient2d_adapt(a, b, c, det_sum),
+    }
+}
+
+#[cold]
+fn orient2d_adapt(a: &Point2, b: &Point2, c: &Point2, det_sum: f64) -> Sign {
+    let ([ax, ay], [bx, by], [cx, cy]) = (a.0, b.0, c.0);
+    let (acx, bcx, acy, bcy) = (ax - cx, bx - cx, ay - cy, by - cy);
+
+    // Stage B: the exact determinant of the rounded differences.
+    let b_exp = cross_diff(acx, bcy, acy, bcx);
+    let det = b_exp.estimate();
+    if certain(det, CCW_B * det_sum) {
+        return Sign::of(det);
+    }
+    let acx_t = two_diff_tail(ax, cx, acx);
+    let bcx_t = two_diff_tail(bx, cx, bcx);
+    let acy_t = two_diff_tail(ay, cy, acy);
+    let bcy_t = two_diff_tail(by, cy, bcy);
+    if acx_t == 0.0 && acy_t == 0.0 && bcx_t == 0.0 && bcy_t == 0.0 {
+        return b_exp.sign();
+    }
+
+    // Stage C: first-order correction for the tails.
+    let errbound = CCW_C * det_sum + RESULT_ERR * det.abs();
+    let det = det + ((acx * bcy_t + bcy * acx_t) - (acy * bcx_t + bcx * acy_t));
+    if certain(det, errbound) {
+        return Sign::of(det);
+    }
+
+    // Stage D: the remaining terms of (acx + acx_t)(bcy + bcy_t) - (acy + acy_t)(bcx + bcx_t).
+    let c1 = Exp::<8>::sum(&b_exp, &cross_diff(acx_t, bcy, acy_t, bcx));
+    let c2 = Exp::<12>::sum(&c1, &cross_diff(acx, bcy_t, acy, bcx_t));
+    Exp::<16>::sum(&c2, &cross_diff(acx_t, bcy_t, acy_t, bcx_t)).sign()
+}
+
+// ---------------------------------------------------------------------------------------------
+// orient3d
+// ---------------------------------------------------------------------------------------------
+
+/// Stage A of [`orient3d`].
+#[inline]
+pub fn orient3d_filter(a: &Point3, b: &Point3, c: &Point3, d: &Point3) -> Option<Sign> {
+    orient3d_a(a, b, c, d).ok()
+}
+
+#[inline(always)]
+fn orient3d_a(a: &Point3, b: &Point3, c: &Point3, d: &Point3) -> Result<Sign, f64> {
+    let ([ax, ay, az], [bx, by, bz], [cx, cy, cz], [dx, dy, dz]) = (a.0, b.0, c.0, d.0);
+    let (adx, bdx, cdx) = (ax - dx, bx - dx, cx - dx);
+    let (ady, bdy, cdy) = (ay - dy, by - dy, cy - dy);
+    let (adz, bdz, cdz) = (az - dz, bz - dz, cz - dz);
+    let (bdxcdy, cdxbdy) = (bdx * cdy, cdx * bdy);
+    let (cdxady, adxcdy) = (cdx * ady, adx * cdy);
+    let (adxbdy, bdxady) = (adx * bdy, bdx * ady);
+    let det = adz * (bdxcdy - cdxbdy) + bdz * (cdxady - adxcdy) + cdz * (adxbdy - bdxady);
+    let permanent = (bdxcdy.abs() + cdxbdy.abs()) * adz.abs() + (cdxady.abs() + adxcdy.abs()) * bdz.abs() + (adxbdy.abs() + bdxady.abs()) * cdz.abs();
+    if certain(det, O3D_A * permanent) { Ok(Sign::of(det)) } else { Err(permanent) }
+}
+
+/// The orientation of `d` relative to the plane through `a, b, c`: [`Sign::Positive`] when `d`
+/// lies below the plane, "above" being the side from which `a, b, c` appear counterclockwise;
+/// [`Sign::Zero`] when the four points are coplanar. Exact. It is the sign of
+/// `det[a - d; b - d; c - d] = -det[b - a; c - a; d - a]` (Shewchuk's convention), so it is
+/// positive when the tetrahedron `a, b, c, d` is negatively oriented in the right-handed sense.
+#[inline]
+pub fn orient3d(a: &Point3, b: &Point3, c: &Point3, d: &Point3) -> Sign {
+    match orient3d_a(a, b, c, d) {
+        Ok(sign) => sign,
+        Err(permanent) => orient3d_adapt(a, b, c, d, permanent),
+    }
+}
+
+#[cold]
+fn orient3d_adapt(a: &Point3, b: &Point3, c: &Point3, d: &Point3, permanent: f64) -> Sign {
+    let ([ax, ay, az], [bx, by, bz], [cx, cy, cz], [dx, dy, dz]) = (a.0, b.0, c.0, d.0);
+    let (adx, bdx, cdx) = (ax - dx, bx - dx, cx - dx);
+    let (ady, bdy, cdy) = (ay - dy, by - dy, cy - dy);
+    let (adz, bdz, cdz) = (az - dz, bz - dz, cz - dz);
+
+    // Stage B.
+    let bc = cross_diff(bdx, cdy, cdx, bdy);
+    let ca = cross_diff(cdx, ady, adx, cdy);
+    let ab = cross_diff(adx, bdy, bdx, ady);
+    let a_det = Exp::<8>::scale(&bc, adz);
+    let b_det = Exp::<8>::scale(&ca, bdz);
+    let c_det = Exp::<8>::scale(&ab, cdz);
+    let fin = Exp::<24>::sum(&Exp::<16>::sum(&a_det, &b_det), &c_det);
+    let det = fin.estimate();
+    if certain(det, O3D_B * permanent) {
+        return Sign::of(det);
+    }
+    let (adx_t, bdx_t, cdx_t) = (two_diff_tail(ax, dx, adx), two_diff_tail(bx, dx, bdx), two_diff_tail(cx, dx, cdx));
+    let (ady_t, bdy_t, cdy_t) = (two_diff_tail(ay, dy, ady), two_diff_tail(by, dy, bdy), two_diff_tail(cy, dy, cdy));
+    let (adz_t, bdz_t, cdz_t) = (two_diff_tail(az, dz, adz), two_diff_tail(bz, dz, bdz), two_diff_tail(cz, dz, cdz));
+    let tails = [adx_t, bdx_t, cdx_t, ady_t, bdy_t, cdy_t, adz_t, bdz_t, cdz_t];
+    if tails.iter().all(|&t| t == 0.0) {
+        return fin.sign();
+    }
+
+    // Stage C.
+    let errbound = O3D_C * permanent + RESULT_ERR * det.abs();
+    let det = det
+        + (adz * ((bdx * cdy_t + cdy * bdx_t) - (bdy * cdx_t + cdx * bdy_t)) + adz_t * (bdx * cdy - bdy * cdx))
+        + (bdz * ((cdx * ady_t + ady * cdx_t) - (cdy * adx_t + adx * cdy_t)) + bdz_t * (cdx * ady - cdy * adx))
+        + (cdz * ((adx * bdy_t + bdy * adx_t) - (ady * bdx_t + bdx * ady_t)) + cdz_t * (adx * bdy - ady * bdx));
+    if certain(det, errbound) {
+        return Sign::of(det);
+    }
+    orient3d_exact(a, b, c, d)
+}
+
+/// `det[a - d; b - d; c - d]` exactly, from the input coordinates.
+fn orient3d_exact(a: &Point3, b: &Point3, c: &Point3, d: &Point3) -> Sign {
+    det3(&[diffs3(a, d), diffs3(b, d), diffs3(c, d)]).sign()
+}
+
+/// The exact coordinate differences `p - o` as expansions.
+fn diffs3(p: &Point3, o: &Point3) -> [Expansion; 3] {
+    let ([px, py, pz], [ox, oy, oz]) = (p.0, o.0);
+    [Expansion::diff(px, ox), Expansion::diff(py, oy), Expansion::diff(pz, oz)]
+}
+
+fn diffs2(p: &Point2, o: &Point2) -> [Expansion; 2] {
+    let ([px, py], [ox, oy]) = (p.0, o.0);
+    [Expansion::diff(px, ox), Expansion::diff(py, oy)]
+}
+
+// ---------------------------------------------------------------------------------------------
+// incircle
+// ---------------------------------------------------------------------------------------------
+
+/// Stage A of [`incircle`].
+#[inline]
+pub fn incircle_filter(a: &Point2, b: &Point2, c: &Point2, d: &Point2) -> Option<Sign> {
+    incircle_a(a, b, c, d).ok()
+}
+
+#[inline(always)]
+fn incircle_a(a: &Point2, b: &Point2, c: &Point2, d: &Point2) -> Result<Sign, f64> {
+    let ([ax, ay], [bx, by], [cx, cy], [dx, dy]) = (a.0, b.0, c.0, d.0);
+    let (adx, bdx, cdx) = (ax - dx, bx - dx, cx - dx);
+    let (ady, bdy, cdy) = (ay - dy, by - dy, cy - dy);
+    let (bdxcdy, cdxbdy) = (bdx * cdy, cdx * bdy);
+    let a_lift = adx * adx + ady * ady;
+    let (cdxady, adxcdy) = (cdx * ady, adx * cdy);
+    let b_lift = bdx * bdx + bdy * bdy;
+    let (adxbdy, bdxady) = (adx * bdy, bdx * ady);
+    let c_lift = cdx * cdx + cdy * cdy;
+    let det = a_lift * (bdxcdy - cdxbdy) + b_lift * (cdxady - adxcdy) + c_lift * (adxbdy - bdxady);
+    let permanent = (bdxcdy.abs() + cdxbdy.abs()) * a_lift + (cdxady.abs() + adxcdy.abs()) * b_lift + (adxbdy.abs() + bdxady.abs()) * c_lift;
+    if certain(det, ICC_A * permanent) { Ok(Sign::of(det)) } else { Err(permanent) }
+}
+
+/// Whether `d` lies inside the circle through `a, b, c`: [`Sign::Positive`] inside,
+/// [`Sign::Negative`] outside, [`Sign::Zero`] on it, for `a, b, c` counterclockwise (the sign is
+/// reversed when they are clockwise). Exact. It is the sign of
+/// `det[a - d, |a - d|²; b - d, |b - d|²; c - d, |c - d|²]`.
+#[inline]
+pub fn incircle(a: &Point2, b: &Point2, c: &Point2, d: &Point2) -> Sign {
+    match incircle_a(a, b, c, d) {
+        Ok(sign) => sign,
+        Err(permanent) => incircle_adapt(a, b, c, d, permanent),
+    }
+}
+
+#[cold]
+fn incircle_adapt(a: &Point2, b: &Point2, c: &Point2, d: &Point2, permanent: f64) -> Sign {
+    let ([ax, ay], [bx, by], [cx, cy], [dx, dy]) = (a.0, b.0, c.0, d.0);
+    let (adx, bdx, cdx) = (ax - dx, bx - dx, cx - dx);
+    let (ady, bdy, cdy) = (ay - dy, by - dy, cy - dy);
+
+    // Stage B.
+    let lifted = |minor: &Exp<4>, x: f64, y: f64| -> Exp<32> {
+        Exp::sum(&Exp::<16>::scale(&Exp::<8>::scale(minor, x), x), &Exp::<16>::scale(&Exp::<8>::scale(minor, y), y))
+    };
+    let a_det = lifted(&cross_diff(bdx, cdy, cdx, bdy), adx, ady);
+    let b_det = lifted(&cross_diff(cdx, ady, adx, cdy), bdx, bdy);
+    let c_det = lifted(&cross_diff(adx, bdy, bdx, ady), cdx, cdy);
+    let fin = Exp::<96>::sum(&Exp::<64>::sum(&a_det, &b_det), &c_det);
+    let det = fin.estimate();
+    if certain(det, ICC_B * permanent) {
+        return Sign::of(det);
+    }
+    let tails = [
+        two_diff_tail(ax, dx, adx),
+        two_diff_tail(bx, dx, bdx),
+        two_diff_tail(cx, dx, cdx),
+        two_diff_tail(ay, dy, ady),
+        two_diff_tail(by, dy, bdy),
+        two_diff_tail(cy, dy, cdy),
+    ];
+    if tails.iter().all(|&t| t == 0.0) {
+        return fin.sign();
+    }
+    incircle_exact(a, b, c, d)
+}
+
+fn incircle_exact(a: &Point2, b: &Point2, c: &Point2, d: &Point2) -> Sign {
+    let row = |p: &Point2| {
+        let [x, y] = diffs2(p, d);
+        let lift = sq2(&x, &y);
+        [x, y, lift]
+    };
+    det3(&[row(a), row(b), row(c)]).sign()
+}
+
+// ---------------------------------------------------------------------------------------------
+// insphere
+// ---------------------------------------------------------------------------------------------
+
+/// Stage A of [`insphere`].
+#[inline]
+pub fn insphere_filter(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) -> Option<Sign> {
+    insphere_a(a, b, c, d, e).ok()
+}
+
+#[inline(always)]
+fn insphere_a(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) -> Result<Sign, f64> {
+    let ([ax, ay, az], [bx, by, bz], [cx, cy, cz], [dx, dy, dz], [ex, ey, ez]) = (a.0, b.0, c.0, d.0, e.0);
+    let (aex, bex, cex, dex) = (ax - ex, bx - ex, cx - ex, dx - ex);
+    let (aey, bey, cey, dey) = (ay - ey, by - ey, cy - ey, dy - ey);
+    let (aez, bez, cez, dez) = (az - ez, bz - ez, cz - ez, dz - ez);
+
+    let (aexbey, bexaey) = (aex * bey, bex * aey);
+    let ab = aexbey - bexaey;
+    let (bexcey, cexbey) = (bex * cey, cex * bey);
+    let bc = bexcey - cexbey;
+    let (cexdey, dexcey) = (cex * dey, dex * cey);
+    let cd = cexdey - dexcey;
+    let (dexaey, aexdey) = (dex * aey, aex * dey);
+    let da = dexaey - aexdey;
+    let (aexcey, cexaey) = (aex * cey, cex * aey);
+    let ac = aexcey - cexaey;
+    let (bexdey, dexbey) = (bex * dey, dex * bey);
+    let bd = bexdey - dexbey;
+
+    let abc = aez * bc - bez * ac + cez * ab;
+    let bcd = bez * cd - cez * bd + dez * bc;
+    let cda = cez * da + dez * ac + aez * cd;
+    let dab = dez * ab + aez * bd + bez * da;
+
+    let a_lift = aex * aex + aey * aey + aez * aez;
+    let b_lift = bex * bex + bey * bey + bez * bez;
+    let c_lift = cex * cex + cey * cey + cez * cez;
+    let d_lift = dex * dex + dey * dey + dez * dez;
+
+    let det = (d_lift * abc - c_lift * dab) + (b_lift * cda - a_lift * bcd);
+
+    let (aez_p, bez_p, cez_p, dez_p) = (aez.abs(), bez.abs(), cez.abs(), dez.abs());
+    let (ab_p, bc_p, cd_p) = (aexbey.abs() + bexaey.abs(), bexcey.abs() + cexbey.abs(), cexdey.abs() + dexcey.abs());
+    let (da_p, ac_p, bd_p) = (dexaey.abs() + aexdey.abs(), aexcey.abs() + cexaey.abs(), bexdey.abs() + dexbey.abs());
+    let permanent = (cd_p * bez_p + bd_p * cez_p + bc_p * dez_p) * a_lift
+        + (da_p * cez_p + ac_p * dez_p + cd_p * aez_p) * b_lift
+        + (ab_p * dez_p + bd_p * aez_p + da_p * bez_p) * c_lift
+        + (bc_p * aez_p + ac_p * bez_p + ab_p * cez_p) * d_lift;
+    if certain(det, ISP_A * permanent) { Ok(Sign::of(det)) } else { Err(permanent) }
+}
+
+/// Whether `e` lies inside the sphere through `a, b, c, d`: [`Sign::Positive`] inside,
+/// [`Sign::Negative`] outside, [`Sign::Zero`] on it, for `orient3d(a, b, c, d) > 0` (the sign is
+/// reversed otherwise). Exact. It is the sign of the 4×4 determinant with rows
+/// `[p - e, |p - e|²]` for `p = a, b, c, d`.
+#[inline]
+pub fn insphere(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) -> Sign {
+    match insphere_a(a, b, c, d, e) {
+        Ok(sign) => sign,
+        Err(permanent) => insphere_adapt(a, b, c, d, e, permanent),
+    }
+}
+
+#[cold]
+fn insphere_adapt(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3, permanent: f64) -> Sign {
+    let ([ax, ay, az], [bx, by, bz], [cx, cy, cz], [dx, dy, dz], [ex, ey, ez]) = (a.0, b.0, c.0, d.0, e.0);
+    let (aex, bex, cex, dex) = (ax - ex, bx - ex, cx - ex, dx - ex);
+    let (aey, bey, cey, dey) = (ay - ey, by - ey, cy - ey, dy - ey);
+    let (aez, bez, cez, dez) = (az - ez, bz - ez, cz - ez, dz - ez);
+
+    // Stage B: the 2×2 minors of the x, y columns, the 3×3 minors, then the lifted terms.
+    let ab = cross_diff(aex, bey, bex, aey);
+    let bc = cross_diff(bex, cey, cex, bey);
+    let cd = cross_diff(cex, dey, dex, cey);
+    let da = cross_diff(dex, aey, aex, dey);
+    let ac = cross_diff(aex, cey, cex, aey);
+    let bd = cross_diff(bex, dey, dex, bey);
+    let minor3 = |m1: &Exp<4>, z1: f64, m2: &Exp<4>, z2: f64, m3: &Exp<4>, z3: f64| -> Exp<24> {
+        Exp::sum(&Exp::<16>::sum(&Exp::<8>::scale(m1, z1), &Exp::<8>::scale(m2, z2)), &Exp::<8>::scale(m3, z3))
+    };
+    let abc = minor3(&bc, aez, &ac, -bez, &ab, cez);
+    let bcd = minor3(&cd, bez, &bd, -cez, &bc, dez);
+    let cda = minor3(&da, cez, &ac, dez, &cd, aez);
+    let dab = minor3(&ab, dez, &bd, aez, &da, bez);
+    let lifted = |m: &Exp<24>, x: f64, y: f64, z: f64| -> Exp<288> {
+        let xx = Exp::<96>::scale(&Exp::<48>::scale(m, x), x);
+        let yy = Exp::<96>::scale(&Exp::<48>::scale(m, y), y);
+        let zz = Exp::<96>::scale(&Exp::<48>::scale(m, z), z);
+        Exp::sum(&Exp::<192>::sum(&xx, &yy), &zz)
+    };
+    let d_term = lifted(&abc, dex, dey, dez);
+    let c_term = lifted(&dab, cex, cey, cez).negated();
+    let b_term = lifted(&cda, bex, bey, bez);
+    let a_term = lifted(&bcd, aex, aey, aez).negated();
+    let fin = Exp::<1152>::sum(&Exp::<576>::sum(&d_term, &c_term), &Exp::<576>::sum(&b_term, &a_term));
+    let det = fin.estimate();
+    if certain(det, ISP_B * permanent) {
+        return Sign::of(det);
+    }
+    let tails = [
+        two_diff_tail(ax, ex, aex),
+        two_diff_tail(bx, ex, bex),
+        two_diff_tail(cx, ex, cex),
+        two_diff_tail(dx, ex, dex),
+        two_diff_tail(ay, ey, aey),
+        two_diff_tail(by, ey, bey),
+        two_diff_tail(cy, ey, cey),
+        two_diff_tail(dy, ey, dey),
+        two_diff_tail(az, ez, aez),
+        two_diff_tail(bz, ez, bez),
+        two_diff_tail(cz, ez, cez),
+        two_diff_tail(dz, ez, dez),
+    ];
+    if tails.iter().all(|&t| t == 0.0) {
+        return fin.sign();
+    }
+    insphere_exact(a, b, c, d, e)
+}
+
+/// The same expansion as stage B (shared 2×2 minors, then the 3×3 minors, then the lifts), on
+/// the exact differences.
+fn insphere_exact(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) -> Sign {
+    let [ax, ay, az] = diffs3(a, e);
+    let [bx, by, bz] = diffs3(b, e);
+    let [cx, cy, cz] = diffs3(c, e);
+    let [dx, dy, dz] = diffs3(d, e);
+    let ab = det2(&ax, &ay, &bx, &by);
+    let bc = det2(&bx, &by, &cx, &cy);
+    let cd = det2(&cx, &cy, &dx, &dy);
+    let da = det2(&dx, &dy, &ax, &ay);
+    let ac = det2(&ax, &ay, &cx, &cy);
+    let bd = det2(&bx, &by, &dx, &dy);
+    let abc = az.clone() * bc.clone() - bz.clone() * ac.clone() + cz.clone() * ab.clone();
+    let bcd = bz.clone() * cd.clone() - cz.clone() * bd.clone() + dz.clone() * bc;
+    let cda = cz.clone() * da.clone() + dz.clone() * ac + az.clone() * cd;
+    let dab = dz.clone() * ab + az.clone() * bd + bz.clone() * da;
+    let lift = sq3(&dx, &dy, &dz) * abc - sq3(&cx, &cy, &cz) * dab + sq3(&bx, &by, &bz) * cda - sq3(&ax, &ay, &az) * bcd;
+    lift.sign()
+}
+
+/// The exact stages alone, for tests that exercise them directly.
+#[cfg(test)]
+pub(crate) mod exact {
+    use super::*;
+
+    pub(crate) fn orient2d(a: &Point2, b: &Point2, c: &Point2) -> Sign {
+        let [ax, ay] = diffs2(a, c);
+        let [bx, by] = diffs2(b, c);
+        det2(&ax, &ay, &bx, &by).sign()
+    }
+    pub(crate) fn orient3d(a: &Point3, b: &Point3, c: &Point3, d: &Point3) -> Sign {
+        orient3d_exact(a, b, c, d)
+    }
+    pub(crate) fn incircle(a: &Point2, b: &Point2, c: &Point2, d: &Point2) -> Sign {
+        incircle_exact(a, b, c, d)
+    }
+    pub(crate) fn insphere(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) -> Sign {
+        insphere_exact(a, b, c, d, e)
+    }
+}

crates/mesh/src/predicates/dyadic.rs

@@ -0,0 +1,253 @@
+//! Exact dyadic rationals `±m · 2^e` (`m` an arbitrary-precision natural number), the exact
+//! stage of the indirect predicates.
+//!
+//! Attene evaluates this stage with floating-point expansions (2020, §5.3). Expansions are exact
+//! only while no partial product underflows, and the polynomials of indirect predicates reach
+//! degree 27 (`orient3d` of four three-plane intersections), far beyond what the coordinate
+//! domain keeps out of underflow. Integer arithmetic has no such limit: every double is a dyadic
+//! rational, and sums and products of dyadic rationals are computed exactly. The stage only runs
+//! when the interval filter cannot decide, so its speed matters little.
+
+use std::cmp::Ordering;
+use std::ops::{Add, Mul, Neg, Sub};
+
+use super::Sign;
+use super::ring::{Ring, Signed};
+
+/// `(-1)^neg · mag · 2^exp`; `mag` holds little-endian base-2³² limbs with neither leading nor
+/// trailing zero limbs, and zero is the empty magnitude.
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub(crate) struct Dyadic {
+    neg: bool,
+    mag: Vec<u32>,
+    exp: i64,
+}
+
+impl Dyadic {
+    fn zero() -> Self {
+        Self { neg: false, mag: Vec::new(), exp: 0 }
+    }
+
+    fn normalized(neg: bool, mut mag: Vec<u32>, mut exp: i64) -> Self {
+        while mag.last() == Some(&0) {
+            mag.pop();
+        }
+        let low_zeros = mag.iter().take_while(|&&l| l == 0).count();
+        if low_zeros > 0 {
+            mag.drain(..low_zeros);
+            exp = exp.saturating_add(32 * low_zeros as i64);
+        }
+        if mag.is_empty() { Self::zero() } else { Self { neg, mag, exp } }
+    }
+
+    fn is_zero(&self) -> bool {
+        self.mag.is_empty()
+    }
+
+    /// The magnitude shifted left by `bits`.
+    fn shifted(&self, bits: u64) -> Vec<u32> {
+        let limbs = (bits / 32) as usize;
+        let bits = (bits % 32) as u32;
+        let mut out = vec![0u32; limbs];
+        out.reserve(self.mag.len() + 1);
+        if bits == 0 {
+            out.extend_from_slice(&self.mag);
+        } else {
+            let mut carry = 0u32;
+            for &l in &self.mag {
+                out.push((l << bits) | carry);
+                carry = l >> (32 - bits);
+            }
+            out.push(carry);
+        }
+        out
+    }
+}
+
+fn cmp_mag(a: &[u32], b: &[u32]) -> Ordering {
+    let (a, b) = (trimmed(a), trimmed(b));
+    a.len().cmp(&b.len()).then_with(|| a.iter().rev().cmp(b.iter().rev()))
+}
+
+fn trimmed(a: &[u32]) -> &[u32] {
+    let len = a.iter().rposition(|&l| l != 0).map_or(0, |i| i + 1);
+    a.get(..len).unwrap_or(a)
+}
+
+fn add_mag(a: &[u32], b: &[u32]) -> Vec<u32> {
+    let (long, short) = if a.len() >= b.len() { (a, b) } else { (b, a) };
+    let mut out = Vec::with_capacity(long.len() + 1);
+    let mut carry = 0u64;
+    for (i, &l) in long.iter().enumerate() {
+        let t = l as u64 + short.get(i).copied().unwrap_or(0) as u64 + carry;
+        out.push(t as u32);
+        carry = t >> 32;
+    }
+    out.push(carry as u32);
+    out
+}
+
+/// `a - b` for `a >= b`.
+fn sub_mag(a: &[u32], b: &[u32]) -> Vec<u32> {
+    let mut out = Vec::with_capacity(a.len());
+    let mut borrow = 0i64;
+    for (i, &l) in a.iter().enumerate() {
+        let mut t = l as i64 - b.get(i).copied().unwrap_or(0) as i64 - borrow;
+        borrow = if t < 0 {
+            t += 1 << 32;
+            1
+        } else {
+            0
+        };
+        out.push(t as u32);
+    }
+    out
+}
+
+fn mul_mag(a: &[u32], b: &[u32]) -> Vec<u32> {
+    let mut out = vec![0u32; a.len() + b.len()];
+    for (i, &x) in a.iter().enumerate() {
+        let mut carry = 0u64;
+        for (slot, &y) in out.iter_mut().skip(i).zip(b) {
+            // x·y + slot + carry ≤ (2³² − 1)² + 2(2³² − 1) = 2⁶⁴ − 1: no overflow.
+            let t = x as u64 * y as u64 + *slot as u64 + carry;
+            *slot = t as u32;
+            carry = t >> 32;
+        }
+        if let Some(slot) = out.get_mut(i + b.len()) {
+            *slot = carry as u32;
+        }
+    }
+    out
+}
+
+impl Ring for Dyadic {
+    fn from_f64(x: f64) -> Self {
+        let bits = x.to_bits();
+        let neg = bits >> 63 == 1;
+        let biased = ((bits >> 52) & 0x7ff) as i64;
+        let fraction = bits & ((1u64 << 52) - 1);
+        let (mant, exp) = match biased {
+            0 => (fraction, -1074),
+            // Not finite: unreachable for admissible coordinates; map to zero to stay total.
+            0x7ff => return Self::zero(),
+            _ => (fraction | (1u64 << 52), biased - 1075),
+        };
+        Self::normalized(neg, vec![mant as u32, (mant >> 32) as u32], exp)
+    }
+}
+
+impl Neg for Dyadic {
+    type Output = Self;
+
+    fn neg(mut self) -> Self {
+        if !self.is_zero() {
+            self.neg = !self.neg;
+        }
+        self
+    }
+}
+
+impl Add for Dyadic {
+    type Output = Self;
+
+    fn add(self, rhs: Self) -> Self {
+        if self.is_zero() {
+            return rhs;
+        }
+        if rhs.is_zero() {
+            return self;
+        }
+        let exp = self.exp.min(rhs.exp);
+        // Exponents stay far from i64 limits (bounded by polynomial degree × 1100), so these
+        // differences are small and non-negative.
+        let a = self.shifted(self.exp.saturating_sub(exp).unsigned_abs());
+        let b = rhs.shifted(rhs.exp.saturating_sub(exp).unsigned_abs());
+        if self.neg == rhs.neg {
+            return Self::normalized(self.neg, add_mag(&a, &b), exp);
+        }
+        match cmp_mag(&a, &b) {
+            Ordering::Equal => Self::zero(),
+            Ordering::Greater => Self::normalized(self.neg, sub_mag(&a, &b), exp),
+            Ordering::Less => Self::normalized(rhs.neg, sub_mag(&b, &a), exp),
+        }
+    }
+}
+
+impl Sub for Dyadic {
+    type Output = Self;
+
+    fn sub(self, rhs: Self) -> Self {
+        self + (-rhs)
+    }
+}
+
+impl Mul for Dyadic {
+    type Output = Self;
+
+    fn mul(self, rhs: Self) -> Self {
+        if self.is_zero() || rhs.is_zero() {
+            return Self::zero();
+        }
+        Self::normalized(self.neg != rhs.neg, mul_mag(&self.mag, &rhs.mag), self.exp.saturating_add(rhs.exp))
+    }
+}
+
+impl Signed for Dyadic {
+    fn sign(&self) -> Option<Sign> {
+        Some(if self.is_zero() {
+            Sign::Zero
+        } else if self.neg {
+            Sign::Negative
+        } else {
+            Sign::Positive
+        })
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use num_bigint::BigInt;
+    use num_rational::BigRational;
+    use num_traits::{FromPrimitive, One};
+
+    use super::*;
+
+    fn to_rational(d: &Dyadic) -> BigRational {
+        let mut m = BigInt::from(0);
+        for &l in d.mag.iter().rev() {
+            m = (m << 32) + BigInt::from(l);
+        }
+        if d.neg {
+            m = -m;
+        }
+        let two = BigRational::from_integer(BigInt::from(2));
+        let r = BigRational::from_integer(m);
+        let mut scale = BigRational::one();
+        for _ in 0..d.exp.unsigned_abs() {
+            scale *= two.clone();
+        }
+        if d.exp >= 0 { r * scale } else { r / scale }
+    }
+
+    fn q(x: f64) -> BigRational {
+        BigRational::from_f64(x).unwrap()
+    }
+
+    #[test]
+    fn matches_big_rationals() {
+        let xs = [0.0, -0.0, 1.0, -3.5, 1e300, -1e-300, 5e-324, 0.1, f64::MAX, 7.0e15 + 1.0, -2.0f64.powi(-1000)];
+        for &a in &xs {
+            assert_eq!(to_rational(&Dyadic::from_f64(a)), q(a), "{a}");
+            for &b in &xs {
+                let (da, db) = (Dyadic::from_f64(a), Dyadic::from_f64(b));
+                assert_eq!(to_rational(&(da.clone() + db.clone())), q(a) + q(b), "{a} + {b}");
+                assert_eq!(to_rational(&(da.clone() - db.clone())), q(a) - q(b), "{a} - {b}");
+                assert_eq!(to_rational(&(da * db)), q(a) * q(b), "{a} * {b}");
+            }
+        }
+        let third = Dyadic::from_f64(1.0 / 3.0);
+        let s = third.clone() + third.clone() + third - Dyadic::from_f64(1.0);
+        assert_eq!(s.sign(), Some(Sign::Negative)); // fl(1/3) < 1/3
+    }
+}

crates/mesh/src/predicates/eft.rs

@@ -0,0 +1,96 @@
+//! Error-free transformations (Shewchuk 1997, §2.3 and §2.5).
+//!
+//! Each function returns `(x, y)` with `x = fl(a ∘ b)` and `y` the exact roundoff, so that
+//! `a ∘ b = x + y` exactly. They assume IEEE 754 binary64 with round-to-nearest-even, no
+//! intermediate extended precision and no fused multiply-add contraction; that is what Rust's
+//! `f64` guarantees on every target ModelCraft builds for (x86_64 with SSE2, aarch64, wasm32):
+//! rustc never contracts `a * b + c` into an FMA and does not use x87. No overflow or underflow
+//! may occur, which the coordinate domain of [`super::Point2`] / [`super::Point3`] guarantees.
+
+/// `2^ceil(53 / 2) + 1`, the splitter for Dekker's product (§2.5, Theorem 17).
+const SPLITTER: f64 = 134_217_729.0;
+
+/// Fast-Two-Sum (Theorem 6). Requires `|a| >= |b|` (or `a == 0`).
+#[inline(always)]
+pub(crate) fn fast_two_sum(a: f64, b: f64) -> (f64, f64) {
+    let x = a + b;
+    let b_virtual = x - a;
+    (x, b - b_virtual)
+}
+
+/// Two-Sum (Theorem 7), no precondition on the magnitudes.
+#[inline(always)]
+pub(crate) fn two_sum(a: f64, b: f64) -> (f64, f64) {
+    let x = a + b;
+    let b_virtual = x - a;
+    let a_virtual = x - b_virtual;
+    let b_round = b - b_virtual;
+    let a_round = a - a_virtual;
+    (x, a_round + b_round)
+}
+
+/// The roundoff of `x = fl(a - b)`, given `x`.
+#[inline(always)]
+pub(crate) fn two_diff_tail(a: f64, b: f64, x: f64) -> f64 {
+    let b_virtual = a - x;
+    let a_virtual = x + b_virtual;
+    let b_round = b_virtual - b;
+    let a_round = a - a_virtual;
+    a_round + b_round
+}
+
+/// Two-Diff: `a - b = x + y` exactly.
+#[inline(always)]
+pub(crate) fn two_diff(a: f64, b: f64) -> (f64, f64) {
+    let x = a - b;
+    (x, two_diff_tail(a, b, x))
+}
+
+/// Dekker's split (Theorem 17): `a = hi + lo` with each half fitting in 26 bits.
+#[inline(always)]
+pub(crate) fn split(a: f64) -> (f64, f64) {
+    let c = SPLITTER * a;
+    let a_big = c - a;
+    let hi = c - a_big;
+    (hi, a - hi)
+}
+
+/// Two-Product with `b` already split (Theorem 18).
+#[inline(always)]
+pub(crate) fn two_product_presplit(a: f64, b: f64, b_hi: f64, b_lo: f64) -> (f64, f64) {
+    let x = a * b;
+    let (a_hi, a_lo) = split(a);
+    let err1 = x - a_hi * b_hi;
+    let err2 = err1 - a_lo * b_hi;
+    let err3 = err2 - a_hi * b_lo;
+    (x, a_lo * b_lo - err3)
+}
+
+/// Two-Product (Theorem 18): `a * b = x + y` exactly.
+#[inline(always)]
+pub(crate) fn two_product(a: f64, b: f64) -> (f64, f64) {
+    let (b_hi, b_lo) = split(b);
+    two_product_presplit(a, b, b_hi, b_lo)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn transforms_are_exact_on_classic_cases() {
+        let (x, y) = two_sum(1.0, 1e-17);
+        assert_eq!((x, y), (1.0, 1e-17));
+        // 2^53 + 1 is a tie that rounds to the even 2^53.
+        let (x, y) = fast_two_sum(2f64.powi(53), 1.0);
+        assert_eq!((x, y), (2f64.powi(53), 1.0));
+        let (x, y) = two_diff(1.0, 1e-17);
+        assert_eq!((x, y), (1.0, -1e-17));
+        // (1 + 2^-30)^2 = 1 + 2^-29 + 2^-60: the last term is the roundoff.
+        let a = 1.0 + 2f64.powi(-30);
+        let (x, y) = two_product(a, a);
+        assert_eq!((x, y), (1.0 + 2f64.powi(-29), 2f64.powi(-60)));
+        let (hi, lo) = split(a);
+        assert_eq!(hi + lo, a);
+    }
+}

crates/mesh/src/predicates/expansion.rs

@@ -0,0 +1,316 @@
+//! Floating-point expansions (Shewchuk 1997, §2.2–§2.6).
+//!
+//! An expansion is a sum of doubles, stored in order of increasing magnitude, whose components
+//! are nonoverlapping (here always nonadjacent or strongly nonoverlapping, as the theorems that
+//! the algorithms below rely on require). Zero components are eliminated, so a zero expansion is
+//! the single component `0.0` and the sign of any expansion is the sign of its last component.
+//!
+//! The two algorithms (Fast-Expansion-Sum with zero elimination, Theorem 13; Scale-Expansion with
+//! zero elimination, Theorem 19) are written once against a sink. [`Exp`] is a fixed-capacity
+//! stack expansion for the adaptive stages, whose capacities are checked at compile time;
+//! [`Expansion`] is a heap expansion with exact `+`, `-` and `*` for the final exact stage.
+
+use std::ops::{Add, Mul, Neg, Sub};
+
+use super::Sign;
+use super::eft::{fast_two_sum, split, two_diff, two_product, two_product_presplit, two_sum};
+use super::ring::{Ring, Signed};
+
+/// Merges two expansions into one stream of components of nondecreasing magnitude.
+struct Merge<'a> {
+    e: &'a [f64],
+    f: &'a [f64],
+}
+
+impl Iterator for Merge<'_> {
+    type Item = f64;
+
+    #[inline(always)]
+    fn next(&mut self) -> Option<f64> {
+        match (self.e.split_first(), self.f.split_first()) {
+            (Some((&a, e_rest)), Some((&b, f_rest))) => {
+                // `(b > a) == (b > -a)` is true when |a| < |b| (Shewchuk's comparison).
+                if (b > a) == (b > -a) {
+                    self.e = e_rest;
+                    Some(a)
+                } else {
+                    self.f = f_rest;
+                    Some(b)
+                }
+            }
+            (Some((&a, rest)), None) => {
+                self.e = rest;
+                Some(a)
+            }
+            (None, Some((&b, rest))) => {
+                self.f = rest;
+                Some(b)
+            }
+            (None, None) => None,
+        }
+    }
+}
+
+/// Fast-Expansion-Sum with zero elimination: emits `e + f`, at most `e.len() + f.len()`
+/// components.
+#[inline]
+fn sum_into(e: &[f64], f: &[f64], mut emit: impl FnMut(f64)) {
+    let mut merged = Merge { e, f };
+    let Some(mut q) = merged.next() else {
+        emit(0.0);
+        return;
+    };
+    let mut emitted = false;
+    if let Some(x) = merged.next() {
+        let (q_new, h) = fast_two_sum(x, q);
+        q = q_new;
+        if h != 0.0 {
+            emit(h);
+            emitted = true;
+        }
+        for x in merged {
+            let (q_new, h) = two_sum(q, x);
+            q = q_new;
+            if h != 0.0 {
+                emit(h);
+                emitted = true;
+            }
+        }
+    }
+    if q != 0.0 || !emitted {
+        emit(q);
+    }
+}
+
+/// Scale-Expansion with zero elimination: emits `e * b`, at most `2 * e.len()` components.
+#[inline]
+fn scale_into(e: &[f64], b: f64, mut emit: impl FnMut(f64)) {
+    let (b_hi, b_lo) = split(b);
+    let mut components = e.iter().copied();
+    let Some(first) = components.next() else {
+        emit(0.0);
+        return;
+    };
+    let mut emitted = false;
+    let (mut q, h) = two_product_presplit(first, b, b_hi, b_lo);
+    if h != 0.0 {
+        emit(h);
+        emitted = true;
+    }
+    for component in components {
+        let (p1, p0) = two_product_presplit(component, b, b_hi, b_lo);
+        let (s, h) = two_sum(q, p0);
+        if h != 0.0 {
+            emit(h);
+            emitted = true;
+        }
+        let (q_new, h) = fast_two_sum(p1, s);
+        q = q_new;
+        if h != 0.0 {
+            emit(h);
+            emitted = true;
+        }
+    }
+    if q != 0.0 || !emitted {
+        emit(q);
+    }
+}
+
+/// The sign of a zero-eliminated expansion: the sign of its most significant component.
+fn sign_of(components: &[f64]) -> Sign {
+    components.last().map_or(Sign::Zero, |&x| Sign::of(x))
+}
+
+/// A stack expansion of at most `N` components. Constructors that combine expansions check at
+/// compile time that `N` holds every component the operation can produce, so `push` never runs
+/// out of room.
+#[derive(Clone, Copy)]
+pub(crate) struct Exp<const N: usize> {
+    c: [f64; N],
+    n: usize,
+}
+
+impl<const N: usize> Exp<N> {
+    #[inline(always)]
+    fn empty() -> Self {
+        Self { c: [0.0; N], n: 0 }
+    }
+
+    #[inline(always)]
+    fn push(&mut self, x: f64) {
+        // In bounds by the compile-time capacity checks of the constructors.
+        if let Some(slot) = self.c.get_mut(self.n) {
+            *slot = x;
+            self.n += 1;
+        }
+    }
+
+    #[inline(always)]
+    pub(crate) fn as_slice(&self) -> &[f64] {
+        self.c.get(..self.n).unwrap_or(&[])
+    }
+
+    /// The expansion `hi + lo` produced by an error-free transformation (`lo` the roundoff).
+    #[inline(always)]
+    pub(crate) fn from_pair((hi, lo): (f64, f64)) -> Self {
+        const { assert!(N >= 2) };
+        let mut out = Self::empty();
+        if lo != 0.0 {
+            out.push(lo);
+        }
+        if hi != 0.0 || out.n == 0 {
+            out.push(hi);
+        }
+        out
+    }
+
+    /// `e + f`.
+    #[inline]
+    pub(crate) fn sum<const A: usize, const B: usize>(e: &Exp<A>, f: &Exp<B>) -> Self {
+        const { assert!(N >= A + B) };
+        let mut out = Self::empty();
+        sum_into(e.as_slice(), f.as_slice(), |x| out.push(x));
+        out
+    }
+
+    /// `e * b`.
+    #[inline]
+    pub(crate) fn scale<const A: usize>(e: &Exp<A>, b: f64) -> Self {
+        const { assert!(N >= 2 * A) };
+        let mut out = Self::empty();
+        scale_into(e.as_slice(), b, |x| out.push(x));
+        out
+    }
+
+    /// `-e`, exactly.
+    #[inline]
+    pub(crate) fn negated(mut self) -> Self {
+        if let Some(live) = self.c.get_mut(..self.n) {
+            live.iter_mut().for_each(|x| *x = -*x);
+        }
+        self
+    }
+
+    /// An approximation of the value (Shewchuk's `estimate`): the sum of the components.
+    #[inline]
+    pub(crate) fn estimate(&self) -> f64 {
+        self.as_slice().iter().sum()
+    }
+
+    pub(crate) fn sign(&self) -> Sign {
+        sign_of(self.as_slice())
+    }
+}
+
+/// The exact two-by-two product difference `a * b - c * d` as a four-component expansion
+/// (Shewchuk's Two-Product followed by Two-Two-Diff).
+#[inline]
+pub(crate) fn cross_diff(a: f64, b: f64, c: f64, d: f64) -> Exp<4> {
+    let ab = Exp::<2>::from_pair(two_product(a, b));
+    let cd = Exp::<2>::from_pair(two_product(c, d)).negated();
+    Exp::sum(&ab, &cd)
+}
+
+/// A heap expansion with exact ring operations, for the final exact stage of a predicate.
+#[derive(Clone, Debug)]
+pub(crate) struct Expansion(Vec<f64>);
+
+impl Expansion {
+    /// The exact difference `a - b`.
+    pub(crate) fn diff(a: f64, b: f64) -> Self {
+        Self(Exp::<2>::from_pair(two_diff(a, b)).as_slice().to_vec())
+    }
+
+    fn scaled(&self, b: f64) -> Self {
+        let mut v = Vec::with_capacity(2 * self.0.len());
+        scale_into(&self.0, b, |x| v.push(x));
+        Self(v)
+    }
+
+    pub(crate) fn sign(&self) -> Sign {
+        sign_of(&self.0)
+    }
+}
+
+impl Add for Expansion {
+    type Output = Self;
+
+    fn add(self, rhs: Self) -> Self {
+        let mut v = Vec::with_capacity(self.0.len() + rhs.0.len());
+        sum_into(&self.0, &rhs.0, |x| v.push(x));
+        Self(v)
+    }
+}
+
+impl Neg for Expansion {
+    type Output = Self;
+
+    fn neg(mut self) -> Self {
+        for x in self.0.iter_mut() {
+            *x = -*x;
+        }
+        self
+    }
+}
+
+impl Sub for Expansion {
+    type Output = Self;
+
+    fn sub(self, rhs: Self) -> Self {
+        self + (-rhs)
+    }
+}
+
+impl Mul for Expansion {
+    type Output = Self;
+
+    /// The sum of `self` scaled by each component of `rhs` (distributivity; every partial
+    /// product is exact).
+    #[allow(clippy::suspicious_arithmetic_impl)] // the product is a sum of scaled expansions
+    fn mul(self, rhs: Self) -> Self {
+        let (long, short) = if self.0.len() >= rhs.0.len() { (self, rhs) } else { (rhs, self) };
+        let mut terms = short.0.iter().map(|&b| long.scaled(b));
+        let first = terms.next().unwrap_or_else(|| Self(vec![0.0]));
+        terms.fold(first, |acc, term| acc + term)
+    }
+}
+
+impl Ring for Expansion {
+    fn from_f64(x: f64) -> Self {
+        Self(vec![x])
+    }
+}
+
+impl Signed for Expansion {
+    fn sign(&self) -> Option<Sign> {
+        Some(sign_of(&self.0))
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn sums_and_products_are_exact() {
+        let big = Expansion::from_f64(1e300);
+        let tiny = Expansion::from_f64(1e-300);
+        let s = big.clone() + tiny.clone();
+        assert_eq!(s.0, vec![1e-300, 1e300]);
+        assert_eq!((s - big).0, vec![1e-300]);
+        let a = Expansion::diff(1.0, 1e-20); // 1 - 1e-20, two components
+        let sq = a.clone() * a;
+        // (1 - x)^2 = 1 - 2x + x^2 with x = 1e-20 (as doubles: the x^2 term is exact too).
+        let expect = Expansion::from_f64(1.0) - Expansion::from_f64(2e-20) + Expansion::from_f64(1e-20) * Expansion::from_f64(1e-20);
+        assert_eq!((sq - expect).sign(), Sign::Zero);
+    }
+
+    #[test]
+    fn stack_and_heap_agree() {
+        let e = cross_diff(3.0, 1e17, 1.0, 3e17 - 64.0);
+        assert_eq!(e.estimate(), 64.0);
+        let s = Exp::<8>::scale(&e, 1.0 + f64::EPSILON);
+        assert_eq!(s.sign(), Sign::Positive);
+        assert_eq!(Exp::<4>::from_pair((0.0, 0.0)).as_slice(), &[0.0]);
+    }
+}

crates/mesh/src/predicates/indirect.rs

@@ -0,0 +1,305 @@
+//! Indirect predicates on implicit points (Attene 2020, §4–§5), from the paper only.
+//!
+//! An implicit point is kept as the explicit points that define it and is never rounded: a
+//! line–plane intersection ([`Lpi`], §4.2) or a three-plane intersection ([`Tpi`]). Each has
+//! homogeneous coordinates `(X, W)` that are polynomials in the defining coordinates, so a
+//! predicate on implicit points is the sign of a polynomial times the signs of the `W`s
+//! (Attene §4.3–§4.6, §5.1: "count how many of them have a negative sign").
+//!
+//! Evaluation follows Attene §5: the interval filter (§5.2, with the defining points' `(X, W)`
+//! cached as intervals, §5.4), then exact evaluation. Two deliberate differences: there is no
+//! semi-static floating-point stage yet (it needs a per-formula error bound; the interval stage
+//! is the dynamic equivalent), and the exact stage uses [`super::dyadic::Dyadic`] integers
+//! instead of expansions, which could underflow at these polynomial degrees. All-explicit
+//! inputs go straight to Shewchuk's predicates.
+//!
+//! A construction is accepted only when its `W` is exactly nonzero (the line is not parallel to
+//! the plane, the three planes meet in one point), so the predicates are total.
+
+use std::cmp::Ordering;
+
+use super::dyadic::Dyadic;
+use super::interval::Interval;
+use super::ring::{Ring, Signed, det2, det3};
+use super::{Point2, Point3, Sign, orient2d, orient3d};
+use crate::{Error, Result};
+
+/// A coordinate axis: the axis compared by [`generic_compare`], or dropped by
+/// [`generic_orient2d`].
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
+pub enum Axis {
+    X,
+    Y,
+    Z,
+}
+
+impl Axis {
+    pub(crate) fn index(self) -> usize {
+        match self {
+            Axis::X => 0,
+            Axis::Y => 1,
+            Axis::Z => 2,
+        }
+    }
+
+    /// The two remaining axes in cyclic order (`X → (y, z)`, `Y → (z, x)`, `Z → (x, y)`), so
+    /// the projection is seen from the positive side of the dropped axis.
+    pub(crate) fn projected(self) -> (usize, usize) {
+        match self {
+            Axis::X => (1, 2),
+            Axis::Y => (2, 0),
+            Axis::Z => (0, 1),
+        }
+    }
+}
+
+/// Homogeneous coordinates `X / W`; `w: None` means `W = 1` (an explicit point), which spares
+/// the filter needless multiplications.
+#[derive(Debug, Clone)]
+struct Hom<T> {
+    x: [T; 3],
+    w: Option<T>,
+}
+
+fn lift<T: Ring>(p: &[f64; 3]) -> [T; 3] {
+    p.map(T::from_f64)
+}
+
+fn sub3<T: Ring>(a: &[T; 3], b: &[T; 3]) -> [T; 3] {
+    let ([a0, a1, a2], [b0, b1, b2]) = (a, b);
+    [a0.clone() - b0.clone(), a1.clone() - b1.clone(), a2.clone() - b2.clone()]
+}
+
+fn coord<T: Clone>(v: &[T; 3], k: usize) -> T {
+    let [x, y, z] = v;
+    match k {
+        0 => x.clone(),
+        1 => y.clone(),
+        _ => z.clone(),
+    }
+}
+
+/// `v · w`, or `v` when `w = 1`.
+fn mulw<T: Ring>(v: T, w: &Option<T>) -> T {
+    match w {
+        None => v,
+        Some(w) => v * w.clone(),
+    }
+}
+
+/// The sign of `W` (positive for explicit points), `None` if uncertain.
+fn w_sign<T: Signed>(w: &Option<T>) -> Option<Sign> {
+    w.as_ref().map_or(Some(Sign::Positive), Signed::sign)
+}
+
+/// Line `q1 q2` meets plane `r s t` (Attene §4.2): with `d = det[q1 - q2; s - r; t - r]` and
+/// `n = det[q1 - r; s - r; t - r]`, the point is `(d·q1 + n·(q2 - q1)) / d`.
+fn lpi_hom<T: Ring>([q1, q2, r, s, t]: &[[f64; 3]; 5]) -> Hom<T> {
+    let (q1, q2, r, s, t) = (lift::<T>(q1), lift::<T>(q2), lift::<T>(r), lift::<T>(s), lift::<T>(t));
+    let (sr, tr) = (sub3(&s, &r), sub3(&t, &r));
+    let d = det3(&[sub3(&q1, &q2), sr.clone(), tr.clone()]);
+    let n = det3(&[sub3(&q1, &r), sr, tr]);
+    let dir = sub3(&q2, &q1);
+    let x = std::array::from_fn(|k| d.clone() * coord(&q1, k) + n.clone() * coord(&dir, k));
+    Hom { x, w: Some(d) }
+}
+
+/// Three planes `(r, s, t)` meet (Cramer's rule): normals `N_i = (s_i - r_i) × (t_i - r_i)`,
+/// offsets `c_i = N_i · r_i`, `W = det[N_1; N_2; N_3]`, `X_k` = the same determinant with
+/// column `k` replaced by `c`.
+fn tpi_hom<T: Ring>(planes: &[[[f64; 3]; 3]; 3]) -> Hom<T> {
+    let plane = |[r, s, t]: &[[f64; 3]; 3]| -> ([T; 3], T) {
+        let r = lift::<T>(r);
+        let u = sub3(&lift::<T>(s), &r);
+        let v = sub3(&lift::<T>(t), &r);
+        let ([ux, uy, uz], [vx, vy, vz]) = (&u, &v);
+        let n = [det2(uy, uz, vy, vz), det2(uz, ux, vz, vx), det2(ux, uy, vx, vy)];
+        let [nx, ny, nz] = &n;
+        let [rx, ry, rz] = &r;
+        let c = nx.clone() * rx.clone() + ny.clone() * ry.clone() + nz.clone() * rz.clone();
+        (n, c)
+    };
+    let [p1, p2, p3] = planes;
+    let [(n1, c1), (n2, c2), (n3, c3)] = [plane(p1), plane(p2), plane(p3)];
+    // Row `n` with column `k` replaced by `c`.
+    let with = |n: &[T; 3], c: &T, k: usize| -> [T; 3] { std::array::from_fn(|j| if j == k { c.clone() } else { coord(n, j) }) };
+    let x = std::array::from_fn(|k| det3(&[with(&n1, &c1, k), with(&n2, &c2, k), with(&n3, &c3, k)]));
+    Hom { x, w: Some(det3(&[n1, n2, n3])) }
+}
+
+/// Decides the sign of `W` exactly: filter first, integers if needed.
+fn exact_w_sign(filter: &Hom<Interval>, exact: impl FnOnce() -> Hom<Dyadic>) -> Sign {
+    w_sign(&filter.w).or_else(|| w_sign(&exact().w)).unwrap_or(Sign::Zero)
+}
+
+/// The intersection of the line through `q1, q2` with the plane through `r, s, t`.
+#[derive(Debug, Clone)]
+pub struct Lpi {
+    points: [[f64; 3]; 5],
+    filter: Hom<Interval>,
+}
+
+impl Lpi {
+    /// Fails when the point is undefined: the line is parallel to the plane, `q1 = q2`, or
+    /// `r, s, t` are collinear (decided exactly).
+    pub fn new(q1: &Point3, q2: &Point3, r: &Point3, s: &Point3, t: &Point3) -> Result<Self> {
+        let points = [q1.0, q2.0, r.0, s.0, t.0];
+        let filter = lpi_hom::<Interval>(&points);
+        if exact_w_sign(&filter, || lpi_hom::<Dyadic>(&points)).is_zero() {
+            return Err(Error::Invalid("line-plane intersection is undefined (parallel or degenerate)".into()));
+        }
+        Ok(Self { points, filter })
+    }
+}
+
+/// The intersection point of three planes, each through three points.
+#[derive(Debug, Clone)]
+pub struct Tpi {
+    planes: [[[f64; 3]; 3]; 3],
+    filter: Hom<Interval>,
+}
+
+impl Tpi {
+    /// Fails when the planes do not meet in exactly one point, or a plane's points are collinear
+    /// (decided exactly).
+    pub fn new(planes: [[&Point3; 3]; 3]) -> Result<Self> {
+        let planes = planes.map(|p| p.map(|q| q.0));
+        let filter = tpi_hom::<Interval>(&planes);
+        if exact_w_sign(&filter, || tpi_hom::<Dyadic>(&planes)).is_zero() {
+            return Err(Error::Invalid("three-plane intersection is undefined (planes not independent)".into()));
+        }
+        Ok(Self { planes, filter })
+    }
+}
+
+/// An explicit point or an implicit one, for the `generic_*` predicates.
+#[derive(Debug, Clone)]
+pub enum GenericPoint {
+    Explicit(Point3),
+    Lpi(Lpi),
+    Tpi(Tpi),
+}
+
+impl From<Point3> for GenericPoint {
+    fn from(p: Point3) -> Self {
+        GenericPoint::Explicit(p)
+    }
+}
+
+impl GenericPoint {
+    fn explicit(&self) -> Option<&Point3> {
+        match self {
+            GenericPoint::Explicit(p) => Some(p),
+            _ => None,
+        }
+    }
+
+    fn filter_hom(&self) -> Hom<Interval> {
+        match self {
+            GenericPoint::Explicit(p) => Hom { x: lift(&p.0), w: None },
+            GenericPoint::Lpi(l) => l.filter.clone(),
+            GenericPoint::Tpi(t) => t.filter.clone(),
+        }
+    }
+
+    fn exact_hom(&self) -> Hom<Dyadic> {
+        match self {
+            GenericPoint::Explicit(p) => Hom { x: lift(&p.0), w: None },
+            GenericPoint::Lpi(l) => lpi_hom(&l.points),
+            GenericPoint::Tpi(t) => tpi_hom(&t.planes),
+        }
+    }
+
+    /// Approximate Cartesian coordinates, for output and display only: never base a decision on
+    /// them.
+    pub fn approx(&self) -> [f64; 3] {
+        match self {
+            GenericPoint::Explicit(p) => p.0,
+            _ => {
+                let h = self.filter_hom();
+                let w = h.w.map_or(1.0, |w| w.mid());
+                h.x.map(|x| x.mid() / w)
+            }
+        }
+    }
+}
+
+/// Runs a formula through the interval filter, then exactly.
+fn decide<const N: usize>(
+    points: [&GenericPoint; N],
+    filter: impl FnOnce([Hom<Interval>; N]) -> Option<Sign>,
+    exact: impl FnOnce([Hom<Dyadic>; N]) -> Option<Sign>,
+) -> Sign {
+    filter(points.map(GenericPoint::filter_hom)).or_else(|| exact(points.map(GenericPoint::exact_hom))).unwrap_or(Sign::Zero)
+}
+
+/// Moves an explicit point (if any) to the last position, the pivot of the formulas, and
+/// reports whether that took a transposition (which flips an orientation's sign).
+fn pivot_last<const N: usize>(mut points: [&GenericPoint; N]) -> ([&GenericPoint; N], bool) {
+    let last = N.saturating_sub(1);
+    match points.iter().position(|p| p.explicit().is_some()) {
+        Some(k) if k != last => {
+            points.swap(k, last);
+            (points, true)
+        }
+        _ => (points, false),
+    }
+}
+
+/// `det[a - d; b - d; c - d] = det[R_a; R_b; R_c] / (W_a W_b W_c W_d³)` with
+/// `R_i = W_d X_i - W_i X_d`.
+fn orient3d_formula<T: Ring + Signed>([a, b, c, d]: [Hom<T>; 4]) -> Option<Sign> {
+    let row = |p: &Hom<T>| -> [T; 3] { std::array::from_fn(|k| mulw(coord(&p.x, k), &d.w) - mulw(coord(&d.x, k), &p.w)) };
+    let det = det3(&[row(&a), row(&b), row(&c)]).sign()?;
+    Some(det * w_sign(&a.w)? * w_sign(&b.w)? * w_sign(&c.w)? * w_sign(&d.w)?)
+}
+
+/// [`orient3d`] on explicit or implicit points, same sign convention. Exact.
+pub fn generic_orient3d(a: &GenericPoint, b: &GenericPoint, c: &GenericPoint, d: &GenericPoint) -> Sign {
+    if let (Some(a), Some(b), Some(c), Some(d)) = (a.explicit(), b.explicit(), c.explicit(), d.explicit()) {
+        return orient3d(a, b, c, d);
+    }
+    let (points, flipped) = pivot_last([a, b, c, d]);
+    let sign = decide(points, orient3d_formula, orient3d_formula);
+    if flipped { -sign } else { sign }
+}
+
+/// `det[a' - c'; b' - c']` of the projections: `det[R_a; R_b] / (W_a W_b W_c²)`.
+fn orient2d_formula<T: Ring + Signed>([a, b, c]: [Hom<T>; 3], (u, v): (usize, usize)) -> Option<Sign> {
+    let row = |p: &Hom<T>, k: usize| mulw(coord(&p.x, k), &c.w) - mulw(coord(&c.x, k), &p.w);
+    let det = det2(&row(&a, u), &row(&a, v), &row(&b, u), &row(&b, v)).sign()?;
+    Some(det * w_sign(&a.w)? * w_sign(&b.w)?)
+}
+
+/// [`orient2d`] of the projections of `a, b, c` onto the coordinate plane that drops `drop`,
+/// seen from the positive side of that axis (`Z` → the `xy` plane, `X` → `yz`, `Y` → `zx`).
+/// Exact. Attene's `orient2d3d` (§4.5) is this with the axis of the largest normal component.
+pub fn generic_orient2d(a: &GenericPoint, b: &GenericPoint, c: &GenericPoint, drop: Axis) -> Sign {
+    let (u, v) = drop.projected();
+    if let (Some(a), Some(b), Some(c)) = (a.explicit(), b.explicit(), c.explicit()) {
+        let proj = |p: &Point3| Point2([coord(&p.0, u), coord(&p.0, v)]);
+        return orient2d(&proj(a), &proj(b), &proj(c));
+    }
+    let (points, flipped) = pivot_last([a, b, c]);
+    let sign = decide(points, |h| orient2d_formula(h, (u, v)), |h| orient2d_formula(h, (u, v)));
+    if flipped { -sign } else { sign }
+}
+
+/// `X_a / W_a - X_b / W_b = (W_b X_a - W_a X_b) / (W_a W_b)` along one axis.
+fn compare_formula<T: Ring + Signed>([a, b]: [Hom<T>; 2], k: usize) -> Option<Sign> {
+    let diff = mulw(coord(&a.x, k), &b.w) - mulw(coord(&b.x, k), &a.w);
+    Some(diff.sign()? * w_sign(&a.w)? * w_sign(&b.w)?)
+}
+
+/// Compares the `axis` coordinates of two points exactly (Attene's `lessThan`).
+pub fn generic_compare(a: &GenericPoint, b: &GenericPoint, axis: Axis) -> Ordering {
+    let k = axis.index();
+    if let (Some(a), Some(b)) = (a.explicit(), b.explicit()) {
+        return coord(&a.0, k).partial_cmp(&coord(&b.0, k)).unwrap_or(Ordering::Equal);
+    }
+    match decide([a, b], |h| compare_formula(h, k), |h| compare_formula(h, k)) {
+        Sign::Negative => Ordering::Less,
+        Sign::Zero => Ordering::Equal,
+        Sign::Positive => Ordering::Greater,
+    }
+}

crates/mesh/src/predicates/interval.rs

@@ -0,0 +1,106 @@
+//! Interval arithmetic, the dynamic filter of the indirect predicates (Attene 2020, §5.2;
+//! Brönnimann, Burnikel & Pion 1998).
+//!
+//! Attene switches the FPU to directed rounding; Rust has no portable control of the rounding
+//! mode (and wasm none at all), so each bound is computed with the default round-to-nearest and
+//! then moved outward by one unit in the last place (`next_down` / `next_up`). A round-to-nearest
+//! result is within half an ulp of the exact value, so the widened interval always contains it,
+//! subnormals included. Overflow gives infinite bounds and `inf - inf` gives NaN; both make the
+//! sign uncertain, which sends the predicate on to the exact stage.
+
+use std::ops::{Add, Mul, Sub};
+
+use super::Sign;
+use super::ring::{Ring, Signed};
+
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub(crate) struct Interval {
+    lo: f64,
+    hi: f64,
+}
+
+impl Interval {
+    const NAN: Interval = Interval { lo: f64::NAN, hi: f64::NAN };
+
+    #[inline(always)]
+    fn widened(lo: f64, hi: f64) -> Self {
+        if lo.is_nan() || hi.is_nan() { Self::NAN } else { Self { lo: lo.next_down(), hi: hi.next_up() } }
+    }
+
+    /// The midpoint, an approximation of the value.
+    pub(crate) fn mid(&self) -> f64 {
+        0.5 * self.lo + 0.5 * self.hi
+    }
+}
+
+impl Add for Interval {
+    type Output = Self;
+
+    #[inline(always)]
+    fn add(self, rhs: Self) -> Self {
+        Self::widened(self.lo + rhs.lo, self.hi + rhs.hi)
+    }
+}
+
+impl Sub for Interval {
+    type Output = Self;
+
+    #[inline(always)]
+    fn sub(self, rhs: Self) -> Self {
+        Self::widened(self.lo - rhs.hi, self.hi - rhs.lo)
+    }
+}
+
+impl Mul for Interval {
+    type Output = Self;
+
+    #[inline(always)]
+    fn mul(self, rhs: Self) -> Self {
+        let p = [self.lo * rhs.lo, self.lo * rhs.hi, self.hi * rhs.lo, self.hi * rhs.hi];
+        // `f64::min`/`max` skip NaN (from `0 * inf`), which would drop a bound: check first.
+        if p.iter().any(|x| x.is_nan()) {
+            return Self::NAN;
+        }
+        let lo = p.iter().copied().fold(f64::INFINITY, f64::min);
+        let hi = p.iter().copied().fold(f64::NEG_INFINITY, f64::max);
+        Self::widened(lo, hi)
+    }
+}
+
+impl Ring for Interval {
+    #[inline(always)]
+    fn from_f64(x: f64) -> Self {
+        Self { lo: x, hi: x }
+    }
+}
+
+impl Signed for Interval {
+    #[inline(always)]
+    fn sign(&self) -> Option<Sign> {
+        if self.lo > 0.0 {
+            Some(Sign::Positive)
+        } else if self.hi < 0.0 {
+            Some(Sign::Negative)
+        } else if self.lo == 0.0 && self.hi == 0.0 {
+            Some(Sign::Zero)
+        } else {
+            None
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn intervals_enclose_and_flag_uncertainty() {
+        let third = Interval::from_f64(1.0) * Interval::from_f64(1.0 / 3.0);
+        let x = third + third + third - Interval::from_f64(1.0);
+        assert!(x.lo <= 0.0 && x.hi >= 0.0);
+        assert_eq!(x.sign(), None);
+        let big = Interval::from_f64(f64::MAX) * Interval::from_f64(2.0);
+        assert_eq!((big - big).sign(), None);
+        assert_eq!((Interval::from_f64(2.0) * Interval::from_f64(3.0)).sign(), Some(Sign::Positive));
+    }
+}

crates/mesh/src/predicates/mod.rs

@@ -0,0 +1,169 @@
+//! Exact geometric predicates: the sign of `orient2d`, `orient3d`, `incircle` and `insphere`
+//! for every admissible input, symbolic perturbation for the Delaunay tie-breaks, and indirect
+//! predicates on points constructed as line–plane and three-plane intersections.
+//!
+//! # Sources
+//! - J. R. Shewchuk, *Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric
+//!   Predicates*, DCG 18:305–363, 1997: expansions (§2), adaptive evaluation (§3), the four
+//!   predicates and their error bounds (§4). [`orient2d`], [`orient3d`], [`incircle`],
+//!   [`insphere`] and their `*_filter` stage-A variants.
+//! - O. Devillers, M. Teillaud, *Perturbations for Delaunay and weighted Delaunay 3D
+//!   triangulations*, CGTA 44(3):160–168, 2011; H. Edelsbrunner, E. P. Mücke, *Simulation of
+//!   Simplicity*, ACM TOG 9(1), 1990: [`insphere_sos`], [`incircle_sos`].
+//! - M. Attene, *Indirect Predicates for Geometric Constructions*, CAD 126:102856, 2020 (paper
+//!   only; his code is LGPL and is not read): [`Lpi`], [`Tpi`], [`GenericPoint`] and the
+//!   `generic_*` predicates.
+//!
+//! # Arithmetic assumptions
+//! IEEE 754 binary64, round-to-nearest-even, no extended intermediate precision, no FMA
+//! contraction. Rust's `f64` gives exactly that on x86_64 (SSE2), aarch64 and
+//! wasm32-unknown-unknown: rustc never fuses `a * b + c` and never uses x87 on these targets, and
+//! the code uses no platform intrinsics, so every predicate returns the same sign on every
+//! platform. The results are deterministic functions of the inputs.
+//!
+//! # Domain
+//! Shewchuk's error bounds and the exactness of expansion arithmetic need no overflow and no
+//! underflow (§5). Every coordinate of a [`Point2`] / [`Point3`] is therefore either zero or has
+//! magnitude in [[`MIN_MAGNITUDE`], [`MAX_MAGNITUDE`]] = [2⁻¹²⁰, 2¹⁶⁰]: then every intermediate
+//! value of the degree-5 `insphere` is either zero or above 2⁻⁸⁶⁰ and below 2⁸¹², all normal
+//! doubles. The constructors reject anything else (NaN, infinities, out-of-range magnitudes) with
+//! [`crate::Error::Invalid`]; the predicates themselves are total and never fail. Callers snap
+//! tessellation noise (e.g. `1e-40`) to zero before building points.
+
+mod adaptive;
+mod dyadic;
+mod eft;
+mod expansion;
+mod indirect;
+mod interval;
+mod ring;
+mod sos;
+
+pub use adaptive::{incircle, incircle_filter, insphere, insphere_filter, orient2d, orient2d_filter, orient3d, orient3d_filter};
+pub use indirect::{Axis, GenericPoint, Lpi, Tpi, generic_compare, generic_orient2d, generic_orient3d};
+pub use sos::{incircle_sos, insphere_sos};
+
+use crate::{Error, Result};
+
+/// The smallest nonzero coordinate magnitude a point may have: 2⁻¹²⁰ ≈ 7.5e-37.
+pub const MIN_MAGNITUDE: f64 = f64::from_bits((1023 - 120) << 52);
+/// The largest coordinate magnitude a point may have: 2¹⁶⁰ ≈ 1.5e48.
+pub const MAX_MAGNITUDE: f64 = f64::from_bits((1023 + 160) << 52);
+
+/// The sign of a predicate.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, PartialOrd, Ord)]
+pub enum Sign {
+    Negative = -1,
+    Zero = 0,
+    Positive = 1,
+}
+
+impl Sign {
+    /// The sign of a double (`Zero` for both zeros; NaN, which the domain excludes, maps to
+    /// `Zero`).
+    #[inline(always)]
+    pub fn of(x: f64) -> Sign {
+        if x > 0.0 {
+            Sign::Positive
+        } else if x < 0.0 {
+            Sign::Negative
+        } else {
+            Sign::Zero
+        }
+    }
+
+    #[inline(always)]
+    pub fn is_zero(self) -> bool {
+        self == Sign::Zero
+    }
+}
+
+impl std::ops::Neg for Sign {
+    type Output = Sign;
+
+    #[inline(always)]
+    fn neg(self) -> Sign {
+        match self {
+            Sign::Negative => Sign::Positive,
+            Sign::Zero => Sign::Zero,
+            Sign::Positive => Sign::Negative,
+        }
+    }
+}
+
+impl std::ops::Mul for Sign {
+    type Output = Sign;
+
+    #[inline(always)]
+    fn mul(self, rhs: Sign) -> Sign {
+        match (self, rhs) {
+            (Sign::Zero, _) | (_, Sign::Zero) => Sign::Zero,
+            (a, b) if a == b => Sign::Positive,
+            _ => Sign::Negative,
+        }
+    }
+}
+
+/// Whether `x` may be a point coordinate (see the module docs).
+#[inline]
+pub fn is_admissible(x: f64) -> bool {
+    x == 0.0 || (x.abs() >= MIN_MAGNITUDE && x.abs() <= MAX_MAGNITUDE)
+}
+
+fn check<const N: usize>(c: [f64; N]) -> Result<[f64; N]> {
+    match c.iter().find(|x| !is_admissible(**x)) {
+        None => Ok(c),
+        Some(x) => {
+            Err(Error::Invalid(format!("point coordinate {x:e} is outside the exact-predicate domain (zero, or magnitude in [2^-120, 2^160])")))
+        }
+    }
+}
+
+/// A 2D point whose coordinates are admissible for the exact predicates.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Point2([f64; 2]);
+
+impl Point2 {
+    pub fn new(x: f64, y: f64) -> Result<Self> {
+        check([x, y]).map(Self)
+    }
+
+    #[inline(always)]
+    pub fn coords(&self) -> [f64; 2] {
+        self.0
+    }
+}
+
+impl TryFrom<[f64; 2]> for Point2 {
+    type Error = Error;
+
+    fn try_from(c: [f64; 2]) -> Result<Self> {
+        check(c).map(Self)
+    }
+}
+
+/// A 3D point whose coordinates are admissible for the exact predicates.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Point3([f64; 3]);
+
+impl Point3 {
+    pub fn new(x: f64, y: f64, z: f64) -> Result<Self> {
+        check([x, y, z]).map(Self)
+    }
+
+    #[inline(always)]
+    pub fn coords(&self) -> [f64; 3] {
+        self.0
+    }
+}
+
+impl TryFrom<[f64; 3]> for Point3 {
+    type Error = Error;
+
+    fn try_from(c: [f64; 3]) -> Result<Self> {
+        check(c).map(Self)
+    }
+}
+
+#[cfg(test)]
+mod tests;

crates/mesh/src/predicates/ring.rs

@@ -0,0 +1,46 @@
+//! Predicate formulas written once over exact or certified number types.
+//!
+//! A [`Ring`] is a number type with `+`, `-`, `*`; a [`Signed`] one can report a sign it is sure
+//! of. The formulas below are instantiated with [`super::expansion::Expansion`] (the exact stage
+//! of the explicit predicates), [`super::interval::Interval`] (the dynamic filter of the indirect
+//! predicates) and [`super::dyadic::Dyadic`] (their exact stage). Writing each polynomial once
+//! and letting the compiler monomorphize it plays the role of the code generators of Lévy's PCK
+//! (2016) and Attene (2020, §6): one source formula, several arithmetic models.
+
+use std::ops::{Add, Mul, Sub};
+
+use super::Sign;
+
+pub(crate) trait Ring: Clone + Add<Output = Self> + Sub<Output = Self> + Mul<Output = Self> {
+    /// The exact value of a double.
+    fn from_f64(x: f64) -> Self;
+}
+
+pub(crate) trait Signed {
+    /// The sign, when the representation is sure of it (always, for exact types).
+    fn sign(&self) -> Option<Sign>;
+}
+
+/// `a * d - b * c`.
+#[inline]
+pub(crate) fn det2<T: Ring>(a: &T, b: &T, c: &T, d: &T) -> T {
+    a.clone() * d.clone() - b.clone() * c.clone()
+}
+
+/// The determinant of the three rows, expanded along the last column.
+#[inline]
+pub(crate) fn det3<T: Ring>(rows: &[[T; 3]; 3]) -> T {
+    let [[a0, a1, a2], [b0, b1, b2], [c0, c1, c2]] = rows;
+    a2.clone() * det2(b0, b1, c0, c1) - b2.clone() * det2(a0, a1, c0, c1) + c2.clone() * det2(a0, a1, b0, b1)
+}
+
+/// `x² + y²` and `x² + y² + z²`.
+#[inline]
+pub(crate) fn sq2<T: Ring>(x: &T, y: &T) -> T {
+    x.clone() * x.clone() + y.clone() * y.clone()
+}
+
+#[inline]
+pub(crate) fn sq3<T: Ring>(x: &T, y: &T, z: &T) -> T {
+    sq2(x, y) + z.clone() * z.clone()
+}

crates/mesh/src/predicates/sos.rs

@@ -0,0 +1,94 @@
+//! Symbolic perturbation of the in-sphere and in-circle tests (Devillers & Teillaud 2011, after
+//! Edelsbrunner & Mücke's Simulation of Simplicity).
+//!
+//! `insphere(a, b, c, d, e)` is the 5×5 determinant with rows `[p, |p|², 1]` (translate by `e`
+//! and expand along the last column to see it). Lift every point by an infinitesimal
+//! `|p|² + ε^k(p)`, where `k(p)` is the rank of `p` in decreasing lexicographic order
+//! (`(x, y, z)`; the largest point gets the largest perturbation `ε¹`). The determinant is linear
+//! in the lift column, so the perturbed value is
+//!
+//! `insphere + Σ_p ε^k(p) · C_p`,  `C_p = (-1)^(row(p) + 4) · orient3d(the other four points)`
+//!
+//! (the cofactor of `p`'s lift is a determinant with rows `[q, 1]`, which is `orient3d` of the
+//! other points in their order). Its sign is the sign of the first nonzero term in order of
+//! decreasing rank. Lifting by these weights is a regular triangulation of the same points, so
+//! the decisions are consistent across calls, and a Bowyer–Watson cavity stays a star-shaped
+//! ball. When `a, b, c, d` is a proper tetrahedron the term of `e` (`-orient3d(a, b, c, d)`) is
+//! nonzero, so the result is never zero.
+//!
+//! The 2D version replaces `orient3d` by `orient2d` and the 5×5 by the 4×4 determinant.
+
+use std::cmp::Ordering;
+
+use super::{Point2, Point3, Sign, incircle, insphere, orient2d, orient3d};
+
+/// Lexicographic order on coordinates (total: `f64::total_cmp`).
+fn lex<const N: usize>(p: &[f64; N], q: &[f64; N]) -> Ordering {
+    p.iter().zip(q).map(|(a, b)| a.total_cmp(b)).find(|o| o.is_ne()).unwrap_or(Ordering::Equal)
+}
+
+/// Argument positions sorted by decreasing lexicographic order of their points, ties (equal
+/// points, which a triangulation never passes) by position.
+fn ranked<const N: usize, const M: usize>(points: &[[f64; N]; M]) -> [usize; M] {
+    let mut order: [usize; M] = std::array::from_fn(|i| i);
+    order.sort_by(|&i, &j| match (points.get(i), points.get(j)) {
+        (Some(p), Some(q)) => lex(q, p).then(i.cmp(&j)),
+        _ => i.cmp(&j),
+    });
+    order
+}
+
+/// [`insphere`] with symbolic perturbation: never [`Sign::Zero`] when `a, b, c, d` is not
+/// coplanar. Points must be pairwise distinct for the perturbation to be consistent (a
+/// Delaunay triangulation never passes duplicates); with all five points coplanar it returns
+/// `Zero`.
+#[inline]
+pub fn insphere_sos(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) -> Sign {
+    let exact = insphere(a, b, c, d, e);
+    if exact.is_zero() { insphere_perturbed(a, b, c, d, e) } else { exact }
+}
+
+#[cold]
+fn insphere_perturbed(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) -> Sign {
+    for i in ranked(&[a.0, b.0, c.0, d.0, e.0]) {
+        let [p, q, r, s] = match i {
+            0 => [b, c, d, e],
+            1 => [a, c, d, e],
+            2 => [a, b, d, e],
+            3 => [a, b, c, e],
+            _ => [a, b, c, d],
+        };
+        let minor = orient3d(p, q, r, s);
+        if !minor.is_zero() {
+            // Row i + 1 (1-based), lift column 4: cofactor sign (-1)^(i + 5).
+            return if i % 2 == 1 { minor } else { -minor };
+        }
+    }
+    Sign::Zero
+}
+
+/// [`incircle`] with symbolic perturbation: never [`Sign::Zero`] when `a, b, c` is not
+/// collinear. Same conventions as [`insphere_sos`].
+#[inline]
+pub fn incircle_sos(a: &Point2, b: &Point2, c: &Point2, d: &Point2) -> Sign {
+    let exact = incircle(a, b, c, d);
+    if exact.is_zero() { incircle_perturbed(a, b, c, d) } else { exact }
+}
+
+#[cold]
+fn incircle_perturbed(a: &Point2, b: &Point2, c: &Point2, d: &Point2) -> Sign {
+    for i in ranked(&[a.0, b.0, c.0, d.0]) {
+        let [p, q, r] = match i {
+            0 => [b, c, d],
+            1 => [a, c, d],
+            2 => [a, b, d],
+            _ => [a, b, c],
+        };
+        let minor = orient2d(p, q, r);
+        if !minor.is_zero() {
+            // Row i + 1 (1-based), lift column 3: cofactor sign (-1)^(i + 4).
+            return if i % 2 == 0 { minor } else { -minor };
+        }
+    }
+    Sign::Zero
+}

crates/mesh/src/predicates/tests.rs

@@ -0,0 +1,607 @@
+//! Differential tests against independent oracles: exact rationals (`num-rational`, Leibniz
+//! determinants written here from the definitions), `i128` on integer grids, and georust's
+//! `robust` (a port of Shewchuk's public-domain `predicates.c`).
+
+use std::cmp::Ordering;
+
+use num_bigint::BigInt;
+use num_rational::BigRational;
+use num_traits::{FromPrimitive, Signed as _, ToPrimitive, Zero};
+
+use super::adaptive::exact;
+use super::*;
+
+// ---------------------------------------------------------------------------------------------
+// Oracles and generators
+// ---------------------------------------------------------------------------------------------
+
+fn q(x: f64) -> BigRational {
+    BigRational::from_f64(x).unwrap()
+}
+
+fn rsign(x: &BigRational) -> Sign {
+    if x.is_zero() {
+        Sign::Zero
+    } else if x.is_positive() {
+        Sign::Positive
+    } else {
+        Sign::Negative
+    }
+}
+
+/// Laplace expansion along the first row.
+fn rdet(m: &[Vec<BigRational>]) -> BigRational {
+    if m.len() == 1 {
+        return m[0][0].clone();
+    }
+    let mut total = BigRational::zero();
+    for j in 0..m.len() {
+        let minor: Vec<Vec<BigRational>> =
+            m[1..].iter().map(|row| row.iter().enumerate().filter(|(k, _)| *k != j).map(|(_, v)| v.clone()).collect()).collect();
+        let term = m[0][j].clone() * rdet(&minor);
+        if j % 2 == 0 { total += term } else { total -= term }
+    }
+    total
+}
+
+fn rdiff(p: &[f64], o: &[f64]) -> Vec<BigRational> {
+    p.iter().zip(o).map(|(a, b)| q(*a) - q(*b)).collect()
+}
+
+fn lifted(mut row: Vec<BigRational>) -> Vec<BigRational> {
+    let lift = row.iter().fold(BigRational::zero(), |acc, v| acc + v.clone() * v.clone());
+    row.push(lift);
+    row
+}
+
+fn o2(a: &Point2, b: &Point2, c: &Point2) -> Sign {
+    rsign(&rdet(&[rdiff(&a.0, &c.0), rdiff(&b.0, &c.0)]))
+}
+fn o3(a: &Point3, b: &Point3, c: &Point3, d: &Point3) -> Sign {
+    rsign(&rdet(&[rdiff(&a.0, &d.0), rdiff(&b.0, &d.0), rdiff(&c.0, &d.0)]))
+}
+fn ic(a: &Point2, b: &Point2, c: &Point2, d: &Point2) -> Sign {
+    rsign(&rdet(&[lifted(rdiff(&a.0, &d.0)), lifted(rdiff(&b.0, &d.0)), lifted(rdiff(&c.0, &d.0))]))
+}
+fn is(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) -> Sign {
+    rsign(&rdet(&[lifted(rdiff(&a.0, &e.0)), lifted(rdiff(&b.0, &e.0)), lifted(rdiff(&c.0, &e.0)), lifted(rdiff(&d.0, &e.0))]))
+}
+
+/// SplitMix64: a deterministic generator for reproducible tests.
+struct Rng(u64);
+
+impl Rng {
+    fn next_u64(&mut self) -> u64 {
+        self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
+        let mut z = self.0;
+        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
+        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
+        z ^ (z >> 31)
+    }
+    /// Uniform in [-1, 1).
+    fn unit(&mut self) -> f64 {
+        (self.next_u64() >> 11) as f64 / (1u64 << 52) as f64 - 1.0
+    }
+    fn int(&mut self, n: i64) -> i64 {
+        (self.next_u64() % (2 * n as u64 + 1)) as i64 - n
+    }
+    /// `x` moved by up to `k` ulps either way.
+    fn nudge(&mut self, mut x: f64, k: i64) -> f64 {
+        let steps = self.int(k);
+        for _ in 0..steps.abs() {
+            x = if steps > 0 { x.next_up() } else { x.next_down() };
+        }
+        x
+    }
+}
+
+fn p2(x: f64, y: f64) -> Point2 {
+    Point2::new(x, y).unwrap()
+}
+fn p3(x: f64, y: f64, z: f64) -> Point3 {
+    Point3::new(x, y, z).unwrap()
+}
+
+fn rc2(p: &Point2) -> robust::Coord<f64> {
+    robust::Coord { x: p.0[0], y: p.0[1] }
+}
+fn rc3(p: &Point3) -> robust::Coord3D<f64> {
+    robust::Coord3D { x: p.0[0], y: p.0[1], z: p.0[2] }
+}
+
+/// Every stage agrees with every oracle, and the filter, when it answers, agrees too.
+fn check2(a: &Point2, b: &Point2, c: &Point2, d: &Point2) {
+    let want = o2(a, b, c);
+    assert_eq!(orient2d(a, b, c), want, "orient2d {a:?} {b:?} {c:?}");
+    assert_eq!(exact::orient2d(a, b, c), want);
+    assert_eq!(Sign::of(robust::orient2d(rc2(a), rc2(b), rc2(c))), want);
+    assert!(orient2d_filter(a, b, c).is_none_or(|s| s == want));
+    let want = ic(a, b, c, d);
+    assert_eq!(incircle(a, b, c, d), want, "incircle {a:?} {b:?} {c:?} {d:?}");
+    assert_eq!(exact::incircle(a, b, c, d), want);
+    assert_eq!(Sign::of(robust::incircle(rc2(a), rc2(b), rc2(c), rc2(d))), want);
+    assert!(incircle_filter(a, b, c, d).is_none_or(|s| s == want));
+}
+
+fn check3(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) {
+    let want = o3(a, b, c, d);
+    assert_eq!(orient3d(a, b, c, d), want, "orient3d {a:?} {b:?} {c:?} {d:?}");
+    assert_eq!(exact::orient3d(a, b, c, d), want);
+    assert_eq!(Sign::of(robust::orient3d(rc3(a), rc3(b), rc3(c), rc3(d))), want);
+    assert!(orient3d_filter(a, b, c, d).is_none_or(|s| s == want));
+    let want = is(a, b, c, d, e);
+    assert_eq!(insphere(a, b, c, d, e), want, "insphere {a:?} {b:?} {c:?} {d:?} {e:?}");
+    assert_eq!(exact::insphere(a, b, c, d, e), want);
+    assert_eq!(Sign::of(robust::insphere(rc3(a), rc3(b), rc3(c), rc3(d), rc3(e))), want);
+    assert!(insphere_filter(a, b, c, d, e).is_none_or(|s| s == want));
+}
+
+// ---------------------------------------------------------------------------------------------
+// Explicit predicates
+// ---------------------------------------------------------------------------------------------
+
+#[test]
+fn random_points_match_the_oracles() {
+    let mut r = Rng(1);
+    for _ in 0..300 {
+        let pt2 = |r: &mut Rng| p2(r.unit(), r.unit());
+        let (a, b, c, d) = (pt2(&mut r), pt2(&mut r), pt2(&mut r), pt2(&mut r));
+        check2(&a, &b, &c, &d);
+        let pt3 = |r: &mut Rng| p3(r.unit(), r.unit(), r.unit());
+        let (a, b, c, d, e) = (pt3(&mut r), pt3(&mut r), pt3(&mut r), pt3(&mut r), pt3(&mut r));
+        check3(&a, &b, &c, &d, &e);
+    }
+}
+
+#[test]
+fn near_degenerate_points_match_the_oracles() {
+    // Points on a line, plane, circle or sphere, rounded and then moved by a few ulps: the
+    // filter fails and every adaptive stage is exercised.
+    let mut r = Rng(2);
+    for i in 0..400 {
+        let k = if i % 4 == 0 { 0 } else { 2 };
+        let scale = [1.0, 1e-3, 1e5, 3.0e-20][i % 4];
+        let (ax, ay, bx, by) = (r.unit() * scale, r.unit() * scale, r.unit() * scale, r.unit() * scale);
+        let t = r.unit();
+        let c = p2(r.nudge(ax + t * (bx - ax), k), r.nudge(ay + t * (by - ay), k));
+        let (a, b) = (p2(ax, ay), p2(bx, by));
+        let on_circle = |r: &mut Rng| {
+            let th = r.unit() * 3.2;
+            p2(r.nudge(scale * th.cos() + scale * 0.25, k), r.nudge(scale * th.sin() - scale * 0.5, k))
+        };
+        let (ca, cb, cc, cd) = (on_circle(&mut r), on_circle(&mut r), on_circle(&mut r), on_circle(&mut r));
+        check2(&a, &b, &c, &c);
+        check2(&ca, &cb, &cc, &cd);
+
+        let pt = |r: &mut Rng| [r.unit() * scale, r.unit() * scale, r.unit() * scale];
+        let (a, b, c) = (pt(&mut r), pt(&mut r), pt(&mut r));
+        let (s, t) = (r.unit(), r.unit());
+        let d: Vec<f64> = (0..3).map(|j| r.nudge(a[j] + s * (b[j] - a[j]) + t * (c[j] - a[j]), k)).collect();
+        let on_sphere = |r: &mut Rng| {
+            let (th, ph) = (r.unit() * 3.2, r.unit() * 1.6);
+            p3(r.nudge(scale * ph.cos() * th.cos() + scale, k), r.nudge(scale * ph.cos() * th.sin(), k), r.nudge(scale * ph.sin() - scale * 0.5, k))
+        };
+        let (a, b, c, d) = (p3(a[0], a[1], a[2]), p3(b[0], b[1], b[2]), p3(c[0], c[1], c[2]), p3(d[0], d[1], d[2]));
+        let s: Vec<Point3> = (0..5).map(|_| on_sphere(&mut r)).collect();
+        check3(&a, &b, &c, &d, &s[0]);
+        check3(&s[0], &s[1], &s[2], &s[3], &s[4]);
+    }
+}
+
+#[test]
+fn ulp_grid_near_a_line_matches_the_oracle() {
+    // The classic failure of naive orient2d (Kettner et al., "Classroom examples of robustness
+    // problems"): p on a 64×64 grid of ulps around (0.5, 0.5), q = (12, 12), r = (24, 24).
+    // Variants with points off the diagonal make the coordinate differences inexact.
+    let ulp = 2f64.powi(-53);
+    for (q2, r2) in [(p2(12.0, 12.0), p2(24.0, 24.0)), (p2(0.1, 0.3), p2(17.3, 0.7)), (p2(1e-9, 3.0), p2(0.9, 1e-7))] {
+        for i in 0..64 {
+            for j in 0..64 {
+                let p = p2(0.5 + i as f64 * ulp, 0.5 + j as f64 * ulp);
+                let want = o2(&p, &q2, &r2);
+                assert_eq!(orient2d(&p, &q2, &r2), want);
+                assert_eq!(orient2d(&q2, &r2, &p), want);
+                assert_eq!(orient2d(&r2, &p, &q2), want);
+            }
+        }
+    }
+}
+
+fn i128_det(m: &[Vec<i128>]) -> i128 {
+    if m.len() == 1 {
+        return m[0][0];
+    }
+    (0..m.len())
+        .map(|j| {
+            let minor: Vec<Vec<i128>> = m[1..].iter().map(|row| row.iter().enumerate().filter(|(k, _)| *k != j).map(|(_, v)| *v).collect()).collect();
+            let t = m[0][j] * i128_det(&minor);
+            if j % 2 == 0 { t } else { -t }
+        })
+        .sum()
+}
+
+#[test]
+fn integer_grids_match_i128() {
+    // Small grids are full of exact degeneracies (collinear, coplanar, cocircular,
+    // cospherical); large coordinates (2^30 for orient, 2^14 for in-sphere) stress the filter.
+    let mut r = Rng(3);
+    for &(n, big) in &[(2i64, 1i64), (3, 1 << 30), (3, 1 << 14)] {
+        for _ in 0..400 {
+            let v: Vec<[i64; 3]> = (0..5).map(|_| [r.int(n) * big, r.int(n) * big, r.int(n) * big]).collect();
+            let pt = |p: &[i64; 3]| p3(p[0] as f64, p[1] as f64, p[2] as f64);
+            let d = |p: &[i64; 3], o: &[i64; 3]| (0..3).map(|k| (p[k] - o[k]) as i128).collect::<Vec<_>>();
+            let o3 = i128_det(&[d(&v[0], &v[3]), d(&v[1], &v[3]), d(&v[2], &v[3])]);
+            assert_eq!(orient3d(&pt(&v[0]), &pt(&v[1]), &pt(&v[2]), &pt(&v[3])), Sign::of(o3 as f64));
+            let o2v = |p: &[i64; 3]| p2(p[0] as f64, p[1] as f64);
+            let d2 = |p: &[i64; 3], o: &[i64; 3]| d(p, o)[..2].to_vec();
+            let want = i128_det(&[d2(&v[0], &v[2]), d2(&v[1], &v[2])]);
+            assert_eq!(orient2d(&o2v(&v[0]), &o2v(&v[1]), &o2v(&v[2])), Sign::of(want as f64));
+            if big <= 1 << 14 || n == 2 {
+                let lift = |mut row: Vec<i128>| {
+                    row.push(row.iter().map(|x| x * x).sum());
+                    row
+                };
+                let s = i128_det(&[lift(d(&v[0], &v[4])), lift(d(&v[1], &v[4])), lift(d(&v[2], &v[4])), lift(d(&v[3], &v[4]))]);
+                assert_eq!(insphere(&pt(&v[0]), &pt(&v[1]), &pt(&v[2]), &pt(&v[3]), &pt(&v[4])), Sign::of(s as f64));
+                let c = i128_det(&[lift(d2(&v[0], &v[3])), lift(d2(&v[1], &v[3])), lift(d2(&v[2], &v[3]))]);
+                assert_eq!(incircle(&o2v(&v[0]), &o2v(&v[1]), &o2v(&v[2]), &o2v(&v[3])), Sign::of(c as f64));
+            }
+        }
+    }
+}
+
+#[test]
+fn extreme_magnitudes_inside_the_domain_stay_exact() {
+    // Coordinates at both ends of the domain, mixed, near-degenerate: no overflow or underflow
+    // may corrupt a sign.
+    let mut r = Rng(4);
+    let tiny = MIN_MAGNITUDE;
+    let huge = MAX_MAGNITUDE;
+    for i in 0..150 {
+        let s = [tiny, huge, 1.0, tiny * 1e10][i % 4];
+        let pick = |r: &mut Rng| {
+            let m = r.unit();
+            if m.abs() < 0.2 { 0.0 } else { (m * s).clamp(-huge, huge) }
+        };
+        let fix = |x: f64| if is_admissible(x) { x } else { 0.0 };
+        let v: Vec<Point3> = (0..5)
+            .map(|_| {
+                let (x, y, z) = (pick(&mut r), pick(&mut r), pick(&mut r));
+                p3(fix(r.nudge(x, 1)), fix(y), fix(z))
+            })
+            .collect();
+        check3(&v[0], &v[1], &v[2], &v[3], &v[4]);
+        // Collinear and coplanar with huge and tiny offsets mixed.
+        let base = p3(huge / 2.0, tiny, 0.0);
+        let shifted = p3(huge / 2.0, tiny * 2.0, tiny);
+        check3(&base, &shifted, &v[2], &v[3], &v[4]);
+    }
+}
+
+#[test]
+fn permutations_follow_the_determinant() {
+    let mut r = Rng(5);
+    for _ in 0..200 {
+        // Grid points so that zero results also occur.
+        let pt = |r: &mut Rng| p3(r.int(2) as f64, r.int(2) as f64, r.int(2) as f64 * 0.5);
+        let (a, b, c, d, e) = (pt(&mut r), pt(&mut r), pt(&mut r), pt(&mut r), pt(&mut r));
+        let s = orient3d(&a, &b, &c, &d);
+        assert_eq!(orient3d(&b, &a, &c, &d), -s);
+        assert_eq!(orient3d(&a, &b, &d, &c), -s);
+        assert_eq!(orient3d(&b, &c, &a, &d), s);
+        assert_eq!(orient3d(&b, &c, &d, &a), -s); // a 4-cycle is odd
+        let s = insphere(&a, &b, &c, &d, &e);
+        assert_eq!(insphere(&b, &a, &c, &d, &e), -s);
+        assert_eq!(insphere(&a, &b, &c, &e, &d), -s);
+        assert_eq!(insphere(&c, &a, &b, &d, &e), s);
+        let (a2, b2, c2, d2) = (p2(a.0[0], a.0[1]), p2(b.0[0], b.0[1]), p2(c.0[0], c.0[1]), p2(d.0[0], d.0[2]));
+        let s = orient2d(&a2, &b2, &c2);
+        assert_eq!(orient2d(&b2, &a2, &c2), -s);
+        assert_eq!(orient2d(&b2, &c2, &a2), s);
+        let s = incircle(&a2, &b2, &c2, &d2);
+        assert_eq!(incircle(&b2, &a2, &c2, &d2), -s);
+        assert_eq!(incircle(&a2, &b2, &d2, &c2), -s);
+    }
+}
+
+#[test]
+fn conventions_match_the_documentation() {
+    let (o, x, y, z) = (p3(0.0, 0.0, 0.0), p3(1.0, 0.0, 0.0), p3(0.0, 1.0, 0.0), p3(0.0, 0.0, 1.0));
+    // o, x, y counterclockwise seen from +z; a point below gives Positive.
+    assert_eq!(orient3d(&o, &x, &y, &p3(0.0, 0.0, -1.0)), Sign::Positive);
+    assert_eq!(orient3d(&o, &x, &y, &z), Sign::Negative);
+    let tet = [o, x, y, p3(0.0, 0.0, -1.0)];
+    assert_eq!(insphere(&tet[0], &tet[1], &tet[2], &tet[3], &p3(0.1, 0.1, -0.1)), Sign::Positive);
+    assert_eq!(insphere(&tet[0], &tet[1], &tet[2], &tet[3], &p3(5.0, 5.0, 5.0)), Sign::Negative);
+    assert_eq!(insphere(&tet[0], &tet[1], &tet[2], &tet[3], &p3(1.0, 1.0, 0.0)), Sign::Zero);
+    let (a, b, c) = (p2(0.0, 0.0), p2(1.0, 0.0), p2(0.0, 1.0));
+    assert_eq!(orient2d(&a, &b, &c), Sign::Positive);
+    assert_eq!(incircle(&a, &b, &c, &p2(0.5, 0.5)), Sign::Positive);
+    assert_eq!(incircle(&a, &b, &c, &p2(1.0, 1.0)), Sign::Zero);
+}
+
+#[test]
+fn the_domain_is_enforced_at_construction() {
+    for bad in [f64::NAN, f64::INFINITY, -f64::INFINITY, 1e-40, -1e-300, 5e-324, 1e49, f64::MAX] {
+        assert!(Point3::new(0.0, bad, 1.0).is_err(), "{bad}");
+        assert!(Point2::try_from([bad, 0.0]).is_err(), "{bad}");
+    }
+    for good in [0.0, -0.0, MIN_MAGNITUDE, -MAX_MAGNITUDE, 1e-36, 1e48] {
+        assert!(Point3::new(good, good, good).is_ok(), "{good}");
+    }
+}
+
+#[test]
+fn hostile_bit_patterns_never_panic() {
+    let mut r = Rng(6);
+    let mut accepted = Vec::new();
+    for _ in 0..20_000 {
+        let c = [f64::from_bits(r.next_u64()), f64::from_bits(r.next_u64()), f64::from_bits(r.next_u64())];
+        if let Ok(p) = Point3::try_from(c) {
+            accepted.push(p);
+        }
+    }
+    accepted.push(p3(0.0, 0.0, 0.0));
+    for w in accepted.windows(5).take(200) {
+        let s = orient3d(&w[0], &w[1], &w[2], &w[3]);
+        assert_eq!(s, o3(&w[0], &w[1], &w[2], &w[3]));
+        let s = insphere(&w[0], &w[1], &w[2], &w[3], &w[4]);
+        assert_eq!(s, is(&w[0], &w[1], &w[2], &w[3], &w[4]));
+        let _ = insphere_sos(&w[0], &w[1], &w[2], &w[3], &w[4]);
+    }
+}
+
+// ---------------------------------------------------------------------------------------------
+// Symbolic perturbation
+// ---------------------------------------------------------------------------------------------
+
+/// The perturbed determinant evaluated with an actual tiny ε, lifts `|p|² + ε^rank` (rank 1 for
+/// the lexicographically largest point): an independent check of the cofactor signs.
+fn perturbed_sign(points: &[Vec<f64>]) -> Sign {
+    let mut order: Vec<usize> = (0..points.len()).collect();
+    order.sort_by(|&i, &j| points[j].partial_cmp(&points[i]).unwrap_or(Ordering::Equal));
+    let eps = BigRational::new(BigInt::from(1), BigInt::from(1u64 << 40));
+    let rows: Vec<Vec<BigRational>> = points
+        .iter()
+        .enumerate()
+        .map(|(i, p)| {
+            let rank = order.iter().position(|&k| k == i).unwrap() + 1;
+            let mut row = lifted(p.iter().map(|x| q(*x)).collect());
+            let mut e = BigRational::from_integer(BigInt::from(1));
+            for _ in 0..rank {
+                e *= eps.clone();
+            }
+            if let Some(lift) = row.last_mut() {
+                *lift += e;
+            }
+            row.push(BigRational::from_integer(BigInt::from(1)));
+            row
+        })
+        .collect();
+    rsign(&rdet(&rows))
+}
+
+#[test]
+fn perturbed_insphere_matches_an_infinitesimal_lift() {
+    let mut r = Rng(7);
+    let mut zeros_resolved = 0;
+    for _ in 0..600 {
+        let v: Vec<Vec<f64>> = (0..5).map(|_| vec![r.int(1) as f64, r.int(1) as f64, r.int(1) as f64]).collect();
+        if (0..5).any(|i| (0..i).any(|j| v[i] == v[j])) {
+            continue;
+        }
+        let p: Vec<Point3> = v.iter().map(|c| p3(c[0], c[1], c[2])).collect();
+        if orient3d(&p[0], &p[1], &p[2], &p[3]).is_zero() {
+            continue;
+        }
+        let got = insphere_sos(&p[0], &p[1], &p[2], &p[3], &p[4]);
+        assert_ne!(got, Sign::Zero);
+        assert_eq!(got, perturbed_sign(&v), "{v:?}");
+        zeros_resolved += insphere(&p[0], &p[1], &p[2], &p[3], &p[4]).is_zero() as usize;
+    }
+    assert!(zeros_resolved > 20, "the grid must produce cospherical cases ({zeros_resolved})");
+}
+
+#[test]
+fn perturbed_incircle_matches_an_infinitesimal_lift() {
+    let mut r = Rng(8);
+    let mut zeros_resolved = 0;
+    for _ in 0..600 {
+        let v: Vec<Vec<f64>> = (0..4).map(|_| vec![r.int(1) as f64, r.int(1) as f64]).collect();
+        if (0..4).any(|i| (0..i).any(|j| v[i] == v[j])) {
+            continue;
+        }
+        let p: Vec<Point2> = v.iter().map(|c| p2(c[0], c[1])).collect();
+        if orient2d(&p[0], &p[1], &p[2]).is_zero() {
+            continue;
+        }
+        let got = incircle_sos(&p[0], &p[1], &p[2], &p[3]);
+        assert_ne!(got, Sign::Zero);
+        assert_eq!(got, perturbed_sign(&v), "{v:?}");
+        zeros_resolved += incircle(&p[0], &p[1], &p[2], &p[3]).is_zero() as usize;
+    }
+    assert!(zeros_resolved > 20, "{zeros_resolved}");
+}
+
+// ---------------------------------------------------------------------------------------------
+// Indirect predicates
+// ---------------------------------------------------------------------------------------------
+
+type RPoint = [BigRational; 3];
+
+fn rp(p: &Point3) -> RPoint {
+    p.0.map(q)
+}
+fn rsub(a: &RPoint, b: &RPoint) -> RPoint {
+    std::array::from_fn(|k| a[k].clone() - b[k].clone())
+}
+fn rdot(a: &RPoint, b: &RPoint) -> BigRational {
+    (0..3).fold(BigRational::zero(), |acc, k| acc + a[k].clone() * b[k].clone())
+}
+fn rcross(a: &RPoint, b: &RPoint) -> RPoint {
+    [
+        a[1].clone() * b[2].clone() - a[2].clone() * b[1].clone(),
+        a[2].clone() * b[0].clone() - a[0].clone() * b[2].clone(),
+        a[0].clone() * b[1].clone() - a[1].clone() * b[0].clone(),
+    ]
+}
+
+/// Line–plane intersection by the parametric form (independent of the determinant formula).
+fn r_lpi(q1: &Point3, q2: &Point3, r: &Point3, s: &Point3, t: &Point3) -> RPoint {
+    let (q1, q2, r, s, t) = (rp(q1), rp(q2), rp(r), rp(s), rp(t));
+    let n = rcross(&rsub(&s, &r), &rsub(&t, &r));
+    let dir = rsub(&q2, &q1);
+    let u = rdot(&n, &rsub(&r, &q1)) / rdot(&n, &dir);
+    std::array::from_fn(|k| q1[k].clone() + u.clone() * dir[k].clone())
+}
+
+/// Three-plane intersection by Gaussian elimination.
+fn r_tpi(planes: &[[Point3; 3]; 3]) -> RPoint {
+    let mut m: Vec<Vec<BigRational>> = planes
+        .iter()
+        .map(|[r, s, t]| {
+            let (r, s, t) = (rp(r), rp(s), rp(t));
+            let n = rcross(&rsub(&s, &r), &rsub(&t, &r));
+            let c = rdot(&n, &r);
+            vec![n[0].clone(), n[1].clone(), n[2].clone(), c]
+        })
+        .collect();
+    for col in 0..3 {
+        let piv = (col..3).find(|&i| !m[i][col].is_zero()).unwrap();
+        m.swap(col, piv);
+        for i in 0..3 {
+            if i != col {
+                let f = m[i][col].clone() / m[col][col].clone();
+                let pivot_row = m[col].clone();
+                for (cell, p) in m[i].iter_mut().zip(&pivot_row) {
+                    *cell -= p.clone() * f.clone();
+                }
+            }
+        }
+    }
+    std::array::from_fn(|k| m[k][3].clone() / m[k][k].clone())
+}
+
+fn r_orient3d(p: &[RPoint; 4]) -> Sign {
+    let rows: Vec<Vec<BigRational>> = (0..3).map(|i| rsub(&p[i], &p[3]).to_vec()).collect();
+    rsign(&rdet(&rows))
+}
+
+fn rand_pt(r: &mut Rng, scale: f64) -> Point3 {
+    p3(r.unit() * scale, r.unit() * scale, r.unit() * scale)
+}
+
+#[test]
+fn indirect_orient3d_matches_rational_constructions() {
+    let mut r = Rng(9);
+    let mut checked = 0;
+    for i in 0..150 {
+        let pts: Vec<Point3> = (0..14).map(|_| rand_pt(&mut r, 1.0)).collect();
+        let Ok(lpi) = Lpi::new(&pts[0], &pts[1], &pts[2], &pts[3], &pts[4]) else { continue };
+        let Ok(tpi) = Tpi::new([[&pts[5], &pts[6], &pts[7]], [&pts[8], &pts[9], &pts[10]], [&pts[11], &pts[12], &pts[13]]]) else {
+            continue;
+        };
+        let rl = r_lpi(&pts[0], &pts[1], &pts[2], &pts[3], &pts[4]);
+        let rt = r_tpi(&[[pts[5], pts[6], pts[7]], [pts[8], pts[9], pts[10]], [pts[11], pts[12], pts[13]]]);
+        let g = [GenericPoint::Lpi(lpi.clone()), GenericPoint::Tpi(tpi.clone())];
+        let rg = [rl.clone(), rt.clone()];
+        let ex: Vec<Point3> = (0..4).map(|_| rand_pt(&mut r, 1.0)).collect();
+        // Every mix of implicit and explicit points in every position.
+        for mask in 0u32..81 {
+            let mut gp = Vec::new();
+            let mut rpv = Vec::new();
+            let mut m = mask;
+            for e in &ex {
+                match m % 3 {
+                    0 => {
+                        gp.push(GenericPoint::Explicit(*e));
+                        rpv.push(rp(e));
+                    }
+                    k => {
+                        gp.push(g[k as usize - 1].clone());
+                        rpv.push(rg[k as usize - 1].clone());
+                    }
+                }
+                m /= 3;
+            }
+            if !(i + mask as usize).is_multiple_of(7) {
+                continue; // a deterministic sample of the masks keeps the test fast
+            }
+            let want = r_orient3d(&[rpv[0].clone(), rpv[1].clone(), rpv[2].clone(), rpv[3].clone()]);
+            assert_eq!(generic_orient3d(&gp[0], &gp[1], &gp[2], &gp[3]), want, "mask {mask}");
+            checked += 1;
+        }
+        // The constructed points lie exactly on their defining planes.
+        let lp = GenericPoint::Lpi(lpi);
+        let tp = GenericPoint::Tpi(tpi);
+        let e = |p: &Point3| GenericPoint::Explicit(*p);
+        assert_eq!(generic_orient3d(&e(&pts[2]), &e(&pts[3]), &e(&pts[4]), &lp), Sign::Zero);
+        assert_eq!(generic_orient3d(&lp, &e(&pts[2]), &e(&pts[3]), &e(&pts[4])), Sign::Zero);
+        for k in 0..3 {
+            let pl = &pts[5 + 3 * k..8 + 3 * k];
+            assert_eq!(generic_orient3d(&e(&pl[0]), &e(&pl[1]), &tp, &e(&pl[2])), Sign::Zero);
+        }
+        // ... and on their line.
+        for axis in [Axis::X, Axis::Y, Axis::Z] {
+            assert_eq!(generic_orient2d(&e(&pts[0]), &e(&pts[1]), &lp, axis), Sign::Zero);
+        }
+        // Two different constructions of the same point compare equal on every axis.
+        let again = GenericPoint::Lpi(Lpi::new(&pts[1], &pts[0], &pts[3], &pts[4], &pts[2]).unwrap());
+        for axis in [Axis::X, Axis::Y, Axis::Z] {
+            assert_eq!(generic_compare(&lp, &again, axis), Ordering::Equal);
+            let want = rl[axis.index()].cmp(&rt[axis.index()]);
+            assert_eq!(generic_compare(&lp, &tp, axis), want);
+            let want = rl[axis.index()].cmp(&q(ex[0].0[axis.index()]));
+            assert_eq!(generic_compare(&lp, &e(&ex[0]), axis), want);
+        }
+        let approx = lp.approx();
+        for k in 0..3 {
+            assert!((approx[k] - rl[k].to_f64().unwrap()).abs() <= 1e-9 * (1.0 + approx[k].abs()));
+        }
+    }
+    assert!(checked > 500, "{checked}");
+}
+
+#[test]
+fn indirect_orient2d_matches_rational_constructions() {
+    let mut r = Rng(10);
+    for _ in 0..200 {
+        let pts: Vec<Point3> = (0..8).map(|_| rand_pt(&mut r, 10.0)).collect();
+        let Ok(lpi) = Lpi::new(&pts[0], &pts[1], &pts[2], &pts[3], &pts[4]) else { continue };
+        let rl = r_lpi(&pts[0], &pts[1], &pts[2], &pts[3], &pts[4]);
+        for axis in [Axis::X, Axis::Y, Axis::Z] {
+            let (u, v) = axis.projected();
+            for slot in 0..3 {
+                let mut g: Vec<GenericPoint> = pts[5..8].iter().map(|p| GenericPoint::Explicit(*p)).collect();
+                let mut rr: Vec<RPoint> = pts[5..8].iter().map(rp).collect();
+                g[slot] = GenericPoint::Lpi(lpi.clone());
+                rr[slot] = rl.clone();
+                let d = |p: &RPoint, o: &RPoint| vec![p[u].clone() - o[u].clone(), p[v].clone() - o[v].clone()];
+                let want = rsign(&rdet(&[d(&rr[0], &rr[2]), d(&rr[1], &rr[2])]));
+                assert_eq!(generic_orient2d(&g[0], &g[1], &g[2], axis), want);
+            }
+        }
+    }
+}
+
+#[test]
+fn undefined_constructions_are_rejected() {
+    let o = p3(0.0, 0.0, 0.0);
+    let (x, y, z) = (p3(1.0, 0.0, 0.0), p3(0.0, 1.0, 0.0), p3(0.0, 0.0, 1.0));
+    // Line along x, plane z = 0 containing it: parallel.
+    assert!(Lpi::new(&o, &x, &o, &x, &y).is_err());
+    // Degenerate line.
+    assert!(Lpi::new(&z, &z, &o, &x, &y).is_err());
+    // Collinear plane points.
+    assert!(Lpi::new(&o, &z, &o, &x, &p3(2.0, 0.0, 0.0)).is_err());
+    // A nearly parallel line is fine: z rises by one ulp over the segment.
+    assert!(Lpi::new(&p3(0.0, 0.0, -1.0), &p3(1.0, 0.0, -1.0 + f64::EPSILON), &o, &x, &y).is_ok());
+    // Two parallel planes.
+    let lifted = |p: &Point3| p3(p.0[0], p.0[1], 1.0);
+    assert!(Tpi::new([[&o, &x, &y], [&lifted(&o), &lifted(&x), &lifted(&y)], [&o, &x, &z]]).is_err());
+    let t = Tpi::new([[&o, &x, &y], [&o, &y, &z], [&o, &z, &x]]).unwrap();
+    let t = GenericPoint::Tpi(t);
+    for axis in [Axis::X, Axis::Y, Axis::Z] {
+        assert_eq!(generic_compare(&t, &GenericPoint::Explicit(o), axis), Ordering::Equal);
+    }
+}

crates/mesh/tests/predicates_bench.rs

@@ -0,0 +1,182 @@
+//! Throughput and filter hit rate of the exact predicates. Ignored by default; run with
+//! `cargo test --release -p modelcraft-mesh --test predicates_bench -- --ignored --nocapture`.
+
+use std::hint::black_box;
+use std::time::Instant;
+
+use modelcraft_mesh::predicates::{
+    GenericPoint, Lpi, Point3, Sign, generic_orient3d, insphere, insphere_filter, insphere_sos, orient3d, orient3d_filter,
+};
+
+struct Rng(u64);
+
+impl Rng {
+    fn next_u64(&mut self) -> u64 {
+        self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
+        let mut z = self.0;
+        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
+        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
+        z ^ (z >> 31)
+    }
+    fn unit(&mut self) -> f64 {
+        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
+    }
+    fn index(&mut self, n: usize) -> usize {
+        (self.next_u64() % n as u64) as usize
+    }
+}
+
+fn k3(p: &Point3) -> robust::Coord3D<f64> {
+    let [x, y, z] = p.coords();
+    robust::Coord3D { x, y, z }
+}
+
+fn p3(x: f64, y: f64, z: f64) -> Point3 {
+    Point3::new(x, y, z).unwrap()
+}
+
+/// Uniform points in the unit cube.
+fn random_cloud(r: &mut Rng, n: usize) -> Vec<Point3> {
+    (0..n).map(|_| p3(r.unit(), r.unit(), r.unit())).collect()
+}
+
+/// Vertices of a tessellated CAD part: a 100 × 60 × 40 mm box with 5 mm grid faces and a
+/// Ø30 cylinder boss (64 segments, 8 rings) on top. Coplanar and cospherical subsets abound.
+fn cad_cloud() -> Vec<Point3> {
+    let mut v = Vec::new();
+    let (lx, ly, lz, h) = (100.0, 60.0, 40.0, 5.0);
+    let steps = |l: f64| (0..=(l / h) as usize).map(move |i| i as f64 * h);
+    for x in steps(lx) {
+        for y in steps(ly) {
+            v.push(p3(x, y, 0.0));
+            v.push(p3(x, y, lz));
+        }
+        for z in steps(lz) {
+            v.push(p3(x, 0.0, z));
+            v.push(p3(x, ly, z));
+        }
+    }
+    for y in steps(ly) {
+        for z in steps(lz) {
+            v.push(p3(0.0, y, z));
+            v.push(p3(lx, y, z));
+        }
+    }
+    for ring in 0..8 {
+        for k in 0..64 {
+            let t = k as f64 * std::f64::consts::TAU / 64.0;
+            v.push(p3(50.0 + 15.0 * t.cos(), 30.0 + 15.0 * t.sin(), lz + ring as f64 * 2.5));
+        }
+    }
+    v
+}
+
+fn queries<const N: usize>(r: &mut Rng, cloud: &[Point3], n: usize) -> Vec<[Point3; N]> {
+    (0..n).map(|_| std::array::from_fn(|_| cloud[r.index(cloud.len())])).collect()
+}
+
+fn time<T>(label: &str, items: &[T], mut f: impl FnMut(&T) -> Sign) {
+    let start = Instant::now();
+    let mut acc = 0i64;
+    for _ in 0..5 {
+        for it in items {
+            acc += f(it) as i64;
+        }
+    }
+    let ns = start.elapsed().as_nanos() as f64 / (5 * items.len()) as f64;
+    println!("  {label:<34} {ns:7.1} ns/call");
+    black_box(acc);
+}
+
+fn hit_rate<T>(items: &[T], f: impl Fn(&T) -> Option<Sign>) -> f64 {
+    100.0 * items.iter().filter(|it| f(it).is_some()).count() as f64 / items.len() as f64
+}
+
+#[test]
+#[ignore = "benchmark: run in release with --ignored --nocapture"]
+fn predicate_throughput_and_filter_hit_rate() {
+    let mut r = Rng(42);
+    let n = 200_000;
+    for (name, cloud) in [("random unit cube", random_cloud(&mut r, 10_000)), ("CAD tessellation", cad_cloud())] {
+        let o: Vec<[Point3; 4]> = queries(&mut r, &cloud, n);
+        let s: Vec<[Point3; 5]> = queries(&mut r, &cloud, n);
+        println!("{name} ({} points, {n} queries)", cloud.len());
+        println!(
+            "  filter hit rate: orient3d {:.2}%, insphere {:.2}%",
+            hit_rate(&o, |[a, b, c, d]| orient3d_filter(a, b, c, d)),
+            hit_rate(&s, |[a, b, c, d, e]| insphere_filter(a, b, c, d, e))
+        );
+        time("orient3d", &o, |[a, b, c, d]| orient3d(a, b, c, d));
+        time("orient3d (robust crate)", &o, |[a, b, c, d]| {
+            let k = |p: &Point3| {
+                let [x, y, z] = p.coords();
+                robust::Coord3D { x, y, z }
+            };
+            Sign::of(robust::orient3d(k(a), k(b), k(c), k(d)))
+        });
+        time("orient3d naive f64 (inexact)", &o, |[a, b, c, d]| {
+            let ([ax, ay, az], [bx, by, bz], [cx, cy, cz], [dx, dy, dz]) = (a.coords(), b.coords(), c.coords(), d.coords());
+            let (adx, ady, adz, bdx, bdy, bdz, cdx, cdy, cdz) = (ax - dx, ay - dy, az - dz, bx - dx, by - dy, bz - dz, cx - dx, cy - dy, cz - dz);
+            Sign::of(adz * (bdx * cdy - cdx * bdy) + bdz * (cdx * ady - adx * cdy) + cdz * (adx * bdy - bdx * ady))
+        });
+        time("insphere", &s, |[a, b, c, d, e]| insphere(a, b, c, d, e));
+        time("insphere (robust crate)", &s, |[a, b, c, d, e]| {
+            let k = |p: &Point3| {
+                let [x, y, z] = p.coords();
+                robust::Coord3D { x, y, z }
+            };
+            Sign::of(robust::insphere(k(a), k(b), k(c), k(d), k(e)))
+        });
+        time("insphere_sos", &s, |[a, b, c, d, e]| insphere_sos(a, b, c, d, e));
+    }
+
+    // Fully degenerate CAD queries: four vertices of one planar face (orient3d = 0), and five
+    // grid vertices of one face (insphere = 0 when not all coplanar is impossible: use the
+    // cylinder rings, which are cospherical in groups of points on two rings).
+    let face: Vec<Point3> = (0..=20).flat_map(|i| (0..=12).map(move |j| p3(i as f64 * 5.0, j as f64 * 5.0, 40.0))).collect();
+    let o: Vec<[Point3; 4]> = queries(&mut r, &face, n);
+    let ring = |z: f64| -> Vec<Point3> {
+        (0..64)
+            .map(|k| {
+                let t = k as f64 * std::f64::consts::TAU / 64.0;
+                p3(50.0 + 15.0 * t.cos(), 30.0 + 15.0 * t.sin(), z)
+            })
+            .collect()
+    };
+    let (r0, r1) = (ring(40.0), ring(42.5));
+    let s: Vec<[Point3; 5]> = (0..n)
+        .map(|i| {
+            let mut q: [Point3; 5] = std::array::from_fn(|k| if k < 3 { r0[r.index(64)] } else { r1[r.index(64)] });
+            // Same angle on both rings: symmetric configurations that are exactly cospherical.
+            if i % 2 == 0 {
+                let k = r.index(64);
+                q[3] = r1[k];
+                q[4] = r1[(k + 32) % 64];
+            }
+            q
+        })
+        .collect();
+    println!("degenerate CAD queries ({n} each)");
+    println!(
+        "  filter hit rate: coplanar orient3d {:.2}%, ring insphere {:.2}%",
+        hit_rate(&o, |[a, b, c, d]| orient3d_filter(a, b, c, d)),
+        hit_rate(&s, |[a, b, c, d, e]| insphere_filter(a, b, c, d, e))
+    );
+    time("orient3d (coplanar face)", &o, |[a, b, c, d]| orient3d(a, b, c, d));
+    time("insphere (two cylinder rings)", &s, |[a, b, c, d, e]| insphere(a, b, c, d, e));
+    time("insphere_sos (two cylinder rings)", &s, |[a, b, c, d, e]| insphere_sos(a, b, c, d, e));
+    time("insphere (robust crate, rings)", &s, |[a, b, c, d, e]| Sign::of(robust::insphere(k3(a), k3(b), k3(c), k3(d), k3(e))));
+    time("orient3d (robust crate, face)", &o, |[a, b, c, d]| Sign::of(robust::orient3d(k3(a), k3(b), k3(c), k3(d))));
+
+    // Indirect: orient3d with one line-plane intersection and three explicit points.
+    let cloud = random_cloud(&mut r, 1000);
+    let mut items = Vec::new();
+    while items.len() < 20_000 {
+        let q: [Point3; 8] = std::array::from_fn(|_| cloud[r.index(cloud.len())]);
+        if let Ok(l) = Lpi::new(&q[0], &q[1], &q[2], &q[3], &q[4]) {
+            items.push([GenericPoint::Lpi(l), q[5].into(), q[6].into(), q[7].into()]);
+        }
+    }
+    println!("indirect (20000 queries)");
+    time("generic_orient3d (LPI + 3 explicit)", &items, |[a, b, c, d]| generic_orient3d(a, b, c, d));
+}

docs/adr/0001-libraries.md

@@ -19,7 +19,12 @@ through the same crate. Only `modelcraft-geometry` may depend on SolidCraft.
 **Meshing: our own**, implemented from published papers (`docs/meshing/literature.md`).
 Robust predicates follow Shewchuk's adaptive-precision arithmetic (the paper; his reference
 `predicates.c` is public domain). A permissive crate may be adopted instead after review; record
-it here.
+it here. *Decided 2026-10-13:* our own implementation, `modelcraft_mesh::predicates`, written from
+Shewchuk 1997, Devillers–Teillaud 2011 and Attene 2020 (no `unsafe`, deterministic, wasm-safe, a
+validated coordinate domain instead of Shewchuk's overflow/underflow caveat, and indirect
+predicates that `robust` lacks). `robust` (georust, MIT OR Apache-2.0) and `num-bigint` /
+`num-rational` / `num-traits` (MIT OR Apache-2.0) are dev-dependencies only: differential test
+oracles, never linked into ModelCraft.
 
 **Linear algebra:** to be decided in the solver task, among permissively licensed pure-Rust
 crates (for example `faer`, MIT). Record the choice here.