From 4667655e08a1fba38de492612ef7783be2c023bd Mon Sep 17 00:00:00 2001 From: Irakli Kourtis Date: Tue, 13 Oct 2026 15:04:56 -0400 Subject: [PATCH] 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 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 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 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm --------- Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com> Co-authored-by: Claude Opus 5.5 --- ATTRIBUTION.md | 9 + Cargo.lock | 34 ++ ROADMAP.md | 2 +- crates/mesh/Cargo.toml | 8 + crates/mesh/src/lib.rs | 2 + crates/mesh/src/predicates/adaptive.rs | 449 ++++++++++++++++++ crates/mesh/src/predicates/dyadic.rs | 253 ++++++++++ crates/mesh/src/predicates/eft.rs | 96 ++++ crates/mesh/src/predicates/expansion.rs | 316 ++++++++++++ crates/mesh/src/predicates/indirect.rs | 305 ++++++++++++ crates/mesh/src/predicates/interval.rs | 106 +++++ crates/mesh/src/predicates/mod.rs | 169 +++++++ crates/mesh/src/predicates/ring.rs | 46 ++ crates/mesh/src/predicates/sos.rs | 94 ++++ crates/mesh/src/predicates/tests.rs | 607 ++++++++++++++++++++++++ crates/mesh/tests/predicates_bench.rs | 182 +++++++ docs/adr/0001-libraries.md | 7 +- 17 files changed, 2683 insertions(+), 2 deletions(-) create mode 100644 crates/mesh/src/predicates/adaptive.rs create mode 100644 crates/mesh/src/predicates/dyadic.rs create mode 100644 crates/mesh/src/predicates/eft.rs create mode 100644 crates/mesh/src/predicates/expansion.rs create mode 100644 crates/mesh/src/predicates/indirect.rs create mode 100644 crates/mesh/src/predicates/interval.rs create mode 100644 crates/mesh/src/predicates/mod.rs create mode 100644 crates/mesh/src/predicates/ring.rs create mode 100644 crates/mesh/src/predicates/sos.rs create mode 100644 crates/mesh/src/predicates/tests.rs create mode 100644 crates/mesh/tests/predicates_bench.rs diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index 5861894..2e87f8d 100644 --- a/ATTRIBUTION.md +++ b/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 | diff --git a/Cargo.lock b/Cargo.lock index 3100976..4ceb38a 100644 --- a/Cargo.lock +++ b/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" diff --git a/ROADMAP.md b/ROADMAP.md index a2bf6cf..5215267 100644 --- a/ROADMAP.md +++ b/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 | diff --git a/crates/mesh/Cargo.toml b/crates/mesh/Cargo.toml index 00c27a8..d1df9b1 100644 --- a/crates/mesh/Cargo.toml +++ b/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" diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs index 4fb0843..d802175 100644 --- a/crates/mesh/src/lib.rs +++ b/crates/mesh/src/lib.rs @@ -20,3 +20,5 @@ pub enum Error { } pub type Result = std::result::Result; + +pub mod predicates; diff --git a/crates/mesh/src/predicates/adaptive.rs b/crates/mesh/src/predicates/adaptive.rs new file mode 100644 index 0000000..7f6739c --- /dev/null +++ b/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 { + 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 { + 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 { + orient3d_a(a, b, c, d).ok() +} + +#[inline(always)] +fn orient3d_a(a: &Point3, b: &Point3, c: &Point3, d: &Point3) -> Result { + 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 { + incircle_a(a, b, c, d).ok() +} + +#[inline(always)] +fn incircle_a(a: &Point2, b: &Point2, c: &Point2, d: &Point2) -> Result { + 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 { + insphere_a(a, b, c, d, e).ok() +} + +#[inline(always)] +fn insphere_a(a: &Point3, b: &Point3, c: &Point3, d: &Point3, e: &Point3) -> Result { + 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) + } +} diff --git a/crates/mesh/src/predicates/dyadic.rs b/crates/mesh/src/predicates/dyadic.rs new file mode 100644 index 0000000..4737b73 --- /dev/null +++ b/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, + exp: i64, +} + +impl Dyadic { + fn zero() -> Self { + Self { neg: false, mag: Vec::new(), exp: 0 } + } + + fn normalized(neg: bool, mut mag: Vec, 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 { + 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 { + 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 { + 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 { + 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 { + 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 + } +} diff --git a/crates/mesh/src/predicates/eft.rs b/crates/mesh/src/predicates/eft.rs new file mode 100644 index 0000000..857dd95 --- /dev/null +++ b/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); + } +} diff --git a/crates/mesh/src/predicates/expansion.rs b/crates/mesh/src/predicates/expansion.rs new file mode 100644 index 0000000..1f17290 --- /dev/null +++ b/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 { + 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 { + c: [f64; N], + n: usize, +} + +impl Exp { + #[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(e: &Exp, f: &Exp) -> 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(e: &Exp, 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); + +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 { + 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]); + } +} diff --git a/crates/mesh/src/predicates/indirect.rs b/crates/mesh/src/predicates/indirect.rs new file mode 100644 index 0000000..766c882 --- /dev/null +++ b/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 { + x: [T; 3], + w: Option, +} + +fn lift(p: &[f64; 3]) -> [T; 3] { + p.map(T::from_f64) +} + +fn sub3(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(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(v: T, w: &Option) -> T { + match w { + None => v, + Some(w) => v * w.clone(), + } +} + +/// The sign of `W` (positive for explicit points), `None` if uncertain. +fn w_sign(w: &Option) -> Option { + 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([q1, q2, r, s, t]: &[[f64; 3]; 5]) -> Hom { + let (q1, q2, r, s, t) = (lift::(q1), lift::(q2), lift::(r), lift::(s), lift::(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(planes: &[[[f64; 3]; 3]; 3]) -> Hom { + let plane = |[r, s, t]: &[[f64; 3]; 3]| -> ([T; 3], T) { + let r = lift::(r); + let u = sub3(&lift::(s), &r); + let v = sub3(&lift::(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, exact: impl FnOnce() -> Hom) -> 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, +} + +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 { + let points = [q1.0, q2.0, r.0, s.0, t.0]; + let filter = lpi_hom::(&points); + if exact_w_sign(&filter, || lpi_hom::(&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, +} + +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 { + let planes = planes.map(|p| p.map(|q| q.0)); + let filter = tpi_hom::(&planes); + if exact_w_sign(&filter, || tpi_hom::(&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 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 { + 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 { + 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( + points: [&GenericPoint; N], + filter: impl FnOnce([Hom; N]) -> Option, + exact: impl FnOnce([Hom; N]) -> Option, +) -> 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(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([a, b, c, d]: [Hom; 4]) -> Option { + let row = |p: &Hom| -> [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([a, b, c]: [Hom; 3], (u, v): (usize, usize)) -> Option { + let row = |p: &Hom, 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([a, b]: [Hom; 2], k: usize) -> Option { + 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, + } +} diff --git a/crates/mesh/src/predicates/interval.rs b/crates/mesh/src/predicates/interval.rs new file mode 100644 index 0000000..665f575 --- /dev/null +++ b/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 { + 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)); + } +} diff --git a/crates/mesh/src/predicates/mod.rs b/crates/mesh/src/predicates/mod.rs new file mode 100644 index 0000000..d3d4f56 --- /dev/null +++ b/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(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 { + 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 { + 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 { + 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 { + check(c).map(Self) + } +} + +#[cfg(test)] +mod tests; diff --git a/crates/mesh/src/predicates/ring.rs b/crates/mesh/src/predicates/ring.rs new file mode 100644 index 0000000..6467c7b --- /dev/null +++ b/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 + Sub + Mul { + /// 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; +} + +/// `a * d - b * c`. +#[inline] +pub(crate) fn det2(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(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(x: &T, y: &T) -> T { + x.clone() * x.clone() + y.clone() * y.clone() +} + +#[inline] +pub(crate) fn sq3(x: &T, y: &T, z: &T) -> T { + sq2(x, y) + z.clone() * z.clone() +} diff --git a/crates/mesh/src/predicates/sos.rs b/crates/mesh/src/predicates/sos.rs new file mode 100644 index 0000000..91d6923 --- /dev/null +++ b/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(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(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 +} diff --git a/crates/mesh/src/predicates/tests.rs b/crates/mesh/src/predicates/tests.rs new file mode 100644 index 0000000..c9afd53 --- /dev/null +++ b/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 { + if m.len() == 1 { + return m[0][0].clone(); + } + let mut total = BigRational::zero(); + for j in 0..m.len() { + let minor: Vec> = + 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 { + p.iter().zip(o).map(|(a, b)| q(*a) - q(*b)).collect() +} + +fn lifted(mut row: Vec) -> Vec { + 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 { + robust::Coord { x: p.0[0], y: p.0[1] } +} +fn rc3(p: &Point3) -> robust::Coord3D { + 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 = (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 = (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 { + if m.len() == 1 { + return m[0][0]; + } + (0..m.len()) + .map(|j| { + let minor: Vec> = 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::>(); + 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| { + 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 = (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]) -> Sign { + let mut order: Vec = (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> = 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> = (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 = 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> = (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 = 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> = 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> = (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 = (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 = (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 = (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 = pts[5..8].iter().map(|p| GenericPoint::Explicit(*p)).collect(); + let mut rr: Vec = 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); + } +} diff --git a/crates/mesh/tests/predicates_bench.rs b/crates/mesh/tests/predicates_bench.rs new file mode 100644 index 0000000..469d226 --- /dev/null +++ b/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 { + 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 { + (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 { + 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(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(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(items: &[T], f: impl Fn(&T) -> Option) -> 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 = (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 { + (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)); +} diff --git a/docs/adr/0001-libraries.md b/docs/adr/0001-libraries.md index 64519eb..341186d 100644 --- a/docs/adr/0001-libraries.md +++ b/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.