因为懒得写博客所以大部分是搬运的。

参考文献

二维计算几何基础 - OI Wiki (oi-wiki.org)

计算几何初步 - 洛谷专栏 (luogu.com)

判断两个线段是否相交 - balabala已被注册 - 博客园 (cnblogs.com)

计算几何三维入门 - 洛谷专栏 (luogu.com)

【笔记】二维凸包 - 洛谷专栏 (luogu.com)

平面计算几何基础

线段平行

设线段端点分别为$P_1,P_2$和$P_3,P_4$

则$$\overrightarrow{P_1 P_2} \times \overrightarrow{P_3 P_4} = 0$$

线段共线

设线段端点分别为$P_1,P_2$和$P_3,P_4$

则$$\overrightarrow{P_1 P_2} \times \overrightarrow{P_2 P_3} = 0$$ $$\overrightarrow{P_1 P_2} \times \overrightarrow{P_2 P_4} = 0$$

判断一个点在线段的哪边

假设有向线段为$\overrightarrow{AB}$ ，点为$C$，首先计算外积$\overrightarrow{AC} × \overrightarrow{AB}=( \overrightarrow{AC}_{x} ⋅\overrightarrow{AB}_{y}-\overrightarrow{AC}_{y} ⋅\overrightarrow{AB}_{x})\overrightarrow{k}$（因为有$\overrightarrow{AB}_{z}=0$，$\overrightarrow{AB}_{z}=0$）。根据右手螺旋定则，如果$\overrightarrow{k}$的系数为正数，说明点$C$在线段$AB$的右侧；如果为负数，说明点$C$在线段$AB$的左侧；如果为0，说明点$C$在线段$AB$所在的直线上。

//*************************************************************************  
// \brief: 计算两个向量的外积（叉乘）。可以根据结果的符号判断三个点的位置关系。  
// \Param: Point A 两个向量的公共起点。  
// \Param: Point B 第一个向量的终点。  
// \Param: Point C 第二个向量的终点。  
// \Returns: double 向量AC与向量AB的外积。
//如果结果为正数，表明点C在直线AB（直线方向为从A到B）的右侧；  
//如果结果为负数，表明点C在直线AB（直线方向为从A到B）的左侧；
//如果结果为0，表明点C在直线AB上。  
//*************************************************************************
double cross(Point A, Point B, Point C) { 
	double cross1 = (C.x - A.x) * (B.y - A.y);
    double cross2 = (C.y - A.y) * (B.x - A.x);
    return (cross1 - cross2);
}

线段相交

跨立实验

T1 = cross(A1, A2, B1);
T2 = cross(A1, A2, B2);
T3 = cross(B1, B2, A1);
T4 = cross(B1, B2, A2);

快速排斥实验

(1). 如果$(T1 * T2) > 0$ || $(T3 * T4) > 0$，说明一条线段的两个端点在另一条线段的同侧，这两条线段肯定不相交。

(2). 如果$T1 == 0$ && $T2 == 0$，说明两条线段共线，是否相交还需要进一步判断。这时可以通过判断两条线段张成的矩形是否相交来判断，而两个矩形是否相交可以通过快速排斥实验来判断。

//*************************************************************************
// \brief: 	快速排斥实验，判断两个线段张成的矩形区域是否相交。
// \Param: 	Point S1 第一条线段的起点。
// \Param: 	Point E1 第一条线段的终点。
// \Param: 	Point S2 第二条线段的起点。
// \Param: 	Point E2 第二条线段的终点。
// \Returns: 	bool 两个线段张成的矩形区域是否相交。
//具有对称性，即交换两条线段（参数S1与S2交换、E1与E2交换），结果不变。
//*************************************************************************  
bool rectsIntersect(Point S1, Point E1, Point S2, Point E2) {  
　　if ( min(S1.y, E1.y) <= max(S2.y, E2.y) && max(S1.y, E1.y) >= min(S2.y, E2.y) 
	&& min(S1.x, E1.x) <= max(S2.x, E2.x) && max(S1.x, E1.x) >= min(S2.x, E2.x)) {
            return true; } return false;
}

(3). 其他情况，两个线段一定相交。

代码实现：

