The Problem, Restated Precisely
Part 1 showed that an epsilon threshold can't fix a floating-point orientation test: rounding error in (b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x) scales with the magnitude of the inputs, so no single $\varepsilon$ is simultaneously tight enough to distinguish genuinely distinct nearly-collinear points and loose enough to absorb error on inputs with large coordinates. Two independent problems need two independent fixes: getting the sign of a predicate exactly right even when it's numerically close to zero, and having a well-defined policy for input that is exactly degenerate (truly collinear, truly cocircular). This deep dive covers both.
Fix 1: Exact Arithmetic
The most direct fix is to never round in the first place. If input coordinates are integers (or can be scaled to integers, e.g. fixed-point GPS or CAD coordinates), the orientation determinant's exact value fits in a wider integer type: two int32 coordinates produce a determinant that fits safely in int64, and the equivalent bound for the in-circle determinant (a degree-2 polynomial in the coordinates) requires roughly double the input bit width again. This is why production geometry kernels document the exact input coordinate range they support before guaranteeing predicate correctness — the promise is only as strong as the integer width backing it.
When coordinates aren't naturally integral, exact rational arithmetic (arbitrary-precision numerator/denominator pairs) gives the same exactness guarantee at the cost of speed: every rational operation can grow the numerator and denominator, and after a long chain of computations (e.g. many Delaunay flips) the numbers involved can become enormous. Libraries like CGAL and LEDA made "exact when needed, fast otherwise" a first-class design goal in the 1990s specifically because neither pure floating point nor pure exact rational arithmetic alone was acceptable for production use — one is unreliable, the other is too slow to use unconditionally.
Fix 2: Adaptive-Precision Floating Point
Jonathan Shewchuk's 1996 paper "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates" is the paper that made exact predicates fast enough for everyday use, and its reference implementation (predicates.c) still ships inside widely used mesh generators (including Shewchuk's own Triangle and, downstream, TetGen) decades later. The core idea:
- Compute the predicate (e.g. the orientation determinant) using ordinary double-precision floating point — this is fast, and correct the overwhelming majority of the time.
- Alongside that computation, track a cheap error bound: an upper bound on how far the floating-point result could possibly be from the true value, derivable purely from the magnitudes of the inputs and the number of floating-point operations performed.
- If the computed result's magnitude exceeds the error bound, its sign is certainly correct — return immediately. This is the fast path and it is taken almost every time in practice.
- Only when the result falls within the error bound (meaning the true sign is genuinely in doubt) does the algorithm fall back to a slower, exact computation — using a technique for representing floating-point sums as an exact expansion of multiple floating-point components, only as precise as the specific input actually requires ("adaptive" precision, rather than always using maximum precision).
Core Insight: the entire performance story of adaptive-precision predicates rests on one observation — the slow, exact path is only needed for inputs that are close to degenerate, which in practice is a small fraction of all predicate evaluations. Paying for exactness only when the fast path's own error bound says it's needed is what makes "always exact, usually fast" achievable at all.
Where This Shows Up in Practice
Early mesh-generation and CAD systems in the 1980s were notorious for crashing, hanging, or silently producing self-intersecting geometry on input that looked perfectly ordinary to a human — three points that looked collinear on screen but weren't quite, or four points that looked like they were on a common circle but weren't. Shewchuk's adaptive predicates are the reason mesh generators built after the mid-1990s stopped exhibiting this class of bug on ordinary CAD input, without becoming too slow to use.
Fix 3: Simulation of Simplicity
Exact and adaptive-precision arithmetic solve the numerical half of the problem: getting a predicate's sign exactly right. They do nothing about the combinatorial half — what an algorithm should actually do when a predicate's true, exact answer really is zero (three points really are exactly collinear; four points really do lie on one exact circle). Handling every such case with bespoke logic scattered through an algorithm is exactly the kind of special-casing this series avoids.
Herbert Edelsbrunner and Ernst Mücke's 1990 paper "Simulation of Simplicity: A Technique to Cope with Degenerate Cases in Geometric Algorithms" gives a general recipe instead: conceptually perturb every input point by an infinitesimal amount, using a fixed, deterministic rule based on the point's index (not a random perturbation, and never actually changing the stored coordinates) so that no predicate, applied to the perturbed points, can ever return exactly zero. Concretely, whenever a predicate's leading (exact) term evaluates to zero, a fully determined tie-breaking rule — a polynomial in a symbolic infinitesimal $\varepsilon$, evaluated using the point's index as a tiebreaker — decides the sign instead. The practical effect: an algorithm that was written assuming "generic position" (no three points collinear, no four points cocircular) can run unmodified on genuinely degenerate input, and will always make the same, well-defined choice.
Choosing an Approach
| Situation | Recommended approach |
|---|---|
| Coordinates are naturally integral, bounded range known | Exact integer arithmetic, sized for the worst-case determinant |
| Coordinates are floating point, performance matters | Adaptive-precision predicates (Shewchuk-style) |
| Algorithm assumes "general position" but real input isn't | Simulation of simplicity as a tie-breaking layer on top of either of the above |
| Prototyping, visualization, no correctness guarantee needed | Plain floating point is fine |
Every implementation-oriented part later in this series states which of these it assumes, rather than leaving robustness as an unstated assumption — the single most common way a "textbook-correct" geometric algorithm fails on real input.