Orientation tests and exact arithmetic

Why this matters when you build with AI

Ask an AI tool to sort points by angle and it usually calls atan2 and compares the results with a tolerance such as 1e-9. That code passes the sample. On a hidden test, two points with almost the same direction get the same angle, the tolerance treats them as equal, and the output order is wrong. The same tool also tends to write "is this point left of the line" with double arithmetic, which can say "on the line" for a point that is one unit away.

The fix is a rule you can check by hand: every geometric decision in this lesson is a sign of an integer expression. When the integers are exact, the sign is exact. This lesson gives that rule, a real sort and a real failing example, and a way to pick the integer type from the input size.

Two terms. The orientation of three points a, b, c tells whether walking from a to b and then to c makes a left turn, a right turn, or a straight line. A predicate is a yes-or-no or sign-valued test used inside an algorithm. Orientation is the basic predicate of computational geometry, and most other predicates in this course are built from it.

The idea: the orientation predicate

The orientation of (a, b, c) is the sign of the cross product of the two arrows from a: cross(b - a, c - a). From lesson 01, a positive cross product means the second arrow is counterclockwise from the first. So the sign tells the turn.

Worked example, with a = (0, 0), b = (4, 0), c = (1, 3). The arrows are b - a = (4, 0) and c - a = (1, 3). The cross product is 4 * 3 - 0 * 1 = 12, which is positive. So the three points make a left turn: going from a to b and then turning to c is a turn to the left. Swap b and c and the sign flips, because (1, 3) is now the first arrow: cross((1, 3), (4, 0)) = 1 * 0 - 3 * 4 = -12.

The same test can be written as a 3 by 3 determinant of the rows (ax, ay, 1), (bx, by, 1) and (cx, cy, 1). That form is useful as a checker, because it is a different formula that must give the same sign. The random test later in this lesson uses it.

Three rules follow from the sign:

• Positive: the points make a counterclockwise (left) turn.
• Negative: the points make a clockwise (right) turn.
• Zero: the three points lie on one line, in any order.

Zero is the case that floating point gets wrong. A point just to the left of a long line has a very small positive value in exact arithmetic, and a rounded calculation can return zero for it. With integers the value is never rounded, so zero means exactly collinear.

Sorting points by polar angle

The polar angle of a point is the angle of the arrow from the origin to that point, measured counterclockwise from the positive x axis. Sorting by polar angle is a standard step in many geometry algorithms, such as convex hulls and some sweep methods.

The standard exact method splits the plane into two halves, then compares within each half with the cross product.

• Half 0 holds the angles from 0 up to, but not including, pi. A point is in half 0 when y > 0, or when y = 0 and x > 0.
• Half 1 holds the angles from pi up to, but not including, 2 pi. Every other point is in half 1.

Within one half, every two arrows differ by less than pi, so the sign of their cross product gives the order. If cross(u, v) > 0, then v is counterclockwise from u, so u comes first. A cross product of 0 means the two points are in the same direction from the origin, and the tie is broken by the distance, nearest first.

Worked example, with four points (-1, -1), (2, 0), (0, 3) and (-3, 0). Their angles are 225, 0, 90 and 180 degrees. The sorted order is (2, 0) at 0, (0, 3) at 90, (-3, 0) at 180, then (-1, -1) at 225. The program below gives that order, and the first-run output is shown after the code.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
 
struct P {
    ll x, y;
};
 
// Sign of the cross product of (b - a) and (c - a): +1 counterclockwise, -1 clockwise, 0 collinear.
// Coordinates up to 1e9 in size keep every product below 4e18, so long long is enough.
int orient(P a, P b, P c) {
    ll v = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
    return (v > 0) - (v < 0);
}
 
// Half-plane index: 0 for angles in [0, pi), 1 for angles in [pi, 2*pi).
int half(P p) {
    return (p.y > 0 || (p.y == 0 && p.x > 0)) ? 0 : 1;
}
 
