计算几何学习与总结

Start on Mar. 5th 2026

学习一下计算几何，听说 icpc 每场都有。

陆续更新中...

参考：
https://oi-wiki.org/geometry/
https://www.luogu.com.cn/article/rnbm18tt
https://www.cnblogs.com/lotus-pjl/p/19968637

计算几何基础

Pick 定理

在格点上，对于一个简单多边形，令多边形内部的点数为 \(a\)，其边界上的点数为 \(b\)，则有：
\(\displaystyle S = a + \frac{b}{2} - 1\)

常用公式与代码

\(\pi = acos(-1.0)\)

正弦、余弦定理

在 \(\triangle \text{ABC}\) 中：\(\frac{a}{\sin A}=\frac{b}{\sin B}=\frac{c}{\sin C}=2R\)。其中，\(R\) 为 \(\triangle \text{ABC}\) 的外接圆半径。

在 \(\triangle \text{ABC}\) 中：\(\cos A = \frac{b^2+c^2-a^2}{2bc}\)。

海伦公式

\(let \ p = \frac{a + b + c}{2}, \\ S = \sqrt{p(p - a)(p - b)(p - c)}\)

点乘与叉乘

\(when\ \ \vec a = (x_1, y_1),\ \vec b = (x_2, y_2):\\ \vec a \cdot \vec b = |a||b| \cos \theta = x_1 x_2 + y_1 y_2 \\ \vec a \times \vec b = |a||b| \sin \theta = x_1 y_2 - y_1 x_2 \)

int operator * (const NODE &rhs) {
	return x * rhs.x + y * rhs.y;
}
int operator ^ (const NODE &rhs) {
	return x * rhs.y - y * rhs.x;
}
求向量的交点 inter()
NODE inter (LINE a, LINE b) {
    double t = cross(b.v, a.x - b.x) / cross(a.y, b.y);
    return NODE {a.x.x + a.y.x * t, a.x.y + a.y.y * t};
}
判断一个点在直线的哪边

求 \(C\) 在 \(\overrightarrow{AB}\) 的左还是右，计算 $\overrightarrow{AB} \times \overrightarrow{AC} $ ：
若为负，则在右侧；若为正，则在左侧；若为 \(0\)，则 \(A, B, C\) 共线。

向量逆时针旋转

将 \(\vec a = (x, y)\) 逆时针旋转 \(r\) 度，变为 \(\vec b = (x\cos r - y \sin r， x \sin r + y \cos r)\)。顺时针取反即可。

注意事项

  1. 使用 acos(x) asin(x) 等函数时，判断 x 会不会因精度等问题超过定义域。若超过，可能会导致 Runtime Error。
  2. 注意 EPS 的取值。

关于double...

from https://www.cnblogs.com/oyking/p/3959905.html

1、在竞赛中，可能存在给一个double多次计算后，非常接近0的情况，但是，它可能是-0.000000000000000001（一下子举不出实际的例子……），这时输出会变成-0.00，在一部分精确比较的题目上可能会出错，解决方案为输出的时候加上一个 EPS（当然不能排除出题人自己煞笔的情况……），即 printf("%f", ans + EPS)。

现在考虑一种情况,题目要求输出保留两位小数。有个case的正确答案的精确值是0.005,按理应该输出0.01,但你的结果可能是0.005000000001(对)，也有可能是0.004999999999(错),如果按照printf(“%.2lf”, a)输出，第二种答案就会错误。解决办法是，如果a为正，则输出a+eps, 否则输出a-eps

ICPC题目输出有个不成文的规定(有时也成文)，不要输出: -0.000那我们首先要弄清，什么时候按printf(“%.3lf\n”, a)输出会出现这个结果。直接给出结果好了：a∈(-0.000499999……, -0.000……1)所以，如果你发现a落在这个范围内，请直接输出0.000。更保险的做法是用sprintf直接判断输出结果是不是-0.000再予处理。

2、如果一个double，要把一个数组里的浮点数累加起来（即 double sum = accumlate(arr.begin(), arr.end(), 0)），如果数组里的数相差可能会比较大，应该从小到大累加。否则可能会发生加了一个很大的浮点数之后，再加一个很小的浮点数，因为精度的原因，跟没加一样（比如一个极端的例子，1e100 + 1e-100还是等于1e100）。因为比较小的数多了，还是会影响到答案的，并不是可以简简单单被忽略的。

