XCPC wiki - 计算几何

平面计算几何

using T = double;
const double eps = 1e-9, pi = acos(-1.0);
int sgn(T x) {
    return (fabs(x) <= eps) ? 0 : (x < 0 ? -1 : 1);
}

struct P {
	T x, y;
	P() : x(0), y(0) {}
	P(T x, T y) : x(x), y(y) {}
	// P(T r, double alpha, int _) : x(r * cos(alpha)), y(r * sin(alpha)) {}
	bool operator == (P o) { return sgn(x - o.x) == 0 && sgn(y - o.y) == 0; }
    P operator + (P o) { return P(x + o.x, y + o.y); }
    P operator - (P o) { return P(x - o.x, y - o.y); }
    T operator * (P o) { return x * o.y - y * o.x; }    // 叉乘 
    T operator ^ (P o) { return x * o.x + y * o.y; }    // 点乘 
    friend T abs(P o) { return sqrt(o.x * o.x + o.y * o.y); }
    // T angle() { T t = atan2(y, x); if (sgn(t) < 0) t += 2 * pi; return t; }
    // T quanrant() { return (y < 0) << 1 | (x < 0) ^ (y < 0); }
    // bool cmp1(P o) { return sgn(x - o.x) == 0 ? y < o.y : x < o.x; }
    // bool cmp2(P o) { return sgn(y - o.y) == 0 ? x < o.x : y < o.y; }
    // bool cmp3(P o) { return sgn(angle() - o.angle()) ? angle() < o.angle() : x < o.x; }
    // bool operator<(P o) { return cmp1(o); }
};

// 特别提醒：如果 T 是整型 不能写 /2
T Area(P a, P b, P c) {
    return (b - a) * (c - a) / 2;
}

struct L {
    P p1, p2, v;
    L(P P1, P P2) { p1 = P1, p2 = P2, v = P2 - P1; }
    friend T dis(P p, L l) { return (l.v * (p - l.p1)) / abs(l.v); } // 点到直线的距离 
    friend bool Int1(L L1, L L2) {
        if (sgn(L1.v * L2.v) == 0)
            return sgn(dis(L1.p1, L2)) == 0;
        return 1;
    } // 直线相交 
    friend bool Int2(L L1, L L2) {
        if (sgn(L1.v * L2.v) == 0 && sgn(dis(L1.p1, L2)) == 0)
            if (max(L1.p1.x, L1.p2.x) < min(L2.p1.x, L2.p2.x) || max(L2.p1.x, L2.p2.x) < min(L1.p1.x, L1.p2.x))
                return 0;
            else return 1;
        if (dis(L1.p1, L2) * dis(L1.p2, L2) <= 0 && dis(L2.p1, L1) * dis(L2.p2, L1) <= 0)
            return 1;
        return 0;
    } // 朴素线段相交 
    friend bool Int3(L L1, L L2) {
        if (dis(L1.p1, L2) == 0 || dis(L1.p2, L2) == 0 || dis(L2.p1, L1) == 0 || dis(L2.p2, L1) == 0)
            return 0;
        if (dis(L1.p1, L2) * dis(L1.p2, L2) < 0 && dis(L2.p1, L1) * dis(L2.p2, L1) < 0)
            return 1;
        return 0;
    } // 线段规范相交 
};

pair<vector<int>, vector<P>> Andrew(vector<P>& pos) {
    const int n = pos.size();
    sort(pos.begin(), pos.end());
    vector<int> stk(n + 1), vis(n);
    int top = -1;
    stk[++top] = 0;
    for (int i = 1; i < n; i++) {
        while (top > 0 && sgn((pos[stk[top]] - pos[stk[top - 1]]) * (pos[i] - pos[stk[top]])) <= 0) {
            vis[stk[top--]] = 0;
        }
        vis[stk[++top] = i] = 1;
    }
    int tmp = top;
    for (int i = n - 2; i >= 0; i--) {
        if (vis[i]) continue;
        while (top > tmp && sgn((pos[stk[top]] - pos[stk[top - 1]]) * (pos[i] - pos[stk[top]])) <= 0) {
            vis[stk[top--]] = 0;
        }
        vis[stk[++top] = i] = 1;
    }
    vector<int> inConvex(n);
    vector<P> Convex;
    for (int i = 0; i <= top; i++) {
        inConvex[stk[i]] = 1;
        Convex.push_back(pos[stk[i]]);
    }
    return pair(inConvex, Convex);
}

