带无穷判断的半平面交

网上的半平面交一般不能区分半平面交不存在和面积无穷的情况。这里给出一个方法。

为了叙述的简洁，先假设所有直线都不共线，并且一些边界情况不予讨论，详情请见代码。

如果半平面交面积无穷，那么一定存在一个点，使得从这个点朝某个方向射去能够无限延伸。假设存在一条有向直线，其方向向量为v，极角为$\theta$，那么从任意一个点，以极角$(\theta ,\theta +\pi )$射出去后一定会被这条直线挡住。因此所有直线的这样的区间并起来就是不能无限延伸的方向的集合。然而区间的并并不好求。注意到区间的并的补集为这些区间的补集的交，所以$(\theta -\pi ,\theta )$的交即为可以无限延伸的方向的集合。 然而求极角需要用到atan2，可能会有精度问题，所以考虑用向量叉积代替极角来判断向量之间的关系。 一个区间用两个向量s, t来表示，原区间即为s转到t的部分。然后对于当前的有向直线，设其方向向量为v。如果v在s的左侧，则s = v，然后判断一下v是否在t左侧，如果是，说明交为空，直接返回false。如果v在s的右侧，则讨论-v，如果-v在t的右侧，则t = -v，否则直线的区间包含了当前区间，继续处理下一条直线。 如果最后交出来的区间不为空，那么说明半平面交面积是无穷的，否则是有穷的。

直线共线的情况比较复杂，需要特判。

网上没有专门的题来测试，就把这个无穷判断放到半平面交的函数里，随便找了个题测。

bzoj 2732 代码（含C++11特性，只能到darkbzoj上交）