// Exact polar order around the origin, counterclockwise from angle 0.
// Points in the same direction are ordered by distance from the origin.
bool polar_less(P u, P v) {
    int hu = half(u), hv = half(v);
    if (hu != hv) return hu < hv;
    ll cr = u.x * v.y - u.y * v.x;     // > 0 means v is counterclockwise from u
    if (cr != 0) return cr > 0;
    return u.x * u.x + u.y * u.y < v.x * v.x + v.y * v.y;
}
 
int main() {
    int n;
    if (!(cin >> n)) return 0;
    vector<P> p(n);
    for (auto &q : p) cin >> q.x >> q.y;
    sort(p.begin(), p.end(), polar_less);
    for (auto &q : p) cout << q.x << " " << q.y << "\n";
    return 0;
}

The function orient is the predicate from the previous section. It is written with long long, and the comment states the bound, because the bound is the only thing that makes it correct. The function polar_less is a strict weak ordering: for any two points at most one of polar_less(u, v) and polar_less(v, u) is true, and equal directions compare as equivalent after the distance tie break. std::sort needs that property, and a comparator that breaks it is undefined behaviour.

Real output, compiled with g++ -std=c++17 -O2 -Wall -Wextra, on the six points (1, 0), (0, -1), (-1, 0), (2, 2), (1, 1), (0, 1) in that input order:

input:
6
1 0
0 -1
-1 0
2 2
1 1
0 1
 
output:
1 0
1 1
2 2
0 1
-1 0
0 -1

The points (1, 1) and (2, 2) have the same direction, and the nearer one, (1, 1), comes first. The points (1, 0) at 0 degrees and (0, -1) at 270 degrees show the wrap-around: the half-plane split puts the bottom point last.

The epsilon trap

Here is the failing version that many first attempts use. It computes the angle with atan2 and calls two angles equal when they differ by less than 1e-9. The program below takes two points that differ in direction by about 5 * 10^-19 radians, which is far smaller than 1e-9, and prints the output of both comparisons.

The two points are u = (10^9, 10^9 - 1) and v = (10^9 - 1, 10^9 - 2). The exact cross product is cross(u, v) = -1, so v is clockwise from u, and v must come first in the polar order.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
 
struct P {
    ll x, y;
};
 
// Floating-point angle, then an epsilon comparison, the way many first attempts are written.
double angle_of(P p) { return atan2((double)p.y, (double)p.x); }
bool angle_less_eps(P u, P v) {
    double a = angle_of(u), b = angle_of(v);
    const double eps = 1e-9;
    if (fabs(a - b) < eps) return false;   // treated as equal
    return a < b;
}
 
// Exact comparison from the previous section: the cross product sign.
bool exact_less(P u, P v) {
    ll cr = u.x * v.y - u.y * v.x;
    return cr > 0;   // assumes u and v lie in the same half-plane
}
 
int main() {
    P u{1000000000LL, 999999999LL};
    P v{999999999LL, 999999998LL};
    printf("atan2(u) = %.17g\n", angle_of(u));
    printf("atan2(v) = %.17g\n", angle_of(v));
    printf("difference of angles = %.3g\n", angle_of(u) - angle_of(v));
    ll cr = u.x * v.y - u.y * v.x;
    printf("exact cross(u, v) = %lld, so v is clockwise from u\n", cr);
    printf("double cross(u, v) = %.1f\n", (double)u.x * (double)v.y - (double)u.y * (double)v.x);
 
    vector<P> in = {u, v};
    vector<P> by_eps = in, by_exact = in;
    stable_sort(by_eps.begin(), by_eps.end(), angle_less_eps);
    stable_sort(by_exact.begin(), by_exact.end(), exact_less);
    printf("epsilon sort (input order u, v):");
    for (auto &q : by_eps) printf(" (%lld,%lld)", q.x, q.y);
    printf("\nexact sort   (input order u, v):");
    for (auto &q : by_exact) printf(" (%lld,%lld)", q.x, q.y);
    printf("\n");
    return 0;
}

