⚙️AlgorithmAdvanced

Half-Plane Intersection

Key Points

  • •
    Half-plane intersection (HPI) computes the common region that satisfies many linear side-of-line constraints in the plane.
  • •
    The standard O(n n) algorithm sorts half-planes by angle and maintains a deque of active constraints, pruning ends when they are cut off by the new constraint.
  • •
    Represent each half-plane as an oriented line with its feasible side on the left; a point is inside if cross(v, x − p) ≥ 0.
  • •
    Handle parallel half-planes by keeping the most restrictive one (largest offset along the left normal) and use an EPS to fight floating-point noise.
  • •
    The intersection, if non-empty and bounded, is a convex polygon; if empty or unbounded, detect and report accordingly.
  • •
    HPI underlies 2D linear programming, convex polygon clipping, and Voronoi cell computation within a bounding box.
  • •
    Numerical robustness is crucial: normalize directions, compare angles with tolerance, and re-check endpoints after inserting all half-planes.
  • •
    The resulting polygon has at most n vertices, and evaluating linear objectives on it reduces to checking its vertices.

Prerequisites

  • →2D vectors and operations — HPI relies on dot and cross products, norms, and orientation tests.
  • →Sorting and custom comparators — Half-planes must be sorted by angle with careful tie-breaking for parallels.
  • →Deque data structure — The algorithm maintains a double-ended queue of active constraints to prune ends efficiently.
  • →Lines and intersections — Computing intersection points between lines is central to updating polygon corners.
  • →Floating-point arithmetic and EPS — Robustness requires tolerances, normalization, and guarding against near-parallel cases.
  • →Convex sets and polygons — Understanding that intersections of half-planes are convex underpins the deque trimming logic.
  • →Linear programming basics (2D) — HPI directly solves 2D LP and motivates the need for half-plane intersections.

Detailed Explanation

Tap terms for definitions

01Overview

Imagine many straight lines drawn on the plane, and for each line you choose one side as allowed (for example, “above this line” or “to the left of this line”). The half-plane intersection (HPI) problem asks: which points satisfy all these side-of-line constraints at once? The answer is either empty (no common point), unbounded (extends to infinity), or a bounded convex polygon. The classical algorithm solves HPI in O(n \log n) time by sorting constraints by angle and maintaining a double-ended queue (deque) of the “currently active” half-planes that still matter. When a new half-plane arrives, you remove from the back (or front) any planes whose pairwise intersection points lie outside the new constraint, thus trimming the feasible region from ends. This method exploits convexity: every cut keeps the feasible set convex, and each half-plane can remove only a contiguous portion from the current boundary. HPI is widely used in 2D linear programming (maximizing a linear function under linear inequalities), polygon clipping (e.g., Sutherland–Hodgman generalization), and computational geometry tasks like computing a site’s Voronoi cell within a bounding box. Implementations require careful numeric handling (EPS tolerances, parallel lines, and final cleanup).

02Intuition & Analogies

Think of each half-plane as a huge sheet of paper that covers the allowed side of its defining line. If you stack all these sheets together, the overlap is the set of points that every sheet covers—that’s the feasible region. Another analogy is slicing a block of clay with knives: every knife makes a straight cut and you keep only the clay on one chosen side. After all cuts, you either have no clay (empty), a slab that stretches infinitely (unbounded), or a neat convex polyhedron in 2D—a convex polygon. The algorithm’s deque acts like a smart way to keep track of the “outer rim” of the remaining clay. When a new cut arrives, only the two ends (front and back) of this rim can become invalid first, because the region stays convex and only the newest constraint can invalidate corners that stick out too far on either end. Sorting by angle is like organizing cuts by their directions so that parallel or near-parallel cuts are processed consistently. Handling parallel cuts is like deciding, among two knives with the same direction, which one actually matters: the one that bites deeper into the clay (more restrictive) determines the final surface. EPS tolerances are the safety gloves: without them, tiny numerical wiggles can make you misjudge whether a corner is inside or outside, leading to wrong conclusions about emptiness or polygon shape.

