计算几何

计算几何模板

计算几何

前言

模板有两种使用方式，一种是直接用 double，一种是用 Frac<i128>

Frac<i128>

template<class T>
struct Frac {
    T num;
    T den;
    Frac(T num_, T den_) : num(num_), den(den_) { if (den < 0) { den = -den, num = -num; } }
    Frac() : Frac(0, 1) {}
    Frac(T num_) : Frac(num_, 1) {}
    explicit operator double() const { return 1. * num / den; }
    Frac &operator+=(const Frac &rhs) {
        num = num * rhs.den + rhs.num * den;
        den *= rhs.den;
        return *this;
    }
    Frac &operator-=(const Frac &rhs) {
        num = num * rhs.den - rhs.num * den;
        den *= rhs.den;
        return *this;
    }
    Frac &operator*=(const Frac &rhs) {
        num *= rhs.num;
        den *= rhs.den;
        return *this;
    }
    Frac &operator/=(const Frac &rhs) {
        num *= rhs.den;
        den *= rhs.num;
        if (den < 0) { num = -num, den = -den; }
        return *this;
    }
    friend Frac operator+(Frac lhs, const Frac &rhs) { return lhs += rhs; }
    friend Frac operator-(Frac lhs, const Frac &rhs) { return lhs -= rhs; }
    friend Frac operator*(Frac lhs, const Frac &rhs) { return lhs *= rhs; }
    friend Frac operator/(Frac lhs, const Frac &rhs) { return lhs /= rhs; }
    friend Frac operator-(const Frac &a) { return Frac(-a.num, a.den); }
    friend bool operator==(const Frac &lhs, const Frac &rhs) { return lhs.num * rhs.den == rhs.num * lhs.den; }
    friend bool operator!=(const Frac &lhs, const Frac &rhs) { return lhs.num * rhs.den != rhs.num * lhs.den; }
    friend bool operator<(const Frac &lhs, const Frac &rhs) { return lhs.num * rhs.den < rhs.num * lhs.den; }
    friend bool operator>(const Frac &lhs, const Frac &rhs) { return lhs.num * rhs.den > rhs.num * lhs.den; }
    friend bool operator<=(const Frac &lhs, const Frac &rhs) { return lhs.num * rhs.den <= rhs.num * lhs.den; }
    friend bool operator>=(const Frac &lhs, const Frac &rhs) { return lhs.num * rhs.den >= rhs.num * lhs.den; }
    friend ostream &operator<<(ostream &os, Frac x) {
        T g = gcd(x.num, x.den);
        if (x.den == g) { return os << x.num / g; }
        else { return os << x.num / g << "/" << x.den / g; }
    }
};

精度控制和常用常量

const double pi = 3.141592653589793238;
const double eps = 1e-8;


template <class T>
bool eq(T a, T b) { return abs(a - b) < eps; }  // ==
template <class T>
bool gt(T a, T b) { return a - b > eps; }  // >
template <class T>
bool lt(T a, T b) { return a - b < -eps; }  // <
template <class T>
bool ge(T a, T b) { return a - b > -eps; }  // >=
template <class T>
bool le(T a, T b) { return a - b < eps; }  // <=


template <class T>
const Point<T> O = { 0, 0 };

二维计算几何

