Exact Geometric Predicates and Adaptive Precision

Geometric algorithms do not consume coordinates; they consume decisions. Which side of a line a point falls on, whether three points are collinear, whether a circumcircle contains a fourth point. Each decision is a sign, and the whole correctness argument of a convex hull or a Delaunay triangulation rests on those signs being right. Floating-point evaluation of the corresponding determinants gets them right most of the time, and “most of the time” is not a hypothesis these algorithms can tolerate: a single wrong sign can make an output that is structurally invalid rather than merely slightly displaced.

The standard remedy is not to compute in exact arithmetic everywhere, which would be correct but slow, nor to add epsilon tolerances, which move the inconsistency rather than remove it. It is to compute a cheap floating-point approximation together with a bound on its error, accept the sign when the approximation is far enough from zero that the bound cannot flip it, and fall back to more precise arithmetic only in the rare cases where it can.

Where the error comes from

An orientation determinant computed as (ax-cx)*(by-cy) - (ay-cy)*(bx-cx) involves four subtractions, two multiplications, and one subtraction of the products. The subtractions of coordinates are where catastrophic cancellation begins: if ax and cx are large and nearly equal, their difference is exact only when the exponents are close, and the relative error of the difference can be enormous even though the absolute error is small.

The final subtraction is worse. When three points are nearly collinear, the two products are nearly equal, and their difference loses every leading bit they share. The absolute error of each product is on the order of the product’s magnitude times the unit roundoff, so the absolute error of the difference is on the same order, while the true value may be many orders of magnitude smaller. The computed sign then depends on rounding noise.

With double precision, the unit roundoff is 2^-53, about 1.1×10^-16. Coordinates near 10^6 therefore give products near 10^12 and an absolute error in the difference near 10^-4. Any true determinant smaller than that in magnitude has an unreliable sign. That is not a pathological corner: points sampled along a smooth curve with coordinates in the millions routinely produce determinants at that scale, which is why naive hull code fails on data that looks perfectly ordinary.

An error bound that makes the fast path safe

Shewchuk’s construction attaches a static error bound to each predicate. For the orientation determinant evaluated in floating point, the computed value det satisfies

|computed - exact| <= errbound * (|term1| + |term2|)

where the two terms are the products before subtraction and errbound is a small constant derived from the unit roundoff and the number of operations, roughly 3.3×10^-16 for doubles with the standard evaluation order. The bound uses the magnitudes actually computed, so it adapts to the data instead of relying on a fixed epsilon.

The decision procedure is then: if |det| > errbound * (|term1| + |term2|), the sign of det equals the sign of the exact value, and the fast path is done. Otherwise the result is uncertain and a more precise evaluation is required. On random or well-separated input the uncertain case essentially never fires, so the average cost is that of the plain floating-point expression plus a few fabs and one comparison.

/* Two-product: exact product of doubles as hi + lo, using FMA when available. */
static inline void two_product(double a, double b, double *hi, double *lo) {
    *hi = a * b;
    *lo = fma(a, b, -*hi);   /* lo is the rounding error of the product */
}

This is the primitive on which the exact expansions are built. A two_sum plays the same role for addition. Composing them yields a representation of the exact determinant as a sum of nonoverlapping doubles, whose leading term determines the sign with certainty. The expansions grow with each operation, which is why the fallback is staged: evaluate at the next level of precision only if the current level’s bound is inconclusive, and stop as soon as a term dominates the sum of the rest.

What the fallback costs and when it runs

Exact evaluation of the orientation determinant for double inputs can be done with a handful of two_product and two_sum operations, giving an exact sign at a cost of perhaps an order of magnitude more arithmetic than the approximate evaluation. That cost is paid only on the uncertain cases. For a Delaunay triangulation of points in general position, the proportion of uncertain cases is negligible; for points on a regular grid, it is large, because grid points are exactly collinear in vast numbers and every such test lands precisely on zero.

This is the real reason adaptive predicates matter. They make the cost proportional to the difficulty of the input rather than to its size. A grid input does more work on the exact path than random input of the same size, but it produces a correct triangulation instead of a contradictory one.

An alternative is to store coordinates as integers and evaluate the determinant in exact integer arithmetic, an approach used by libraries that convert inputs to a fixed grid. With 32-bit coordinates each difference fits in 33 bits and each product in 64, which a signed 64-bit word cannot hold but a 128-bit accumulator can; with 64-bit coordinates the products need close to 128 bits, so the accumulator must still leave headroom for the subtraction. The sign is then exact by construction with no error analysis at all. The cost is a wider multiply, and the constraint is that the input must be representable on the grid without losing distinctions that the algorithm depends on. Conversions that round coordinates to the nearest grid point can merge distinct points, which changes the input rather than the arithmetic, and that is a different problem with a different fix.

The tolerance trap

The intuitive alternative to exact predicates is to treat values within some epsilon of zero as zero, and to treat distances within epsilon as equal. This fails for two reasons. First, a fixed epsilon is not scale-invariant: 10^-9 is a meaningful separation at unit scale and meaningless at scale 10^6, so any fixed threshold is wrong for some part of the input domain. Second, and more damaging, a tolerance relation is not transitive. If a is within epsilon of b and b is within epsilon of c, a may be far from c. Algorithms that rely on equality being an equivalence relation, which includes every sorting and grouping step in geometry, then produce output whose parts contradict each other.

The practical symptom is an inconsistency that survives any amount of debugging because it is not a bug in one place. A point is classified as left of an edge by the hull routine and right of the same edge by the triangulation routine, and both routines are locally correct. The fix is to make the predicate exact and to let the topology follow from it, rather than to tune thresholds until the current dataset stops complaining.

Topology from predicates, not from coordinates

The deeper principle is that the output structure should be derived from the signs alone. If the algorithm’s branches are all taken on exact predicate results, the output is a valid combinatorial structure for the input, even if the coordinates are floating point. The vertices may be positioned with rounding, but the edges will not cross, no triangle will be inverted, and no vertex will be lost.

This separation is what makes it possible to store coordinates as doubles and still guarantee a correct triangulation. It also explains why algorithms are designed around predicates rather than around computed intersection points: computing the intersection of two lines introduces new rounding, while a predicate about existing points does not. Where an intersection point is genuinely needed, the Bentley–Ottmann sweep must handle it, and the usual technique is to report the crossing pair without materialising the coordinates, or to represent the intersection symbolically.

Where the same problem reappears

The in-circle predicate for Delaunay triangulation is a 4×4 determinant with terms of degree four in the coordinates, so its floating-point error is larger and its exact evaluation is correspondingly more involved. The same adaptive structure applies, with a different bound and a deeper expansion. Delaunay triangulation by empty circumcircles depends on that predicate being exact, because a single wrong in-circle result flips an edge into a configuration that is not a triangulation.

Point location and containment tests have the same character, as described in point in polygon: crossing number and winding number: the parity of crossings is a topological fact that a floating-point comparison can silently invert when the query point is nearly on an edge. The remedy is identical, and it is worth adopting as a house rule: any function whose return value is a discrete classification of geometric data should be evaluated by an exact or provably-correct-sign predicate, and any function that returns a coordinate should be understood to return an approximation.