T Diameter(vector<P>& pos) { // 旋转卡壳 
    const int n = pos.size();
    if (n == 3) return abs(pos[1] - pos[0]);
    T res = 0;
    for (int l = 0, r = 0; l < n - 1; l++) {
        L line(pos[l], pos[l + 1]);
        while (sgn(fabs(dis(pos[r], line)) - fabs(dis(pos[(r + 1) % n], line))) <= 0) {
            r = (r + 1) % n;
        }
        res = max(res, max(abs(pos[l] - pos[r]), abs(pos[l + 1] - pos[r])));
    }
    return res;
}

2024 牛客多校 9H - Two Convex Polygons

[[…/contests/2024 牛客多校 /2024-08-13：2024 牛客暑期多校训练营 9|2024-08-13：2024 牛客暑期多校训练营 9]]

long double res = 0;

int m;
cin >> m;
vector<P> a(m);
for (int i = 0; i < m; i++) {
	cin >> a[i].x >> a[i].y;
}
for (int i = 0; i < m; i++) {
	res += abs(a[i] - a[(i + 1) % m]);
}

int n;
cin >> n;
vector<P> b(n);
for (int i = 0; i < n; i++) {
	cin >> b[i].x >> b[i].y;
}

b = Andrew(b);
double v = Diameter(b);
res += 2 * v * pi;

cout << res << endl;

2023SDCPC - M. Computational Geometry

给定一个有$n$ 个顶点的凸多边形$P$，选择$P$ 的三个顶点$a$，$b$ 和$c$（按逆时针顺序）。要求在$b$ 沿逆时针方向到$c$ 之间恰有$k$ 条边（也就是说，$a$ 不是这$k$ 条边的端点）。将由线段$ab$，$ac$，以及$b$ 和$c$ 之间的$k$ 条边围成的$(k + 2)$ 边形记作$Q$。求$Q$ 可能的最大面积。

int n, k;
cin >> n >> k;

vector<P> pos(n);
for (int i = 0; i < n; ++i) {
	cin >> pos[i].x >> pos[i].y;
}

double area = 0, ans = 0;
for (int i = 0; i < k; ++i) {
	area += Area(pos[0], pos[i], pos[(i + 1) % n]);
}

for (int i = 0; i <= n; ++i) {
	int a = i + k, b = i + n;
	auto check = [&](int mid) {
		return Area(pos[a % n], pos[b % n], pos[mid % n]);
	};
	int lt = a, rt = b;
	while (lt < rt) {
		int lm = lt + (rt - lt) / 3, rm = lt + (rt - lt) * 2 / 3;
		if (check(lm) < check(rm)) lt = lm + 1;
		else rt = rm;
	}
	ans = max(ans, area + check(lt));
	area += Area(pos[(i + k + 1) % n], pos[i % n], pos[(i + k) % n]);
	area -= Area(pos[(i + k + 1) % n], pos[i % n], pos[(i + 1) % n]);
}

cout << ans << '\n';

Space Ant

一只蚂蚁，只会向前或左转，现在给出平面上很多个点，从某个$(0,y)$ 开始，求解一种走法，能使得蚂蚁能经过的点最多，每个顶点该蚂蚁只能经过一次，且所行走的路线不能发生交叉。

int n;
cin >> n;

deque<P> pos(n);
for (int i = 0; i < n; i++) {
	int id;
	cin >> id;
	cin >> pos[i].x >> pos[i].y;
}
auto old = pos;
sort(pos.begin(), pos.end(), [&](P a, P b) {
	return a.y < b.y;
});

vector<P> res(n);
for (int t = 0; t < n; ++t) {
	res[t] = pos[0];
	pos.pop_front();
	sort(pos.begin(), pos.end(), [&](P a, P b) {
		return (a - res[t]) * (b - res[t]) > 0;
	});
}

cout << n << ' ';
for (int i = 0; i < n; ++i) {
	for (int j = 0; j < n; ++j) {
		if (res[i] == old[j]) {
			cout << j + 1 << " ";
			break;
		}
	}
}
cout << '\n';

2023 ICPC 济南 M. Almost Convex

求凹一次的凹包数量。$n \leqslant 2\times 10^{3}$。(Link, [[…/contests/VP/2024-10-07-Autumn6-2023ICPC 济南 |2024-10-07-Autumn6-2023ICPC 济南]])

WA：排序应为极角排序，而不是按横坐标排序。

int n;
cin >> n;

vector<P> pos(n);
for (int i = 0; i < n; ++i) {
	cin >> pos[i].x >> pos[i].y;
}

sort(pos.begin(), pos.end(), [&](const P& x, const P& y) {
	return x.x == y.x ? x.y < y.y : x.x < y.x;
});