定义
点（向量）
template <class T>
struct Point {
    T x, y;
    Point(T x_ = 0, T y_ = 0) : x(x_), y(y_) {}
    template <class U>
    operator Point<U>() { return Point<U>(U(x), U(y)); }
    Point& operator+=(Point p) & { x += p.x, y += p.y;return *this; }
    Point& operator-=(Point p) & { x -= p.x, y -= p.y;return *this; }
    Point& operator*=(T v) & { x *= v, y *= v;return *this; }
    Point operator-() const { return Point(-x, -y); }
    friend Point operator+(Point a, Point b) { return a += b; }
    friend Point operator-(Point a, Point b) { return a -= b; }
    friend Point operator*(Point a, T b) { return a *= b; }
    friend Point operator*(T a, Point b) { return b *= a; }
    friend bool operator==(Point a, Point b) { return a.x == b.x && a.y == b.y; }
    friend bool operator!=(Point a, Point b) { return !(a == b); }
    friend bool operator<(Point a, Point b) {
        if(eq(a.x, b.x)) return lt(a.y, b.y);
        return lt(a.x, b.x);
    }
    friend istream &operator>>(istream &is, Point& p) { return is >> p.x >> p.y; }
    friend ostream &operator<<(ostream &os, Point p) { return os << "(" << p.x << ", " << p.y << ")"; }
};
template <class T>
using Vec = Point<T>;
直线（两点式）
template <class T>
struct Line {
    Point<T> a, b;
    Line(Point<T> a_ = Point<T>(), Point<T> b_ = Point<T>()) : a(a_), b(b_) {}
};
template <class T>
using Seg = Line<T>;
圆
template <class T>
struct Circle {
    Point<T> p;
    T r;
    Circle() {}
    Circle(Point<T> p, T r) : p(p), r(r) {}
};
向量操作
// 点乘
template <class T>
T dot(Vec<T> a, Vec<T> b) { return a.x * b.x + a.y * b.y; }


// 模长平方
template <class T>
T square(Vec<T> p) { return dot(p, p); }


// 向量长度
// square
template <class T>
double len(Vec<T> p) { return sqrt(double(square(p))); }


// 叉乘
template <class T>
T cross(Vec<T> a, Vec<T> b) { return a.x * b.y - a.y * b.x; }


// 向量顺时针旋转90度
template <class T>
Vec<T> r90c(Vec<T> v) { return { v.y, -v.x }; }


// 向量逆时针旋转90度
template <class T>
Vec<T> r90a(Vec<T> v) { return { -v.y, v.x }; }


// 点绕某点旋转 theta
template <class T>
Point<T> rot(Point<T> p, Point<T> o, double theta) {
    auto v = p - o;
    auto c = cos(theta), s = sin(theta);
    return {o.x + v.x * c - v.y * s, o.y + v.x * s + v.y * c};
}


// 向量夹角
// dot, len
template <class T>
double cos_t(Vec<T> a, Vec<T> b) { return dot(a, b) / len(a) / len(b); }


// 方向向量
// len
template <class T>
Vec<T> norm(Vec<T> v) { return { v.x / len(v), v.y / len(v) }; }


// 斜率
template <class T>
double slope(Vec<T> v) { return v.y / v.x; }


// 求夹角
double rotangle(Vec a, Vec b) {
    return atan2(cross(a, b), dot(a, b));
}
// 返回 a 逆时针转动多少度会到达 b，返回负值表示顺时针转动，可以通过如果为负就 += 2Pi 全部变为逆时针转动
距离
// 点到点距离
template <class T>
double dis(Point<T> a, Point<T> b) { return hypot(a.x - b.x, a.y - b.y); }


// 点到点距离平方
// 这里用 T 因为整数点为整数
// square
template <class T>
T dis2(Point<T> a, Point<T> b) { return square(a - b); }


// 点到直线距离
template <class T>
double dis(Point<T> p, Line<T> l) { return abs(cross(p - l.a, p - l.b)) / dis(l.a, l.b); }


// 点到直线距离的平方
template <class T>
double dis2(Point<T> p, Line<T> l) { return cross(p - l.a, p - l.b) * cross(p - l.a, p - l.b) / dis2(l.a, l.b); }


// 点到线段最近点
// cross,dis(p,p),inter(l,l)
template <class T>
Point<T> closet_seg(Point<T> p, Seg<T> l) {
    auto t = p;
    t.x += l.a.y - l.b.y;
    t.y += l.b.x - l.a.x;
    if (gt(cross(l.a - p, t - p) * cross(l.b - p, t - p), (T)0)) {
        return lt(dis(p, l.a), dis(p, l.b)) ? l.a : l.b;
    }
    return inter(Line<T>(p, t), l);
}