03Formal Definition

A half-plane in \(\) is the set \(H = \{ : a x + b y c\}\) (closed version), bounded by the line \(a x + b y = c\). Equivalently, with an oriented line given by point \(\) and direction vector \(\), define its left normal \( = (-, )\). The associated left half-plane is \(\{ : (, - ) 0\}\), which equals \(\{ : \}\). Given a set of half-planes \(\{\}_{i=1}^{n}\), the half-plane intersection is \(P = \). The classical HPI algorithm: (1) sort half-planes by their line angles \( = (, )\); (2) scan in this order while maintaining a deque of candidate half-planes; (3) before pushing a new half-plane \(H\), remove from the back/front any half-planes whose intersection points with their neighbors lie outside \(H\); (4) after all insertions, perform a final cleanup to ensure the first and last constraints are mutually feasible; (5) if at least three constraints remain, their pairwise intersections in order form the convex polygon boundary of \(P\). The method relies on convexity and angle order to ensure amortized O(1) deque pops per constraint and thus O(n) after sorting.

04When to Use

Use HPI when you need the exact feasible region of many linear inequalities in 2D and can tolerate floating-point computation: (1) 2D linear programming—maximize or minimize a linear functional (\mathbf{c}\cdot\mathbf{x}) subject to (A\mathbf{x} \le \mathbf{b}). The optimum, if bounded, lies on a vertex of the HPI polygon. (2) Voronoi cells within a bounding box—each competitor site contributes a half-plane defined by a perpendicular bisector; intersecting them yields the cell polygon. (3) Convex polygon clipping—intersect an initial convex polygon (or large bounding box) with a set of half-planes to model visibility, shadow volumes, or collision regions. (4) Feasibility checks—quickly detect emptiness of constraints in 2D. (5) Robust geometric preprocessing—trim search spaces for further algorithms (e.g., rotating calipers, support mapping) by enforcing linear bounds. Do not use HPI if constraints are non-linear, in high-dimensional spaces (2D-specific algorithm), or if exact rational arithmetic is required (floating-point HPI needs careful tolerance). For 3D and higher, consider half-space intersection using different data structures or linear programming solvers.

⚠️Common Mistakes

  • Inconsistent orientation: mixing “left-of-line” and “right-of-line” conventions causes silent errors. Pick one convention (e.g., inside if cross(v, x−p) ≥ 0) and stick to it for all half-planes.
- Mishandling parallel lines: if two lines have the same angle, keep only the most restrictive half-plane (largest offset along the left normal). Keeping both can make the deque logic fail or lead to incorrect unbounded/empty detection.
- No EPS or too-tight EPS: floating-point comparisons like cross(...) ≥ 0 must use a tolerance (e.g., 1e−12 with long double). Too small EPS fails under noise; too large EPS distorts geometry.
- Forgetting final cleanup: after processing all half-planes, you must re-check the ends of the deque against the first/last constraints; otherwise, stale corners might remain.
- Wrong intersection formula or not guarding near-parallel denominators: always check |cross(v1, v2)| against EPS.
- Assuming boundedness: the intersection can be unbounded; if you need a polygon, first clip with a big bounding rectangle (or add constraints) to bound it.
- Angle tie-breaking by atan2 only: equal angles need extra handling; normalizing direction vectors reduces surprises.
- Returning degenerate outputs: if fewer than three constraints survive, the intersection is empty or a line/point—treat as empty polygon for many applications.

Key Formulas

2D Cross Product

Explanation: This scalar measures the signed area of the parallelogram spanned by a and b. Its sign tells if b is to the left (positive) or right (negative) of a.

Left Normal

Explanation: Rotating a direction vector v by +90° gives the left normal n. The half-plane left of v satisfies n·.

Line-Line Intersection Parameter

Explanation: For two non-parallel lines, the parameter t locates the intersection on the first line. The denominator being zero indicates parallel lines.

Half-plane Equivalence

