Half-plane intersection

Why this matters when you build with AI

Half-plane intersection looks small: you have a list of lines, and you want the region that lies on the correct side of every line. AI tools often write it in two broken ways.

The first broken way skips the bounding box. Without a box, an unbounded region can look finished too early, and the code returns the wrong answer. I ran a box-free version on the region x >= 0 and y >= 0 (the top-right quadrant). It printed EMPTY. The right answer is UNBOUNDED.

The second broken way decides "is this corner on the line?" with a floating point tolerance, written eps. A tolerance is a small number, such as 1e-9, used to treat nearly equal values as equal. Floating point numbers are approximate, so the tolerance can give a wrong answer for nearly parallel lines. A cleaner rule is to use exact integer arithmetic for every sign test. Lesson 02 explains why floating point signs can be wrong. As a small reminder, I computed the orientation of the points (0, 0), (10^9 + 1, 10^9) and (10^9, 10^9 - 1) in two ways, and the two results have different signs:

exact -1 float 0.0

How to catch these errors: keep four small test files in your own folder (a triangle, an unbounded region, an empty region and a region with zero area). Then compare any version you write with a slow exact brute force on thousands of random small inputs. The steps are in the hands-on section below.

Terms used here

A directed line is a point p and a direction vector d with integer parts. The left side of that line is every point q with cross(d, q - p) >= 0. Here cross(u, v) is u.x * v.y - u.y * v.x, as defined in lesson 01.

A half-plane is the set of points on one side of a line, including the line itself. A constraint is a half-plane written as an inequality, such as x + y <= 4. The feasible region is the set of points that satisfy every constraint.

A region is bounded when it fits inside some circle. It is unbounded when it extends forever in at least one direction. In this lesson a region with zero area (a single point, a segment or a whole line) counts as empty, because the answer asks for an area.

A bounding box is four very large half-planes that enclose everything. They cut any unbounded region into a finite polygon, so the algorithm can treat every region the same way. A deque is a double-ended queue: you can add and remove items at both the front and the back.

The idea: one half-plane is one linear inequality

Start with a concrete example. Take three constraints: x >= 0, y >= 0 and x + y <= 4. Here is each one written as a directed line.

For x >= 0 use p = (0, 0) and d = (0, -1). Then cross(d, q - p) = 0 * y - (-1) * x = x, so the left side is x >= 0.

For y >= 0 use p = (0, 0) and d = (1, 0). Then cross(d, q - p) = 1 * y - 0 * x = y.

For x + y <= 4 use p = (4, 0) and d = (-1, 1). Then cross(d, q - p) = (-1) * y - 1 * (x - 4) = 4 - x - y, which is at least zero exactly when x + y <= 4.

The feasible region is the triangle with corners (0, 0), (4, 0) and (0, 4). Its area is 8.

Two facts drive the algorithm.

The first fact: walk around the boundary of a convex polygon counterclockwise, and the direction of each edge turns to the left a little at each corner. So the lines that form the edges appear in increasing order of their direction angle. This is why the algorithm sorts by angle.

The second fact: a line can only own an edge if the corner points at both of its ends lie inside the other half-planes. A line whose corners fall outside a newer line has no edge left, so it can be removed.

Sorting by angle needs no atan2. Put the directions into two halves: the upper half (angle from 0 up to, not including, pi) and the lower half. Within a half, one direction comes before another when the cross product of the first with the second is positive. The comparison is exact because it only uses integers.

The standard technique: sort by angle and keep a deque

The method works in five steps.

1. Add the four box lines. The code uses B = 10^7, which is far beyond any corner of a bounded answer.
2. Sort all lines by direction angle. If two lines point the same way, keep the more restrictive one. Line b is more restrictive than line a when the point of b lies strictly on the left of a.
3. Go through the lines in sorted order. Before adding line L, remove the last deque line while the corner of the last two lines is strictly outside L. Then remove the first deque line while the corner of the first two lines is strictly outside L. If the last line in the deque is parallel to L (only the opposite direction can remain, because same-direction lines were merged), the answer is empty. Otherwise add L to the back.
4. After all lines are processed, clean up the two ends again. Remove the last line while its corner with the line before it is outside the first line. Remove the first line while its corner with the line after it is outside the last line.
5. If fewer than three lines are left, the answer is empty. Otherwise the corners of each pair of neighbouring lines, including the pair (last, first), are the vertices of the answer polygon.