3、在计算一个数减去一组数的时候（即double sum = x - accumlate(arr.begin(), arr.end(), 0)），应该先把数组里的数全加起来，再用那个数来减。否则跟上面一样，可能会出现减去了一个很小的数，跟没减一样。

4、二分的时候，EPS设置不当，可能会出现奇怪的问题（特别是有多次二分而对着两次二分要求的精度不一样的时候），这时可以采取另一种二分固定次数的写法：

double binary_search (double l, double r) {
    int t = 100;
    while (t--) {
        double mid = (l + r) / 2.0;
        if (check(mid)) l = mid;
        else r = mid;
    }
    return l;
}

循环次数可按需设置，一般100都够，对时间要求高可以设小一点。

5、在连乘的时候，浮点数可能会丢比较多的精度，此时可以利用公式 x1 * x2 * …… * xn = exp(ln(x1 * x2 * …… * xn)) = exp(ln(x1) + ln(x2) + ... + ln(xn))，取他们的对数相加，再求次幂。

比如在计算阶乘的时候，可以：

double Factorial(int n) {
    double res = 0;
    for(int i = 1; i <= n; ++i)
        res += log(i);
    return exp(res);
}

当然有必要的时候（其实是大多时候）我们可以先不exp()先返回，作完后续运算再exp()

------- 以下 算法部分 -------

凸包

Graham 算法

参考: https://www.luogu.com.cn/article/td3ah746

从最下面的那个点（称其为 1 号点）开始，其他点按照他们与 1 号点的角度排序，进行枚举，开一个栈，当前点进栈，若可以，同时移除栈顶点。

学习到: 使用 \(atan2(y, x) \in (-π, π]\) 而不是 \(atan(x) \in [−\frac{π}{2}, \frac{π}{2}]\)

P2742 【模板】二维凸包 / [USACO5.1] 圈奶牛Fencing the Cows

#include <stdio.h>
#include <algorithm>
#include <math.h>

int n, s[100003], snt;

struct NODE {
	double x, y, z;
} a[100003];

inline double dis (const NODE &x, const NODE &y) {
	double p = x.x - y.x, q = x.y - y.y;
	return sqrt(p * p + q * q);
}

inline double turn (const NODE &x, const NODE &y, const NODE &z) { // p x q 
	return (y.x - x.x) * (z.y - x.y) - (y.y - x.y) * (z.x - x.x); 
} 

inline void kagari () {
	scanf("%d", &n);
	for (int i = 1; i <= n; ++i) scanf("%lf %lf", &a[i].x, &a[i].y);
	for (int i = 2; i <= n; ++i) 
		if (a[i].y < a[1].y || a[i].y == a[1].y && a[i].x < a[1].x) 
			std:: swap(a[1], a[i]);
	
	for (int i = 2; i <= n; ++i) a[i].z = atan2(a[i].y - a[1].y, a[i].x - a[1].x);
	std:: sort(a + 2, a + n + 1, [](const NODE &x, const NODE &y) 
				{ return x.z < y.z || x.z == y.z && dis(a[1], x) < dis(a[1], y); });
	
	snt = 0; s[++snt] = 1, s[++snt] = 2;
	for (int i = 3; i <= n; ++i) {
		while (snt >= 2 && turn(a[s[snt - 1]], a[s[snt]], a[i]) <= 0.0) --snt;
		s[++snt] = i;
	}
	
	double ans = dis(a[s[1]], a[s[snt]]);
	for (int i = 1; i < snt; ++i) ans += dis(a[s[i]], a[s[i + 1]]);
	printf("%.2lf\n", ans);
	return;
}
int main () {
	kagari();
	return 0;
}

半平面交

给一些凸多边形，求他们的交集的面积。

首先要明确 POINT VECTOR LINE 结构体的编写方式，以及利用叉乘判断某点是否在某向量左边 onleft(), 两个向量的交点 inter() 函数的编写方式。

大致步骤如下：

  1. 将各个凸多边形转化为从节点出发的向量
  2. 所有向量依照角度排序，并去重（同角度用叉乘留下靠左的）
  3. 开 deque 存储点与线。依次枚举所有线，对一个新的线，若最后一个点不在这条线左边，即这个点一定不在答案里，则将队尾的点与线移除队列；第一个点同理。然后将 这个线与队尾线的交点 与 这个线 入队
  4. 枚举结束后，判断队列中队尾点是否在队首线的左边，若不在，继续出列。最后将队首与队尾线的交点入队
  5. 计算面积，得出答案

P4196 【模板】半平面交 / [CQOI2006] 凸多边形

