带无穷判断的半平面交

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

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

如果半平面交面积无穷，那么一定存在一个点，使得从这个点朝某个方向射去能够无限延伸。假设存在一条有向直线，其方向向量为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;
}