[笔记] 计算几何

其他的咕咕咕了，不喜欢计算几何 qaq

基本定义

向量

• (overrightarrow{PQ}=Q-P)

叉积

二维 a.x * b.y-a.y * b.x

三维

sum += a.x * b.y * c.z + a.y * b.z * c.x + a.z * b.x * c.y

sum -= a.z * b.y * c.x + a.y * b.x * c.z + a.x * b.z * c.y

基本操作

将一个向量旋转一定角度

Vec rotate(Vec a, LD b){ 	LD s = sin(b), c = cos(b); 	return {a.x * c - a.y * s, a.x * s + a.y * c}; } 

判断一个点在直线的哪边

有直线上的一点 (P)，直线的方向向量 (v)，想知道 (Q) 在直线哪边：

利用叉积的性质，若 (overrightarrow{PQ}times v >0)，则 (Q) 在直线逆时针方向，否则在顺时针方向。

bool Left(Vec a, Line b){ return cross(a - b.p , b.v) > 0; } 

判断两圆之间的关系

返回的是公切线个数

int Pos(Circle a, Circle b){ 	LD dis = Dis(a.O, b.O); if(a.r < b.r) swap(a, b); 	if(dis > a.r + b.r) return 4; 	if(dis == a.r + b.r) return 3; 	if(dis > a.r - b.r) return 2; 	if(dis == a.r - b.r) return 1; 	return 0; } 

求两条直线的交点

有直线 (AB,CD)，求交点 (E)：

首先确定是否只有一个交点，然后因为记录的是直线上的一个点和直线的方向向量，所以只需要知道这个点与交点的距离 (l) ，再将这个点沿方向向量平移 (l) 个单位长度即可。利用正弦定理求解

由上图可知，(|mathbf{atimes b|=|a||b|sinbeta,|utimes b|=|u||b|sintheta})。

作商得：

交点即点 (B+Tmathbf{a})。

Vec Inter(Line a, Line b){ return a.p + cross(b.v, (b.p - a.p)) / cross(b.v, a.v) * a.v;}  

求两圆的交点

注意要先判断有没有交点

pair <Vec, Vec> Inter(Circle a, Circle b){ 	LD x = a.r, y = b.r, z = Dis(a.O, b.O); 	LD tar1 = acos((x * x + z * z - y * y) / (2 * x * z)); 	Vec i1 = a.O + rotate(Line(a.O, b.O), tar1) / z * x; 	LD tar2 = acos((y * y + z * z - x * x) / (2 * y * z)); 	Vec i2 = b.O + rotate(Line(b.O, a.O), tar2) / z * y; 	return {i1, i2}; } 

求多边形周长

double PloygonDis(Vec *a, int n){ 	double sum = 0; 	lfor(i, 1, n) sum += Dis(a[i], a[i % n + 1]); 	return sum; } 

求多边形面积

double PloygonArea(Vec *a, int n){ 	double sum = 0; 	lfor(i, 3, n) sum += cross(a[i - 1] - a[1], a[i] - a[1]); 	return sum / 2; } 

极角排序

实现

1. 使用 atan2(y, x) 函数，返回值的范围是 ([-pi,pi])；
2. 使用叉积大于 (0) 的性质，注意因为叉积无法判 180 度以上，所以可能要结合象限排序；

凸包

二维凸包

特殊结论

• 二维平面四个点求凸包面积 (rightarrow) 任选三个点面积之和 / 2

  二维平面三个点面积 (rightarrow) 二个二维向量行列式值的绝对值 / 2

• 三维空间五个点求凸包体积 (rightarrow) 任选四个点体积之和 / 2

  三维空间四个点体积 (rightarrow) 三个三维向量行列式值的绝对值 / 6

半平面交

定义

有若干条有向直线，要求保留每条直线其中一侧，求最后保留的范围。

离线 (O(nlog n))

用有向直线（一个点和一个方向向量）表示半平面，以下默认半平面在有向直线的左侧。
对有向直线按方向向量的极角排序，维护一个双端队列，存储当前构成半平面的直线以及相邻两直线的交点。
每次加入一条有向直线，如果队首 / 队尾的交点在直线右侧（用叉积判）则弹掉队首 / 队尾的直线。
需要注意的细节：

