Numerical model

The guarantee#

For every finite input, each predicate returns the sign of the exact determinant of its input points: the value the formula would have if it were evaluated in real arithmetic, with no rounding at any step. There is no range of magnitudes, no spread between the smallest and largest coordinate, and no degree of near-degeneracy beyond which the answer becomes approximate.

The mechanism is certification. A predicate evaluates its determinant in a ladder of stages. Each stage either proves that its result has the sign of the exact determinant, or passes the input to the next stage. No stage ever returns a sign it has not proved.

Why floating point fails#

Each predicate is the sign of a determinant. For orient2d:

Near a degenerate configuration the two products are almost equal, and their difference is smaller than the rounding error of computing them. The rounded result can then be zero when the true value is not, or have the opposite sign. Downstream code that branches on that sign — a triangulation, a polygon overlay, a point-in-polygon test — builds a structure that is topologically impossible, with no indication that anything went wrong.

use arithmetica::{Sign, orient2d};

// The textbook formula, evaluated in plain f64.
let naive = |a: [f64; 2], b: [f64; 2], c: [f64; 2]| {
    (a[0] - c[0]) * (b[1] - c[1]) - (a[1] - c[1]) * (b[0] - c[0])
};

// The third point is 2^-48 above the line y = x through the other two.
let (a, b, c) = ([0.5, 0.5], [12.0, 12.0], [24.0, 24.0 + 2f64.powi(-48)]);
assert_eq!(naive(a, b, c), 0.0);                   // rounding erased the answer
assert_eq!(orient2d(a, b, c), Ok(Sign::Positive)); // the exact sign

Shewchuk’s adaptive predicates fix this with exact floating-point arithmetic, but his proofs assume an unbounded exponent range. Binary64 has a bounded one: a product can overflow to infinity, or underflow to a subnormal with fewer than 53 significant bits or to zero. At the extremes of the range, the stages can then certify a sign that is wrong. Arithmetica closes that gap; the next two sections describe how.

The four stages#

Stage 1: the filter

The determinant d̃ is computed in plain f64, alongside a magnitude M built from the same terms in the same order, with the smallest normal number u_N = 2⁻¹⁰²² added after every product. The sign is accepted only if

where A is the predicate’s error-bound constant. The two absolute terms cover the error of products that land among the subnormals, where relative error bounds do not hold. The comparison is strict and the bound is never zero, so the filter never certifies an exactly degenerate input. Most ordinary input is decided here, for the cost of a few dozen floating-point operations.

Stages 2 and 3: adaptive expansions

An expansion is a sum of f64 values whose significands do not overlap; it can represent a sum or product of f64 values exactly. Stage 2 builds the determinant from exact products of the rounded coordinate differences and tests it against a tighter bound. Stage 3 adds the rounding errors of the differences themselves, giving the exact determinant. (insphere has no Stage 3: input that Stage 2 cannot decide goes to Stage 4.) These are Shewchuk’s algorithms; they run only on input inside the exactness window, below.

Stage 4: exact integers

Every finite f64 is an integer multiple of 2⁻¹⁰⁷⁴, the smallest subnormal. Scaling every coordinate by 2¹⁰⁷⁴ therefore turns the determinant into an integer polynomial whose sign is the sign of the original:

Stage 4 evaluates that polynomial in fixed-width two’s-complement integers, sized in advance so that no intermediate value can overflow:

Integer arithmetic has no rounding error, so this stage needs no error bound. It costs about two orders of magnitude more than the filter, but input that is neither exactly degenerate nor spread across more bits than the window never reaches it.

The exactness window#

For one call, let E be the largest exponent of a leading bit among the non-zero coordinates, and q the smallest exponent of a lowest set bit, so every coordinate is an integer multiple of 2^q and smaller than 2^E+1. For a predicate of degree d, the exactness window is the set of inputs with

Inside the window, every operation of Stages 2 and 3 returns the same value in binary64 as it would with an unbounded exponent range, so Shewchuk’s theorems hold unchanged. The lower limit keeps every exact value on the 2⁻¹⁰⁷⁴ grid and every error-bound product in the normal range, with a margin of 19 bits; the upper limit keeps every value below the overflow threshold, with a margin of 26 bits.

Outside the window, if the input’s spread E − q fits within the window’s width, every coordinate is multiplied by a single power of two. That multiplication is exact, and it multiplies the determinant by a positive number, so the sign is unchanged. If the spread is wider than the window, no scale factor works, and the call goes straight to Stage 4.

For ordinary data the window is far away: any orient2d input whose non-zero coordinates lie between about 2⁻³⁹⁸ and 2⁴⁷⁸ is inside it as given. The full statements and proofs are in docs/exactness_window.md in the repository.

Non-finite input#

A NaN or infinite coordinate has no geometric meaning, so there is no correct sign to return. Every predicate and construction checks its input first and returns Err(NonFinite). Returning Sign::Zero instead would be indistinguishable from a genuinely degenerate configuration.

A non-finite intermediate is different. Finite coordinates near f64::MAX can overflow a difference or a product to infinity. That is not an error: the filter declines to certify such a call, and the exact stages answer it like any other finite input.

Verification#

Random inputs almost never land where a predicate can go wrong: a predicate can only fail where the true determinant is smaller than its rounding error, next to the surface on which the sign changes. The test suite is built around inputs that are there.

• Boundary search. For each predicate and five magnitude profiles — ordinary, tiny, huge, mixed and the window edges — one coordinate is bisected to the pair of adjacent floats on either side of a sign change, and the predicate is checked on those floats and their neighbors. A deep run makes about 18 million such calls with no wrong sign and no false zero.
• Exact oracles. Every expected sign is computed in arbitrary-precision integer or rational arithmetic, never in floating point.
• Fixed inputs at the extremes of the range: subnormal products, determinants far below the smallest subnormal, huge and tiny coordinates in one call.
• Fuzzing with cargo-fuzz against the exact oracles, and a stable-toolchain replay of the stored corpus on every test run.
• Property tests with proptest.

Tests find defects; they do not prove their absence. The proof is the argument above and in docs/exactness_window.md. The tests check the one part of it that rests on inspection rather than derivation: that the code implements Shewchuk’s algorithms without an omitted term.

References#

• Shewchuk, J. R. (1997). Adaptive precision floating-point arithmetic and fast robust geometric predicates. Discrete & Computational Geometry, 18(3), 305–363. https://doi.org/10.1007/PL00009321. The expansion arithmetic, the adaptive predicates and their error bounds.
• Shewchuk, J. R. (1996). Routines for arbitrary precision floating-point arithmetic and fast robust geometric predicates (predicates.c). Public domain. https://www.cs.cmu.edu/~quake/robust.html. The reference implementation the adaptive stages are ported from.
• Bartels, T., Fisikopoulos, V., & Weiser, M. (2023). Fast floating-point filters for robust predicates. BIT Numerical Mathematics, 63, 31. https://doi.org/10.1007/s10543-023-00975-x. The underflow-protected filter (Theorem 1(2)) and the zero filter (Definition 6).