vector<int> stk(n + 1), vis(n);
int top = -1;
stk[++top] = 0;
for (int i = 1; i < n; i++) {
	while (top > 0 && sgn((pos[stk[top]] - pos[stk[top - 1]]) * (pos[i] - pos[stk[top]])) <= 0) {
		vis[stk[top--]] = 0;
	}
	vis[stk[++top] = i] = 1;
}
int tmp = top;
for (int i = n - 2; i >= 0; i--) {
	if (vis[i]) continue;
	while (top > tmp && sgn((pos[stk[top]] - pos[stk[top - 1]]) * (pos[i] - pos[stk[top]])) <= 0) {
		vis[stk[top--]] = 0;
	}
	vis[stk[++top] = i] = 1;
}

vector<int> inConvex(n);
vector<P> Convex;
for (int i = 0; i <= top; i++) {
	inConvex[stk[i]] = 1;
	Convex.push_back(pos[stk[i]]);
}

vector<P> npos;
for (int i = 0; i < n; ++i) {
	if (!inConvex[i]) {
		npos.push_back(pos[i]);
	}
}
pos = npos;
n = pos.size();

ll ans = 1;
for (int i = 1; i < Convex.size(); ++i) {
	P S = Convex[i], T = Convex[i - 1];
	sort(pos.begin(), pos.end(), [&](const P& x, const P& y) {
		return Area(S, x, y) > 0;
	});
	vector<int> stk(n + 1);
	int top = -1;
	for (int i = 0; i < n; i++) {
		while (top >= 0 && Area(pos[i], pos[stk[top]], T) > 0) {
			top--;
		}
		if (top == -1 || Area(pos[i], pos[stk[top]], S) < 0) {
			stk[++top] = i;
		}
	}
	ans += top + 1;
}

cout << ans << "\n";

2024 ICPC 网络赛 (I) K. Minimum Euclidean Distance

给定一个凸包，多次询问，每次给出一个圆，在凸包内（包括边界）找到一个点，使得这个点到圆内所有点的期望平方距离最小。

答案是$\cfrac{1}{2}r^{2}+d^{2}$，其中$r$ 为圆的半径，$d$ 为答案点到圆心的距离。问题转化为求凸包内一个点到圆心距离最近。如果圆心在凸包内$d=0$，否则计算圆心到凸包每一条边的距离。

int n, q;
cin >> n >> q;

vector<P> pos(n);
for (auto& [x, y] : pos) {
    cin >> x >> y;
}

while (q--) {
    double x1, y1, x2, y2;
    cin >> x1 >> y1 >> x2 >> y2;

    P cir((x1 + x2) / 2, (y1 + y2) / 2);
    double r = abs(P((x1 - x2) / 2, (y1 - y2) / 2));

    bool in = 1;
    for (int i = 0; i < n; ++i) {
        if (Area(cir, pos[i], pos[(i + 1) % n]) < 0) {
            in = 0;
        }
    }

    double minn = 8e18;
    if (in) {
        minn = 0;
    } else for (int i = 0; i < n; ++i) {
        auto P = pos[i], Q = pos[(i + 1) % n];
        if (((P - Q) ^ (P - cir)) >= 0 && ((Q - P) ^ (Q - cir)) >= 0) {
            minn = min(minn, abs(2 * Area(cir, P, Q) / abs(P - Q)));
        }
        minn = min(minn, abs(abs(P - cir)));
    }

    cout << fixed << setprecision(10) << (minn * minn + r * r / 2) << '\n';
}

2024 ICPC 南京 M. Ordainer of Inexorable Judgment

思路 临界状态时$2n$ 条切线，算出这些切线的倾斜角，找到多边形所在半平面内的最大最小的倾斜角。

一些 Tip：

• 如何求切线：可以按高中解析几何思路暴力解方程，也可以直接用角度关系得到切点的倾斜角。

• 初始方向的处理：将所有切线的倾斜角减去初始倾斜角。

• 如何找所在半平面：排序后看相邻两项，如果超过$\pi$，这两项就是分界线。

• 半平面跨过$x=0$ 如何处理：取反，求打不到的区间，这个区间是不会跨过$x=0$ 的。

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

using T = double;
const double eps = 1e-9, pi = acos(-1.0);
int sgn(T x) {
	return (fabs(x) <= eps) ? 0 : (x < 0 ? -1 : 1);
}

