Robust Orientation Predicate, in depth: exact signs from floating-point inputs with error-free transformations and an adaptive filter

What the predicate decides, and why one wrong sign matters

orient2d(a, b, c) is the sign of a 2x2 determinant, twice the signed area of triangle abc:

det = (ax - cx) * (by - cy) - (ay - cy) * (bx - cx)
sign(det) > 0   ->  a, b, c counter-clockwise (c left of a->b)
sign(det) < 0   ->  clockwise (c right of a->b)
sign(det) == 0  ->  collinear

Two facts drive everything else. The caller needs only the sign, never the exact value. And every double is an exact rational number, so the six stored values either are collinear or are not. Robustness means answering that correctly for the numbers actually stored.

A wrong answer is damaging because algorithms rely on consistency between calls. A predicate that says left for one ordering of three points and right for a rotation of them violates the axioms the algorithm's correctness proof uses. Kettner, Mehlhorn, Pion, Schirra and Yap's paper Classroom Examples of Robustness Problems in Geometric Computations shows hull and triangulation algorithms failing this way on innocent-looking input.

How the textbook formula fails, measured

To see the failure concretely, fix b = (12, 12) and c = (24, 24) and move a through a 64 by 64 grid of the smallest steps doubles allow near (0.5, 0.5): a = (0.5 + i·u, 0.5 + j·u) with u = 2^-53, which is exactly one unit in the last place at 0.5. All 4,096 points are distinct doubles. Exactly 64 of them, the diagonal i = j, lie on the line y = x; every other point is strictly on one side. The textbook formula (bx - ax)*(cy - ay) - (by - ay)*(cx - ax) pivots on a, the perturbed point. Checked against exact fractions.Fraction arithmetic:

False zeros drop points that belong on a hull, and flipped signs create inconsistent triples. The input is not exotic: points a few ulps off a line are routine in snapped GPS tracks, CAD edges and shared mesh vertices.

Two roundings cause it. The differences such as 12 - ax land where the ulp is 16u, so small offsets vanish before any multiplication. Then each product is about 270.25 with an ulp near 5.7e-14, while the true det, 12(j - i)u, is between 1.3e-15 and 8.4e-14 in size. An epsilon threshold does not fix this. It turns some wrong signs into zeros that disagree with exact collinearity, and it must be tuned to the coordinate scale.

Error-free transformations and exact expansions

The tool that fixes it is the error-free transformation. With IEEE 754 round-to-nearest arithmetic, and no overflow or underflow, the rounding error of a sum or a product is itself exactly representable as a double, and a few more flops compute it:

import math

def two_sum(a, b):
    """Return (s, e) with s = fl(a + b) and a + b == s + e exactly (Knuth)."""
    s = a + b
    bb = s - a
    e = (a - (s - bb)) + (b - bb)
    return s, e

def two_product(a, b):
    """Return (p, e) with p = fl(a * b) and a * b == p + e exactly.
    math.fma needs Python 3.13+; older code uses Dekker splitting instead."""
    p = a * b
    return p, math.fma(a, b, -p)

two_sum takes six operations and no branches. two_product uses a fused multiply-add, which computes a·b - p with one rounding, so the residual is exact. Without FMA, Dekker's algorithm splits each operand into 26-bit halves with the constant 2²⁷ + 1. That costs about 17 flops and can overflow near the top of the double range.

An expansion is a list of doubles whose exact sum is the value. If the components are nonoverlapping and sorted by increasing magnitude, the expansion's sign is that of its largest component. grow_expansion keeps that invariant with a chain of two_sums. The determinant's six products become twelve components, summed exactly.

def grow_expansion(e, b):
    """Add double b to expansion e (nonoverlapping, increasing magnitude). Zeros are dropped."""
    out, q = [], b
    for x in e:
        q, h = two_sum(q, x)
        if h != 0.0:
            out.append(h)
    if q != 0.0:
        out.append(q)
    return out

def orient2d_exact(a, b, c):
    """Exact sign of det, computed from the six products with no input subtraction."""
    (ax, ay), (bx, by), (cx, cy) = a, b, c
    e = []
    for u, v, s in ((ax, by, 1.0), (ax, cy, -1.0), (ay, bx, -1.0),
                    (ay, cx, 1.0), (bx, cy, 1.0), (by, cx, -1.0)):
        hi, lo = two_product(u, v)
        e = grow_expansion(e, s * hi)
        e = grow_expansion(e, s * lo)
    return 0 if not e else (1 if e[-1] > 0 else -1)

Expanding the determinant before multiplying avoids rounding the differences ax - cx, which are not exact in general. Multiplying by ±1 is exact. The result is exact for any finite inputs whose products neither overflow nor underflow into subnormals, because the two_product residual of a tiny product may not be representable.

The adaptive filter

The exact version is correct but slow. The trick that makes robust predicates practical is a floating-point filter: compute the determinant in plain floats, compute a rigorous bound on how far rounding could have moved it, and return the float sign whenever the result is farther from zero than the bound. Shewchuk's stage A is the first such filter:

EPS = 2.0 ** -53                          # unit roundoff for doubles
CCW_ERRBOUND_A = (3.0 + 16.0 * EPS) * EPS

def sign(x):
    return (x > 0) - (x < 0)

def orient2d(a, b, c):
    (ax, ay), (bx, by), (cx, cy) = a, b, c
    detleft = (ax - cx) * (by - cy)
    detright = (ay - cy) * (bx - cx)
    det = detleft - detright
    if detleft > 0.0:
        if detright <= 0.0:
            return sign(det)              # terms have opposite signs: no cancellation
        detsum = detleft + detright
    elif detleft < 0.0:
        if detright >= 0.0:
            return sign(det)
        detsum = -detleft - detright
    else:
        return sign(det)                  # detleft is exactly zero
    errbound = CCW_ERRBOUND_A * detsum
    if det >= errbound or -det >= errbound:
        return sign(det)                  # rounding cannot have crossed zero
    return orient2d_exact(a, b, c)        # in doubt: pay for exactness