#include <stdio.h>
#include <algorithm>
#include <math.h>
#include <queue>
#define ll long long

const double EPS = 1e-7;
inline bool dbe (double x, double y) {
	return fabs(x - y) <= EPS;
}

int n, m;
struct NODE {
	double x, y;
	NODE operator - (const NODE &b) const { return NODE { x - b.x, y - b.y }; }
	bool operator == (const NODE &b) const { return dbe(x, b.x) && dbe(y, b.y); }
	bool operator < (const NODE &b) const { return dbe(x, b.x) ? y < b.y : x < b.x; }
};
struct LINE {
	NODE p, v;
	double z;
	LINE () { }
	LINE (const NODE &x, const NODE &y) { 
		p = x, v = y - x;
		z = atan2(v.y, v.x);
	}
	
} a[1003];
inline double cross (const NODE &x, const NODE &y) { return x.x * y.y - x.y * y.x; }
inline bool onleft (const LINE &x, const NODE &y) { return cross(x.v, y - x.p) >= 0.0; }
NODE inter (LINE a, LINE b) {
    double t = cross(b.v, a.p - b.p) / cross(a.v, b.v);
    return NODE {a.p.x + a.v.x * t, a.p.y + a.v.y * t};
}

inline void kagari () {
	scanf("%d", &m);
	for (int i = 1; i <= m; ++i) {
		int t; double x0, y0, x, y, x00, y00; scanf("%d %lf %lf", &t, &x00, &y00); x0 = x00, y0 = y00;
		for (int j = 2; j <= t; ++j) {
			scanf("%lf %lf", &x, &y);
			a[++n] = LINE(NODE{x0, y0}, NODE{x, y});
			x0 = x, y0 = y;
		}
		a[++n] = LINE(NODE{x, y}, NODE{x00, y00});
	}
	
	std:: sort(a + 1, a + n + 1, [](const LINE &x, const LINE &y) {
		return dbe(x.z, y.z) ? cross(x.v, y.p - x.p) < 0.0 : x.z < y.z;
	});
	int nn = 1;
	for (int i = 2; i <= n; ++i) 
		if (!dbe(a[i].z, a[i - 1].z)) a[++nn] = a[i];
	n = nn;
	
	std:: deque < NODE > p;
	std:: deque < LINE > q;
	for (int i = 1; i <= n; ++i) {
		while (p.size() && !onleft(a[i], p.back())) p.pop_back(), q.pop_back();
		while (p.size() && !onleft(a[i], p.front())) p.pop_front(), q.pop_front();
		if (q.size()) p.push_back(inter(a[i], q.back()));
		q.push_back(a[i]);
	}
	while (p.size() && !onleft(q.front(), p.back())) p.pop_back(), q.pop_back();
	p.push_back(inter(q.front(), q.back()));
	
	if (p.size() <= 2) { puts("0.000"); return; }
	double ans = 0.0;
	for (int i = 1; i < p.size() - 1; ++i) ans += cross(p[i] - p[0], p[i + 1] - p[0]) * 0.5;
	printf("%.3f\n", ans);
	return;
}
int main () {
	kagari();
	return 0;
}

旋转卡壳

09/03/2026

主要是答案单调的思想。
注意这个单调性在何时成立（如第二题的单调性不恒成立，在计算 u 的时候要先 $ u:=q $ 以跳过单调性相反的区间）。

P1452 【模板】旋转卡壳 / [USACO03FALL] Beauty Contest G

题意：平面 n 个点中，求其凸包中两点距离最大值。

做法：枚举凸包上的每条边 \(\overrightarrow{AB}\)，对这条边，找到距其最远的那个点 \(C\)，所有这样的\(|AC|, |BC|\) 的最大值，就是凸包中两点距离的最大。而这个 \(C\) 点的确定由凸包的凸性与答案单调性，时间复杂度和枚举边一同为 \(O(n)\)。

启发：有时候枚举点是错的，我们要转变思路，从凸包上的边入手。

#include <stdio.h>
#include <algorithm>
#include <vector>
#include <stack>

int n;
struct NODE {
	double x, y;
	NODE operator - (const NODE &t) const { return { x - t.x, y - t.y }; }
} a[50003];
inline int dis2 (NODE x, NODE y) {
	return (x.x - y.x) * (x.x - y.x) + (x.y - y.y) * (x.y - y.y);
}
inline int cross (NODE x, NODE y) {
	return x.x * y.y - x.y * y.x;
}