Explanation: Using v = (b, −a) and choosing p so that a + b makes the algebraic half-plane equal to the geometric left-of-line condition.

Line Angle

Explanation: The angle of a direction vector is used to sort lines. Equal or nearly equal angles imply parallelism that needs special handling.

Shoelace Area

Explanation: Given the polygon vertices in cyclic order (with ), this computes the polygon’s area. Useful to detect emptiness or degenerate outputs.

HPI Time Complexity

Explanation: Sorting the n half-planes by angle costs O(n log n). The deque phase is linear because each half-plane is inserted once and popped at most once from each end.

HPI Space Complexity

Explanation: We store the sorted half-planes, a deque of active half-planes, and O(n) intersection points, all linear in n.

Voronoi Half-plane

Explanation: Being closer to s than to t can be written as a linear inequality, yielding a half-plane whose boundary is the perpendicular bisector of s and t.

Support Function and LP Optimum

Explanation: Maximizing a linear objective over the convex polygon P equals its support function in direction c and is attained at a vertex of P.

Complexity Analysis

The half-plane intersection algorithm runs in O(n log n) time and O(n) space. The dominant cost is sorting the n half-planes by their line angles, which is O(n log n). After sorting, the deque maintenance is linear: each half-plane is inserted once and can be removed at most once from the front and once from the back, leading to O(n) total pops. Computing an intersection between two lines is O(1), and we do it only when pushing/pruning neighbors, so intersection computations are also O(n) in total. Constructing the final polygon from the deque—computing pairwise intersections of adjacent lines, including the wrap-around—is O(k), where is the number of surviving half-planes (and thus the number of polygon vertices). The space usage is linear: we store the input list, possibly a filtered list to remove duplicates in angle, the deque of active half-planes (size O(n)), and at the end the list of polygon vertices (O(n)). In floating-point implementations, using long double modestly increases constant factors but does not change asymptotic complexity. If you add a bounding box (four extra constraints) to guarantee boundedness, the complexity stays O(n log n). Note that post-processing steps such as polygon area computation or evaluating a linear objective over vertices are O(k) and usually negligible compared to the sorting step. In adversarial numeric conditions, extra robustness checks may slightly increase constant factors but preserve the same asymptotic bounds.

Code Examples