// 线段间最近距离
// closet_seg
template <class T>
T dis_seg(Seg<T> s1, Seg<T> s2) {
    return min({closet_seg(s1.a, s2), closet_seg(s1.b, s2),
                closet_seg(s2.a, s1), closet_seg(s2.b, s1)});
}
对称
// 点a关于点p对称
template <class T>
Point<T> reflect(Point<T> a, Point<T> p) { return { p.x * 2 - a.x, p.y * 2 - a.y }; }


// l关于p对称
// reflect(p,p)
template <class T>
Line<T> reflect(Line<T> l, Point<T> p) { return { reflect(l.a, p), reflect(l.b, p) }; }


// 关于直线对称
// NOTE 向量和点在这里的表现不同，求向量关于某直线的对称向量需要用reflect_v
// reflect(p,p),pedal
template <class T>
Point<T> reflect(Point<T> a, Line<T> ax) { return reflect(a, pedal(a, ax)); }


// 向量关于直线对称
// reflect(p,l)
template <class T>
Vec<T> reflect_v(Vec<T> v, Line<T> ax) { return reflect(v, ax) - reflect(O, ax); }


// 直线关于直线对称
// reflect(p,l)
template <class T>
Line<T> reflect(Line<T> l, Line<T> ax) { return { reflect(l.a, ax), reflect(l.b, ax) }; }
关系判定
// 判定点是否在线段上
// cross
template <class T>
bool on_seg(Point<T> p, Line<T> l) {
    bool ok[5];
    ok[0] = eq(cross(p - l.a, l.b - l.a), (T)0);
    ok[1] = le(min(l.a.x, l.b.x), p.x);
    ok[2] = le(p.x, max(l.a.x, l.b.x));
    ok[3] = le(min(l.a.y, l.b.y), p.y);
    ok[4] = le(p.y, max(l.a.y, l.b.y));
    return ok[0] & ok[1] & ok[2] & ok[3] & ok[4];
}


// 判定 q 是否在 p 左侧
// -1: 右边
//  0: 共线
//  1: 左边
// cross
template <class T>
int isleft(Vec<T> p, Vec<T> q) {
    auto t = cross(p, q);
    if(eq(t, (T)0)) return 0;
    if(gt(t, (T)0)) return 1;
    else return -1;
}


// 判定是否在直线左侧
// -1: 右侧
//  0: 线上
//  1: 左侧
// cross
template <class T>
int isleft(Point<T> p, Line<T> l) {
    auto t = cross(l.b - l.a, p - l.a);
    if(eq(t, (T)0)) return 0;
    if(gt(t, (T)0)) return 1;
    else return -1;
}


// 判定点是否在多边形内，使用绕行数计算
// 特判在边上    {true, 1}
// 内部     {true, 0}
// 外部         {false, res}
// on_seg,isleft(p,l)
template <class T>
pair<bool, int> in(Point<T> p, vector<Point<T>>& P) {
    int n = P.size();
    int res = 0;
    for(int i = 0; i < n; i++) {
        Line<T> l(P[i] ,P[(i + 1) % n]);
        if(on_seg(p, l)) return {true, 1};
        if(eq(l.a.y, l.b.y)) continue;
        if(le(l.a.y, p.y)) {
            if(ge(l.b.y, p.y) && isleft(p, l) == 1) ++res;
        } else {
            if(le(l.b.y, p.y) && isleft(p, l) == 0) --res;
        }
    }
    if(!res) return {true, 0};
    return {false, res};
}


// -1: 重合
//  0: 不重合
//  1: 端点重合
//  2: 严格相交
// isleft(p,l)
template <class T>
int inter_lineseg(Line<T> l, Seg<T> s) {
    if(!isleft(s.a, l) || !isleft(s.b, l)) {
        if(isleft(s.a, l) == isleft(s.b, l)) return -1;
        else return 1;
    }
    if(isleft(s.a, l) != isleft(s.b, l)) return 2;
    else return 0;
}