bool segmentsIntersect(Point A1, Point A2, Point B1, Point B2) { 
    double T1 = cross(A1, A2, B1);
    double T2 = cross(A1, A2, B2); 
    double T3 = cross(B1, B2, A1);
    double T4 = cross(B1, B2, A2);　　
    if (((T1 * T2) > 0) || ((T3 * T4) > 0)) {    
    // 一条线段的两个端点在另一条线段的同侧，不相交。
    //（可能需要额外处理以防止乘法溢出，视具体情况而定。）
        return false;
    } else if(T1 == 0 && T2 == 0) {             
    // 两条线段共线，利用快速排斥实验进一步判断。此时必有 T3 == 0 && T4 == 0。
        return rectsIntersect(A1, A2, B1, B2);
    } else {                                    // 其它情况，两条线段相交。
        return true;  
 }
}

可以看到，这种方法不需要对线段的起终点重合（线段退化为一个点）做特殊判断，也不需要对线段平行（除了共线的情况）做特殊判断。纯几何方法，逻辑更简洁。

线段交点

参考文献：计算几何之两条线段的交点 - Huntto - 博客园 (cnblogs.com)

设线段端点分别为$P_1,P_2$和$P_3,P_4$，交点为$P_0$，$$P_0 = P_1+t * \overrightarrow{P_1 P_2} $$ $$P_0 = P_3+s * \overrightarrow{P_3 P_4} $$

将点坐标代入公式$$\begin{cases}x_1+(x_2-x_1)t=x_3+(x_4-x_3)s \\ \\ y_1+(y_2-y_1)t=y_3+(y_4-y_3)s \end{cases}$$

解方程得$$\begin{cases}t=\frac{(x_3-x_1)(y_4-y_3)-(y_3-y_1)(x_4-x_3)}{(x_2-x_1)(y_4-y_3)-(y_2-y_1)(x_4-x_3)} \\ \\ s=\frac{(x_3-x_1)(y_2-y_1)-(y_3-y_1)(x_2-x_1)}{(x_2-x_1)(y_4-y_3)-(y_2-y_1)(x_4-x_3)} \end{cases}$$

写作向量形式为 $$\begin{align*}t=\frac{\overrightarrow{P_3 P_1} \times \overrightarrow{P_4 P_3}}{\overrightarrow{P_2 P_1} \times \overrightarrow{P_4 P_3}} \\s=\frac{\overrightarrow{P_3 P_1} \times \overrightarrow{P_2 P_1}}{\overrightarrow{P_2 P_1} \times \overrightarrow{P_4 P_3}} \end{align*}$$

因此$P_0$的坐标为$$\begin{cases}x_0 = x_1+t * (x_2-x_1) \\y_0 = y_1+t * (y_2-y_1) \end{cases}$$或$$\begin{cases}x_0 = x_3+s * (x_4-x_3) \\y_0 = y_3+s * (y_4-y_3) \end{cases}$$

向量旋转

设$\overrightarrow{a} = (x,y)$，逆时针旋转$\alpha$角得到

$$\overrightarrow{b} =(x\cos\alpha −y\sin\alpha ,y\cos\alpha+x\sin\alpha )$$

可以用复数的乘法证明。

三角剖分求面积

把相邻每两个顶点与原点构成的向量的叉积的数值的一半依次累加起来，就能得到多边形的面积。

$$ S_{ABCDEF} = \frac{\overrightarrow{OA} × \overrightarrow{OB} + \overrightarrow{OB} × \overrightarrow{OC} +⋯+ \overrightarrow{OF} × \overrightarrow{OA}}{2} $$

凸包

首先把所有点以横坐标为第一关键字，纵坐标为第二关键字排序。

显然排序后最小的元素和最大的元素一定在凸包上。而且因为是凸多边形，我们如果从一个点出发逆时针走，轨迹总是 “左拐” 的，一旦出现右拐，就说明这一段不在凸包上。因此我们可以用一个单调栈来维护上下凸壳。

因为从左向右看，上下凸壳所旋转的方向不同，为了让单调栈起作用，我们首先升序枚举求出下凸壳，然后降序求出上凸壳。