Half-Plane Intersection Core: Compute Intersection Polygon (with robustness)
1#include <bits/stdc++.h>
2using namespace std;
3
4struct Point {
5    long double x, y;
6    Point() : x(0), y(0) {}
7    Point(long double _x, long double _y) : x(_x), y(_y) {}
8    Point operator+(const Point& o) const { return Point(x + o.x, y + o.y); }
9    Point operator-(const Point& o) const { return Point(x - o.x, y - o.y); }
10    Point operator*(long double k) const { return Point(x * k, y * k); }
11};
12
13static const long double EPS = 1e-12L;
14
15long double dot(const Point& a, const Point& b) { return a.x * b.x + a.y * b.y; }
16long double cross(const Point& a, const Point& b) { return a.x * b.y - a.y * b.x; }
17long double norm(const Point& a) { return sqrtl(dot(a, a)); }
18Point rotCCW90(const Point& a) { return Point(-a.y, a.x); }
19Point rotCW90(const Point& a) { return Point(a.y, -a.x); }
20
21struct HalfPlane {
22    Point p;      // a point on the line
23    Point v;      // direction (unit length)
24    long double ang; // atan2(v.y, v.x)
25    HalfPlane() {}
26    HalfPlane(Point _p, Point _v) {
27        long double len = norm(_v);
28        if (len == 0) throw runtime_error("Zero direction in HalfPlane");
29        v = Point(_v.x / len, _v.y / len);
30        p = _p;
31        ang = atan2l(v.y, v.x);
32    }
33    // Inside check: point x is inside if it is to the left of oriented line (p, v)
34    bool contains(const Point& x) const {
35        return cross(v, x - p) >= -EPS;
36    }
37    // For ax + by <= c, produce left-of-line half-plane with v = rotateCW90(n)
38    static HalfPlane from_le(long double a, long double b, long double c) {
39        Point n(a, b);
40        Point v = rotCW90(n); // so that left normal equals n
41        long double n2 = dot(n, n);
42        if (n2 == 0) throw runtime_error("Invalid inequality normal");
43        // choose p on the line n·p = c
44        Point p = n * (c / n2);
45        return HalfPlane(p, v);
46    }
47};
48
49// Intersection point of lines of two half-planes (assumed non-parallel)
50Point line_intersection(const HalfPlane& A, const HalfPlane& B) {
51    // Solve A.p + A.v * t = B.p + B.v * s
52    Point u = B.p - A.p;
53    long double den = cross(B.v, A.v);
54    // Caller should ensure |den| > EPS
55    long double t = cross(B.v, u) / den;
56    return A.p + A.v * t;
57}
58
59bool halfplane_cmp(const HalfPlane& a, const HalfPlane& b) {
60    if (fabsl(a.ang - b.ang) > 1e-15L) return a.ang < b.ang;
61    // For equal angles, keep the most restrictive: larger offset along left normal n = rotCCW90(v)
62    Point na = rotCCW90(a.v);
63    long double da = dot(na, a.p);
64    long double db = dot(na, b.p); // same normal because angles equal
65    return da > db; // sort descending so first is most restrictive
66}
67
68// Core HPI: returns true with polygon if non-empty and bounded by surviving half-planes; false if empty.
69bool halfplane_intersection(vector<HalfPlane> hp, vector<Point>& poly) {
70    // 1) Sort by angle; for equal angles, most restrictive first
71    sort(hp.begin(), hp.end(), halfplane_cmp);
72
73    // 2) Remove exact duplicates of angle keeping the first (most restrictive)
74    vector<HalfPlane> filt; filt.reserve(hp.size());
75    for (size_t i = 0; i < hp.size(); ++i) {
76        if (i == 0 || fabsl(hp[i].ang - hp[i-1].ang) > 1e-15L) filt.push_back(hp[i]);
77        // else skip less restrictive parallels
78    }
79
80    deque<HalfPlane> dq;
81    deque<Point> pts; // intersection points between consecutive half-planes in dq
82
83    auto is_outside = [&](const HalfPlane& H, const Point& X) {
84        return !H.contains(X);
85    };
86
87    auto push_halfplane = [&](const HalfPlane& H) {
88        // Prune from back while last corner is outside H
89        while (!dq.empty() && !pts.empty() && is_outside(H, pts.back())) {
90            dq.pop_back();
91            pts.pop_back();
92        }
93        // Prune from front while first corner is outside H
94        while (!dq.empty() && !pts.empty() && is_outside(H, pts.front())) {
95            dq.pop_front();
96            pts.pop_front();
97        }
98        // If last is parallel to H, keep the more restrictive (H is already most restrictive per sort)
99        if (!dq.empty() && fabsl(cross(dq.back().v, H.v)) < EPS) {
100            // Replace last with H (since H is more restrictive due to sorting)
101            dq.pop_back();
102            if (!pts.empty()) pts.pop_back();
103            // Also check parallel with front if dq size==1 and equal angle
104            if (!dq.empty() && fabsl(cross(dq.front().v, H.v)) < EPS) {
105                dq.pop_front();
106                if (!pts.empty()) pts.pop_front();
107            }
108        }
109        if (!dq.empty()) {
110            // Compute intersection with previous last
111            if (fabsl(cross(dq.back().v, H.v)) < EPS) {
112                // Parallel after clean-up shouldn't happen, but guard
113                if (!H.contains(dq.back().p)) return; // will be emptied later
114                // Otherwise, ignore as duplicate
115            } else {
116                Point ip = line_intersection(dq.back(), H);
117                pts.push_back(ip);
118            }
119        }
120        dq.push_back(H);
121    };
122
123    for (const auto& H : filt) push_halfplane(H);
124
125    // Final cleanup: ensure ends are consistent
126    while (dq.size() >= 2 && !pts.empty() && !dq.front().contains(pts.back())) {
127        dq.pop_back();
128        pts.pop_back();
129    }
130    while (dq.size() >= 2 && !pts.empty() && !dq.back().contains(pts.front())) {
131        dq.pop_front();
132        pts.pop_front();
133    }
134
135    if (dq.size() < 3) {
136        poly.clear();
137        return false; // empty or unbounded in our representation
138    }
139
140    // Recompute all intersection points in order (including wrap-around)
141    vector<HalfPlane> hs(dq.begin(), dq.end());
142    poly.clear();
143    for (size_t i = 0; i < hs.size(); ++i) {
144        size_t j = (i + 1) % hs.size();
145        long double den = fabsl(cross(hs[i].v, hs[j].v));
146        if (den < EPS) return false; // nearly parallel neighbors => treat as empty/degenerate
147        poly.push_back(line_intersection(hs[i], hs[j]));
148    }
149    return true;
150}
151
152long double polygon_area(const vector<Point>& P) {
153    long double s = 0;
154    int n = (int)P.size();
155    for (int i = 0; i < n; ++i) {
156        int j = (i + 1) % n;
157        s += P[i].x * P[j].y - P[i].y * P[j].x;
158    }
159    return fabsl(s) * 0.5L;
160}
161
162int main() {
163    ios::sync_with_stdio(false);
164    cin.tie(nullptr);
165
166    // Example: Intersect five half-planes forming a pentagon-like region
167    // Constraints in ax + by <= c form
168    vector<HalfPlane> hp;
169    hp.push_back(HalfPlane::from_le( 1, 0, 5));  //  x <= 5
170    hp.push_back(HalfPlane::from_le(-1, 0, 3));  // -x <= 3  => x >= -3
171    hp.push_back(HalfPlane::from_le( 0, 1, 4));  //  y <= 4
172    hp.push_back(HalfPlane::from_le( 0,-1, 2));  // -y <= 2  => y >= -2
173    hp.push_back(HalfPlane::from_le( 1, 1, 6));  //  x + y <= 6
174
175    vector<Point> poly;
176    bool ok = halfplane_intersection(hp, poly);
177
178    if (!ok) {
179        cout << "Empty or unbounded intersection.\n";
180        return 0;
181    }
182
183    cout.setf(std::ios::fixed); cout << setprecision(10);
184    cout << "Vertices (counter-clockwise):\n";
185    for (auto &pt : poly) cout << pt.x << " " << pt.y << "\n";
186    cout << "Area: " << polygon_area(poly) << "\n";
187    return 0;
188}
189