inline auto gettb() {
	for (int i = 2; i <= n; ++i) if (a[i].x < a[1].x || a[i].x == a[1].x && a[i].y < a[1].y) std:: swap(a[1], a[i]);
	std:: sort(a + 2, a + n + 1, [](const NODE &x, const NODE &y) {
		return cross(x - a[1], y - a[1]) == 0 ? (x.x == y.x ? x.y < y.y : x.x < y.x) : cross(x - a[1], y - a[1]) > 0; }
	);
	for (int i = 1; i <= n; ++i) printf("~%lf %lf\n", a[i].x, a[i].y);
	std:: stack < int > s;
	s.push(1); s.push(2);
	for (int i = 3; i <= n; ++i) {
		while (s.size() >= 2) {
			int p = s.top(); s.pop();
			int q = s.top(); 
			if (cross(a[p] - a[q], a[i] - a[p]) > 0) { s.push(p); break; }
		}
		s.push(i);
	}
	std:: vector < NODE > tb;
	while (!s.empty()) tb.push_back(a[s.top()]), s.pop();
	std:: reverse(tb.begin(), tb.end());
	return tb;
}

inline int getdia (std:: vector < NODE > &tb) {
	int m = tb.size(), res = 0, r = 1;
	for (int i = 0; i < m; ++i) {
		while (abs(cross(tb[(i+1) % m] - tb[i], tb[(r + 1) % m] - tb[i])) > 
			   abs(cross(tb[(i+1) % m] - tb[i], tb[r % m] - tb[i]))) ++r;
		res = std:: max(res, std:: max(dis2(tb[i], tb[r % m]), dis2(tb[(i + 1) % m], tb[r % m])));
	}
	return res;
} 

inline void kagari () {
	scanf("%d", &n);
	for (int i = 1; i <= n; ++i) scanf("%lf %lf", &a[i].x, &a[i].y);
	auto tb = gettb();
	int ans = getdia(tb);
	printf("%d\n", ans);
	return;
}
int main () {
	kagari();
	return 0;
}

P3187 [HNOI2007] 最小矩形覆盖

题意：平面内 n 个点内，求能够覆盖所有点的最小面积的矩形的面积和四个顶点坐标。

旋转卡壳模板题，有一些向量的计算。

要确定这个矩形，也就是需要确定一条边与三个点 \(p, q, u\)，并注意保持 \(p\leq q\leq u\)。顺时针依次枚举边，同时确定点，即可得出答案。

#include <stdio.h>
#include <algorithm>
#include <math.h>

int n, m;
struct NODE {
	double x, y;
	NODE operator - () const { return NODE { -x, -y }; }
	NODE operator + (const NODE t) const { return NODE { x + t.x, y + t.y }; }
	NODE operator - (const NODE t) const { return NODE { x - t.x, y - t.y }; }
	NODE operator * (const double z) const { return NODE { x * z, y * z }; }
} a[50003], b[50003];

inline double dot (const NODE x, const NODE y) { return x.x * y.x + x.y * y.y; }
inline double cross (const NODE x, const NODE y) { return x.x * y.y - x.y * y.x; }
inline double dis2 (const NODE x) { return x.x * x.x + x.y * x.y; }
inline double dis (const NODE x) { return sqrt(dis2(x)); }