求凸壳时，一旦发现即将进栈的点（$P$）和栈顶的两个点（$S1​,S2​$，其中 $S1​$ 为栈顶 ）行进的方向向右旋转，即叉积小于 00：$\overrightarrow{S_{2}​S_{1}}​​×\overrightarrow{S_{1}P}​<0$，则弹出栈顶，回到上一步，继续检测，直到 $\overrightarrow{S_{2}​S_{1}}​​×\overrightarrow{S_{1}P}​\ge 0$​或者栈内仅剩一个元素为止。

通常情况下不需要保留位于凸包边上的点，因此上面一段中 $\overrightarrow{S_{2}​S_{1}}​​×\overrightarrow{S_{1}P}​<0$ 这个条件中的 “$<$” 可以视情况改为 $\le$，同时后面一个条件应改为 $>$。

最后不要忘了把最小的元素与栈顶进行比较，以保证最后一段也是凸壳。

时间复杂度 $O(n\log_{}{n} )$。

sort(p+1,p+1+n); 
stk[++tp]=1; 
for(int i=2;i<=n;++i) { 
	while(tp>1&&(p[stk[tp]]-p[stk[tp-1]])*(p[i]-p[stk[tp]])<=0) 
		[stk[tp--]]=0; 
	used[i]=1; 
	stk[++tp]=i; 
} 
int qaq=tp; 
for(int i=n-1;i>0;--i) 
	if(!used[i]) { 
	while(tp>qaq&&(p[stk[tp]]-p[stk[tp-1]])*(p[i]-p[stk[tp]])<=0) 
		used[stk[tp--]]=0; 
	used[i]=1; 
	stk[++tp]=i; 
} 
for(int i=1;i<=tp;++i) 
	h[i]=p[stk[i]];

eps误差处理

$$a==0 : fabs(a)<eps$$ $$a<0 : a<−eps$$ $$a>0 : a>eps$$ $$a<=0 : a<eps$$ $$a>=0 : a>-eps$$

在一些卡精度的题（如UVA10173），允许的情况下建议使用long double

注意long double 对应 %Lf

旋转卡壳

旋转卡壳是用来求凸包直径的算法。

对踵点

如果过凸多边形上两点作一对平行线，使得整个多边形都在这两条线之间，那么这两个点被称为一对对踵点。

凸多边形的直径

即凸多边形上任意两个点之间距离的最大值。直径一定会在对踵点中产生。

算法思路

我们以凸包上的每条边为对象考虑，求出每一条边的对踵点。

边的对踵点示意图

点到边的距离可以用三角形面积*2除以底边长度来得到，而三角形面积可以通过叉积求解，又注意到对踵点在逆时针枚举边的过程中也逆时针移动，所以对于所有边的对踵点，我们使用双指针法求解。

难点在于编码细节极其多。

struct vec{
	ll x,y;
	friend vec operator -(const vec &a,const vec &b)
		{return vec{a.x-b.x,a.y-b.y};}
	friend ll operator *(const vec &a,const vec &b)
		{return (a.x*b.y-a.y*b.x);}
	ll dis_(){return x*x+y*y;}
	friend ll dis(const vec &a,const vec &b,const vec &x)
		{return abs((x-a)*(x-b));}
	friend bool operator <(const vec &a,const vec &b)
		{return (a.x<b.x)||(a.x==b.x&&a.y<b.y);}
}p[100005]; 
int n,st[100005],top=0;
char used[100005];
inline void convexhull(){	//求凸包 
	sort(p+1,p+1+n);
	for(int i=1;i<=n;i++){
		while(top>1&&(p[st[top]]-p[st[top-1]])*(p[i]-p[st[top]])<=0)
			used[st[top--]]=0;
		st[++top]=i;
		used[i]=1;
	}
	int qaq=top;
	for(int i=n-1;i>1;i--){
		if(!used[i]){
			while(top>qaq&&(p[st[top]]-p[st[top-1]])*(p[i]-p[st[top]])<=0)
				used[st[top--]]=0;
			st[++top]=i;
			used[i]=1;
		}	
	}
	while(top>qaq&&(p[st[top]]-p[st[top-1]])*(p[1]-p[st[top]])<=0)
		used[st[top--]]=0;
	st[++top]=1;
}
inline ll rotating(){
	if(top<4)		//注意点1会在开头和结尾各自入栈一次，所以栈的大小为凸包点数+1 
		return (p[st[1]]-p[st[2]]).dis_();	//如果凸包只有两个点直接输出两点距离 
	int t=3;
	ll mx=0;
	for(int i=1;i< top;i++){
	    while(dis(p[st[i]],p[st[i+1]],p[st[t]])<=dis(p[st[i]],p[st[i+1]],p[st[t%top+1]]))
			t=t%top+1;
	    mx=max(mx,max((p[st[t]]-p[st[i]]).dis_(),(p[st[t]]-p[st[i+1]]).dis_()));
	}
	return mx;
}