Here is a trace on the triangle, without the box. The box is not needed for this input, and the code's answer is the same with or without it.

The angles are: y >= 0 at 0 degrees, x + y <= 4 at 135 degrees, and x >= 0 at 270 degrees. Call them Y, S and X.

Push Y. The deque is [Y].

Push S. The deque has one line, so no corner test runs. The deque is [Y, S].

Add X. The corner of Y and S is (4, 0). Its x value is 4, which is not outside x >= 0, so nothing is removed. The deque is [Y, S, X].

Clean up. The corner of S and X is (0, 4). Its y value is 4, which is not outside y >= 0. The corner of Y and S is (4, 0), which is inside X. Nothing is removed.

The three corners are (4, 0), (0, 4) and (0, 0) (the last one is the corner of X and Y). The shoelace sum over these corners is 16, so the area is 8.

Now the empty case, again without the box: x <= 0 and x >= 1. The two lines point in opposite directions, (0, 1) at 90 degrees and (0, -1) at 270 degrees. The first one goes into the deque. When the second one is added, the back of the deque is parallel to it, and the code returns EMPTY. The reason is that two opposite lines can only be neighbours in the deque when every line between them was removed, which means the strip between them has no room left.

The unbounded case is the quadrant from earlier. Without the box, the deque ends with two lines, so the answer looks empty. With the box, the two box lines on the far side (x = B and y = B) survive with edges of positive length, which means the region reaches the edge of the box. That is the test for UNBOUNDED in the code: the answer has positive area, and one of the box lines is still in the final polygon.

Exact sign tests instead of a tolerance

The core question is: is the corner of lines A and B strictly outside the half-plane of C? The corner is P = pA + dA * t, where t = N / D, D = cross(dB, dA) and N = cross(dB, pB - pA). The test is cross(dC, P - pC) < 0.

Substitute the corner into the test and multiply by D. If D is made positive first, the test becomes c0 * D + N * c1 < 0, where c0 = cross(dC, pA - pC) and c1 = cross(dC, dA). Only integers are used, so there is no tolerance.

Worked example. Line A has p = (0, 0) and d = (1, 1), which is y >= x. Line B has p = (2, 0) and d = (0, 1), which is x <= 2. Their corner is P = (2, 2). Line C has p = (0, 3) and d = (1, 0), which is y >= 3.

Compute D = cross((0, 1), (1, 1)) = -1 and N = cross((0, 1), (2, 0)) = -2. D is negative, so flip both: D = 1, N = 2. Then c0 = cross((1, 0), (0, -3)) = -3 and c1 = cross((1, 0), (1, 1)) = 1. The value is c0 * D + N * c1 = -3 + 2 = -1, which is below zero, so the corner is outside. Check: the corner has y = 2, and C needs y >= 3. The answer agrees.

Overflow: with input coordinates and directions up to 100 and the box at 10^7, every product in the test stays below about 1.6 * 10^14, which fits in a 64-bit integer. If the inputs go up to 10^9, the products need 128-bit integers (__int128), which lesson 02 covers.

The area is the one place where the code uses floating point. The corners are converted to long double and the shoelace formula is applied. This works because the problem promises that a positive area is at least 1, so the code can treat any area below 0.5 as zero. A tiny positive area below 1 cannot be told from zero this way. An exact area needs rational arithmetic, and the problem statement avoids that case.

The code

This is a complete C++17 program. It reads n, then n lines px py dx dy, and prints EMPTY, UNBOUNDED or the area with six decimals. Save it as hpi.cpp, then compile with g++ -O2 -std=c++17 -Wall -Wextra -o hpi hpi.cpp.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
 
// Coordinates of input points and directions are at most 100 in size.
// The bounding box is at distance BIG, far beyond any vertex of a bounded answer.
const ll BIG = 10000000LL;
 
