Exact geometric predicates · onemath cookbook

Decide orientation, in-circle and segment intersection with OneMath’s exact predicates, and see how often the naive floating-point determinant gets the sign wrong.

The classic test, after Shewchuk: a line through (12, 12) and (24, 24), and a grid of 256 × 256 points around (0.5, 0.5), one unit in the last place apart. Each of them lies on, above or below the line; the naive formula is asked which.

using RecursiveAI.OneMath.Predicates;
 
// Shewchuk's example: is p left of, right of, or on the line through a and b?
// p sweeps a 256 x 256 grid of the doubles nearest (0.5, 0.5), one ulp apart.
double ax = 12, ay = 12, bx = 24, by = 24;
double ulp = Math.Pow(2, -53);
int wrong = 0;
for (int i = 0; i < 256; i++)
for (int j = 0; j < 256; j++)
{
    double px = 0.5 + i * ulp, py = 0.5 + j * ulp;
    double naive = (ax - px) * (by - py) - (ay - py) * (bx - px);
    if (Math.Sign(naive) != ExactPredicates.Orient2D(ax, ay, bx, by, px, py)) wrong++;
}
Console.WriteLine($"naive determinant has the wrong sign for {wrong} of {256 * 256} points");
 
double qx = 0.5 + 3 * ulp, qy = 0.5 + 1 * ulp;
Console.WriteLine($"one of them: naive {Math.Sign((ax - qx) * (by - qy) - (ay - qy) * (bx - qx))}, exact {ExactPredicates.Orient2D(ax, ay, bx, by, qx, qy)}");
 
// In-circle test for Delaunay triangulations: +1 inside, -1 outside, 0 on the circle.
Console.WriteLine($"InCircle (1,1): {ExactPredicates.InCircle(0, 0, 1, 0, 0, 1, 1, 1)}");
Console.WriteLine($"InCircle (0.5,0.5): {ExactPredicates.InCircle(0, 0, 1, 0, 0, 1, 0.5, 0.5)}");
 
// Segment classification.
Console.WriteLine(SegmentPredicates.Classify(0, 0, 2, 2, 0, 2, 2, 0));   // cross in the middle
Console.WriteLine(SegmentPredicates.Classify(0, 0, 2, 2, 2, 2, 3, 0));   // share an endpoint
Console.WriteLine(SegmentPredicates.Classify(0, 0, 2, 2, 1, 1, 3, 3));   // overlap on a line
Console.WriteLine(SegmentPredicates.Classify(0, 0, 1, 0, 0, 1, 1, 1));   // parallel, apart

Text

naive determinant has the wrong sign for 11972 of 65536 points
one of them: naive 0, exact -1
InCircle (1,1): 0
InCircle (0.5,0.5): 1
Proper
Touching
CollinearOverlap
Disjoint

What they return

Orient2D(ax, ay, bx, by, cx, cy): +1 if c is to the left of the directed line a→b (counterclockwise), −1 to the right, 0 on it. Orient3D does the same for a point and a plane.
InCircle(a, b, c, d): +1 if d is inside the circle through a, b and c (taken counterclockwise), −1 outside, 0 on it. InSphere is the 3D form. Here (1, 1) is exactly on the circle through (0, 0), (1, 0) and (0, 1).
SegmentPredicates.Classify returns a SegmentRelation: Disjoint, Proper (they cross at one interior point), Touching (they meet at an endpoint) or CollinearOverlap.

The arithmetic is adaptive: a fast floating-point filter answers when it can prove the sign, and exact expansion arithmetic finishes the rest, so the common case stays cheap.

Good to know

Coordinates must be finite and either zero or of magnitude in [2^-140, 2^200]. Out-of-range input may throw ArgumentOutOfRangeException.
Exact signs are what keep triangulations, convex hulls and polygon clipping consistent: two calls about the same points can never contradict each other.

The algorithm and its error bounds are in exact predicates.