例题：Smallest Bounding Rectangle - UVA 10173

给出平面上的一堆点，找出一个能够覆盖所有点的面积最小的矩形，输出面积及四个顶点的坐标。

首先有结论：这样的矩形一定有一条边与凸包重合。

最小矩形覆盖

那么我们对于每一条边，只需要求出对踵点，最左点和最右点。

最左点和最右点的求法和对踵点类似，我们只需要改用点积来判断即可。

求出这三点后，我们先求出底边的垂线，就可以结合左右两点坐标求出水平长度，再结合边到对踵点的高度求出面积。

这题还是卡精度。

#define double long double
typedef long long ll;
const int MAXN=100001;
int n;
const double eps =1e-8;
struct vec{
	double x,y;
	friend vec operator +(vec x,vec y){return vec{x.x+y.x,x.y+y.y};}
	friend vec operator -(vec x,vec y){return vec{x.x-y.x,x.y-y.y};}
	friend double operator *(vec x,vec y){return (x.x*y.y-x.y*y.x);}
	friend double operator ^(vec x,vec y){return (x.x*y.x+x.y*y.y);}
	friend bool operator <(vec x,vec y)
		{return (x.x-y.x<eps)||(fabs(x.x-y.x)<eps&&x.y-y.y<eps);}
	double dis(){return sqrt(x*x+y*y);}
}p[MAXN],h[MAXN];
int st[MAXN],top;
char used[MAXN];
inline void convexhull(){
	top=0;
	memset(used,0,n+1);
	sort(p+1,p+1+n);
	st[++top]=1;
	for(int i=2;i<=n;i++){
		while(top>1&&(p[st[top]]-p[st[top-1]])*(p[i]-p[st[top]])<eps) 
			used[st[top--]]=0;
		st[++top]=i;
		used[i]=1;
	}
	for(int i=n-1;i>=1;i--){
		if(!used[i]){
			while(top>1&&(p[st[top]]-p[st[top-1]])*(p[i]-p[st[top]])<eps) 
				used[st[top--]]=0;
			st[++top]=i;
			used[i]=1;
		}
	}
	for(int i=1;i<=top;i++) h[i]=p[st[i]];
}
inline double sqr(vec a,vec b,vec x){return (a-x)*(b-x);}
inline double solve(){
	if(top<4) return 0;
	int t=3;
	while(sqr(h[1],h[2],h[t])-sqr(h[1],h[2],h[t%top+1])<eps) t=t%top+1;
	int t2=1,t3=t;
	double ans=1e18;
	for(int i=1;i<top;i++){
		while(sqr(h[i],h[i+1],h[t])-sqr(h[i],h[i+1],h[t%top+1])<eps) 
			t=t%top+1;
		while(((h[i+1]-h[i])^(h[t2%top+1]-h[t2]))>-eps) 
			t2=t2%top+1; //记得都要%top再+1 
		while(((h[i+1]-h[i])^(h[t3%top+1]-h[t3]))<eps) 
			t3=t3%top+1;
		vec vert={-(h[i+1].y-h[i].y),h[i+1].x-h[i].x};
		double ar=(sqr(h[i],h[i+1],h[t])/(h[i]-h[i+1]).dis())
					*(sqr(h[t2],h[t2]+vert,h[t3])/(vert).dis());
		ans=min(ans,ar);
	}
	return ans;
}
int main(){
	ios::sync_with_stdio(false);
	cin.tie(0);
	while(cin>>n){
		if(!n) return 0; 
		for(int i=1;i<=n;i++) cin>>p[i].x>>p[i].y;
		convexhull();
		printf("%.4Lf\n",solve());
	}
	return 0;
}

