计算几何 模板

#include <algorithm>
#include <cmath>
#include <iostream>
#include <vector>

typedef double db;
const db eps = 1e-9;
const db pi = acos(-1);

// 判断符号
inline int sign(db ps) { return ps < -eps ? -1 : ps > eps; }
// 比较大小
inline int cmp(db ps, db b) { return sign(ps - b); }

// 点
struct P {
    db x, y;
    P(){};
    P(db _x, db _y) : x(_x), y(_y){};
    // 四则运算
    P operator+(P p) const { return {x + p.x, y + p.y}; };
    P operator-(P p) const { return {x - p.x, y - p.y}; };
    P operator*(db d) const { return {x * d, y * d}; };
    P operator/(db d) const { return {x / d, y / d}; };

    bool operator<(P p) const {
        int c = cmp(x, p.x);
        if (c) return c == -1;
        return cmp(y, p.y) == -1;
    }

    bool operator==(P p) const {
        return cmp(x, p.x) == 0 && cmp(y, p.y) == 0;
    }

    // 两点间距离
    db dis(P p) const { return (*this - p).abs(); }
    // 求与x轴夹角
    db alpha() const { return atan2(y, x); }
    // 模长
    db abs() const { return sqrt(abs2()); }
    db abs2() const { return x * x + y * y; }
    // 逆时针旋转
    P rotate(db phl) const { return {x * cos(phl) - y * sin(phl), x * sin(phl) + y * cos(phl)}; }
    P rotate90() const { return P(-y, x); }
    // 单位化
    P unit() const { return *this / abs(); }
};

// 点积
db dot(P p1, P p2) { return p1.x * p2.x + p1.y * p2.y; }
// 叉积
db cross(P p1, P p2) { return p1.x * p2.y - p1.y * p2.x; }
// 夹角
db rad(P p1, P p2) { return atan2(dot(p1, p2), cross(p1, p2)); }

// 判断直线k1k2和直线k3k4是否平行 true平行
bool checkLL(P k1, P k2, P k3, P k4) {
    return cmp(cross(k3 - k1, k4 - k1), cross(k3 - k2, k4 - k2)) == 0;
}
// 求直线k1k2和直线k3k4的交点
P getLL(P k1, P k2, P k3, P k4) {
    db w1 = cross(k1 - k3, k4 - k3), w2 = cross(k4 - k3, k2 - k3);
    return (k1 * w2 + k2 * w1) / (w1 + w2);
}

// 判断区间[l1, r1]和区间[l2, r2]是否相交
bool intersect(db l1, db r1, db l2, db r2) {
    if (l1 > r1) std::swap(l1, r1);
    if (l2 > r2) std::swap(l2, r2);
    return cmp(r1, l2) != -1 && cmp(r2, l1) != -1;
}
// 线段相交
bool checkSS(P k1, P k2, P k3, P k4) {
    return intersect(k1.x, k2.x, k3.x, k4.x) && intersect(k1.y, k2.y, k3.y, k4.y) &&
           sign(cross(k3 - k1, k4 - k1)) * sign(cross(k3 - k2, k4 - k2)) <= 0 &&
           sign(cross(k1 - k3, k2 - k3)) * sign(cross(k1 - k4, k2 - k4)) <= 0;
}
// 线段严格相交
bool checkSS_strict(P k1, P k2, P k3, P k4) {
    return sign(cross(k3 - k1, k4 - k1)) * sign(cross(k3 - k2, k4 - k2)) < 0 &&
           sign(cross(k1 - k3, k2 - k3)) * sign(cross(k1 - k4, k2 - k4)) < 0;
}

// m 在 [k1,k2] 内
bool isMiddle(db k1, db m, db k2) {
    return sign(k1 - m) * sign(k2 - m) <= 0;
}
// 点m在点k1和点k2中间
bool isMiddle(P k1, P m, P k2) {
    return isMiddle(k1.x, m.x, k2.x) && isMiddle(k1.y, m.y, k2.y);
}
// 点p在线段k1k2上
bool onSeg(P k1, P k2, P p) {
    return sign(cross(p - k1, k2 - k1)) == 0 && isMiddle(k1, p, k2);
}
bool onSeg_strict(P k1, P k2, P p) {
    return sign(cross(p - k1, k2 - k1)) == 0 && sign(cross(p - k1, k1 - k2)) * sign(cross(p - k2, k1 - k2)) < 0;
}

