17.17 Floating-Point Robustness
Some geometric problems cannot stay entirely in integer arithmetic.
17.17 Floating-Point Robustness
Problem
Some geometric problems cannot stay entirely in integer arithmetic.
Examples include:
line intersections
circle intersections
rotations
projections
normalization
physics simulation
computer graphics
geographic coordinates
These computations naturally produce real-valued coordinates. Floating-point arithmetic is fast and widely supported, but it introduces rounding, cancellation, overflow, underflow, and comparison errors.
The goal is not to avoid floating point entirely. The goal is to use it deliberately enough that algorithms remain predictable.
Solution
Use floating-point arithmetic for measurements and constructed coordinates, but protect the decisions that control geometry.
A practical design separates:
measurement code
predicate code
construction code
Measurements compute values such as length, angle, and area.
Predicates decide signs, order, containment, and intersection.
Constructions compute new points, such as projections and intersection coordinates.
The most important rule is:
do not let unstable floating-point predicates silently control topology
If a predicate decides whether two segments intersect, whether a point is inside a polygon, or whether a triangle edge should be flipped, its reliability matters more than the exact formatting of a computed coordinate.
Floating-Point Basics
Most programming languages use IEEE 754 double precision for float values.
A double-precision value has about:
15 to 17 decimal digits
of precision.
This means numbers such as:
1.0
and:
1.000000000000001
may be distinguishable, but differences much smaller than that may vanish.
The spacing between representable numbers grows as the numbers get larger.
Near:
1
floating point is dense.
Near:
10^12
the gaps are much larger.
This matters when coordinates have large magnitude and small differences.
Absolute Tolerance
A common comparison helper uses absolute tolerance:
def almost_equal_abs(a, b, eps=1e-9):
return abs(a - b) <= eps
This is useful when values are known to live near a fixed scale.
Example:
if almost_equal_abs(length, 0.0):
...
For coordinates and geometric predicates with unknown scale, absolute tolerance alone is usually insufficient.
Relative Tolerance
Relative tolerance scales with the magnitude of the values.
def almost_equal_rel(a, b, rel=1e-12):
return abs(a - b) <= rel * max(abs(a), abs(b))
This behaves better for large numbers.
However, relative tolerance fails near zero because the right side becomes too small.
Combined Tolerance
Use both absolute and relative tolerance.
def almost_equal(a, b, rel=1e-12, abs_eps=1e-15):
return abs(a - b) <= max(
abs_eps,
rel * max(abs(a), abs(b))
)
This handles both:
values near zero
and:
large-magnitude values
more reliably than either method alone.
Sign with Tolerance
Many geometric predicates reduce to the sign of a value.
def sign(value, eps=1e-12):
if value > eps:
return 1
if value < -eps:
return -1
return 0
This is useful for small and moderate coordinate ranges.
For scale-aware geometric predicates, compute a problem-specific scale.
def sign_scaled(value, scale, rel=1e-12, abs_eps=1e-15):
eps = max(abs_eps, rel * scale)
if value > eps:
return 1
if value < -eps:
return -1
return 0
Robust Orientation Approximation
A floating-point orientation test should account for scale.
def orientation_value(a, b, c):
ux = b.x - a.x
uy = b.y - a.y
vx = c.x - a.x
vy = c.y - a.y
return ux * vy - uy * vx
A scale-aware sign function:
def orientation_sign(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)
return sign_scaled(value, scale)
This still remains an approximation. It is a pragmatic improvement, not a substitute for exact predicates in high-reliability geometry.
Stable Segment Intersection
For floating-point coordinates, use tolerance-aware orientation and bounding checks.
def between(a, b, x, eps=1e-12):
return min(a, b) - eps <= x <= max(a, b) + eps
def on_segment(a, b, p):
if orientation_sign(a, b, p) != 0:
return False
return (
between(a.x, b.x, p.x)
and
between(a.y, b.y, p.y)
)
The segment intersection predicate follows the same structure as the integer version, but uses signs instead of raw orientation values.
def segments_intersect(a, b, c, d):
o1 = orientation_sign(a, b, c)
o2 = orientation_sign(a, b, d)
o3 = orientation_sign(c, d, a)
o4 = orientation_sign(c, d, b)
if o1 * o2 < 0 and o3 * o4 < 0:
return True
if o1 == 0 and on_segment(a, b, c):
return True
if o2 == 0 and on_segment(a, b, d):
return True
if o3 == 0 and on_segment(c, d, a):
return True
if o4 == 0 and on_segment(c, d, b):
return True
return False
This is suitable for many practical tasks, but exact predicates remain preferable for mesh topology, CAD kernels, and planar graph construction.
Avoiding Unnecessary Division
Division often magnifies errors, especially when the denominator is small.
Before computing a line intersection, check whether the lines are parallel or nearly parallel.
def cross(a, b):
return a.x * b.y - a.y * b.x
def line_intersection(a, b, c, d):
r = Vector(b.x - a.x, b.y - a.y)
s = Vector(d.x - c.x, d.y - c.y)
denominator = cross(r, s)
scale = abs(r.x * s.y) + abs(r.y * s.x)
if sign_scaled(denominator, scale) == 0:
return None
ca = Vector(c.x - a.x, c.y - a.y)
t = cross(ca, s) / denominator
return Point(
a.x + t * r.x,
a.y + t * r.y
)
The denominator test prevents unstable division when lines are nearly parallel.
Normalize Carefully
Vector normalization divides by length.
from math import sqrt
def normalize(v):
length = sqrt(v.x * v.x + v.y * v.y)
if almost_equal(length, 0.0):
raise ValueError("cannot normalize zero-length vector")
return Vector(
v.x / length,
v.y / length
)
Always handle zero-length vectors.
They occur in duplicate points, degenerate segments, collapsed polygons, and noisy input.
Prefer Squared Distances for Ordering
When comparing distances, avoid square roots.
def distance_squared(a, b):
dx = a.x - b.x
dy = a.y - b.y
return dx * dx + dy * dy
Use:
if distance_squared(a, b) < distance_squared(c, d):
...
This is faster and removes one source of rounding.
Only compute the square root for final output.
Coordinate Translation
If coordinates are large but local differences are small, translate the data near the origin before computation.
For example, instead of computing with points around:
(1000000000.0001, 1000000000.0002)
subtract a common origin:
def translate_points(points, origin):
return [
Point(
point.x - origin.x,
point.y - origin.y
)
for point in points
]
This can reduce cancellation in local geometric operations.
Translation does not change distances, orientation signs, intersections, or polygon shape.
It only improves numerical conditioning.
Scaling Coordinates
Sometimes input values use inconvenient units.
Example:
latitude and longitude
may need projection into a planar coordinate system before Euclidean geometry is meaningful.
Example:
meters
may be preferable to degrees for local distance calculations.
Scaling can also make tolerances easier to reason about.
A common pattern:
normalize input coordinates into a moderate numeric range
run geometry
transform output back
Sorting Events
Sweep-line algorithms using floating-point event coordinates need deterministic ordering.
Avoid comparing only one computed coordinate when ties are possible.
events.sort(
key=lambda event: (
round(event.x, 12),
event.kind,
event.id
)
)
This is a pragmatic engineering approach. For exact sweep-line algorithms, event ordering should be based on exact predicates rather than rounded keys.
Testing Floating-Point Geometry
Use tests that stress precision.
nearly parallel lines
nearly collinear points
very small triangles
large coordinate offsets
duplicate points
shared endpoints
points just inside and just outside boundaries
A useful strategy is metamorphic testing.
For example, if you translate all points by the same vector, the answer should not change.
def translated(points, dx, dy):
return [
Point(p.x + dx, p.y + dy)
for p in points
]
Test:
convex_hull(points)
and:
convex_hull(translated(points, dx, dy))
should have the same shape after translating back.
Common Pitfalls
Using One Epsilon for Everything
Distance comparisons, orientation tests, and coordinate equality often need different tolerances.
Treating Approximate Collinearity as True Geometry
An epsilon can hide narrow but real features.
Computing Angles Unnecessarily
Angles require trigonometric functions and introduce additional numerical issues.
Many algorithms can use cross products and dot products instead.
Forgetting Zero-Length Cases
Duplicate points create zero vectors and collapsed segments.
Assuming Visual Output Proves Correctness
A polygon may render correctly but still have invalid topology.
Practical Guidelines
Use integer or rational arithmetic when exact predicates are easy.
Use floating point for measurements, rendering, simulation, and approximate constructions.
Use scale-aware comparisons for predicate signs.
Translate and scale coordinates to improve conditioning.
Avoid division, square roots, and trigonometry unless they are necessary.
Test near degeneracies explicitly.
Document the numeric range your implementation expects.
Discussion
Floating-point robustness is a design discipline. The problem is not that floating point is unusable. The problem is that geometric algorithms often use numeric results to make structural decisions. A wrong sign or wrong ordering can invalidate the entire result.
For everyday geometry tasks, tolerance-aware predicates and careful normalization are often enough. For topology-sensitive software, such as mesh generation, CAD, GIS overlay, and planar graph construction, exact or adaptive predicates are the safer foundation.