Real output, compiled with g++ -std=c++17 -O2 -Wall -Wextra:

atan2(u) = 0.78539816289744835
atan2(v) = 0.78539816289744835
difference of angles = 0
exact cross(u, v) = -1, so v is clockwise from u
double cross(u, v) = 0.0
epsilon sort (input order u, v): (1000000000,999999999) (999999999,999999998)
exact sort   (input order u, v): (999999999,999999998) (1000000000,999999999)

The two angles are the same double, so the difference is exactly 0 and the tolerance treats them as equal. The epsilon sort keeps the input order, which puts u first. That order is wrong. The exact sort puts v first, which is right. The last two lines of the first block show the same failure in the cross product itself: the double calculation returns 0.0 where the true value is -1. Both numbers are about 10^18, and a double cannot store the low digits of that product.

The fault appears only when two directions are closer than the tolerance. Such pairs are rare in random data, which is why a sort can pass random tests and still fail a crafted one.

Choosing the number type

Every geometry solution uses one of four number types. Pick by the size of the largest product the code computes, not by the size of the input alone.

• long long holds values up to about 9.2 * 10^18. Use it when every product of two coordinate differences is below that. For differences up to 2 * 10^9, as in the cross product of points in [-10^9, 10^9], the largest cross product is 8 * 10^18, and long long is safe. Lesson 01 gives the derivation.
• __int128 holds values up to about 1.7 * 10^38. It is a GCC and Clang extension, not part of the standard, and it is available on 64-bit targets. Use it when coordinates reach about 10^18. The orientation of three such points needs products up to 4 * 10^36 and a difference up to 8 * 10^36, which fits. The demo below shows a case that long long cannot hold.
• double keeps 53 bits, about 16 decimal digits. It is fine for printing a final answer with a fixed number of decimals when the inputs are small and no sign decision depends on the result. It must not decide the sign of an orientation, or the order of two angles, when the inputs can be large.
• Rationals keep a numerator and a denominator as exact integers. They are the right tool when a point is a division, such as the intersection of two lines, and when two such points must be compared. The Python fractions module gives them directly, and lesson 03 uses them to compute intersection points.

The demo below uses three points with coordinates near 10^18. The true orientation value is -1. It prints the exact value with __int128, and the value that long long gives after the products wrap around.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using i128 = __int128;
 
string to_str(i128 v) {
    if (v == 0) return "0";
    bool neg = v < 0;
    unsigned __int128 u = neg ? -(unsigned __int128)v : (unsigned __int128)v;
    string s;
    while (u > 0) { s += char('0' + (int)(u % 10)); u /= 10; }
    if (neg) s += '-';
    reverse(s.begin(), s.end());
    return s;
}
 
int main() {
    // Three points with coordinates near 1e18. The true cross product of (b - a) and (c - a) is -1.
    ll ax = 0, ay = 0;
    ll bx = 1000000000000000000LL, by = 999999999999999999LL;
    ll cx = 999999999999999999LL, cy = 999999999999999998LL;
    // Products of two differences reach about 1e36, far above the 9.2e18 limit of long long.
    unsigned long long w = (unsigned long long)(bx - ax) * (unsigned long long)(cy - ay)
                         - (unsigned long long)(by - ay) * (unsigned long long)(cx - ax);
    i128 exact = (i128)(bx - ax) * (cy - ay) - (i128)(by - ay) * (cx - ax);
    printf("exact orientation value (__int128) = %s\n", to_str(exact).c_str());
    printf("sign from __int128 = %d\n", (exact > 0) - (exact < 0));
    printf("long long product wrapped (unsigned view) = %llu\n", w);
    return 0;
}

Real output, compiled with g++ -std=c++17 -O2 -Wall -Wextra:

exact orientation value (__int128) = -1
sign from __int128 = -1
long long product wrapped (unsigned view) = 18446744073709551615

The wrapped value 18446744073709551615 is 2^64 - 1, the 64-bit pattern of -1. It matches the exact answer only because this input happens to give -1. Products of this size carry no reliable information in long long, so the same code cannot be trusted on other inputs of this size. Python gives the same exact value, -1, from 10**18 * (10**18 - 2) - (10**18 - 1)**2.

For exact rational work, Python's fractions.Fraction keeps numerator and denominator. The snippet below runs in your browser:

Real output of that snippet:

1/3 > 333333333/1e9 : True
difference as a fraction: 1/3000000000
1/3 + 1/3 + 1/3 == 1: True
0.1 + 0.2 == 0.3 in floats: False
Fraction(1, 10) + Fraction(2, 10) == Fraction(3, 10): True
Fraction(0.1) is the exact binary value: False

The last line matters. Fraction(0.1) is built from the binary value stored in the double, which is not exactly one tenth. Construct fractions from integers, as in Fraction(1, 10), when you want the decimal value.

How this was tested

The polar_less sort was tested against an independent Python version on 1,500 random inputs. Each input had 1 to 12 points with coordinates in [-3, 3], [-10, 10] or [-10^9, 10^9], a point of the same direction was added in half of the trials, and points at the origin were excluded. The Python version uses a different key: (half, 0 or 1, Fraction(-x, y), squared distance), computed with exact fractions. The run had 0 mismatches.

The orient function was tested on 20,000 random triples, with coordinates in [-2, 2], [-5, 5] and [-10^9, 10^9]. The reference value was the sign of the 3 by 3 determinant, computed with Python's unbounded integers. The run had 0 mismatches. The code that ran these tests is shown below, because the test is part of the method. The triple test only checks the sign, not the full sort order, so it does not replace the sort test.

import random
import subprocess
from fractions import Fraction
 
random.seed(7)
 
def key(p):
    x, y = p
    h = 0 if (y > 0 or (y == 0 and x > 0)) else 1
    if y == 0:
        return (h, 0, Fraction(0), x * x + y * y)
    return (h, 1, Fraction(-x, y), x * x + y * y)
 
def run(binary, text):
    return subprocess.run([binary], input=text, capture_output=True, text=True, timeout=20).stdout
 
bad_sort = 0
for trial in range(1500):
    n = random.randint(1, 12)
    R = random.choice([3, 10, 10**9])
    pts = []
    for _ in range(n):
        while True:
            p = (random.randint(-R, R), random.randint(-R, R))
            if p != (0, 0):
                break
        pts.append(p)
    # Add a same-direction copy, to exercise ties.
    if random.random() < 0.5:
        x, y = pts[0]
        pts.append((2 * x, 2 * y))
        if abs(2 * x) > 10**9 or abs(2 * y) > 10**9:
            pts.pop()
    text = f"{len(pts)}\n" + "".join(f"{x} {y}\n" for x, y in pts)
    got = run("./sortpolar", text).split("\n")
    got = [tuple(map(int, line.split())) for line in got if line.strip()]
    want = sorted(pts, key=key)
    if got != want:
        bad_sort += 1
print("polar sort trials 1500, mismatches", bad_sort)
 
bad_orient = 0
triples = []
for _ in range(20000):
    R = random.choice([2, 5, 10**9])
    tri = [(random.randint(-R, R), random.randint(-R, R)) for _ in range(3)]
    triples.append(tri)
text = f"{len(triples)}\n" + "".join(
    f"{a[0]} {a[1]} {b[0]} {b[1]} {c[0]} {c[1]}\n" for a, b, c in triples)
got = [int(v) for v in run("./orient_test", text).split()]
for (a, b, c), g in zip(triples, got):
    # 3x3 determinant of [[ax, ay, 1], [bx, by, 1], [cx, cy, 1]], computed with unbounded integers.
    det = (a[0] * (b[1] - c[1]) - a[1] * (b[0] - c[0]) + (b[0] * c[1] - b[1] * c[0]))
    want = (det > 0) - (det < 0)
    if g != want:
        bad_orient += 1
