Put three points on a table. Do they turn left, turn right, or sit on one line? It is the smallest question in geometry, and your software asks it more often than any other. Every Boolean, convex hull, mesh intersection and Delaunay triangulation is built from millions of these yes/no answers. Computational geometers call them predicates. The orientation test is the simplest one: the sign of a 2×2 determinant, (q_x−p_x)·(r_y−p_y) − (q_y−p_y)·(r_x−p_x). Positive means left, negative means right, zero means collinear.

Here is the surprise. On a computer, that answer can come back wrong. Not slightly off. Wrong. Coordinates are stored as IEEE 754 doubles, a 53-bit significand times a power of two, and every subtraction and multiplication rounds. When the points are almost collinear, the two products are almost equal, and the rounding error can be larger than the real difference. The sign is the whole answer, and it comes out as zero or flips.

The staircase where a line should be

Lutz Kettner, Kurt Mehlhorn, Sylvain Pion, Stefan Schirra and Chee Yap made the problem visible in a paper written for teaching: Classroom Examples of Robustness Problems in Geometric Computations (Computational Geometry: Theory and Applications 40(1), 2008), from the Max Planck Institute for Informatics, INRIA Sophia Antipolis, Otto-von-Guericke-Universität Magdeburg and New York University. They took the diagonal through (12, 12) and (24, 24), tested query points near (0.5, 0.5) in steps of 2⁻⁵³ across a 256 × 256 grid, and coloured each point by the floating-point answer.

Kettner et al. write that they “expected to see a yellow band around the diagonal with nearly straight boundaries.” They found a staircase of blocks instead, with “points that change the side of the line.” The cause: 12 has four binary digits before the point, so the subtraction throws away the last four bits of the step, and the result stays constant for 2⁴ consecutive steps (2⁵ for 24). Nobody ordered stairs.

PAZ reran their setup on 2026-09-28, once in Python with exact fractions and once with the C# below; both gave identical counts. The naive double test disagrees with the exact answer on 11,972 of 65,536 points, about 18%: 11,300 wrongly called collinear, and 672 put on the wrong side outright. Only 256 grid points truly lie on the line. These are PAZ’s figures, not numbers from the paper.

One wrong sign is one bad answer. The trouble is that algorithms combine many answers and assume they agree. When they don’t, the result is topologically impossible. The classroom paper shows double-precision convex hulls that leave out an extreme point or intersect themselves, and implementations that “may even run forever.” According to the CGAL kernel manual, predicates “encapsulate decisions”, and inaccurate arithmetic “may lead to inconsistent decisions, causing unexpected failures for some correct input data.”

The fix: cheap first, exact only when it matters

Jonathan Richard Shewchuk’s answer, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates (Discrete & Computational Geometry 18(3):305–363, 1997, written at Carnegie Mellon; he is now a professor at UC Berkeley), fits in one sentence. Do the cheap calculation, work out how wrong it could possibly be, and only do more work when that error could flip the sign.

Two ideas carry it. An exact number can be stored as an expansion, a sum of non-overlapping doubles: the rounding error of a + b is itself exactly a double, recoverable in a few operations (the TWO-SUM lineage of Dekker, Knuth and Priest). And the computation runs in four stages. Stage A is the plain determinant plus an error-bound check; if the result beats the bound, the sign is guaranteed and work stops. Stages B and C add corrections, stage D computes the exact value, each reusing the last.

Shewchuk measured the payoff on a 2D Delaunay triangulation of a million-point tilted grid, one of his nastiest inputs: about 9.4 million orientation calls, 9,318,610 settled at stage A, 121,081 at B, 118 at C, and only 3 needing full exact arithmetic. According to Shewchuk, the cost is about 8% on random points and up to 30% on hard ones, “precisely for the point sets that are most likely to cause difficulties.” The trade-off is real: in 3D, robust predicates cost about 35% on random input and a factor of eleven on points near a sphere. Yet his 3D mesher Pyramid, run on the tilted grid without robust arithmetic, “failed to terminate”, caught in an infinite loop by a roundoff inconsistency. His verdict: “Robust arithmetic is not always slower after all.”

Look at who gets the credit. Shewchuk’s acknowledgements name Steven Fortune, Douglas Priest and Christopher Van Wyk as the foundations, and say Fortune “unwittingly sparked this research in mid-1994 with a few brief email responses.” A few emails became code that now sits, mostly unseen, under a great deal of geometry software.

←TODAY: Shewchuk’s 1997 predicates.c is public domain and still ships, directly or as ports, including Vladimir Agafonkin’s JavaScript robust-predicates. →3012: Zurich’s rotated grids and aligned façades still hang on the same four-stage sign decision. Fulcrum: code anyone may open and recompile outlives any vendor.

Why a building is the worst possible input

Buildings are full of degenerate cases on purpose: aligned walls, coplanar slabs, repeated modules, parts that meet exactly. A structural grid rotated to follow a plot boundary is the tilted grid from both papers. Rhino and Grasshopper handle this with the document’s absolute tolerance: agree in advance what distance counts as “the same point.” Tolerance manages constructions, where a point lands. Robust predicates keep decisions consistent, which side it lands on. Two layers of one problem, and your own scripts need both.