Each early return can be justified in a sentence. Rounding a subtraction never changes its sign, and it yields zero only when the operands are equal, so detleft and detright carry the true signs of the true products. If those signs differ, or one is zero, the true difference cannot cancel and its sign is known. Otherwise both products contribute and could cancel. Shewchuk's error analysis shows that the computed det is within (3 + 16ε)ε times |detleft| + |detright| of the true value. Past that distance, the sign is certain.

Shewchuk's code adds stages B and C between the float test and the full expansion. Each stage reuses the previous one's work. Use his published code, or a port, in production.

Testing against an exact oracle

Test against an oracle that is obviously right. Fraction converts a double to its exact rational value. Half the cases are uniform triples. The other half are interpolated onto a random segment, so they are collinear up to one rounding, which is the case that matters.

import random
from fractions import Fraction

def orient2d_reference(a, b, c):
    ax, ay, bx, by, cx, cy = map(Fraction, (*a, *b, *c))
    return sign((bx - ax) * (cy - ay) - (by - ay) * (cx - ax))

rnd = random.Random(1)
for k in range(200_000):
    a = (rnd.uniform(-1, 1), rnd.uniform(-1, 1))
    b = (rnd.uniform(-1, 1), rnd.uniform(-1, 1))
    if k % 2:
        c = (rnd.uniform(-1, 1), rnd.uniform(-1, 1))
    else:
        t = rnd.random()
        c = (a[0] + t * (b[0] - a[0]), a[1] + t * (b[1] - a[1]))
    want = orient2d_reference(a, b, c)
    assert orient2d(a, b, c) == want
    assert orient2d_exact(a, b, c) == want
    # consistency under rotation and reflection of the argument order
    assert orient2d(b, c, a) == want and orient2d(b, a, c) == -want

On CPython 3.13 both functions agreed with Fraction on all 200,000 cases. In a separate count, none of 100,000 uniform triples reached the slow path, but about 64,000 of 100,000 near-collinear ones did. The rotation assertions matter as much as the value checks, because consistency is what the algorithms consume.

Worked example: three ulps off the line

Trace one grid point through the filter: a = (0.5 + 3u, 0.5), b = (12, 12), c = (24, 24). The point is three ulps right of the line y = x, so the true orientation is clockwise, -1.

1. The naive formula pivots on a. Its two products round to the same double, so it returns 0.0, a false collinear.
2. The filter pivots on c. detleft = (ax - 24)(12 - 24) and detright = (0.5 - 24)(12 - 24) are both positive, about 282. They can cancel, so detsum is about 564 and errbound is about 564 × 3.3e-16 ≈ 1.9e-13.
3. Both products round to exactly 282.0, so the float det is 0.0, inside the bound. The filter declares itself in doubt. That is the honest answer: the true det is -36u ≈ -4e-15.
4. orient2d_exact expands the six products into twelve components. The large terms cancel exactly in the two_sum chain, and the surviving top component is negative. It returns -1, which matches the Fraction reference.

Timing per call on this machine under CPython 3.13 was about 210 ns for the naive formula, 340 ns for the filtered predicate on a non-degenerate triple, 9.1 µs for the pure expansion and 20 µs for Fraction. Those are interpreter numbers. Measure compiled code separately. The lesson carries over: the filter costs a small factor over the naive formula, and the exact path should run only where the filter sends you.

What silently breaks a robust predicate

• FMA contraction. C and C++ compilers may fuse a*b - c*d into an FMA when -ffp-contract allows it, and some do so by default on FMA-capable targets. That changes the rounding the error bound was proven for, and it breaks two_sum and the splitting version of two_product outright. Compile predicate files with -ffp-contract=off, or the equivalent, and test the build you ship.
• Fast-math flags. -ffast-math lets the compiler simplify (a - (s - bb)) + (b - bb) algebraically to zero, deleting the error term.
• Extended precision. 32-bit x87 code keeps 80-bit intermediates, so values are double-rounded and the transformations stop being error-free. SSE2 avoids this.
• Overflow and underflow. The bounds assume neither happens. Coordinates near 1e154 overflow the products, and tiny coordinates produce subnormal residuals. Normalize or reject such input at the boundary.
• Exact predicates, rounded constructions. A float intersection point lies on neither segment, so feeding it back into orient2d asks an exact question about a wrong point.
• Already-rounded input. The predicate is exact for the doubles it receives. If parsing or a transform already rounded away collinearity, that is a data question, not a predicate bug.

Trade-offs and alternatives

Symbolic perturbation (Edelsbrunner and Mücke's Simulation of Simplicity) breaks exact zeros consistently and needs an exact predicate underneath. Snap rounding handles constructions by rounding intersection points to a grid while preserving topology.

What to do next

1. Run the test harness, and add the 64 by 64 grid as a regression test against Fraction.
2. Find every orientation or cross-product sign in your codebase, including the inline ones, and route all of them through a single orient2d function.
3. Delete epsilon comparisons used for geometric decisions.
4. In C or C++, use Shewchuk's predicates.c or a maintained port, compile it with -ffp-contract=off and without fast-math, and run the same Fraction-checked tests against the compiled build.
5. Log how often the slow path fires. A sudden rise means your data has become more degenerate.
6. Keep learning: the 2D geometry primitives kernel, Andrew's monotone chain hull, Graham scan and convex hull algorithms compared, all of which stand or fall on this predicate.