自适应辛普森法

这玩意是计算几何吗

好像还真是

Simpson公式：用抛物线来拟合原函数。$$\int\limits_{a}^{b}f(x)dx=\frac{1}{6}(b-a)[f(a)+f(b)+4f(\frac{a+b}{2})]$$

通过不断将区间二分直到满足需要的精度。

inline double simpson(double l,double r) {
	return (r-l)*(f(l)+f(r)+4*f((l+r)/2))/6; //f(x)是原函数
}
double asr(double l,double r,double eps,double ans) {
	double mid=(l+r)/2;
	double l_=simpson(l,mid),r_=simpson(mid,r);
	if(fabs(l_+r_-ans)<=15*eps) return l_+r_+(l_+r_-ans)/15;
	return asr(l,mid,eps/2,l_)+asr(mid,r,eps/2,r_);
}
inline double asr(double l,double r,double eps) {
	return asr(l,r,eps,simpson(l,r));
}

为什么要乘15参考：自适应辛普森法详解 - 知乎 (zhihu.com)

题外话：原来卡西欧是这么算定积分的

例题：BZOJ-2178 圆的面积并

用自适应辛普森法求解，$f(x_{0})$即直线$x=x_{0}$与所有圆的相交的$y$坐标区间的长度。

trick：注意这道题看似对精度要求不高，实则测试点13会卡精度，所以为了避免被特意构造的数据卡掉，我们将坐标系旋转一个随机角度。这是一个常用技巧，在旋转卡壳中也可以使用。

#define double long double
int n;
double x[1005],y[1005],r[1005];
typedef pair<double,double> P;
P q[1005];
inline double f(double v){
	int tot=0;
	for(int i=1;i<=n;i++){
		double dx=x[i]-v;
		if(r[i]*r[i]-dx*dx>1e-8){
			double dy=sqrt(r[i]*r[i]-dx*dx);
			q[++tot]=P(y[i]-dy,y[i]+dy);
		}
	}
	if(!tot) return 0;
	sort(q+1,q+1+tot);
	//以下将x=v穿过的所有y区间进行合并 
	double ans=0,st=q[1].first,ed=q[1].second;
	for(int i=2;i<=tot;i++){
		if(q[i].first<=ed)
			ed=max(ed,q[i].second);
		else{
			ans+=ed-st;
			st=q[i].first;
			ed=q[i].second;
		}
	}
	return ans+ed-st;
}
inline double simpson(double l,double r){
	return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6;
}
inline double asr(double l,double r,double eps,double ans){
	double mid=(l+r)/2;
	double ansl=simpson(l,mid),ansr=simpson(mid,r);
	if(fabs(ansl+ansr-ans)<15*eps) return ansl+ansr+(ansl+ansr-ans)/15;
	return asr(l,mid,eps/2,ansl)+asr(mid,r,eps/2,ansr);
}
void _rotate(double &x, double &y,double angle)
{
	double x_=x;
	double y_=y;
	x=x_ * cos(angle) + y_ * sin(angle);y=-x_ * sin(angle) + y_ * cos(angle);
}
int main(){
	ios::sync_with_stdio(0);cin.tie(0);
	cin>>n>>x[1]>>y[1]>>r[1];
	int rot=1;	//坐标系旋转一个随机角度 
	_rotate(x[1],y[1],rot);
	double lf=x[1]-r[1],rf=x[1]+r[1];
	for(int i=2;i<=n;i++){
		cin>>x[i]>>y[i]>>r[i];
		_rotate(x[i],y[i],rot);
		lf=min(x[i]-r[i],lf);
		rf=max(x[i]+r[i],rf);
	}
	printf("%.3Lf",asr(lf-10,rf+10,1e-5,simpson(lf-10,rf+10)));
	return 0;
}

由于自适应辛普森法是一种近似方法，因此本题还有另一种解法是使用格林公式：[BZOJ2178]圆的面积并(格林公式) - coder66 - 博客园 (cnblogs.com)

UPD：不要在XCPC中使用自适应辛普森法，因为这种解法极有可能被出题人卡掉。参考2022 CCPC Guilin Site F的题解。

2022 CCPC Guilin Site-tutorial