// 判定向量垂直
// dot
template <class T>
bool verti(Vec<T> a, Vec<T> b) { return eq(dot(a, b), (T)0); }


// 判定向量平行
// cross
template <class T>
bool paral(Vec<T> a, Vec<T> b) { return eq(cross(a, b), (T)0); }


// 判定直线关系
// -1: 重合
//  0: 平行
//  1: 垂直
//  2: 相交
// paral,verti
template <class T>
int line_relation(Line<T> l1, Line<T> l2) {
    if(paral(l1.b - l1.a, l2.b - l2.a)) return paral(l1.b - l2.b, l1.b - l1.a) ? 0 : -1;
    if(verti(l1.b - l1.a, l2.b - l2.a)) return 1;
    return 2;
}
直线操作
// 点到直线的垂足
// dot,square
template <class T>
Point<T> pedal(Point<T> p, Line<T> l) {
    Vec<T> ab = l.b - l.a;  // 向量 AB
    Vec<T> ap = p - l.a;    // 向量 AP
    // 投影长度比例
    double t = dot(ap, ab) / double(square(ab));
    // 垂足点 F
    Point<T> f = l.a + ab * t;
    return f;
}


// 计算过点的直线垂线
// norm
template <class T>
Line<T> perpline(Line<T> l, Point<T> p) {
    Point<T> dir = l.b - l.a;
    Point<T> norm(-dir.y, dir.x);
    return Line<T>(p, p + norm);
}
相交
// 直线相交
// cross
template <class T>
Point<T> inter(Line<T> l1, Line<T> l2) {
    return l1.a + (l1.b - l1.a) * (cross(l2.b - l2.a, l1.a - l2.a) / cross(l2.b - l2.a, l1.a - l1.b));
}


// 直线与圆交点
// pedal,len,norm
template <class T>
vector<Point<T>> inter(Line<T> l, Circle<T> c) {
    Point<T> P = pedal(c.p, l);
    double h = len(P - c.p);
    if (gt(h, c.r)) return {};
    if (eq(h, c.r)) return { P, P };
    double d = sqrt(c.r * c.r - h * h);
    Vec<T> vec = d * norm(l.b - l.a);
    return { P - vec, P + vec };
}


// 圆与圆的交点
// 一个用法是判定两圆关系
// 返回长度为 0: 外离或内含
// 返回长度为 1: 外切或内切
// 返回长度为 2: 相交
// r90c,len
template <class T>
vector<Point<T>> inter(Circle<T> c1, Circle<T> c2) {
    auto v1 = c2.p - c1.p, v2 = r90c(v1);
    double d = len(v1);
    if (gt(d, c1.r + c2.r) || gt(abs(c1.r - c2.r), d)) return {};
    if (eq(d, c1.r + c2.r) || eq(d, abs(c1.r - c2.r))) return { c1.p + c1.r / d * v1 };
    double a = ((c1.r * c1.r - c2.r * c2.r) / d + d) / 2;
    double h = sqrt(c1.r * c1.r - a * a);
    auto av = a / len(v1) * v1, hv = h / len(v2) * v2;
    return { c1.p + av + hv, c1.p + av - hv };
}