inline void kagari () {
	scanf("%d", &n);
	for (int i = 1; i <= n; ++i) scanf("%lf %lf", &a[i].x, &a[i].y);
	
	for (int i = 2; i <= n; ++i) if (a[i].x < a[1].x || a[i].x == a[1].x && a[i].y < a[1].y) std:: swap(a[1], a[i]);
	std:: sort(a + 2, a + n + 1, [](const NODE &x, const NODE &y) {
		return cross(x - a[1], y - a[1]) == 0.0 ? (x.x == y.x ? x.y < y.y : x.x < y.x) : cross(x - a[1], y - a[1]) >= 0.0;
	});
	
	b[m++] = a[1]; b[m++] = a[2];
	for (int i = 3; i <= n; ++i) {
		while (m >= 1 && cross(b[m - 1] - b[m - 2], a[i] - b[m - 1]) <= 0.0) --m;
		b[m++] = a[i];
	}
	
	int p = 1, q = 1, u = 1;
	double ans = 10000.0; NODE apos[4];
	b[m] = b[0];
	for (int i = 0; i < m; ++i) {
		NODE t = b[i + 1] - b[i];
		while (dot(t, b[(p + 1) % m] - b[i]) > dot(t, b[p % m] - b[i])) p++;
		if (q < p) q = p;
		while (abs(cross(t, b[(q + 1) % m] - b[i])) > abs(cross(t, b[q % m] - b[i]))) q++;
		if (u < q) u = q;
		while (dot(-t, b[(u + 1) % m] - b[i + 1]) > dot(-t, b[u % m] - b[i + 1])) u++;
		
		double res = (abs(dot(t, b[u % m] - b[i + 1])) + abs(dot(t, b[p % m] - b[i])) - dis2(t)) * abs(cross(t, b[q % m] - b[i])) / dis2(t);
		if (ans > res) {
			ans = res;
			NODE vx1 = t * (dot(t, b[p % m] - b[i]) / dis2(t));
			NODE vx2 = t * (dot(-t, b[u % m] - b[i + 1]) / dis2(t));
			NODE vy = NODE { -t.y, t.x } * (abs(cross(t, b[q % m] - b[i])) / dis2(t));
			apos[0] = b[i] + vx1;
			apos[1] = b[i] + vy + vx1;
			apos[2] = b[i + 1] + vy - vx2;
			apos[3] = b[i + 1] - vx2;
		}
	}
	
	printf("%.5f\n%.5f %.5f\n%.5f %.5f\n%.5f %.5f\n%.5f %.5f\n", 
		ans, apos[0].x, apos[0].y, apos[1].x, apos[1].y, apos[2].x, apos[2].y, apos[3].x, apos[3].y);
	return;
}

int main () {
	kagari();
	return 0;
}

平面最近点对

P1429 平面最近点对（加强版）

采用分治的方法，先按照 x 排序，分为两块后再按照 y 排序，时间复杂度 \(O(n\log^2n)\) （块内排序由 sort 改为归并排序则为 \(O(n\log n)\)）。

NOTE
  1. 由于 \(sqrt()\) 速度较慢，可以计算距离的平方（但是记得在比较距离的时候要平方/开方）。
  2. STL 归并排序 std:: inplace_merge() 的用法：当对数组 \(a[l] \sim a[mid]\ \ \&\ \ a[mid+1]\sim a[r]\) 进行排序（这两部分必须已经有序）时：
    std:: inplace_merge(a + l, a + mid + 1, a + r + 1, cmp);
#include <stdio.h>
#include <algorithm>
#include <math.h>

int n;
struct NODE {
	double x, y;
} a[200003], b[200003];
double dis (const double x, const double y) {
	return sqrt(x * x + y * y);
}
double dis (const NODE x, const NODE y) {
	return dis(x.x - y.x, x.y - y.y);
}

int v[200003], vnt; 
inline double merge (int l, int r) {
	if (l >= r) return 1e18;
	if (l + 1 == r) {
		if (a[l].y > a[r].y) std:: swap(a[l], a[r]); 
		return dis(a[l], a[r]);
	}
	
	// 1. 处理子块 
	int mid = l + r >> 1; NODE amid = a[mid]; 
	double d = merge(l, mid);
	double d2 = merge(mid + 1, r);
	if (d2 < d) d = d2;
	
	
	// 2. 归并排序 
	/* std:: inplace_merge(a + l, a + mid + 1, a + r + 1, [](const NODE &x, const NODE &y) {
		return x.y < y.y;
	}); */
	int p = l, q = mid + 1;
	for (int i = l; i <= r; ++i) 
		if (p <= mid && a[p].y <= a[q].y || q > r) b[i] = a[p++];
		else b[i] = a[q++];
	for (int i = l; i <= r; ++i) a[i] = b[i];
	
	// 3. 合并 
	vnt = 0; 
	for (int i = l; i <= r; ++i) if (fabs(a[i].x - amid.x) < d) v[++vnt] = i;
	for (int i = 1; i < vnt; ++i)
		for (int j = i + 1; j <= vnt; ++j) {
			if (fabs(a[v[i]].y - a[v[j]].y) >= d) break;
			d = std:: min(d, dis(a[v[i]], a[v[j]]));
		}
	return d;
}

inline void kagari () {
	scanf("%d", &n);
	for (int i = 1; i <= n; ++i) scanf("%lf %lf", &a[i].x, &a[i].y);
	std:: sort(a + 1, a + n + 1, [](const NODE &x, const NODE &y) {
		return x.x == y.x ? x.y < y.y : x.x < y.x;
	});
	double ans = merge(1, n);
	printf("%.4f\n", ans);
	return;
}
int main () {
	kagari();
	return 0;
}
posted @ 2026-03-09 19:36  dbg_8  阅读(14)  评论(0)    收藏  举报