计算几何全家桶

计算几何模板。

• Pick 定理：给定顶点均为整点的简单多边形，其面积 $S$ 和内部格点数目 $a$ 、边上格点数目 $b$ 的关系：$S = a + \frac b 2 - 1$ 。
• 欧拉公式：多面体的点数 $V$ 与边数 $E$ 与面数 $F$ 的关系：$V - E + F = 2$ 。

查看代码

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const int N = 1e5 + 10;
const double pi = acos(-1), eps = 1e-10;
int cmp(double a, double b)
{
    if (a - b > eps) return 1;
    if (b - a > eps) return -1;
    return 0;
}
int sign(double a) { return cmp(a, 0); }
struct Point // 一个点或者一个向量
{
    double x, y;
    void input() { scanf("%lf%lf", &x, &y); }
    void output() { printf("%.5lf %.5lf\n", x + eps, y + eps); } // 防止输出 -0.00
    bool operator<(const Point &_) const
    {
        int t = cmp(x, _.x);
        if (t == 1) return false;
        if (t == -1) return true;
        return y < _.y;
    }
    Point operator+(const Point &_) const { return (Point){x + _.x, y + _.y}; }
    Point operator-(const Point &_) const { return (Point){x - _.x, y - _.y}; }
    Point operator*(const double &_) const { return (Point){x * _, y * _}; }
    Point operator/(const double &_) const { return (Point){x / _, y / _}; }
    double operator^(const Point &_) const { return x * _.x + y * _.y; } // 点乘
    double operator&(const Point &_) const { return x * _.y - _.x * y; } // 叉乘
    double Len() { return sqrt(*this ^ *this); }
    Point unit() { return *this / Len(); }
};
double area(Point a, Point b, Point c) { return (b - a) & (c - a); }  // 向量 ab 与 向量 ac 的面积
double project(Point a, Point b, Point c) { return ((b - a) ^ (c - a)) / (b - a).Len(); } // 向量 ac 在 向量 ab 上的投影
double Angle(Point a, Point b) { return acos((a ^ b) / a.Len() / b.Len()); } // 两个向量夹角
double Distance(Point a, Point b)
{
    double dx = a.x - b.x, dy = a.y - b.y;
    return sqrt(dx * dx + dy * dy);
}
Point rotate(Point a, double theta) // 逆时针旋转
    { return (Point){a.x * cos(theta) - a.y * sin(theta), a.x * sin(theta) + a.y * cos(theta)}; }
struct Segment // 起点+方向表示或者起点+终点表示
{
    Point x, y;
    void input() { x.input(), y.input(); }
    double Angle() { return atan2(y.y - x.y, y.x - x.x); }
    bool operator<(Segment &_) // 按照极角排序
    {
        double a = Angle(), b = _.Angle();
        if (!cmp(a, b)) return sign(area(x, y, _.y)) == -1;
        return a < b;
    }
    int onSegment(Point p) { return -sign(area(x, y, p)); }  // 点相对线段的位置 (-1左，0上，1右)
};
Segment Vertical(Point a, Point b) // 中垂线 注意得到的答案是一个起点+方向的线段
    { return (Segment){(a + b) / 2, rotate(b - a, pi / 2)}; }
/*Point Intersection(Segment a, Segment b) // 注意两条线段是起点+方向的表示
{
    Point u = a.x - b.x;
    double t = ((b.y) & u) / (a.y & b.y);
    return a.x + a.y * t;
}*/
Point Intersection(Segment a, Segment b) // 注意两条线段是起点+终点的表示
{
    Point u = a.x - b.x;
    double t = ((b.y - b.x) & u) / ((a.y - a.x) & (b.y - b.x));
    return a.x + (a.y - a.x) * t;
}
bool ExistIntersection(Segment a, Segment b)
    { return sign((a.y - a.x) & (b.x - a.x)) * sign((a.y - a.x) & (b.y - a.x)) <= 0 && sign((b.y - b.x) & (a.x - b.x)) * sign((b.y - b.x) & (a.y - b.x)) <= 0; }