// 线段相交
// 0 : 不相交
// 1 : 严格相交
// 2 : 重合
// 3 : 端点处相交
// on_seg,cross,inter(l,l)
template <class T>
tuple<int, Point<T>, Point<T>> inter_seg(Line<T> l1, Line<T> l2) {
    if (lt(max(l1.a.x, l1.b.x), min(l2.a.x, l2.b.x))) {
        return { 0, Point<T>(), Point<T>() };
    }
    if (gt(min(l1.a.x, l1.b.x), max(l2.a.x, l2.b.x))) {
        return { 0, Point<T>(), Point<T>() };
    }
    if (lt(max(l1.a.y, l1.b.y), min(l2.a.y, l2.b.y))) {
        return { 0, Point<T>(), Point<T>() };
    }
    if (gt(min(l1.a.y, l1.b.y), max(l2.a.y, l2.b.y))) {
        return { 0, Point<T>(), Point<T>() };
    }
    if (eq(cross(l1.b - l1.a, l2.b - l2.a), (T)0)) {
        if (eq(cross(l1.b - l1.a, l2.a - l1.a), (T)0)) {
            auto maxx1 = max(l1.a.x, l1.b.x);
            auto minx1 = min(l1.a.x, l1.b.x);
            auto maxy1 = max(l1.a.y, l1.b.y);
            auto miny1 = min(l1.a.y, l1.b.y);
            auto maxx2 = max(l2.a.x, l2.b.x);
            auto minx2 = min(l2.a.x, l2.b.x);
            auto maxy2 = max(l2.a.y, l2.b.y);
            auto miny2 = min(l2.a.y, l2.b.y);
            Point<T> p1(max(minx1, minx2), max(miny1, miny2));
            Point<T> p2(min(maxx1, maxx2), min(maxy1, maxy2));
            if (!on_seg(p1, l1)) {
                swap(p1.y, p2.y);
            }
            if (p1 == p2) {
                return { 3, p1, p2 };
            } else {
                return { 2, p1, p2 };
            }
        } else {
            return { 0, Point<T>(), Point<T>() };
        }
    }
    auto cp1 = cross(l2.a - l1.a, l2.b - l1.a);
    auto cp2 = cross(l2.a - l1.b, l2.b - l1.b);
    auto cp3 = cross(l1.a - l2.a, l1.b - l2.a);
    auto cp4 = cross(l1.a - l2.b, l1.b - l2.b);


    if ((gt(cp1, (T)0) && gt(cp2, (T)0)) || (lt(cp1, (T)0) && lt(cp2, (T)0)) ||
        (gt(cp3, (T)0) && gt(cp4, (T)0)) || (lt(cp3, (T)0) && lt(cp4, (T)0))) {
        return { 0, Point<T>(), Point<T>() };
    }


    Point<T> p = inter(l1, l2);
    if (eq(cp1, (T)0) || eq(cp2, (T)0) || eq(cp3, (T)0) || eq(cp4, (T)0)) {
        return { 3, p, p };
    } else {
        return { 1, p, p };
    }
}
凸包
// 判断凸性
// cross
template <class T>
bool is_convex(vector<Point<T>>& P) {
    int n = P.size();
    bool ok1, ok2;
    ok1 = ok2 = false;
    for (int i = 0; i < n; ++i) {
        auto v1 = P[(i + 1) % n] - P[i];
        auto v2 = P[(i + 2) % n] - P[(i + 1) % n];
        if (gt(cross(v1, v2), (T)0)) ok1 = true;
        else if (lt(cross(v1, v2), (T)0)) ok2 = true;
        if (ok1 && ok2) return false;
    }
    return true;
}


// 求凸包
// cross
template <class T>
vector<Point<T>> hull(vector<Point<T>> p) {
    vector<Point<T>> h, l;
    sort(p.begin(), p.end(),[&](auto a, auto b) {
        if (a.x != b.x) {
            return a.x < b.x;
        } else {
            return a.y < b.y;
        }
    });
    p.erase(unique(p.begin(), p.end()), p.end());
    if (p.size() <= 1) return p;
    for (auto a : p) {
        while (h.size() > 1 && le(cross(a - h.back(), a - h[h.size() - 2]), (T)0)) {
            h.pop_back();
        }
        while (l.size() > 1 && ge(cross(a - l.back(), a - l[l.size() - 2]), (T)0)) {
            l.pop_back();
        }
        l.push_back(a);
        h.push_back(a);
    }
    l.pop_back();
    reverse(h.begin(), h.end());
    h.pop_back();
    l.insert(l.end(), h.begin(), h.end());
    return l;
}


