17.15 Precision Errors

Problem

Geometric algorithms often look exact on paper but fail in code because numeric computations are approximate.

A segment intersection test may say two segments do not intersect when they barely touch. A convex hull algorithm may remove the wrong point because a nearly collinear turn is misclassified. A polygon area calculation may lose meaningful digits when coordinates are large and differences are small.

The difficulty is that computational geometry depends heavily on predicates:

is this point left of this line?
are these points collinear?
do these segments intersect?
is this point inside this circle?

When these predicates are wrong, the algorithm’s control flow becomes wrong.

Solution

Treat numeric precision as part of the algorithm design.

For robust geometry, separate two kinds of computation:

predicates: decisions about signs and ordering
constructions: computed geometric objects

Predicates decide branches. They must be as reliable as possible.

Constructions produce coordinates, such as intersection points. They may tolerate more approximation depending on the application.

A practical rule:

use exact arithmetic for predicates when possible
use floating-point arithmetic carefully for measurements and output

The Classic Failure

Consider the orientation test:

def orientation(a, b, c):
    return (
        (b.x - a.x) * (c.y - a.y)
        -
        (b.y - a.y) * (c.x - a.x)
    )

Mathematically, the sign is exact.

In floating-point arithmetic, the computed result may be:

0.000000000000000003

or:

-0.000000000000000002

for points that are intended to be collinear.

If the algorithm interprets these tiny values as real left or right turns, the result can be topologically invalid.

Example: Nearly Collinear Points

Suppose:

A = (0, 0)
B = (1000000000, 1000000000)
C = (2000000000, 2000000000.000001)

The points are almost collinear.

The orientation result comes from subtracting two large products:

(b.x - a.x) × (c.y - a.y)
-
(b.y - a.y) × (c.x - a.x)

Each product is around:

10^18

The useful difference is tiny by comparison.

This is called cancellation.

When cancellation occurs, floating-point arithmetic may lose the information needed to classify the turn correctly.

Epsilon Comparisons

A common first defense is to use an epsilon threshold.

EPS = 1e-9

def sign(value):
    if value > EPS:
        return 1

    if value < -EPS:
        return -1

    return 0

Then:

def orientation_sign(a, b, c):
    return sign(orientation(a, b, c))

This is better than direct equality checks, but it is not a complete solution.

The right epsilon depends on coordinate scale.

For coordinates around 1, 1e-9 may be reasonable.

For coordinates around 1e12, it may be far too small.

For coordinates around 1e-12, it may be far too large.

Scale-Aware Epsilon

One improvement is to scale the threshold with the magnitude of the inputs.

def sign_scaled(value, scale):
    eps = 1e-12 * scale

    if value > eps:
        return 1

    if value < -eps:
        return -1

    return 0

For orientation, a rough scale can be based on coordinate magnitudes:

def orientation_sign_scaled(a, b, c):
    ux = b.x - a.x
    uy = b.y - a.y
    vx = c.x - a.x
    vy = c.y - a.y

    value = ux * vy - uy * vx
    scale = (
        abs(ux * vy)
        +
        abs(uy * vx)
        +
        1.0
    )

    return sign_scaled(value, scale)

This is still heuristic, but it behaves better across different coordinate ranges.

Prefer Integer Coordinates When Possible

Many geometric problems use integer coordinates:

pixels
grid positions
map tiles
competitive programming input
database coordinates after scaling

For integer coordinates, orientation and area computations can be exact if the numeric type is large enough.

In Python, integers have arbitrary precision:

def orientation_int(a, b, c):
    return (
        (b.x - a.x) * (c.y - a.y)
        -
        (b.y - a.y) * (c.x - a.x)
    )

No rounding occurs.

In C, C++, Java, Go, or Rust, use 64-bit or wider integer types when products can be large.

If coordinates may be up to:

10^9

then products may reach:

10^18

which requires 64-bit signed integers.

If coordinates may be larger, use arbitrary-precision integers or carefully bounded arithmetic.

Exact Predicates

For high-reliability geometry, use exact predicates.

An exact predicate computes only the sign of a determinant, such as:

orientation
in-circle

without necessarily computing all resulting geometry exactly.

This gives the algorithm reliable branching decisions.