// The half-plane is the set of points on the left of the directed line
// from (px, py) along (dx, dy), including the line itself.
struct Line {
    ll px, py, dx, dy;
    int id;  // index in the input, or >= n for the four box lines
};
 
ll cross2(ll ax, ll ay, ll bx, ll by) { return ax * by - ay * bx; }
 
// True when the intersection point of lines A and B lies strictly outside C.
// A and B must not be parallel.
// The intersection is P = pA + dA * N / D with D = cross(dB, dA) and N = cross(dB, pB - pA).
// Writing E = cross(dC, P - pC) = c0 + (N / D) * c1, we get E * D = c0 * D + N * c1.
// After making D positive, E < 0 exactly when c0 * D + N * c1 < 0.
bool outsideInter(const Line& A, const Line& B, const Line& C) {
    ll D = cross2(B.dx, B.dy, A.dx, A.dy);
    ll N = cross2(B.dx, B.dy, B.px - A.px, B.py - A.py);
    if (D < 0) { D = -D; N = -N; }
    ll c0 = cross2(C.dx, C.dy, A.px - C.px, A.py - C.py);
    ll c1 = cross2(C.dx, C.dy, A.dx, A.dy);
    return c0 * D + N * c1 < 0;
}
 
// Upper half (angle in [0, pi)) is 0, lower half (angle in [pi, 2pi)) is 1.
int halfOf(const Line& L) { return (L.dy < 0 || (L.dy == 0 && L.dx < 0)) ? 1 : 0; }
 
// Strict angle order without atan2: compare halves first, then the cross product.
bool byAngle(const Line& a, const Line& b) {
    if (halfOf(a) != halfOf(b)) return halfOf(a) < halfOf(b);
    return cross2(a.dx, a.dy, b.dx, b.dy) > 0;
}
 
// Same direction as a? Two directions are the same when they are in the same half and cross to zero.
bool sameDir(const Line& a, const Line& b) {
    return halfOf(a) == halfOf(b) && cross2(a.dx, a.dy, b.dx, b.dy) == 0;
}
 
int main() {
    int n;
    cin >> n;
    vector<Line> in(n + 4);
    for (int i = 0; i < n; i++) {
        cin >> in[i].px >> in[i].py >> in[i].dx >> in[i].dy;
        in[i].id = i;
    }
    // The four box lines, each with the box on its left.
    in[n + 0] = {-BIG, -BIG, 1, 0, n + 0};
    in[n + 1] = {BIG, -BIG, 0, 1, n + 1};
    in[n + 2] = {BIG, BIG, -1, 0, n + 2};
    in[n + 3] = {-BIG, BIG, 0, -1, n + 3};
 
    sort(in.begin(), in.end(), byAngle);
 
    // Keep one line per direction: the most restrictive one.
    // Line b is more restrictive than a when b's point lies strictly left of a.
    vector<Line> L;
    for (const Line& l : in) {
        if (!L.empty() && sameDir(L.back(), l)) {
            Line& a = L.back();
            if (cross2(a.dx, a.dy, l.px - a.px, l.py - a.py) > 0) a = l;
        } else {
            L.push_back(l);
        }
    }
 
    int m = L.size();
    vector<Line> dq(m);
    int h = 0, t = 0;  // the deque is dq[h..t-1]
    bool empty = false;
    for (int i = 0; i < m && !empty; i++) {
        const Line& l = L[i];
        while (t - h >= 2 && outsideInter(dq[t - 2], dq[t - 1], l)) t--;
        while (t - h >= 2 && outsideInter(dq[h], dq[h + 1], l)) h++;
        if (t - h >= 1 && cross2(l.dx, l.dy, dq[t - 1].dx, dq[t - 1].dy) == 0) {
            // Adjacent lines are parallel: they point the same way only after dedupe, so they are opposite.
            empty = true;
            break;
        }
        dq[t++] = l;
    }
    while (!empty && t - h >= 3 && outsideInter(dq[t - 2], dq[t - 1], dq[h])) t--;
    while (!empty && t - h >= 3 && outsideInter(dq[h], dq[h + 1], dq[t - 1])) h++;
    if (empty || t - h < 3) { cout << "EMPTY\n"; return 0; }
 
    int k = t - h;
    vector<Line> poly(dq.begin() + h, dq.begin() + t);
    // Vertices: vertex i is the intersection of poly[i] and poly[i+1] (cyclic).
    vector<long double> X(k), Y(k);
    for (int i = 0; i < k; i++) {
        const Line& A = poly[i];
        const Line& B = poly[(i + 1) % k];
        ll D = cross2(B.dx, B.dy, A.dx, A.dy);
        if (D == 0) { cout << "EMPTY\n"; return 0; }  // wrap-around pair is parallel: degenerate
        ll N = cross2(B.dx, B.dy, B.px - A.px, B.py - A.py);
        long double s = (long double)N / (long double)D;
        X[i] = A.px + A.dx * s;
        Y[i] = A.py + A.dy * s;
    }
    long double twiceArea = 0;
    for (int i = 0; i < k; i++) {
        int j = (i + 1) % k;
        twiceArea += X[i] * Y[j] - X[j] * Y[i];
    }
    long double area = twiceArea / 2;
    // Positive area is promised to be at least 1, so anything below 0.5 is zero area.
    if (area < 0.5L) { cout << "EMPTY\n"; return 0; }
    for (const Line& l : poly)
        if (l.id >= n) { cout << "UNBOUNDED\n"; return 0; }
    cout << fixed << setprecision(6) << (double)area << "\n";
    return 0;
}