// 凸包直径
// dis,cross
template <class T>
T hullDiameter(vector<Point<T>>& P) {
    int n = P.size();
    if (n < 2) return 0;
    if (n == 2) return dis(P[0], P[1]);
    T mx = 0;
    int j = 0;
    for (int i = 0; i < n; i++) {
        auto u = P[i], v = P[(i + 1) % n];
        while(le(cross(u - P[j], v - P[j]), cross(u - P[(j + 1) % n], v - P[(j + 1) % n]))) j = (j + 1) % n;
        mx = max({mx, dis(P[(i + 1) % n], P[j]), dis(P[i], P[j])});
    }
    return mx;
}


// 旋转卡壳求两凸包间最大三角形
// cross,area_t
template <class T>
i64 RotateCalipers(vector<Point<T>> P1, vector<Point<T>> P2) {
    int n = P1.size(), m = P2.size();
    int j = 0;
    i64 res = LONG_LONG_MAX;
    for(int i = 0; i < n; i++) {
        auto u = P1[i], v = P1[(i + 1) % n];
        while(abs(cross(u - P2[(j + 1) % m], v - P2[(j + 1) % m])) <
              abs(cross(u - P2[j], v - P2[j]))) j = (j + 1) % m;
        res = min(res, area_t(P2[j], {u, v}));
    }
    return res;
}
圆相关
// 圆上取点
template <class T>
Point<T> point(Circle<T> c, double theta) {
    return Point<T>(c.p.x + cos(theta) * c.r, c.p.y + sin(theta) * c.r);
}


// 点作圆切线交点
// dis,r90c
template <class T>
vector<Point<T>> tangent(Point<T> p, Circle<T> c) {
    auto t = dis(p, c.p);
    if (eq(dis(p, c.p), c.r)) {
        return { p, p };
    }
    auto a = c.r * c.r / t, b = sqrt(c.r * c.r - a * a);
    auto e1 = (1 / t) * (p - c.p);
    auto e2 = r90c(e1);
    auto p1 = c.p + a * e1 + b * e2, p2 = c.p + a * e1 - b * e2;
    if (gt(p1.x, p2.x))
        swap(p1, p2);
    else if (eq(p1.x, p2.x) && gt(p1.y, p2.y))
        swap(p1, p2);
    return { p1, p2 };
}


// 两圆公切线
// inter(c,c),dis,point
template <class T>
vector<Line<T>> tangent(Circle<T> a, Circle<T> b) {
    if (a.r < b.r) {
        swap(a, b);
    }
    vector<Point<T>> P = inter(a, b);
    auto d = dis(a.p, b.p);
    // 包含
    if (gt(a.r - b.r, d)) {
        return {};
    }
    vector<Line<T>> cur;
    if (P.size() == 1) {
        auto p = P[0];
        cur.push_back(Line<T>(p, p));
        if (eq(d, a.r - b.r)) return cur;
    }
    auto base = atan2(b.p.y - a.p.y, b.p.x - a.p.x);
    auto ang = acos((a.r - b.r) / d);
    auto p1 = point(a, base + ang), p2 = point(b, base + ang);
    cur.push_back(Line<T>(p1, p2));
    p1 = point(a, base - ang), p2 = point(b, base - ang);
    cur.push_back(Line<T>(p1, p2));
    if (P.size() == 2) return cur;
    ang = acos((a.r + b.r) / d);
    p1 = point(a, base + ang);
    p2 = point(b, base + ang + pi);
    cur.push_back(Line<T>(p1, p2));
    p1 = point(a, base - ang);
    p2 = point(b, base - ang + pi);
    cur.push_back(Line<T>(p1, p2));
    return cur;
}
网格
// 多边形上的网格点个数
template <class T>
int grid_onedge(int n, vector<Point<T>> P) {
    int n = P.size(), ret = 0;
    for (int i = 0; i < n; i++)
        ret += gcd(abs(P[i].x - P[(i + 1) % n].x), abs(P[i].y - P[(i + 1) % n].y));
    return ret;
}


// 多边形内的网格点个数
// grid_onedge
int grid_inside(int n, vector<Point<T>> P) {
    int n = P.size(), ret = 0;
    for (int i = 0; i < n; i++)
        ret += P[(i + 1) % n].y * (P[i].x - P[(i + 2) % n].x);
    return (abs(ret) - grid_onedge(n, P)) / 2 + 1;
}