struct P {
	T x, y;
	P() : x(0), y(0) {}
	P(T x, T y) : x(x), y(y) {}
	P(T r, double alpha, int _) : x(r * cos(alpha)), y(r * sin(alpha)) {}
	bool operator==(P o) { return sgn(x - o.x) == 0 && sgn(y - o.y) == 0; }
	P operator+(P o) { return P(x + o.x, y + o.y); }
	P operator-(P o) { return P(x - o.x, y - o.y); }
	T operator*(P o) { return x * o.y - y * o.x; }
	T operator^(P o) { return x * o.x + y * o.y; }
	friend double abs(P o) { return sqrt(o.x * o.x + o.y * o.y); }
	double angle() { double res = atan2(y, x); if (sgn(res) < 0) res += (2 * pi); return res; }
};

int main() {
	int n;
	cin >> n;

	P P0;
	cin >> P0.x >> P0.y;
	double t0 = P0.angle();

	int rad, t;
	cin >> rad >> t;

	vector<double> alpha;  // 存所有的倾斜角 
	for (int i = 0; i < n; i++) {
		P Q;
		cin >> Q.x >> Q.y;
		double phi = acos(rad / abs(Q));
		alpha.push_back((Q - P(rad, (Q.angle() + phi), 1)).angle());
		alpha.push_back((Q - P(rad, (Q.angle() - phi), 1)).angle());
	}

	for (int i = 0; i < 2 * n; i++) {
		alpha[i] -= t0;
		if (alpha[i] < 0) alpha[i] += 2 * pi;
	}

	double L = -1, R = -1;
	sort(alpha.begin(), alpha.end());
	bool flag = 1;
	if (alpha[0] + pi > alpha.back()) {
		L = alpha[0], R = alpha.back();
		flag = 0;
	} else for (int i = 1; i < 2 * n; i++) {
		if (alpha[i - 1] + pi < alpha[i]) {
			L = alpha[i - 1], R = alpha[i];
		}
	}

	double q = floor(t / (2 * pi));  // quotient
	double r = fmod(t, (2 * pi));    // remainder
	double ans = q * (R - L);  // full circles
	if (r > L) ans += r - L;
	if (r > R) ans -= r - R;   // partial circle

	cout << fixed << setprecision(15);
	if (flag) {
		cout << (t - ans) << endl;
	} else {
		cout << ans << endl;
	}

	return 0;
}

2024 CCPC 哈尔滨 B. Concave Hull

int n;
cin >> n;
vector<P> p(n);
for (int i = 0; i < n; i++) {
    cin >> p[i].x >> p[i].y;
}

auto [inConvex, Convex] = Andrew(p);
Convex.pop_back();
vector<P> tmp;
for (int i = 0; i < p.size(); i++) {
    if (!inConvex[i]) {
        tmp.push_back(p[i]);
    }
}

p = tmp;
if (p.empty()) {
    cout << "-1\n";
    continue;
}
auto [inConvex2, Convex2] = Andrew(p);
if (Convex2.size() > 1) Convex2.pop_back();

int m = Convex.size(), k = Convex2.size();
ll area = 0;

ll res = 8e18;   // WA
for (int i = 0, j = 0; i < m; i++) {
    while (Area(Convex[i], Convex[(i + 1) % m], Convex2[(j + 1) % k]) < Area(Convex[i], Convex[(i + 1) % m], Convex2[j])) {
        j = (j + 1) % k;
    }
    res = min(res, Area(Convex[i], Convex[(i + 1) % m], Convex2[j]));
    area += Area(Convex[i], Convex[(i + 1) % m], Convex[0]);
}
cout << (area - res) << endl;

三维计算几何

模板：2024 杭电多校 4001 - 超维攻坚

求满足如下条件的点集$S$ 的数量：

• $S$ 的凸包表面与内部的所有点都在$S$ 中。

数据范围：$n \leqslant 15$。

#include <iostream>
#include <algorithm>
#include <vector>
using namespace std;

const int eps = 0;
int sgn(int x) {
    return (abs(x) <= eps) ? 0 : (x < 0 ? -1 : 1);
}

struct Point {
    int x, y, z;
    Point(int X = 0, int Y = 0, int Z = 0) { x = X, y = Y, z = Z; }
    Point operator - (const Point& o) const { return Point(x - o.x, y - o.y, z - o.z); }
    Point operator * (const Point& o) const { return Point(y * o.z - z * o.y, z * o.x - x * o.z, x * o.y - y * o.x); }
    int operator ^ (const Point& o) const { return x * o.x + y * o.y + z * o.z; }
    bool operator < (const Point& o) const {
        return x == o.x ? y == o.y ? z < o.z : y < o.y : x < o.x;
    }
};