The program sorts n + 4 lines, so it takes O(n log n) time. The deque loop is O(n) in total, because each line is added once and removed at most once.

Running the three examples

Save the three inputs as separate files, then run ./hpi < file. The outputs below are from the compiled program on Linux with g++ 15.2.0.

Triangle, input 3, then 0 0 0 -1, 0 0 1 0 and 4 0 -1 1:

8.000000

Quadrant, input 2, then 0 0 0 -1 and 0 0 1 0:

UNBOUNDED

Empty, input 2, then 1 0 0 -1 and 0 0 0 1:

EMPTY

If you remove the box from the code (keep only the n input lines, and skip the four box lines), the quadrant prints EMPTY. That is the mistake from the first section, and you can reproduce it with your own copy.

Checking against a brute force

A brute force is a slow program that is easy to trust. This one lists every pair of lines (input lines and the four box lines), computes their corner with exact fractions, keeps the corners that satisfy every constraint, takes the convex hull of those corners and computes the area. If there are no corners the answer is EMPTY. If the area is zero the answer is EMPTY. If the area is positive and a corner lies on the box, the answer is UNBOUNDED.

from fractions import Fraction as F
 
BIG = 10000000
 
def cross(ax, ay, bx, by):
    return ax * by - ay * bx
 
def solve(lines):
    # lines: list of (px, py, dx, dy); region is cross(d, q - p) >= 0
    pts = []
    for i in range(len(lines)):
        for j in range(i + 1, len(lines)):
            pa, da = lines[i][:2], lines[i][2:]
            pb, db = lines[j][:2], lines[j][2:]
            D = cross(db[0], db[1], da[0], da[1])
            if D == 0:
                continue
            N = cross(db[0], db[1], pb[0] - pa[0], pb[1] - pa[1])
            t = F(N, D)
            x = pa[0] + da[0] * t
            y = pa[1] + da[1] * t
            ok = all(cross(l[2], l[3], x - l[0], y - l[1]) >= 0 for l in lines)
            if ok:
                pts.append((x, y))
    pts = sorted(set(pts))
    if not pts:
        return ("EMPTY", None)
    # convex hull (monotone chain) with exact fractions
    def half(o, a, b):
        return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
    if len(pts) < 3:
        return ("EMPTY", None)
    lower, upper = [], []
    for p in pts:
        while len(lower) >= 2 and half(lower[-2], lower[-1], p) <= 0:
            lower.pop()
        lower.append(p)
    for p in reversed(pts):
        while len(upper) >= 2 and half(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        upper.append(p)
    hull = lower[:-1] + upper[:-1]
    if len(hull) < 3:
        return ("EMPTY", None)
    a2 = F(0)
    for i in range(len(hull)):
        x1, y1 = hull[i]
        x2, y2 = hull[(i + 1) % len(hull)]
        a2 += x1 * y2 - x2 * y1
    area = a2 / 2
    if area == 0:
        return ("EMPTY", None)
    box_hit = any(abs(x) == BIG or abs(y) == BIG for (x, y) in pts)
    if box_hit:
        return ("UNBOUNDED", float(area))
    return ("AREA", float(area))

The brute force runs in O(n^3) time, which is fine for up to about 10 lines.

I tested the C++ program against this brute force with two random generators, which I ran in my scratch folder:

• Generic inputs: 1 to 7 lines, points with coordinates from -20 to 20, directions from -5 to 5. This run drew 1500 inputs. Nine were skipped because their positive area was below 1, which the promise excludes. The other 1491 had the outputs UNBOUNDED (823), a positive area (125) or EMPTY (543).
• Bounded and degenerate inputs: four axis-aligned lines with random offsets, plus lines through a few shared points (so many lines meet at one point), parallel lines and repeated directions. This run drew 2500 inputs. 190 had no lines and were discarded, and 38 were skipped for the same reason as above. The other 2272 had the outputs EMPTY (425), a positive area (1191) or UNBOUNDED (656).

The C++ output matched the brute force on every compared case, which is 3763 in total. This is a test, not a proof, and it used small coordinates. Large coordinates are covered by the overflow argument above, which I checked by hand and did not test on random inputs with coordinates near 10^9.

Common mistakes

No bounding box. The quadrant example above is the real symptom: the program prints EMPTY when the answer is UNBOUNDED. Test one unbounded input every time you change the code.

A tolerance in the corner test. A corner that lies exactly on a line can be judged inside or outside depending on the eps value. Use the exact integer test from the previous section, and keep eps out of the sign test.

Parallel lines not merged. If two lines point the same way and both stay in the deque side by side, the corner of those two lines does not exist, because D = 0. The formula then returns a meaningless value. The code above merges same-direction lines before the deque loop.

Printing 0.000000 for an empty region. An empty region or one with zero area should print EMPTY. The code prints EMPTY when the area is below 0.5.

Using a float for the sign test. The example from the first section shows the size of the error: the exact cross product is -1, and the float one is 0.0. A sign test with floats can say a corner is on a line when it is not.

Checklist

• Add a bounding box that is far larger than any corner of a bounded answer.
• Sort by angle with integer comparisons only, and merge lines with the same direction, keeping the most restrictive one.
• Make every sign test an exact integer test, and check the largest intermediate value against the integer size you use.
• Return EMPTY for an empty region or one with zero area, and UNBOUNDED when a box line survives in a region with positive area.
• Test the triangle, the quadrant, the empty pair and one degenerate case on every change.
• Compare with a brute force on random small inputs before you trust the code.

Prompt to use with an AI tool

Use this prompt when you ask an assistant to write the code. A good answer uses the bounding box, integer tests and the four cases above. Check the answer on the triangle, quadrant and empty examples, and run the brute force on random inputs.

Write a C++17 program that reads n, then n lines "px py dx dy". Each line is the
half-plane to the left of the directed line from (px, py) along (dx, dy), boundary
included. Coordinates and direction components are integers with absolute value at
most 100. Print EMPTY when the intersection has zero area, UNBOUNDED when it has
positive area and no finite bound, and otherwise the area with six decimals.
Use a bounding box of size 10^7. Sort by angle with integer comparisons only. Use no
floating point in any sign test. Explain the largest intermediate value you compute.
Show the three examples: the triangle (0,0),(4,0),(0,4) with area 8, the quadrant
(UNBOUNDED) and the pair x <= 0, x >= 1 (EMPTY).

Linear programming and the kernel

Linear programming (LP) asks for a point that maximises a linear function c . x among points that satisfy several linear inequalities. The feasible region is exactly what this lesson computes. The best point of an LP in two dimensions is always at a corner, so once the region is known, one pass over its corners gives the answer. Seidel (1991) gave a randomised method for LP in a fixed number of dimensions with expected running time linear in the number of constraints: R. Seidel, "Small-Dimensional Linear Programming and Convex Hulls Made Easy", Discrete and Computational Geometry 6(5), 1991, pages 423 to 434. This lesson does not implement it.

The kernel of a simple polygon is the set of points from which every point of the polygon is visible without crossing an edge. It equals the intersection of the inner half-planes of the polygon's edges. So this lesson's algorithm, applied to the edges of a polygon, computes its kernel.

Hands-on exercise

Compile hpi.cpp and run the three examples. Your output must match the three blocks above, except that the last digits of the area can differ on a machine where long double has a different width.

Then write a generator of random inputs in the same format and compare your program with the brute force above on at least 1000 inputs. Keep the generator small: between 1 and 7 lines, coordinates from -20 to 20, directions from -5 to 5. Report how many outputs were EMPTY, UNBOUNDED and a number, so that you can see whether every branch was tested. If one branch never appears, change the generator until it does.

Practice

1. Write the half-plane for the line through (1, 2) and (3, 5) where the left side is the side that contains (0, 5).

Show the answer

Use p = (1, 2) and d = (2, 3), which is the direction from (1, 2) to (3, 5). Then cross(d, q - p) = 2 * (y - 2) - 3 * (x - 1) = -3x + 2y + 1. The point (0, 5) gives 11, which is positive, so it is on the left. The half-plane is 3x - 2y - 1 <= 0.

2. Decide by hand whether x >= 0, y >= 0 and x + y >= 1 form a bounded region. Write the third line as a directed line.

Show the answer

Use p = (1, 0) and d = (1, -1), because cross(d, q - p) = x + y - 1. The region is the part of the first quadrant on or above the line x + y = 1. It has positive area and extends forever: the points (t, 0) for every t >= 1 satisfy all three constraints. So the answer is UNBOUNDED.

3. Explain why the code returns EMPTY when the last line in the deque is parallel to the new line.

Show the answer

Two lines that are parallel and point in opposite directions can only be neighbours in the deque after every line between them was removed. A line is removed only when its corner with a neighbour is outside a newer line. If the lines between them are all gone, the strip between the two opposite lines has no room for a corner, so the answer is empty or has zero area, and the code returns EMPTY. With the box, box lines sit between the two opposite lines in angle order, and they must be removed before the two meet. On the random inputs in the test above, no positive-area region reached this return. The full argument that this is always true is not written out here.

4. Change the code so it also prints how many input lines (not box lines) remain in the final polygon. Find an input where the answer is less than the number of lines.

Show the answer

After the cleanup, poly holds the final lines. Count the entries with l.id < n, and print that count before the area. Any input where a line is redundant works. For example, add x <= 10 to the triangle input (10 0 0 1). The triangle's constraint x + y <= 4 already implies x <= 4, so x <= 10 is redundant. The modified program prints 3 input lines remain for the four input lines, then 8.000000.

5. Compute the sign of c0 * D + N * c1 by hand for the lines A: p = (0, 0), d = (1, 1), B: p = (2, 0), d = (0, 1) and C: p = (0, 3), d = (1, 0), without flipping the sign of D. Check the answer against the corner (2, 2).

Show the answer

Here D = -1 and N = -2, and c0 = -3 and c1 = 1. The value is c0 * D + N * c1 = 3 - 2 = 1, which is positive. Without the flip, the test would say the corner is inside. This is wrong, because D is negative. Flipping D and N together gives the value -1, which is correct. This is why the code makes D positive before the test.

Summary

A half-plane is one linear inequality, written as the left side of a directed line. Sorting the lines by angle and keeping a deque gives the region in O(n log n) time. A bounding box makes unbounded regions show up as a finite polygon with a box edge, and exact integer sign tests remove the tolerance problem. Check every version against four cases and a brute force on random inputs.