print("orientation triples 20000, mismatches", bad_orient)

Real output of this test run:

polar sort trials 1500, mismatches 0
orientation triples 20000, mismatches 0

The orientation check used a small C++ program that reads triples and prints each sign, built from the same orient formula. The reference is the 3 by 3 determinant, which is a different formula. Both gave the same sign on every triple.

Complexity

orient takes O(1) time, with a fixed number of multiplications. Sorting n points with polar_less takes O(n log n) comparisons, and each comparison is O(1). The __int128 version takes the same number of steps, so the order of the algorithm does not change.

Common mistakes

Epsilon comparison of angles. The symptom is an order that is wrong only for pairs of points with almost the same direction. The example above shows it: the epsilon sort returns (10^9, 10^9 - 1) before (10^9 - 1, 10^9 - 2), which is the wrong order. Fix: sort with the exact comparator, and never use atan2 in a comparator.

Using doubles for the sign of a cross product. The symptom is a value of 0.0 where the exact value is -1, as in the output above. Fix: keep the cross product in integers. If the coordinates are too big for long long, use __int128.

Overflow in the middle of a formula. The symptom is a wrong sign that depends on the coordinates, with no error message. Fix: compute the size bound of each product, not of the final answer. The products are the problem, even when the difference fits.

A comparator that is not a strict weak ordering. The symptom is an order that changes between runs or a crash, because std::sort needs a consistent comparator. Typical causes are <= in place of <, and a missing tie break when two points share a direction. The sort above returns false for equal directions and uses the distance only as a tie break.

Forgetting the half-plane split. Comparing only by the cross product gives a wrong order whenever two points are more than pi apart in angle. For example, (-1, 0) and (1, 0) have cross product 0, which is a tie, but their angles are 180 and 0 degrees. The split by half is what separates them.

Taking the zero point. The point (0, 0) has no direction. polar_less still places it, in half 1 and before the other points of that half, but that place has no geometric meaning. Remove it before sorting by angle, or handle it as a separate case.

Checklist to keep

• Write every geometric sign as the sign of an integer expression.
• Use orient for turns and collinearity, and never a tolerance.
• Sort by polar angle with the half-plane split, then the cross product, then the distance.
• Check each product against the bound, not only the final answer.
• Use __int128 when coordinates reach about 10^18, and long long for 10^9.
• Use rationals when a point is a division, and compare them by cross multiplication or with Fraction.
• Test every predicate against an independent formula on random small inputs.

Hands-on exercise

Write orient and polar_less on your own machine. Generate 1,000 random inputs of 1 to 10 points, with coordinates in [-5, 5], and check your sorted output against a Python version that uses Fraction keys, as in the test above. Keep the seed fixed.

Then take the epsilon sort and run it on u = (n, n - 1) and v = (n - 1, n - 2) for n = 10^9. Record whether the order is right, and then check your exact sort on the same input. A good answer shows the epsilon sort giving the wrong order and the exact sort giving the right one. If your exact sort gives the epsilon order, check whether you used double anywhere in the comparison.

Prompt to try with an AI tool

Use this prompt, then check the answer with the steps above:

Write a C++17 function that sorts a vector of integer points (coordinates up to
1e9 in size) by counterclockwise polar angle around the origin, starting at
angle 0. Do not use atan2, double, or any tolerance. Ties (same direction) must
be ordered by distance from the origin. Explain in two sentences why the cross
product fits in 64 bits for these coordinates. Then give one input where an
atan2 sort with a 1e-9 tolerance gives the wrong order, with the expected output.