Many robust geometry libraries use adaptive precision: they first try a fast floating-point computation, then fall back to exact arithmetic when the result is too close to call.

The principle is simple:

fast when the answer is obvious
exact when the answer is numerically delicate

Predicate vs Construction

Suppose two segments intersect.

Predicate:

do they intersect?

Construction:

where is the intersection point?

The predicate must be correct because it controls topology.

The constructed point may be approximate if the application only needs drawing or measurement.

For example, a renderer may tolerate a tiny coordinate error.

A planar subdivision or mesh generator may not.

Keep these two responsibilities separate in your design.

Floating-Point Equality

Avoid direct equality:

if x == y:
    ...

Use either exact arithmetic or tolerance-based comparison.

def almost_equal(a, b, eps=1e-9):
    return abs(a - b) <= eps

For scale-aware comparison:

def almost_equal_scaled(a, b, rel=1e-12, abs_eps=1e-15):
    return abs(a - b) <= max(
        abs_eps,
        rel * max(abs(a), abs(b))
    )

The absolute tolerance handles values near zero.

The relative tolerance handles large values.

Stable Distance Comparisons

Do not compute square roots when comparing distances.

Instead of:

if distance(a, b) < distance(c, d):
    ...

use:

if distance_squared(a, b) < distance_squared(c, d):
    ...

Implementation:

def distance_squared(a, b):
    dx = a.x - b.x
    dy = a.y - b.y

    return dx * dx + dy * dy

This avoids an unnecessary rounding step and is faster.

Sorting with Floating-Point Keys

Many geometric algorithms sort by angle, coordinate, or event position.

Floating-point ties are dangerous.

For example:

events.sort(key=lambda event: event.x)

works only if event coordinates are stable enough to sort.

When events are computed from intersections, nearly identical coordinates may appear in inconsistent order.

Use deterministic tie-breaking:

events.sort(
    key=lambda event: (
        round(event.x, 12),
        event.kind,
        event.id
    )
)

This is a pragmatic approach for moderate-precision systems.

For exact computational geometry, sort using exact comparison predicates instead.

Degeneracy and Precision Are Different

Degeneracy means the input has special structure:

three points are collinear
two segments overlap
four points are cocircular
a point lies exactly on a polygon edge

Precision error means arithmetic incorrectly reports the structure.

You need to handle both.

A robust segment intersection algorithm must know what to do with overlapping collinear segments.

A robust numeric implementation must also classify collinearity correctly.

Defensive Testing

Geometry code needs adversarial tests.

Include:

duplicate points
collinear points
near-collinear points
shared endpoints
overlapping segments
very large coordinates
very small coordinates
mixed-scale coordinates
zero-area polygons
nearly parallel lines

A useful pattern is to compare an optimized algorithm against a brute-force or exact implementation on small random inputs.

Example:

def test_against_brute_force(points):
    expected = closest_pair_brute(points)
    actual = closest_pair(points)

    assert expected[1] == actual[1]

For floating-point algorithms, compare within a tolerance, but keep predicate decisions explicit.

Common Pitfalls

Using One Global EPS Everywhere

A single epsilon rarely works across all coordinate scales.

Treating Near-Zero as Always Zero

This can hide real small features in the input.

Ignoring Overflow in Integer Geometry

Exact integer arithmetic only works if intermediate products fit the chosen type.

Computing Intersections Before Testing Existence

First decide whether objects intersect. Then compute the intersection point.

Relying on Visual Correctness

A plot can look correct while the topology is wrong.

Practical Guidelines

Use integer coordinates whenever the problem allows it.

Use squared distances for comparisons.

Use exact or adaptive predicates for production geometry.

Define boundary semantics explicitly.

Keep predicate logic separate from coordinate construction.

Test with degenerate and near-degenerate inputs.

Document the numeric assumptions of each algorithm.

Discussion

Precision errors are not incidental bugs in computational geometry. They are part of the problem domain. Geometric algorithms are built from sign tests, ordering decisions, and boundary classifications. When numeric error changes one of those decisions, the algorithm can fail structurally.

For simple scripts, epsilon comparisons may be enough. For serious geometry engines, exact predicates are the better foundation. The design goal is not to make every coordinate exact. It is to make every decision that affects topology reliable.