计算几何
基础
(1) $pi = acos(-1)$
(2) 余弦定理 $c^2 = a^2 + b^2 - 2abcos(t)$
浮点数的比较
向量
定义的eps太小了在判断直线时会错误或者tle，太大了误差并比较大，，所有需要适量
精度不够使用long double
trick
注意舍入方式(0.5 的舍入方向);
- 防止输出-0，判断是0，输出0
- 几何题注意多测试不对称数据
- 数几何注意 xmult 和 dmult 是否会出界; 符点几何注意 eps 的使用.
- 避免使用斜率;注意除数是否会为 0.
- 公式一定要化简后再代入
- 判断同一个 2*PI 域内两角度差应该是 abs(a1-a2)< beta||abs(a1-a2)>pi+pi-beta;相等应该是 abs(a1-a2)< eps||abs(a1-a2)>pi+pi-eps; //注意周期性
- 需要的话尽量使用 atan2,注意:atan2(0,0)=0
- 误差限缺省使用 1e-8!
- 可以对所有点旋转一个随机的角度，防止毒瘤出题人卡
各种基础函数
int sign(double x){ // 判断一个浮点数的符号
if(fabs(x) < eps) return 0 ;
if(x < 0) return - 1;
return 1 ;
}
int dcmp(double x,double y){ //比较两个浮点数
if(fabs(x - y) < eps) return 0 ;
if(x < y) return - 1;
return 1;
}
pdd operator +(pdd a,pdd b){ //点的重载+
return {a.x+b.x,a.y+b.y} ;
}
pdd operator-(pdd a,pdd b){ //点的重载-
return {a.x-b.x,a.y-b.y} ;
}
pdd operator*(pdd a,double b){ //向量乘实数
return {a.x * b , a.y * b} ;
}
pdd operator/(pdd a,double b){ //向量除以实数
return {a.x / b , a.y / b} ;
}
double operator&(pdd a,pdd b){
return a.x * b.x + a.y * b.y ; //两个向量的点积
}
double operator*(pdd a,pdd b){ //向量的重载叉乘，可以用来求面积，判断点在向量的那边
return a.x * b.y - a.y * b.x ;
}
double area(pdd a,pdd b,pdd c){ //求得面积，从ab转向ac为正方向
return (b - a) * (c - a) ;
}
double get_len(pdd a){ //向量长度，模长
return sqrt(a & a) ;
}
double project(pdd a,pdd b,pdd c){
return ((b - a) & (c - a)) / get_len(b - a) ; //投影长度
}
pdd norm(pdd a){ //单位向量
return a / (double)get_len(a) ;
}
double rand_eps(){
return ((double)rand() / RAND_MAX - 0.5) * eps ; //一个轻微扰动
}
void shake(){
x += rand_eps() , y += rand_eps(), z += rand_eps() ; //一个轻微扰动
}
double dist(pdd a,pdd b){ //求得两个点的距离
double dx = a.x - b.x ;
double dy = a.y - b.y ;
return sqrt(dx * dx + dy * dy) ;
}
pdd rotate(pdd a,double b){ //顺时针旋转b个弧度
return {a.x * cos(b) + a.y * sin(b) , - a.x * sin(b) + a.y * cos(b) } ; //矩阵变换
}
pdd get_line_intersection(pdd a,pdd v,pdd b,pdd w){ //获得两条直线的交点，直线使用点向式来表达
pdd u = a - b ;
double t = w * u / (v * w) ;
return a + v * t ;
}
bool on_segment(pdd a,pdd b,pdd c){ //判断一个点是否在一个线段上
return !sign(area(a,b,c)) && sign(project(a,b,c)) <= 0 ; //当面积是0，代表在直线上，，当投影是负数，代表c在ab中间
}
double get_angle(const node &a){ //使用atan2函数来得到直线与x轴正方向的夹角
return atan2(a.ed.y - a.st.y ,a.ed.x - a.st.x) ;
}
凸包
凸包
数学定义：设S为欧几里得空间$R^n$的任意子集。包含S的最小凸集称为S的凸包，记作conv(S)。
二维凸包
平面凸包
Graham扫描法
将所有的点排序，按照x来排序
然后从第一个点开始，然后每次加入一个点，新加入的点与已有凸包中最后一个点形成直线的斜率与原有凸包中最后一条直线斜率的比较
然后看是否弹出已有凸包中的点，
使用一个栈来记录所有凸包的点（凸包围成的凸多边形是周长最小的可以包含所有点的多边形，但不是面积最小的包含所有点的多边形）
时间复杂度是被排序限制的，所以是$O(nlongn)$
double gradham(){
sort(q,q+n) ; //对所有的点进行排序
for(int i = 0 ; i < n ; i++){ //遍历n个点
while(top >= 2 && sign(area(q[stk[top-2]],q[stk[top-1]],q[i])) <= 0){//如果新形成的直线斜率比最后一条直线小就删除最后一个点
if(sign(area(q[stk[top-2]],q[stk[top-1]],q[i])) < 0) used[stk[--top]] = 0 ;//如果斜率相等，那么就不恢复标记
else top -- ;
}
used[i] = 1; //标记
stk[top++] = i ; //加入凸包栈
}
used[0] = 0 ; //恢复初始点
for(int i = n - 1; i >= 0 ; i--){
if(used[i]) continue ; //如果找下凸包的时候已经用过该点，之间continue
while(top >= 2 && sign(area(q[stk[top-2]],q[stk[top-1]],q[i])) <= 0) top -- ;
stk[top++] = i ;
}
double res = 0 ;
for(int i = 1 ; i < top ; i++){ //遍历所有的点，然后找到周长
res += dist(q[stk[i-1]],q[stk[i]]) ;
}
return res ;
}
三维凸包
多面体欧拉定理
顶点数-棱长数+表面数=2
增量法
增量法
基本思想是把点依次加到凸包中，跟最小覆盖园类似，，，都是先已经有了一个凸包，然后考虑再加入一个点，新加入的点对原先凸包的影响
如果已经在凸包内部，那么就不用处理，如果在外部的话，，考虑投影，新的点可以照射到所有的面将会被更新，
旋转卡壳
在凸包的基础上，找到直线来卡住，以便于找到卡住凸包的一些东西
例如找到最大的两个点的距离，
那么就是找到使用两条直线旋转卡壳，当两条直线旋转一周之后，就会找到答案
那么旋转卡壳就是一种双指针，，在凸包上的双指针
当固定一条直线，，当这条直线旋转的时候，距离这条直线最远的点是单调沿着凸包前进
由于找凸包需要排序，所以时间复杂度是$O(nlogn)$
int rotating_calipers(){ //旋转卡壳
if(top <= 2 ) return dist(q[0],q[n-1]) ; //特例
int ans = 0 ;
for(int i = 1,j = 3 ; i < top ; i++){ //从开始遍历一条边
pii a = q[stk[i]],b = q[stk[i+1]] ;
while(area(a,b,q[stk[j]]) < area(a,b,q[stk[j+1]])){ //找到距离这条边最远的点，单调更新
j ++ ;
if(j == top) j = 1 ;
}
ans = max(ans,dist(q[stk[j]],a)) ; //更新答案
ans = max(ans,dist(q[stk[j]],b)) ;
}
return ans ;
}
半平面交
算法是可以求出多个多边形的交集
前置知识
求一个多边形的面积，可以进行三角剖分，n个点分成n-2个三角性，求出面积，，也可以在多边形外随便找一点，连接该点到多边形的一条边，然后多边形边如果符合从右向左记位面积，否则记位负面积，最后相加(利用的是容斥原理)
然后一个多边形其实就是多个半平面的交集。
记一个条直线的向量的左侧可以记位这个半平面，一个多边形可以转化为半平面的交集，那么多个平面也同样可以，，那么使用一个双端队列来维护半平面的边界，对所有直线排序。
double half_plane_intersection(){ //获得半平面交
sort(line,line+cnt,cmp) ; //按照直线角度排序
int hh = 0 , tt = - 1;
for(int i = 0 ; i < cnt ; i++){
if(i && !dcmp(get_angle(line[i]),get_angle(line[i-1]))) continue ;
while(hh + 1 <= tt && is_right(line[i],line[q[tt-1]],line[q[tt]])) tt -- ; //更新队列头
while(hh + 1 <= tt && is_right(line[i],line[q[hh+1]],line[q[hh]])) hh ++ ; //更新队列尾
q[++tt] = i ;
}
while(hh + 1 <= tt && is_right(line[q[hh]],line[q[tt-1]],line[q[tt]])) tt -- ; //使用队头来更新队尾
while(hh + 1 <= tt && is_right(line[q[tt]],line[q[hh+1]],line[q[hh]])) hh ++ ; //使用队尾来更新队头
q[++tt] = q[hh] ;
int idx = 0 ;
for(int i = hh ; i < tt ; i++){
ans[idx++] = get_line_intersection(line[q[i]],line[q[i+1]]) ; //找到边界上的点
}
double res = 0 ;
for(int i = 1 ; i + 1 < idx; i++){
res += area(ans[0],ans[i],ans[i+1]) ; //使用三角剖分算出面积
}
return res ;
}
最小圆覆盖
三角剖分
凸多边形不好考虑情况，所以可以剖分成多个三角形分别进行考虑
扫描线
把一个复杂图像通过分割，分割乘多个区间，是得每个区间面积容易计算
使用扫描线这种直线进行分割
比如是求多个矩形的面积并，，那么就可以把所有出现过的端点的横坐标当成线来分割
在每个中进行区间合并
自适应辛普森积分
首先辛普森积分就是使用二次函数进行逼近图像的积分
由一段的两个端点和中点，那么二次函数逼近的积分是 (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6 ;
然后说自适应，自适应就是使用二分区间递归的方式来自动适应题目的精度，每次二分左区间和右区间，进行精度比较
double simpon(double l,double r){ //求出l到r使用二次函数逼近面积
double mid = (l + r) / 2 ;
return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6 ; //公式
}
double simpon(double l,double r,double s){
double mid = (l + r) / 2 ;
double left = simpon(l,mid) ,right = simpon(mid,r) ; //分成左区间和右区间
if(fabs(left + right - s) < eps) return left + right ; //查看是否精度足够
return simpon(l,mid,left) + simpon(mid,r,right) ; //递归
}