struct Circle { Point p; double r; };
Circle circle(Point a, Point b, Point c)
{
    Segment u = Vertical(a, b), v = Vertical(a, c);
    Point p = Intersection(u, v);
    return (Circle){p, Distance(p, a)};
}
struct Polygon
{
    int n;
    Point p[N];
    void input()
    {
        scanf("%d", &n);
        for (int i = 0; i < n; i++) p[i].input();
    }
    void insert(Point a) { p[n++] = a; }
    double PolygonArea()
    {
        double res = 0;
        for (int i = 1; i < n - 1; i++)
            res += area(p[0], p[i], p[(i + 1) % n]);
        return res / 2;
    }
    double PolygonCircumference()
    {
        double res = 0;
        for (int i = 0; i < n; i++)
            res += Distance(p[i], p[(i + 1) % n]);
        return res;
    }
    Polygon ConvexHull() // 凸包
    {
        static bool vis[N];
        sort(p, p + n);
        static int top = 0, stk[N];
        for (int i = 0; i < n; i++)
        {
            while (top >= 2)
            {
                int t = (Segment){p[stk[top - 1]], p[stk[top]]}.onSegment(p[i]);
                if (t == 1) vis[stk[top--]] = false;
                else if (t == 0) --top;
                else break;
            }
            vis[stk[++top] = i] = true;
        }
        vis[0] = false;
        for (int i = n - 1; i >= 0; i--)
        {
            if (vis[i]) continue;
            while (top >= 2 && (Segment){p[stk[top - 1]], p[stk[top]]}.onSegment(p[i]) >= 0) top--; // 不保留共线的点
            vis[stk[++top] = i] = true;
        }
        Polygon res;
        res.n = 0;
        for (int i = 1; i < top; i++)
            res.insert(p[stk[i]]);
        return res;
    }
    Circle CircleCover() // 最小圆覆盖
    {
        random_shuffle(p, p + n);
        Circle c;
        c.p = p[0], c.r = 0;
        for (int i = 1; i < n; i++)
        {
            if (cmp(c.r, Distance(c.p, p[i])) >= 0) continue;
            c.p = p[i], c.r = 0;
            for (int j = 0; j < i; j++)
            {
                if (cmp(c.r, Distance(c.p, p[j])) >= 0) continue;
                c.p = (p[i] + p[j]) / 2;
                c.r = Distance(p[i], p[j]) / 2;
                for (int k = 0; k < j; k++)
                {
                    if (cmp(c.r, Distance(c.p, p[k])) >= 0) continue;
                    c = circle(p[i], p[j], p[k]);
                }
            }
        }
        return c;
    }
    Polygon RectangleCover() // 最小矩形覆盖
    {
        p[n] = p[0];
        double mn = 1e20;
        Polygon res;
        res.n = 4;
        for (int i = 0, a = 2, b = 1, c = 2; i < n; i++)
        {
            Point d = p[i], e = p[i + 1];
            while (cmp(area(d, e, p[a]), area(d, e, p[a + 1])) < 0) a = (a + 1) % n;
            while (cmp(project(d, e, p[b]), project(d, e, p[b + 1])) < 0) b = (b + 1) % n;
            if (!i) c = a;
            while (cmp(project(d, e, p[c]), project(d, e, p[c + 1])) > 0) c = (c + 1) % n;
            Point x = p[a], y = p[b], z = p[c];
            double h = area(d, e, x) / (e - d).Len();
            double w = ((y - z) ^ (e - d)) / (e - d).Len();
            if (h * w < mn)
            {
                mn = h * w;
                res.p[0] = d + (e - d).unit() * project(d, e, y);
                res.p[3] = d + (e - d).unit() * project(d, e, z);
                Point u = rotate(e - d, pi / 2).unit();
                res.p[1] = res.p[0] + u * h;
                res.p[2] = res.p[3] + u * h;
            }
        }
        return res;
    }
    double Diameter() // 凸包的直径
    {
        if (n <= 2) return Distance(p[0], p[n - 1]);
        double res = 0;
        p[n] = p[0];
        for (int i = 0, j = 1; i < n; i++)
        {
            Point d = p[i], e = p[i % n];
            while (cmp(area(d, e, p[j]), area(d, e, p[j % n])) == -1) j = (j + 1) % n;
            res = max(res, max(Distance(d, p[j]), Distance(e, p[j])));
        }
        return res;
    }
};
Polygon HalfPlaneIntersection(int idx, Segment *l) // 半平面交
{
    sort(l, l + idx);
    static int hd = 0, tl = -1, q[N];
    for (int i = 0; i < idx; i++)
    {
        if (i && !cmp(l[i].Angle(), l[i - 1].Angle())) continue;
        while (hd + 1 <= tl && l[i].onSegment(Intersection(l[q[tl - 1]], l[q[tl]])) >= 0) tl--;
        while (hd + 1 <= tl && l[i].onSegment(Intersection(l[q[hd + 1]], l[q[hd]])) >= 0) hd++;
        q[++tl] = i;
    }
    while (hd + 1 <= tl && l[q[hd]].onSegment(Intersection(l[q[tl - 1]], l[q[tl]])) >= 0) tl--;
    while (hd + 1 <= tl && l[q[tl]].onSegment(Intersection(l[q[hd + 1]], l[q[hd]])) >= 0) hd++;
    q[++tl] = q[hd];
    Polygon res;
    res.n = 0;
    for (int i = hd; i < tl; i++)
        res.insert(Intersection(l[q[i]], l[q[i + 1]]));
    return res;
}
int main()
{
    return 0;
}