Floating-Point Robustness: Epsilons, Exact Predicates and Filters

In the previous lesson, three points on the line y = 3x, written as 0.1, 0.3 and so on, came out as a left turn in double. The cross product was about 2e-17 instead of 0. In most numerical code that would be harmless: an area off by 2e-17 is still the right area.

A geometric predicate is different. Its result is a sign, and an algorithm branches on it. Report "left" when the truth is "on the line", and a convex hull keeps a point it should drop. Report "left" for one ordering of three points and "right" for another ordering of the same three points, and the algorithm's assumptions stop holding: hulls get dents, triangulations get overlapping triangles, and loops that rely on progress never finish. Kettner, Mehlhorn, Pion, Schirra and Yap collected such failures in a paper whose figures show convex hull programs dropping obvious extreme points and producing non-convex output on small inputs.

This lesson measures the problem first, then looks at the three standard responses: the epsilon (which does not work), integer coordinates (which work when you can use them) and filtered exact predicates (which work for arbitrary double input and cost almost nothing on typical data). From this lesson on, the course uses integer coordinates with stated bounds.

A double on every platform you are likely to use is an IEEE 754 binary64 number: a sign, an 11-bit exponent and a 53-bit significand (52 bits stored, plus an implicit leading 1). The C++ standard does not require this format. std::numeric_limits<double>::is_iec559 tells you whether you have it, and the exact predicate later in this lesson checks it with a static_assert.

Three facts about binary64 carry most of this lesson.

• Each basic operation is correctly rounded. +, -, *, /, sqrt and fma return the exact mathematical result rounded to the nearest representable double (ties go to the even one). The relative error of one operation is at most ε = 2^-53, about 1.1e-16, provided the result does not overflow or fall into the subnormal range.
• Integers are exact up to 2^53. Every integer with absolute value at most 2^53 (about 9.0e15) is a double. A sum, difference or product of such integers is exact if its true result also fits in 53 bits.
• Most decimal fractions are not representable. 0.1 is stored as the nearest binary fraction, which is slightly more than one tenth.

The product of two 26-bit integers needs at most 52 bits, so it is exact. The product of two 30-bit integers needs about 60, so the last seven or so bits are rounded away, and the result is off by 35. The orientation formula multiplies coordinate differences, so this is the line between "double is exact for my data" and "it is not": with integer coordinates below about 2^25 in magnitude, differences stay below 2^26, and the whole determinant is computed exactly in double.