This program implements a robust half-plane intersection. Half-planes are created from linear inequalities ax + by ≤ c using a left-of-line convention. The algorithm sorts by angle, removes less restrictive parallels, and maintains a deque of active constraints while pruning ends whose corner points fall outside the current half-plane. After a final cleanup, consecutive intersections form the output convex polygon. The demo constructs five constraints, computes the intersection polygon, prints vertices, and its area.

Time: O(n log n) due to sorting; deque maintenance and intersection computations are O(n).Space: O(n) for storing constraints, deque, and output polygon.
Voronoi Cell of a Site Clipped by a Bounding Box via HPI
1#include <bits/stdc++.h>
2using namespace std;
3
4struct Point { long double x,y; };
5static const long double EPS = 1e-12L;
6
7long double dot(const Point&a,const Point&b){return a.x*b.x+a.y*b.y;}
8long double cross(const Point&a,const Point&b){return a.x*b.y-a.y*b.x;}
9long double norm(const Point&a){return sqrtl(dot(a,a));}
10Point rotCW90(const Point&a){return Point{a.y,-a.x};}
11Point rotCCW90(const Point&a){return Point{-a.y,a.x};}
12
13struct HalfPlane{
14    Point p,v; long double ang;
15    HalfPlane(){}
16    HalfPlane(Point _p, Point _v){ long double L=norm(_v); v={_v.x/L,_v.y/L}; p=_p; ang=atan2l(v.y,v.x);}    
17    static HalfPlane from_le(long double a,long double b,long double c){
18        Point n{a,b}; Point v=rotCW90(n); long double n2=dot(n,n); Point p{n.x*(c/n2), n.y*(c/n2)}; return HalfPlane(p,v);
19    }
20    bool contains(const Point&x)const{ return cross(v, Point{x.x-p.x, x.y-p.y}) >= -EPS; }
21};
22
23Point isect(const HalfPlane&A,const HalfPlane&B){ Point u{B.p.x-A.p.x, B.p.y-A.p.y}; long double den=cross(B.v, A.v); long double t=cross(B.v,u)/den; return Point{A.p.x + A.v.x*t, A.p.y + A.v.y*t}; }
24
25bool halfplane_cmp(const HalfPlane&a,const HalfPlane&b){ if(fabsl(a.ang-b.ang)>1e-15L) return a.ang<b.ang; Point n=rotCCW90(a.v); long double da=n.x*a.p.x+n.y*a.p.y; long double db=n.x*b.p.x+n.y*b.p.y; return da>db; }
26
27bool HPI(vector<HalfPlane> hp, vector<Point>& poly){
28    sort(hp.begin(), hp.end(), halfplane_cmp);
29    vector<HalfPlane> f; f.reserve(hp.size());
30    for(size_t i=0;i<hp.size();++i) if(i==0 || fabsl(hp[i].ang-hp[i-1].ang)>1e-15L) f.push_back(hp[i]);
31    deque<HalfPlane> dq; deque<Point> pts;
32    auto outside=[&](const HalfPlane&H,const Point&X){return !H.contains(X);} ;
33    auto pushH=[&](const HalfPlane&H){
34        while(!dq.empty() && !pts.empty() && outside(H, pts.back())){ dq.pop_back(); pts.pop_back(); }
35        while(!dq.empty() && !pts.empty() && outside(H, pts.front())){ dq.pop_front(); pts.pop_front(); }
36        if(!dq.empty() && fabsl(cross(dq.back().v,H.v))<EPS){ dq.pop_back(); if(!pts.empty()) pts.pop_back(); }
37        if(!dq.empty() && fabsl(cross(dq.back().v,H.v))>=EPS){ pts.push_back(isect(dq.back(), H)); }
38        dq.push_back(H);
39    };
40    for(auto&H: f) pushH(H);
41    while(dq.size()>=2 && !pts.empty() && !dq.front().contains(pts.back())){ dq.pop_back(); pts.pop_back(); }
42    while(dq.size()>=2 && !pts.empty() && !dq.back().contains(pts.front())){ dq.pop_front(); pts.pop_front(); }
43    if(dq.size()<3){ poly.clear(); return false; }
44    vector<HalfPlane> hs(dq.begin(), dq.end()); poly.clear();
45    for(size_t i=0;i<hs.size();++i){ size_t j=(i+1)%hs.size(); long double den=fabsl(cross(hs[i].v, hs[j].v)); if(den<EPS) return false; poly.push_back(isect(hs[i], hs[j])); }
46    return true;
47}
48
49long double area(const vector<Point>&P){ long double s=0; int n=P.size(); for(int i=0;i<n;++i){ int j=(i+1)%n; s+=P[i].x*P[j].y-P[i].y*P[j].x; } return fabsl(s)*0.5L; }
50
51int main(){
52    ios::sync_with_stdio(false); cin.tie(nullptr);
53    // Input: site s and competitor sites t_i. Build Voronoi cell of s clipped by a bounding box.
54    Point s; int m; cin >> s.x >> s.y >> m; vector<Point> t(m); for(int i=0;i<m;++i) cin >> t[i].x >> t[i].y;
55    long double B = 1e6; // bounding box size (adjust as needed)
56    vector<HalfPlane> hp;
57    // Bounding box: -B <= x <= B, -B <= y <= B
58    hp.push_back(HalfPlane::from_le( 1, 0,  B)); // x <=  B
59    hp.push_back(HalfPlane::from_le(-1, 0,  B)); // -x <= B => x >= -B
60    hp.push_back(HalfPlane::from_le( 0, 1,  B)); // y <=  B
61    hp.push_back(HalfPlane::from_le( 0,-1,  B)); // -y <= B => y >= -B
62
63    // For each competitor t: ||x-s||^2 <= ||x-t||^2 => 2(t-s)·x <= ||t||^2 - ||s||^2
64    long double ss = s.x*s.x + s.y*s.y;
65    for(const auto& u: t){
66        long double tt = u.x*u.x + u.y*u.y;
67        long double A = 2*(u.x - s.x);
68        long double Bc = 2*(u.y - s.y);
69        long double C = tt - ss;
70        hp.push_back(HalfPlane::from_le(A, Bc, C));
71    }
72
73    vector<Point> cell;
74    if(!HPI(hp, cell)){
75        cout << "Empty cell or numerical failure\n";
76        return 0;
77    }
78    cout.setf(std::ios::fixed); cout<<setprecision(8);
79    cout << cell.size() << "\n";
80    for(auto&p:cell) cout << p.x << " " << p.y << "\n";
81    cout << "Area: " << area(cell) << "\n";
82    return 0;
83}
84