1. 加入直线时，先弹队尾，再弹队首。
2. 特判平行直线，在右侧的要弹掉。
3. 最后还要检查队尾交点是否在队首直线的右侧，如果是也要弹掉。
4. 如果题目给出的半平面不一定有限制边界，则应该手动加入一个 INF 边界。

double HPI(Line *a, int n){ 	static deque <Vec> I; while(!I.empty()) I.pop_back();  	static deque <Line> Q; while(!Q.empty()) Q.pop_back(); 	sort(a + 1, a + n + 1); 	Q.push_back(a[1]); 	lfor(i, 2, n){ 		while(!I.empty() && Cross(a[i].v, I.back() - a[i].p) <= 0) I.pop_back(), Q.pop_back(); 		while(!I.empty() && Cross(a[i].v, I.front() - a[i].p) <= 0) I.pop_front(), Q.pop_front(); 		if(a[i].at2 != Q.back().at2) Q.push_back(a[i]); 		else if(Cross(a[i].v, Q.back().p - a[i].p) <= 0){ 			Q.back() = a[i]; if(!I.empty()) I.pop_back(); 		}else continue; 		auto qwq = Q.rbegin(); 		if(Q.size() > 1) I.push_back(Inter(*(++qwq), Q.back())); 	} 	while(!I.empty() && Cross(Q.front().v, I.back() - Q.front().p) <= 0) I.pop_back(), Q.pop_back();  	if(Q.size() > 1) I.push_back(Inter(Q.front(), Q.back())); 	int cnt = 0; static Vec *b = new Vec[I.size() + 1]; 	for(auto x : I) b[++cnt] = x; 	return PloygonDis(b, cnt); } 

积分

对积分的感性理解

• 有一个函数 (f(x))，求其在区间 ([a,b]) 与 (x) 轴围成的面积，(x) 轴上为正，(x) 轴下为负。

• 那么 (int) 类比与 (sum) 符号，同时用 (mathrm{d}x) 表示将 (x) 分成很多很多很小的份。

• 同时也不难意识到，这个空间是封闭的才可以求面积。

如何积分

大部分的函数都是无法精确积分的，于是采用一些公式来逼近。

自适应辛普森积分

公式：

• 二次函数的积分可以精确计算，辛普森积分即是一种拿二次函数来拟合的方式。
• 在二次函数的情况下，该式求出的即准确积分，推导过程。

自适应：

因为 (f(x)) 不是二次函数，那么当然不能直接积分，于是就有了根据误差调整的自适应做法。

LD Ars(LD l, LD r, LD eps, LD val){ 	LD mid = (l + r) / 2; 	LD L = simpson(l, mid), R = simpson(mid, r); 	if(fabs(L + R - val) <= eps) return L + R; 	return Ars(l, mid, eps / 2, L) + Ars(mid, r, eps / 2, R); } 

• 非常直接的想法，如果不满足 eps，就继续细分区间。
• 注意因为要合并两个区间，所以对误差的要求在提高。
• 因为是自适应的，所以复杂度玄学。
• 对于拟合易错的函数，可强制迭代一定层数。

积分在 OI 中的应用

一般用来求各种面积。

[CQOI2005]三角形面积并

没封装的简陋玩意