The last line shows the other face of rounding. Subtracting two nearly equal doubles is exact (this is Sterbenz's lemma: if y/2 <= x <= 2y, then x - y is representable). The subtraction adds no error, but it removes the large equal parts and leaves the earlier errors as the whole answer. The orientation test near a line is exactly this situation: two nearly equal products, subtracted.

Kettner and coauthors give an experiment worth repeating. Take the line through q = (12, 12) and r = (24, 24), which is y = x. Near (0.5, 0.5), consecutive doubles are 2^-53 apart, so the points p = (0.5 + i·2^-53, 0.5 + j·2^-53) for i, j from 0 to 255 form a tiny 256 by 256 grid of exactly representable points. For each one the true answer is known without any arithmetic: p is left of the line from q to r when j > i, on it when j == i, right of it when j < i.

The program evaluates the textbook formula in double with p as the first point, as Kettner et al. do, and counts how often it is wrong. It also evaluates the two rotations (q, r, p) and (r, p, q), which in exact arithmetic always give the same sign.

Two kinds of error appear. The common one is rounding to zero: a point that is strictly off the line gets the answer "collinear". The dangerous one is the 672 sign inversions: a point clearly on one side is reported on the other. And for 13,968 points, more than a fifth of the grid, the answer depends on which of the three points you name first. Any algorithm that asks about the same triple twice, in different orders, can receive contradictory answers.

Here is the same grid as a picture. Light blue is the true left side, white the true right side, and the blue diagonal is the line itself. Dark cells were reported collinear although they are not; red cells were reported on the wrong side. Each cell is one representable point, and the picture is drawn with y pointing up (the Svg helper flips it for the screen).

The errors are not a thin fuzzy band along the line. They form blocks, with arms that reach well into both half-planes, and the wrong-side cells sit right next to the diagonal. The set of points that the formula calls "on the line" is not a line at all. Kettner et al. explain the block structure in detail; the short version is that the rounding of 0.5 + i·2^-53 - 12 happens in steps, because doubles near 11.5 are 16 times further apart than doubles near 0.5.

The next program builds convex hulls with Andrew's monotone chain, which lesson 4 explains in detail. All you need to know now is that it sorts the points and keeps a point only if it makes a left turn with the two before it, so the orientation test is its only geometric operation.

Each test set has the two points (12, 12) and (24, 24) plus 20 points within a few units in the last place of (0.5, 0.5). The hull is built once with the double formula and once with an exact predicate, and each result is checked with the exact predicate: a correct counter-clockwise hull turns strictly left at every corner and has no input point strictly to the right of any edge.

Almost a quarter of the hulls are wrong, and nothing crashes: the program returns a plausible list of points. The algorithm is the same in both runs. Only the predicate changed.

The input was chosen to sit on a line. Real data rarely looks like that everywhere, but often somewhere: points digitized along a straight wall, a mesh with many points on a plane. Those are exactly the inputs where the naive predicate fails.

The usual first fix is a tolerance: treat the cross product as zero if |cross| <= eps. It hides the 2e-17 from the previous lesson. It does not make the predicate correct, for two reasons that the next program shows.

An absolute epsilon depends on units. The same triangle answers "collinear" in meters and "left" in millimeters, because the cross product scales with the square of the coordinates. A relative epsilon, which compares the cross product with the size of the two products, fixes the scale problem. Neither fixes the second problem: "collinear within a tolerance" is not transitive. Exact collinearity is: if a, b, c lie on one line and b, c, d lie on one line (with b != c), then a, b, d lie on that line too. Algorithms rely on this without saying so, for example when they drop the middle point of every collinear triple.

Every consecutive triple on the curve is "collinear", yet the first and last points are not on the line through the first two. A hull algorithm that removes collinear middle points removes all four interior points here, one at a time, each step justified by the tolerance, and loses a real corner. A tolerance large enough to absorb the sign inversions from the grid experiment also turns a band of non-collinear points into "collinear", which is a different wrong answer.

There are legitimate tolerances. "Snap the cursor to a vertex within 5 pixels", "merge survey points closer than 1 mm", "consider two meshes touching if they are within 0.01 mm": these are decisions about the data, made once, outside the predicates, by someone who knows what the units mean. After the data has been cleaned, the predicates should answer exactly about the data they are given.

Every geometric computation does one of two things.

A predicate answers a question with one of a few discrete values: the sign of an orientation, whether a point is inside a circle, which of two distances is smaller. Because the answer is discrete, it can be computed exactly, as long as the inputs are taken as exact numbers.

A construction produces a new geometric object: the intersection point of two segments, the circumcenter of a triangle, a rotated point. The exact result is often not representable (the intersection of two lines through integer points is rational; a rotation usually involves an irrational cosine), so a construction must round.

The algorithms in this course make their decisions with predicates and use constructions only for output. If every decision is exact, the combinatorial result (which points are on the hull, which segments intersect, which triangles exist) is exactly right for the given input, even if a computed coordinate is off in its last bit.

CGAL, the Computational Geometry Algorithms Library, is organized around this split. A CGAL algorithm is written against a kernel, a bundle of number type, point types, predicates and constructions. One of its predefined kernels, CGAL::Exact_predicates_inexact_constructions_kernel, stores double coordinates, and its documentation states that it provides exact predicates while constructions are not exact in general. Its predicates are filtered: CGAL first evaluates them with interval arithmetic and computes exactly only when the interval does not determine the sign. CGAL::Exact_predicates_exact_constructions_kernel makes constructions exact as well, at a higher cost. The rest of this lesson builds the same two ideas, an exact representation and a filter, by hand.

The simplest exact predicate is the one from lesson 1: integer coordinates and integer arithmetic. With |x|, |y| <= 10^9, the orientation determinant is at most 8·10^18 in absolute value, which fits in long long. Every program from lesson 3 on uses this:

Real-valued data gets there by scaling to a grid: multiply by a fixed factor (1000 for meters to millimeters) and round once with std::llround. That rounding moves each point by at most half a grid step and can change the geometry: three nearly collinear points can become exactly collinear. If the application cannot accept that, it needs exact arithmetic on the original doubles, the subject of the next two sections.

Higher-degree predicates need more bits. The in-circle test of lesson 7 multiplies four coordinate differences, so long long overflows quickly. GCC and Clang provide __int128 as an extension (it is not standard C++ and MSVC does not have it), which covers the in-circle determinant for a smaller coordinate bound that lesson 7 works out. Beyond that you need a multiprecision type such as GMP's, or the techniques below.

Signed overflow is undefined behavior in C++, not merely a wrong answer, so the bound is part of the predicate's contract: put it in a named constant and a static_assert, and validate input against it where data enters the program.