This program computes the Voronoi cell of a given site s, clipped to a large bounding box. For each competitor t, the inequality ||x−s||^2 ≤ ||x−t||^2 becomes a linear half-plane 2(t−s)·x ≤ ||t||^2 − ||s||^2. Intersecting these with the bounding box via HPI yields the polygonal cell. The output lists the vertices and area.

Time: O((m+4) log (m+4)) for sorting all half-planes; intersection and polygon assembly are linear.Space: O(m) for constraints, deque, and polygon vertices.
2D Linear Programming via HPI: Maximize c·x subject to Ax ≤ b
1#include <bits/stdc++.h>
2using namespace std;
3
4struct Point { long double x,y; };
5static const long double EPS = 1e-12L;
6
7long double dot(const Point&a,const Point&b){return a.x*b.x+a.y*b.y;}
8long double cross(const Point&a,const Point&b){return a.x*b.y-a.y*b.x;}
9long double norm(const Point&a){return sqrtl(dot(a,a));}
10Point rotCW90(const Point&a){return Point{a.y,-a.x};}
11Point rotCCW90(const Point&a){return Point{-a.y,a.x};}
12
13struct HalfPlane{
14    Point p,v; long double ang;
15    HalfPlane(){}
16    HalfPlane(Point _p, Point _v){ long double L=norm(_v); v={_v.x/L,_v.y/L}; p=_p; ang=atan2l(v.y,v.x);}    
17    static HalfPlane from_le(long double a,long double b,long double c){
18        Point n{a,b}; Point v=rotCW90(n); long double n2=dot(n,n); Point p{n.x*(c/n2), n.y*(c/n2)}; return HalfPlane(p,v);
19    }
20    bool contains(const Point&x)const{ return cross(v, Point{x.x-p.x, x.y-p.y}) >= -EPS; }
21};
22
23Point isect(const HalfPlane&A,const HalfPlane&B){ Point u{B.p.x-A.p.x, B.p.y-A.p.y}; long double den=cross(B.v, A.v); long double t=cross(B.v,u)/den; return Point{A.p.x + A.v.x*t, A.p.y + A.v.y*t}; }
24
25bool halfplane_cmp(const HalfPlane&a,const HalfPlane&b){ if(fabsl(a.ang-b.ang)>1e-15L) return a.ang<b.ang; Point n=rotCCW90(a.v); long double da=n.x*a.p.x+n.y*a.p.y; long double db=n.x*b.p.x+n.y*b.p.y; return da>db; }
26
27bool HPI(vector<HalfPlane> hp, vector<Point>& poly){
28    sort(hp.begin(), hp.end(), halfplane_cmp);
29    vector<HalfPlane> f; f.reserve(hp.size());
30    for(size_t i=0;i<hp.size();++i) if(i==0 || fabsl(hp[i].ang-hp[i-1].ang)>1e-15L) f.push_back(hp[i]);
31    deque<HalfPlane> dq; deque<Point> pts;
32    auto outside=[&](const HalfPlane&H,const Point&X){return !H.contains(X);} ;
33    auto pushH=[&](const HalfPlane&H){
34        while(!dq.empty() && !pts.empty() && outside(H, pts.back())){ dq.pop_back(); pts.pop_back(); }
35        while(!dq.empty() && !pts.empty() && outside(H, pts.front())){ dq.pop_front(); pts.pop_front(); }
36        if(!dq.empty() && fabsl(cross(dq.back().v,H.v))<EPS){ dq.pop_back(); if(!pts.empty()) pts.pop_back(); }
37        if(!dq.empty() && fabsl(cross(dq.back().v,H.v))>=EPS){ pts.push_back(isect(dq.back(), H)); }
38        dq.push_back(H);
39    };
40    for(auto&H: f) pushH(H);
41    while(dq.size()>=2 && !pts.empty() && !dq.front().contains(pts.back())){ dq.pop_back(); pts.pop_back(); }
42    while(dq.size()>=2 && !pts.empty() && !dq.back().contains(pts.front())){ dq.pop_front(); pts.pop_front(); }
43    if(dq.size()<3){ poly.clear(); return false; }
44    vector<HalfPlane> hs(dq.begin(), dq.end()); poly.clear();
45    for(size_t i=0;i<hs.size();++i){ size_t j=(i+1)%hs.size(); long double den=fabsl(cross(hs[i].v, hs[j].v)); if(den<EPS) return false; poly.push_back(isect(hs[i], hs[j])); }
46    return true;
47}
48
49int main(){
50    ios::sync_with_stdio(false); cin.tie(nullptr);
51    // Input: n constraints, objective c, and an optional bounding box to ensure bounded optimum.
52    int n; cin >> n; vector<long double> a(n), b(n), cst(n);
53    for(int i=0;i<n;++i) cin >> a[i] >> b[i] >> cst[i];
54    Point c; cin >> c.x >> c.y; // objective: maximize c·x
55
56    vector<HalfPlane> hp; hp.reserve(n+4);
57    for(int i=0;i<n;++i) hp.push_back(HalfPlane::from_le(a[i], b[i], cst[i]));
58
59    // If the intersection may be unbounded, add a large bounding box (tunable parameter):
60    long double B = 1e6; // adjust based on problem scale
61    hp.push_back(HalfPlane::from_le( 1, 0,  B));
62    hp.push_back(HalfPlane::from_le(-1, 0,  B));
63    hp.push_back(HalfPlane::from_le( 0, 1,  B));
64    hp.push_back(HalfPlane::from_le( 0,-1,  B));
65
66    vector<Point> poly;
67    if(!HPI(hp, poly)){
68        cout << "Infeasible or numerically degenerate\n"; return 0;
69    }
70
71    // Maximize c·x over convex polygon: best at a vertex
72    long double best = -numeric_limits<long double>::infinity();
73    Point arg; for(const auto& p: poly){ long double val = dot(c, p); if(val > best){ best = val; arg = p; } }
74
75    cout.setf(std::ios::fixed); cout<<setprecision(10);
76    cout << "opt_value " << best << "\n";
77    cout << "argmax " << arg.x << " " << arg.y << "\n";
78    return 0;
79}
80

This code solves a 2D linear program by converting all constraints Ax ≤ b into half-planes and intersecting them. To guarantee a bounded feasible region, it adds a large bounding box (common in practice if constraints don’t bound all directions). The optimal value of c·x over a convex polygon is attained at a vertex, so we simply scan all polygon vertices and pick the maximum.

Time: O((n+4) log (n+4)) for HPI; selecting the best vertex is O(k) where k ≤ n.Space: O(n) for constraints, deque, and polygon vertices.