A good answer uses the half-plane split and a cross product, and includes the distance tie break. The bound it states should be near 8 * 10^18 for differences up to 2 * 10^9. Check the second part: the example must have two directions closer than 1e-9 radians. A common failure is a comparator that uses atan2 for the order and the cross product only for ties, which still breaks on the example. Reject any answer that uses a tolerance in the comparator.

Summary

Geometry decisions in integer arithmetic are exact signs. The orientation test is the cross product of two differences, the polar sort uses a half-plane split and that cross product, and the number type follows from the size of the largest product. Doubles and angles from atan2 can hide near-collinear cases, and the epsilon trap shows that they do.

Practice

These take a few minutes each. Try them on paper or in your scratch folder before you open the answers.

Practice 1: orientation by hand

Find the orientation of the triples ((0, 0), (2, 2), (4, 4)), ((0, 0), (2, 0), (2, 2)) and ((0, 0), (2, 2), (2, 0)).

Hint: Subtract the first point from the other two, and take the sign of x1 * y2 - y1 * x2.

Show the answer

The first triple has (2, 2) and (4, 4) as arrows, and 2 * 4 - 2 * 4 = 0, so it is collinear. The second has (2, 0) and (2, 2), and 2 * 2 - 0 * 2 = 4, which is a left turn with sign +1. The third has (2, 2) and (2, 0), and 2 * 0 - 2 * 2 = -4, which is a right turn with sign -1. The program with orient prints 0 1 -1 for these triples.

Practice 2: polar sort of four points

Sort the points (-1, -1), (2, 0), (0, 3) and (-3, 0) by counterclockwise polar angle from angle 0, using polar_less.

Hint: Put each point in half 0 or half 1 first, then order the points within each half by the cross product.

Show the answer

The order is (2, 0), (0, 3), (-3, 0), (-1, -1). The angles are 0, 90, 180 and 225 degrees. Points (0, 3) and (-3, 0) are in half 0 and half 1 respectively, because (-3, 0) has y = 0 and x < 0. The program prints 2 0, 0 3, -3 0, -1 -1 on separate lines.

Practice 3: a near-collinear pair with exact arithmetic

For a positive integer n, let u = (n, n - 1) and v = (n - 1, n - 2). Find cross(u, v) as a polynomial in n. Then find the value for n = 10^9, and say which of u and v comes first in polar order.

Hint: Expand n * (n - 2) - (n - 1) * (n - 1). The n^2 terms cancel.

Show the answer

The cross product is n * (n - 2) - (n - 1)^2 = n^2 - 2n - (n^2 - 2n + 1) = -1. For any n, including n = 10^9, the value is -1. A negative cross product means v is clockwise from u, so v comes first in polar order. The angle between them is about 1 / (2 * 10^18) radians, which a tolerance of 1e-9 cannot separate.

Practice 4: which number type fits

Points have coordinates up to 10^18 in size. Is __int128 enough to compute orient exactly? Give the largest product and the largest sum, and compare them with the __int128 limit of about 1.7 * 10^38.

Hint: The differences are at most 2 * 10^18. Each product of two differences is at most the square of that, and orient subtracts two products.

Show the answer

The differences are at most 2 * 10^18, so each product is at most 4 * 10^36. The difference of two products is at most 8 * 10^36. The __int128 limit is 2^127 - 1 = 170141183460469231731687303715884105727, about 1.7 * 10^38. The sum fits with a wide margin, so __int128 is enough for this input range. In Python, 8 * 10**36 < 2**127 - 1 is True.

Practice 5: exact comparison of two fractions

Decide whether 7/9 is bigger than 0.7778, with the decimal read exactly. Then find the exact difference.

Hint: Write 0.7778 as the fraction 7778/10000, and subtract with Fraction.

Show the answer

7/9 is smaller. The difference 7/9 - 7778/10000 is -1/45000. In Python, Fraction(7, 9) > Fraction(7778, 10000) gives False, and the difference prints as -1/45000. A double comparison gives the same answer here, but the exact result is what you need when the numbers are the output of a division inside a larger calculation.