The Tool: Shewchuk’s adaptive predicates (orient2d, orient3d, incircle, insphere), public-domain C in predicates.c, ported to JavaScript by Vladimir Agafonkin as robust-predicates (Unlicense), with the same filter-first idea in CGAL’s filtered kernels. Today we build a small C# teaching version: a naive Fast test and a brute-force Exact test. Exact is only Shewchuk’s stage D; his real contribution is that you almost never have to run it.

Setup:

dotnet new console -n OrientLab
cd OrientLab
# paste the Orient class and the loop into Program.cs
dotnet run
# expected: 11972 wrong, 672 on the wrong side

using System;
using System.Numerics;

static class Orient
{
    public static int Fast(double px, double py, double qx, double qy, double rx, double ry)
        => Math.Sign((qx - px) * (ry - py) - (qy - py) * (rx - px));

    // Every double is m * 2^e: scale all six to integers, use BigInteger.
    public static int Exact(double px, double py, double qx, double qy, double rx, double ry)
    {
        double[] v = { px, py, qx, qy, rx, ry };
        var m = new BigInteger[6]; var e = new int[6]; int minE = int.MaxValue;
        for (int i = 0; i < 6; i++)
        {
            long bits = BitConverter.DoubleToInt64Bits(v[i]);
            int exp = (int)((bits >> 52) & 0x7FF);
            long man = bits & 0xFFFFFFFFFFFFFL;
            if (exp == 0) exp = 1; else man |= 1L << 52;
            m[i] = bits < 0 ? -man : man; e[i] = exp - 1075;
            if (man != 0) minE = Math.Min(minE, e[i]);
        }
        if (minE == int.MaxValue) return 0;
        for (int i = 0; i < 6; i++) m[i] <<= e[i] - minE;
        BigInteger d = (m[2] - m[0]) * (m[5] - m[1]) - (m[3] - m[1]) * (m[4] - m[0]);
        return d.Sign;
    }
}

First steps:

1. In Rhino 8, drop a C# script component onto the Grasshopper canvas and paste the Orient class below the script class, beside RunScript.
2. Inside RunScript, paste the loop:

   double u = Math.Pow(2, -53);
   int wrong = 0, flipped = 0;
   for (int X = 0; X < 256; X++)
     for (int Y = 0; Y < 256; Y++)
     {
       double px = 0.5 + X * u, py = 0.5 + Y * u;
       int f = Orient.Fast(px, py, 12, 12, 24, 24);
       int x = Orient.Exact(px, py, 12, 12, 24, 24);
       if (f != x) wrong++;
       if (f != 0 && f == -x) flipped++;
     }
   Print($"{wrong} wrong, {flipped} on the wrong side");

3. Read the output: 11972 wrong, 672 on the wrong side.
4. Stretch goal: output a coloured Point3d(X, Y, 0) per cell by Fast answer and bake. Plot grid indices, not real coordinates, which sit 2⁻⁵³ apart, far below any modelling tolerance. The staircase appears.

Atelier: Offices that write their own Grasshopper C# or Python components usually carry a few lines like if (cross == 0) deciding whether a column sits on a grid line or a room edge closes. They pass on the test model and break on the rotated wing. The Monday move: search your script library for zero or sign tests on raw cross products and route each through a filtered predicate like the Hack below before the next tender-stage model runs.

Hack: Gate the brute-force test behind Shewchuk’s stage A so you pay for exactness only when the bound demands it. The constant comes from his public-domain predicates.c, with epsilon 2⁻⁵³. If the cheap determinant beats the worst-case rounding error, its sign is certain.

double l = (qx - px) * (ry - py), r2 = (qy - py) * (rx - px), det = l - r2;
double eps = Math.Pow(2, -53), bound = (3.0 + 16.0 * eps) * eps * (Math.Abs(l) + Math.Abs(r2));
if (Math.Abs(det) > bound) return Math.Sign(det);   // stage A: sign is certain
return Orient.Exact(px, py, qx, qy, rx, ry);        // rare: go exact

From where I sit in the late 2070s, the geometry that kept working was the geometry whose source anyone could still compile. Write office scripts the same way: readable, open, fixable by a 25-year-old. Confidence doesn’t decide the sign.

Learn-it:

• Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, DCG 18(3), 1997 (CMU-CS-96-140R). Section 4 covers the stages.
• Kettner, Mehlhorn, Pion, Schirra, Yap, Classroom Examples of Robustness Problems in Geometric Computations, CGTA 40(1):61–78, 2008. The staircase is in Section 3.
• Vladimir Agafonkin’s robust-predicates: a readable JavaScript port with fast variants.
• CGAL kernel manual, Exact_predicates_inexact_constructions_kernel.
• PAZ note: PAZ Academy has taught Grasshopper 2 and C# scripting since the Alpha. Bring your office’s geometry scripts and we will run them against a tilted grid together.

Paste the class, run the 256 × 256 loop, and watch 672 points land on the wrong side before you trust another == 0.