const LD Pi = acos(-1);  struct Vec{ LD x, y; }; void Out(Vec a){ cerr << a.x << ' ' << a.y << endl; } void In(Vec &a){ scanf("%Lf%Lf", &a.x, &a.y); } bool operator ==(Vec a, Vec b){ return a.x == b.x && a.y == b.y; } bool operator !=(Vec a, Vec b){ return a.x != b.x || a.y != b.y; } bool operator <(Vec a, Vec b){ return atan2(a.y, a.x) < atan2(b.y, b.x); } LD atan2(Vec a){ return atan2(a.y, a.x); } LD Cross(Vec a, Vec b){ return a.x * b.y - a.y * b.x; } LD Dis(Vec a, Vec b){ return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)); } Vec operator -(Vec a, Vec b){ return {a.x - b.x, a.y - b.y}; } Vec operator +(Vec a, Vec b){ return {a.x + b.x, a.y + b.y}; } Vec operator /(Vec a, LD b){ return {a.x / b, a.y / b}; } Vec operator *(Vec a, LD b){ return {a.x * b, a.y * b}; } Vec operator *(LD b, Vec a){ return {a.x * b, a.y * b}; } Vec rotate(Vec a, LD b){ 	LD s = sin(b), c = cos(b); 	return {a.x * c - a.y * s, a.x * s + a.y * c}; }  struct Line{  	Vec p, v; LD at2;  	Line(){} 	Line(Vec a, Vec b, LD c){ p = a, v = b, at2 = c; } 	Line(Vec a, Vec b){ p = a, v = b - a, at2 = atan2(v.y, v.x); }  }; bool operator <(Line a, Line b){ return a.at2 < b.at2; } bool Left(Vec a, Line b){ return Cross(a - b.p , b.v) > 0; } Vec Inter(Line a, Line b){ return a.p + Cross(b.v, (b.p - a.p)) / Cross(b.v, a.v) * a.v;}  Vec operator +(Vec a, Line b){ return {a.x + b.v.x, a.y + b.v.y}; } Line rotate(Line a, LD b){ return {a.p, rotate(a.v, b), a.at2}; } Line operator *(Line a, LD b){ return (Line){a.p, (Vec){a.v.x * b, a.v.y * b}, a.at2}; } Line operator /(Line a, LD b){ return (Line){a.p, (Vec){a.v.x / b, a.v.y / b}, a.at2}; }  struct Circle{ Vec O; LD r; };  void In(Circle &a){ scanf("%Lf", &a.r), In(a.O); } int Pos(Circle a, Circle b){ 	LD dis = Dis(a.O, b.O); if(a.r < b.r) swap(a, b); 	if(dis > a.r + b.r) return 4; 	if(dis == a.r + b.r) return 3; 	if(dis > a.r - b.r) return 2; 	if(dis == a.r - b.r) return 1; 	return 0; } pair <Vec, Vec> Inter(Circle a, Circle b){ 	LD x = a.r, y = b.r, z = Dis(a.O, b.O); 	LD tar1 = acos((x * x + z * z - y * y) / (2 * x * z)); 	Vec i1 = a.O + rotate(Line(a.O, b.O), tar1) / z * x; 	LD tar2 = acos((y * y + z * z - x * x) / (2 * y * z)); 	Vec i2 = b.O + rotate(Line(b.O, a.O), tar2) / z * y; 	return {i1, i2}; }  double PloygonArea(Vec *a, int n){ 	double sum = 0; 	lfor(i, 3, n) sum += Cross(a[i - 1] - a[1], a[i] - a[1]); 	return sum / 2; } double PloygonDis(Vec *a, int n){ 	double sum = 0; 	lfor(i, 1, n) sum += Dis(a[i], a[i % n + 1]); 	return sum; } double HPI(Line *a, int n){ 	static deque <Vec> I; while(!I.empty()) I.pop_back();  	static deque <Line> Q; while(!Q.empty()) Q.pop_back(); 	sort(a + 1, a + n + 1); 	Q.push_back(a[1]); 	lfor(i, 2, n){ 		while(!I.empty() && Cross(a[i].v, I.back() - a[i].p) <= 0) I.pop_back(), Q.pop_back(); 		while(!I.empty() && Cross(a[i].v, I.front() - a[i].p) <= 0) I.pop_front(), Q.pop_front(); 		if(a[i].at2 != Q.back().at2) Q.push_back(a[i]); 		else if(Cross(a[i].v, Q.back().p - a[i].p) <= 0){ 			Q.back() = a[i]; if(!I.empty()) I.pop_back(); 		}else continue; 		auto qwq = Q.rbegin(); 		if(Q.size() > 1) I.push_back(Inter(*(++qwq), Q.back())); 	} 	while(!I.empty() && Cross(Q.front().v, I.back() - Q.front().p) <= 0) I.pop_back(), Q.pop_back();  	if(Q.size() > 1) I.push_back(Inter(Q.front(), Q.back())); 	int cnt = 0; static Vec *b = new Vec[I.size() + 1]; 	for(auto x : I) b[++cnt] = x; 	return PloygonDis(b, cnt); }