// 求点p在直线k1k2上的投影
P project(P k1, P k2, P p) {
    P k = k2 - k1;
    return k1 + k * (dot(p - k1, k) / k.abs2());
}
// 求点p关于直线k1k2的对称点
P reflect(P k1, P k2, P p) {
    return project(k1, k2, p) * 2 - p;
}
// p到线段k1k2的最近距离
db nearest(P k1, P k2, P p) {
    P h = project(k1, k2, p);
    if (isMiddle(k1, h, k2)) return p.dis(h);
    return std::min(k1.dis(p), k2.dis(p));
}
// 线段距离
db disSS(P k1, P k2, P k3, P k4) {
    if (checkSS(k1, k2, k3, k4)) return 0;
    return std::min(std::min(nearest(k1, k2, k3), nearest(k1, k2, k4)), std::min(nearest(k3, k4, k1), nearest(k3, k4, k2)));
}

// 多边形面积
db area(std::vector<P> ps) {
    int n = ps.size();
    db res = 0;
    for (int i = 0; i < n; ++i) res += cross(ps[i], ps[(i + 1) % n]);
    return res / 2;
}

// 点包含 2内部 1边界 0外部
int contain(std::vector<P> ps, P p) {
    int n = ps.size(), res = 0;
    for (int i = 0; i < n; ++i) {
        P u = ps[i], v = ps[(i + 1) % n];
        if (onSeg(u, v, p)) return 1;
        if (cmp(u.y, v.y) <= 0) std::swap(u, v);
        if (cmp(p.y, u.y) > 0 || cmp(p.y, v.y) <= 0) continue;
        res ^= (sign(cross(u - p, v - p)) > 0);
    }
    return res * 2;
}

// 凸包 flag=0 不严格 flag=1 严格
std::vector<P> ConvexHull(std::vector<P> ps, int flag = 1) {
    int n = ps.size(), now = -1;
    std::vector<P> res(n * 2);
    sort(ps.begin(), ps.end());
    for (int i = 0; i < ps.size(); ++i) {
        while (now > 0 && sign(cross(res[now] - res[now - 1], ps[i] - res[now - 1])) < flag) now--;
        res[++now] = ps[i];
    }
    int pre = now;
    for (int i = n - 2; i >= 0; --i) {
        while (now > pre && sign(cross(res[now] - res[now - 1], ps[i] - res[now - 1])) < flag) now--;
        res[++now] = ps[i];
    }
    res.resize(now);
    return res;
}

// 点集直径
db convexDiameter(std::vector<P> ps) {
    int now = 0, n = ps.size();
    db res = 0;
    for (int i = 0; i < ps.size(); ++i) {
        now = std::max(now, i);
        while (true) {
            db k1 = ps[i].dis(ps[now % n]), k2 = ps[i].dis(ps[(now + 1) % n]);
            res = std::max(res, std::max(k1, k2));
            if (k2 > k1)
                now++;
            else
                break;
        }
    }
    return res;
}

// 半径为r的圆o与直线p1p2的交交点
std::vector<P> isCL(P o, db r, P p1, P p2) {
    db dis = cross(o - p1, p2 - p1) / p1.dis(p2);
    if (cmp(abs(dis), r) > 0) return {};
    db x = dot(p1 - o, p2 - p1), y = (p2 - p1).abs2(), d = x * x - y * ((p1 - o).abs2() - r * r);
    d = std::max(d, 0.0);
    P m = p1 - (p2 - p1) * (x / y), dr = (p2 - p1) * sqrt(d) / y;
    return {m - dr, m + dr};
}

// 半径为r1的圆o1 与 半径为r2的圆o2 的交点
std::vector<P> isCC(P o1, db r1, P o2, db r2) {
    db d = o1.dis(o2);
    if (cmp(d, r1 + r2) == 1) return {};
    if (cmp(d, abs(r1 - r2)) == -1) return {};
    d = std::min(d, r1 + r2);
    db y = (r1 * r1 + d * d - r2 * r2) / (2 * d), x = sqrt(r1 * r1 - y * y);
    P dr = (o2 - o1).unit();
    P q1 = o1 + dr * y, q2 = dr.rotate90() * x;
    return {q1 - q2, q1 + q2};
}

// 切线
std::vector<P> tanCP(P o, db r, P p) {
    db x = (p - o).abs2(), d = x - r * r;
    if (sign(d) <= 0) return {};
    P q1 = o + (p - o) * (r * r / x);
    P q2 = (p - o).rotate90() * (r * sqrt(d) / x);
    return {q1 - q2, q1 + q2};
}

(未完待续)