皮克定理：2S=2a+b−2S=a+b2−1a=S−b2+12S = 2a + b- 2 \\ S = a + \frac{b}{2} - 1 \\a = S - \frac{b}{2} + 1

面积
// 求三角形面积
// 如果是整数一般去掉除2
// cross
template <class T>
T area_t(Point<T> p, Seg<T> s) {
    return abs(cross(s.a - p, s.b - p)) / 2.;
};


// 圆和三角形交（三角剖分）
// dot,cross,dis,inter(l,l),inter(l,c),len
template <class T>
T areaofCT(Point<T> a, Point<T> b, Circle<T> c) {
    auto sign = 1.0;
    a = a - c.p;
    b = b - c.p;
    if (eq(cross(a, b), (T)0)) return 0.0;
    if (dis(a, c.p) > dis(b, c.p)) {
        swap(a, b);
        sign = -1.0;
    }
    if (dis(a, c.p) <= c.r) {
        if (dis(b, c.p) <= c.r) {
            return cross(a, b) / 2.0 * sign;
        }
        auto intersections = inter(Line<T>(a, b), c);
        auto p1 = intersections[0], p2 = intersections[1];
        if (len(p1 - b) > len(p2 - b)) swap(p1, p2);
        auto ret1 = fabs(cross(a, p1));
        auto ret2 = acos(dot(p1, b) / len(p1) / len(b)) * c.r * c.r;
        assert(ret2 <= pi * c.r * c.r);
        auto ret = (ret1 + ret2) / 2.0;
        if ((lt(cross(a, b), (T)0) && sign > 0.0) || (gt(cross(a, b), (T)0) && sign < 0.0)) {
            ret = -ret;
        }
        return ret;
    }
    auto ins = closet_seg(c.p, Line<T>(a, b));
    if (gt(dis(c.p, ins), c.r)) {
        auto ret = acos(dot(a, b) / len(a) / len(b)) * c.r * c.r / 2.0;
        assert(ret <= pi * c.r * c.r);
        if ((lt(cross(a, b), (T)0) && sign > 0.0) || (gt(cross(a, b), (T)0) && sign < 0.0)) {
            ret = -ret;
        }
        return ret;
    }
    auto intersections = inter(Line<T>(a, b), c);
    auto p1 = intersections[0], p2 = intersections[1];
    auto cm = c.r / (len(c.p - a) - c.r);
    auto m = Point<T>((c.p.x + cm * a.x) / (1 + cm), (c.p.y + cm * a.y) / (1 + cm));
    auto cn = c.r / (len(c.p - b) - c.r);
    auto n = Point<T>((c.p.x + cn * b.x) / (1 + cn), (c.p.y + cn * b.y) / (1 + cn));
    auto ret1 = acos(dot(m, n) / len(m) / len(n)) * c.r * c.r;
    auto ret2 = acos(dot(p1, p2) / len(p1) / len(p2)) * c.r * c.r - fabs(cross(p1, p2));
    auto ret = (ret1 - ret2) / 2.0;
    if ((lt(cross(a, b), (T)0) && sign > 0.0) || (gt(cross(a, b), (T)0) && sign < 0.0)) {
        ret = -ret;
    }
    return ret;
}


// 多边形面积计算
// 可用于判定顺逆时针，正顺负逆
// cross
template <class T>
double calArea(const vector<Point<T>>& P) {
    int n = P.size();
    double area = 0;
    for (int i = 0; i < n; i++) {
        int j = (i + 1) % n;
        area += cross(P[i], P[j]);
    }
    return area;
}
三角形
// 外心
template <class T>
Point<T> ocenter(Point<T> a, Point<T> b, Point<T> c) {
    T a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1 * a1 + b1 * b1) / 2;
    T a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2 * a2 + b2 * b2) / 2;
    T d = a1 * b2 - a2 * b1;
    return Point<T>(a.x + (c1 * b2 - c2 * b1) / d, a.y + (a1 * c2 - a2 * c1) / d);
}