#include <cstdio> #include <cmath> #include <cstring> #include <algorithm> #include <stack> using namespace std; #define DEBUG 0 #define MAXN 200011 template <typename T> T sq(T x) { return x * x; } const double eps = 1e-12, pi = acos(-1.0); int sgn(double x) { return x < -eps ? -1 : (x > eps ? 1 : 0); } struct VEC { double x, y; bool operator == (const VEC& b) const { return sgn(x - b.x) == 0 && sgn(y - b.y) == 0; } bool operator < (const VEC& b) const { return sgn(y - b.y) == 0 ? sgn(x - b.x) < 0 : y < b.y; } VEC operator - (const VEC& b) const { return VEC{x - b.x, y - b.y}; } VEC operator + (const VEC& b) const { return VEC{x + b.x, y + b.y}; } double len() const { return hypot(x, y); } double dist(const VEC& b) const { return (*this - b).len(); } double len2() const { return sq(x) + sq(y); } double dist2(const VEC& b) const { return (*this - b).len2(); } VEC operator * (double k) const { return VEC{x * k, y * k}; } VEC trunc(double l) const { double ori = len(); if (!sgn(ori)) return *this; return *this * (l / ori); } VEC operator / (double k) const { return VEC{x / k, y / k}; } VEC norm() const { return *this / len(); } double operator ^ (const VEC& b) const { return x * b.y - y * b.x; } double operator * (const VEC& b) const { return x * b.x + y * b.y; } // The angle between *this and b //anticlockwise is plus double rad_di(const VEC& b) const { return atan2(*this ^ b, *this * b); } double rad(const VEC& b) const { return fabs(rad_di(b)); //return asin(fabs(*this ^ b) / (len() * b.len())); } VEC rotleft() const { return VEC{-y, x}; } VEC rotright() const { return VEC{y, -x}; } VEC rot(double a) const { double c = cos(a), s = sin(a); return VEC{x * c - y * s, x * s + y * c}; } //Independent //0 <= slope < pi double slope() const { double ans = atan2(y, x); if (sgn(ans) < 0) ans += pi; else if (sgn(ans - pi) == 0) ans -= pi; return ans; } double cross(const VEC& p1, const VEC& p2) const { return (p1 - *this) ^ (p2 - *this); } //-1: left //0: parallel //1: right int relation(const VEC& b) const { return sgn(b ^ *this); } VEC operator - () const { return VEC{-x, -y}; } bool SameDi(const VEC& b) const { return sgn(x) == sgn(b.x) && sgn(y) == sgn(b.y); } void input() { scanf("%lf%lf", &x, &y); } }; struct LINE { VEC s, v; void set(const VEC& _s, double a) { s = _s; v = {cos(a), sin(a)}; } //ax + by + c = 0 void set(double a, double b, double c) { v = {b, -a}; s = fabs(a) > fabs(b) ? VEC{-c / a, 0} : VEC{0, -c / b}; } //-1: left //0: on it //1: right int relation(const VEC& p) const { return sgn((p - s) ^ v); } //0: parallel //1: coincide //2: intersect int relation(const LINE& b) const { if (v.relation(b.v) == 0) return b.relation(s) == 0; return 2; } //Assume that they intersect VEC CrossPoint(const LINE& b) const { double t = ((b.s - s) ^ b.v) / (v ^ b.v); return s + v * t; } //left is plus double dist_di(const VEC& p) const { return (v ^ (p - s)) / v.len(); } double dist(const VEC& p) const { return fabs(dist_di(p)); } VEC projection(const VEC& p) const { VEC v1 = v.norm(); return v1 * ((p - s) * v1) + s; } VEC symmetry(const VEC& p) const { VEC p0 = projection(p); return p0 * 2 - p; } void SetAngularBisector(const VEC& a, const VEC& b, const VEC& c) { set(a, (atan2(b.y - a.y, b.x - a.x) + atan2(c.y - a.y, c.x - a.x)) / 2); } LINE operator - () const { return LINE{s, -v}; } }; namespace HalfPlane { //No same line //-2: reverse int di(const LINE& a, const LINE& b) { int rel = a.v.relation(b.v); if (0 == rel) { if (a.v.SameDi(b.v)) { return a.relation(b.s); } else { return -2; } } else { return rel; } } bool Infinite(const LINE ls[], int n) { stack<int> sta; LINE s = ls[0]; LINE t = -ls[0]; for (int i = 1; i < n; ++i) { int rel = di(s, ls[i]); if (-2 == rel) { sta.push(i); } else if (rel < 0) { s = ls[i]; if (di(t, s) < 0) return false; } else if (rel > 0) { rel = di(t, -ls[i]); if (rel > 0) { t = -ls[i]; if (di(s, t) > 0) return false; } } } while (!sta.empty()) { int i = sta.top(); sta.pop(); if (0 == ls[i].v.relation(s.v)) { if (s.relation(ls[i].s) <= 0) { s.s = ls[i].s; } else { return false; } } else if (ls[i].v.relation(t.v) || t.relation(ls[i].s) < 0) { return false; } } if (s.v.relation(t.v)) { return true; } else { if (s.v.SameDi(-t.v)) return true; else return s.relation(t.s) <= 0; } } int que[MAXN]; VEC cp[MAXN]; //Cross points int st, ed; //0: The area is 0 //1: Can construct a convex polygon //2: The area is infinite int Intersection(const LINE ls[], int n) { static pair<double, int> ord[MAXN]; for (int i = 0; i < n; ++i) ord[i] = make_pair(atan2(ls[i].v.y, ls[i].v.x), i); sort(ord, ord + n); int m = 0; for (int i = 1; i < n; ++i) if (sgn(ord[i].first - ord[m].first)) ord[++m] = ord[i]; else if (ls[ord[m].second].relation(ls[ord[i].second].s) == -1) ord[m] = ord[i]; n = m + 1; if (Infinite(ls, n)) return 2; que[st=0] = ord[0].second; que[ed=1] = ord[1].second; cp[1] = ls[ord[1].second].CrossPoint(ls[ord[0].second]); for (int i = 2; i < n; ++i) { const LINE &now = ls[ord[i].second]; while (st < ed && now.relation(cp[ed]) == 1) --ed; while (st < ed && now.relation(cp[st + 1]) == 1) ++st; if (now.v.relation(ls[que[ed]].v) == 0) return 0; que[++ed] = ord[i].second; cp[ed] = now.CrossPoint(ls[que[ed - 1]]); } while (st < ed && ls[que[st]].relation(cp[ed]) == 1) --ed; while (st < ed && ls[que[ed]].relation(cp[st + 1]) == 1) ++st; return st + 1 < ed; } void GetConvex(VEC ps[], int& n, const LINE ls[]) { n = ed - st + 1; memcpy(ps, cp + st, n * sizeof(VEC)); } } int main() { int n; static double x[MAXN], y1[MAXN], y2[MAXN]; static LINE ls[MAXN << 1]; scanf("%d", &n); for (int i = 0; i < n; ++i) { scanf("%lf%lf%lf", x + i, y1 + i, y2 + i); } int tot = 0; ls[tot++] = LINE{VEC{0, 0}, VEC{1, 0}}; ls[tot++] = LINE{VEC{0, 0}, VEC{0, 1}}; //ls[tot++] = LINE{VEC{-1e10, 0}, VEC{0, -1}}; //ls[tot++] = LINE{VEC{0, 1e10}, VEC{-1, 0}}; for (int i = 0; i < n; ++i) { ls[tot++] = LINE{VEC{0, y1[i] / x[i]}, VEC{1, -x[i]}}; ls[tot++] = LINE{VEC{0, y2[i] / x[i]}, VEC{-1, x[i]}}; } #if DEBUG for (int i = 0; i < tot; ++i) { printf("%f %f %f %f\n", ls[i].s.x, ls[i].s.y, ls[i].v.x, ls[i].v.y); } #endif int L = 0, R = n-1, ans = -1; while (L <= R) { int mid = (L + R) >> 1; if (HalfPlane::Intersection(ls, (mid + 2) << 1)) { ans = mid; L = mid + 1; } else { R = mid - 1; } } printf("%d\n", ans + 1); return 0; }