// 点到面的距离 
int Point2Plane(const Point& a, const Point& b, const Point& c, const Point& d) {
    return ((b - a) * (c - a)) ^ (d - a);
}

// 三点共线 
bool colinear(const Point& a, const Point& b, const Point& c) {
    Point n = (a - b) * (b - c);   // 法向量 
    return !sgn(n.x) && !sgn(n.y) && !sgn(n.z);
}

struct Plane {
    Point a, b, c;
};

const int inf = ~0U >> 1;

int main() {
    int t;
    cin >> t;

    while (t--) {
	    int n;
	    cin >> n;
	
	    vector<Point> p(n);
	    for (int i = 0; i < n; i++) {
	        cin >> p[i].x >> p[i].y >> p[i].z;
	    }
	
	    vector<int> good(1 << n);
	    for (int bit = 1; bit < 1 << n; bit++) {   // 枚举所有状态 
	        vector<Plane> Convex;
	        bool coplanar = 1;   // 是否所有点共面 
	
	        // 枚举举三个点 ABC 得到所有的面 
	        for (int i = 0; i < n; i++) if (bit >> i & 1) {
	            for (int j = 0; j < i; j++) if (bit >> j & 1) {
	                for (int k = 0; k < j; k++) if (bit >> k & 1) {
	                    auto& A = p[i], & B = p[j], & C = p[k];
	                    if (colinear(A, B, C)) continue;  // 忽略三点共线 
	                    bool ne = 0, po = 0;
	                    for (int o = 0; o < n; o++) if (bit >> o & 1) {
	                        int tmp = Point2Plane(A, B, C, p[o]);
	                        if (tmp < 0) ne = 1, coplanar = 0;  // o 点与 ABC 不共面 
	                        if (tmp > 0) po = 1, coplanar = 0;  // o 点与 ABC 不共面 
	                    }
	                    if (ne && po) continue;  // 有不在同一侧的点 这个面不是凸包的面 
	                    // 所有点都在 ABC 同一侧 
	                    if (ne) Convex.push_back({ A, C, B });
	                    else Convex.push_back({ A, B, C });
	                }
	            }
	        }
	
	        // 找到一个凸包 把凸包中间的点都标记上 
	        int mask = bit;
	        if (Convex.empty()) {  // 凸包退化为一维线段 
	            Point L{ inf, inf, inf }, R{ -inf, -inf, -inf };  // 线段的两个端点 
	            for (int i = 0; i < n; i++) if (bit >> i & 1) {
	                auto& P = p[i];
	                if (P < L) L = P;
	                if (R < P) R = P;
	            }
	            // 判断点 p 是否在 L 和 R 构成的线段上 
	            for (int i = 0; i < n; i++) if (!(bit >> i & 1)) {
	                auto& P = p[i];
	                if (!colinear(L, R, P)) continue;
	                if ((P.x - L.x) * (P.x - R.x) > 0) continue;
	                if ((P.y - L.y) * (P.y - R.y) > 0) continue;
	                if ((P.z - L.z) * (P.z - R.z) > 0) continue;
	                mask |= 1 << i;
	            }
	        }
	        else if (coplanar) {  // 凸包退化为二维平面凸包 
	            for (int i = 0; i < n; i++) if (!(bit >> i & 1)) {
	                auto& P = p[i];
	                // 判断点是否在凸包内 只需判断是否在某个三角形中 
	                bool in = false;
	                for (auto& [A, B, C] : Convex) {
	                    if (sgn(Point2Plane(A, B, C, P))) continue;
	                    if ((((B - A) * (C - A)) ^ ((B - A) * (P - A))) < 0) continue;
	                    if ((((C - B) * (A - B)) ^ ((C - B) * (P - B))) < 0) continue;
	                    if ((((A - C) * (B - C)) ^ ((A - C) * (P - C))) < 0) continue;
	                    in = true;
	                }
	                if (in) mask |= 1 << i;
	            }
	        }
	        else {  // 三维凸包 
	            for (int i = 0; i < n; i++) if (!(bit >> i & 1)) {
	                auto& P = p[i];
	                // 判断点是否在凸包内 
	                bool in = true;
	                for (auto& [A, B, C] : Convex) {
	                    if (sgn(Point2Plane(A, B, C, P)) == -1) in = false;
	                }
	                if (in) mask |= 1 << i;
	            }
	        }
	        good[mask] = 1;
	    }
	
	    int ans = 0;
	    for (int i = 1; i < 1 << n; i++) {
	        ans += good[i];
	    }
	    cout << ans << "\n";
    }

    return 0;
}