// 内心
// dis
template <class T>
Point<T> icenter(Point<T> a, Point<T> b, Point<T> c) {
    T A = dis(b, c), B = dis(a, c), C = dis(a, b);
    return Point<T>((A * a.x + B * b.x + C * c.x) / (A + B + C), (A * a.y + B * b.y + C * c.y) / (A + B + C));
}


// 垂心
// ocenter
template <class T>
Point<T> hcenter(Point<T> a, Point<T> b, Point<T> c) {
    return a + b + c - ocenter(a, b, c) * 2;
}


// 重心
template <class T>
Point<T> gcenter(Point<T> a, Point<T> b, Point<T> c) {
    return 1.0 / 3 * (a + b + c) ;
}
半平面交
// 判定点的象限
// 一二象限或 x 轴正方向
template <class T>
int sgn(Point<T> a) {
    if (gt(a.y, (T)0) || (eq(a.y, (T)0) && gt(a.x, (T)0))) {
        return 1;
    } else {
        return -1;
    }
}


// 判定是否在直线左侧
// cross
template <class T>
bool isleft(Point<T> p, Line<T> l) {
    return gt(cross(l.b - l.a, p - l.a), (T)0);
}


// 半平面交
// sgn,isleft,cross,dot,inter(l,l)
template <class T>
vector<Point<T>> hp(vector<Line<T>> lines) {
    sort(lines.begin(), lines.end(), [&](auto l1, auto l2) {
        auto d1 = l1.b - l1.a;
        auto d2 = l2.b - l2.a;
        if (sgn(d1) != sgn(d2)) {
            return sgn(d1) == 1;
        }
        return gt(cross(d1, d2), (T)0);
    });
    deque<Line<T>> ls;
    deque<Point<T>> ps;
    for (auto l : lines) {
        if (ls.empty()) {
            ls.push_back(l);
            continue;
        }
        while (!ps.empty() && !isleft(ps.back(), l)) {
            ps.pop_back();
            ls.pop_back();
        }
        while (!ps.empty() && !isleft(ps[0], l)) {
            ps.pop_front();
            ls.pop_front();
        }
        if (eq(cross(l.b - l.a, ls.back().b - ls.back().a), (T)0)) {
            if (gt(dot(l.b - l.a, ls.back().b - ls.back().a), (T)0)) {
                if (!isleft(ls.back().a, l)) {
                    assert(ls.size() == 1);
                    ls[0] = l;
                }
                continue;
            }
            return {};
        }
        ps.push_back(inter(ls.back(), l));
        ls.push_back(l);
    }
    while (!ps.empty() && !isleft(ps.back(), ls[0])) {
        ps.pop_back();
        ls.pop_back();
    }
    if (ls.size() <= 2) {
        return {};
    }
    ps.push_back(inter(ls[0], ls.back()));
    return vector(ps.begin(), ps.end());
}
极角排序
// 极角排序
template <class T>
int get_region(Point<T> p) { return lt(p.y, (T)0) ? -1 : gt(p.y, (T)0) | (eq(p.y, (T)0) & lt(p.x, (T)0));}
// -1 下半平面 （不包括x轴） 1 上半平面/x轴负半轴 0 原点/x轴正半轴


// 使用 cmp_arg<T>
// get_region, isleft
template <class T>
bool cmp_arg(Point<T> a, Point<T> b) {
    // 下半平面 < 原点（极角认为是 0） < 正半轴 < 上半平面 < 负半轴
    int p = get_region(a), q = get_region(b);
    return p == q ? isleft(a, b) == 1 : p < q; // 同一区域叉积判断 否则判区间
}
特殊用法
// 方便把点扔进 unordered_set
template<typename T>
struct std::hash<Point<T>> {
    size_t operator()(const Point<T>& p) const noexcept {  // 添加 noexcept
        return std::hash<T>()(p.x) ^ (std::hash<T>()(p.y) << 1);
    }
};