计算几何
计算几何是一个偏冷门的板块（基本只有 SCOI 考/gg），但是基本的了解还是要有的。
平面几何基础
向量与点
由于点和向量的运算和性质等的相似性，我们将其一并定义。
定义向量与点
注意到二者有类似的形式。事实上，我们经常将二者视为同一种东西，因此经常直接将点代做向量进行向量运算。例如点 \(A\)，我们可以直接将其视作 \(\vec A\) 带入运算。
向量的基本运算分为向量与实数运算和向量间运算。
其中，向量间运算主要分为叉积和点积两种。
叉积的几何意义为 \(\vec a\) 与 \(\vec b\) 之间夹成的三角形的面积的两倍，点积的几何意义是一个向量在另一个向量上的投影长度。
若 \(\vec b\) 在 \(\vec a\) 的逆时针方向，叉积为正，反之为负，共线为 0。这一个性质相对关键，比较常用。
而如果两个向量的夹角为钝角，则点积为负，反之为正，恰好垂直则为 0。
在旋转之前，我们先介绍极角。极角是指 \(x\) 轴正半轴与向量逆时针方向夹成的角。这个还是比较好理解的。
向量的旋转可以通过和角公式推导。设旋转角度为 \(\theta\)，\(\vec a\) 的长度为 \(r\)，极角为 \(\varphi\)，有
这就是向量的旋转公式，还比较好推。
直线
由于可能存在竖直的线，因此直线与线段一般都不使用斜率来存，而是用线上的两个点来确定线的形态。
同时，由于计算几何的特殊性，我们是要求直线与线段是有方向性的，其方向的规定与向量类似。
定义一条直线被点 \(s,t\) 确定，方向是 \(s\to t\)。
点 \(B\) 到直线的距离 \(L\) 可以通过常见的面积法得到。由于叉积的几何意义，我们直接有
其中 \(dis_{u,v}\) 表示两点间的距离。
直线交点由相似法得到。假设第一条直线的起终点分别为 \(A,B\)，第二条的为 \(C,D\)，设原点为 \(O'\)，交点为 \(O\)，则
其中分子上是 \(\triangle ADC\) 面积的两倍，分母是四边形 \(ACBD\) 的面积。这个分数相当于是求出了 \(\frac{AO'}{AB}\) 的值。
一个点 \(A\) 向一条直线作垂线的垂足，实际上就是投影向量，可以通过点积快速得到。
设垂足为 \(O\)，有
线段
基本与直线类似。需要注意的是对于线段而言，一个点到线段的距离与直线是不同的。
如果一个点不在线段“对应”的平面内，那么距离就不是直线，而是到最近的端点的距离。所谓“对应”就是指是否这个点与两个端点的两条连线与线段本身的夹角都小于等于九十度。这个东西可以用点积的性质快速判。
多边形
一般我们存储多边形的方式都是通过点集。存储的相邻两个点之间存在连边，最后一个点与第一个点连边，所有边组成一个封闭图形。
简洁起见，我们定义第 \(i\) 个点在多边形上的下一个点为 \(nxt (i)\)，这主要是为了方便最后一个点与第一个点间连边的表示。
在代码实现中，由于使用了 vector 存储，因此点的编号是从 0 开始的。而在算式中，我们为了方便从 1 开始编号。一般将多边形的点积设为 \(P\)，多边形上的第 \(i\) 个点表示为 \(P_i\)，一共有 \(n\) 个点。
我们在这里约定一般的多边形都是点逆时针旋转的方向存储的。
多边形的周长是好算的，直接将两个相邻点的距离加起来即可。
多边形的面积则是一个相对困难的点。由于没有保证多边形是一个凸多边形，因此我们不能采用类似将多边形分成若干个三角形的方式快速计算。我们希望能有一个简洁的公式能够快速计算。事实上，我们有
确实比较简洁，但证明并不显然。对于凸多边形，这个公式还是相对显然的。讨论原点在多边形内部与外部两种情况，可以通过画图比较明显的得出来。
而对于凹多边形，情况就变得复杂了起来。画个图可以发现，由于叉积自带符号的特性，这个式子容斥的意味变得比较深。
其中红色的是负，蓝色的是正。在感性上可能能理解，在理性上感觉不严谨。
叉积性质
\(1\over2\) 叉积即带符号的三角形面积，多边形每条边对应的面积为正/负，称之为 \(+1\) 边 / \(-1\) 边。
对于边 \(AB\)（\(A=A_i,B=A_{(i+1)\bmod n}\)），引射线 \(OA\)，看 \(AB\) 相对于 \(A\) 之后的射线是顺时针还是逆时针转（转动 \(<180^\circ\)），顺时针就是 \(-1\)，逆时针就是 \(+1\)。
射线性质
若 \(O\) 在多边形上，就施以轻微扰动，让它去多边形内/外。以 \(O\) 为起点，引一条射线 \(r\)，要求 \(r\) 不经过多边形的顶点（自然也不与多边形的边重合），称这样的射线为“可用射线”，此时 \(r\) 交的边都“跨过”\(r\)。
由于“叉积性质”，跨过 \(r\) 的边必定是 \(+1\)/\(-1\) 交替出现。
若点 \(O\) 在多边形外，由反证法可得若有边跨过 \(r\)，则第一条跨过 \(r\)（交点离 \(O\) 最近）的边必定是 \(-1\) 边（反证法：保证是第一条，还要形成闭合回路，那么必定会绕成顺时针，与假设矛盾）。同理，\(O\) 在多边形内时，第一条跨过 \(r\) 的边必定是 \(+1\) 边。
贡献系数
要求对平面内每个点 \(B\)，若其在多边形内，则贡献系数为 \(1\)，否则（多边形上/多边形外）贡献系数为 \(0\)。首先忽略多边形上的点、不在可用射线上的点。
从 \(O\) 向 \(B\) 引射线，称射线在 \(B\) 之后的部分为 \(r'\)。跨过 \(r'\) 的边的 \(\pm1\) 之和即该点的贡献系数。
\(r'\) 可看作“\(B\) 逃逸的过程”，每有一条边跨过则必定会转换“是否在多边形内”的状态，否则不转换。因为最终必定在多边形外，所以若 \(B\) 在多边形外则恰有偶数条边跨过 \(r'\)，否则恰有奇数条边跨过 \(r'\)。
于是多边形外的点 \(B\) 贡献系数为 \(0\)。又由于“射线性质”，无论 \(O\) 在何处，多边形内的点 \(B\) 对应的 \(r'\) 都会首先被 \(+1\) 边跨过，故贡献系数是 \(+1\)。
由于我对容斥的理解很差，因此这里给出一个同学的证明qwq。
凸包
求凸包
最先就是求点集的凸包。个人习惯是用 Andrew 算法。大体思路是将凸包分为两部分，分为上凸壳与下凸壳，然后拼起来。这个还是比较符合直觉的。
以求上凸包为例。先直接将点按横坐标排序。然后直接从左到右扫，如果新的点与栈顶的上一个点组成的直线在栈顶的前两个点形成的向量的逆时针方向，那栈顶就可以弹掉了。
下凸包同理，直接反着做一遍即可。
半平面交
然后是重头戏，求半平面交。半平面交是指给定一些平面然后要我们求这些平面的交。实际操作的时候一般是通过给出直线来给定半平面然后求这些半平面的交。
半平面是指一条直线及其某一边的平面。这个时候之前提到的“直线的方向性”的作用就体现于此。具体而言，一般钦定直线左边的平面被称作半平面。更准确地说，是直线的“逆时针方向”的平面。求半平面交的过程之所以比较繁复，就体现在定义的隐性，容易把自己转晕。
然后我们要引入一个叫极角排序的东西。上面在向量的内容中已经粗略介绍了极角的概念。极角是指向量与 \(x\) 正半轴逆时针方向所成的夹角。我们可以通过 c++ 内置的 atan2(y,x) 函数（注意这个函数是纵坐标在前）轻松求出极角的大小。不过需要注意的是，如果极角大小超过了 \(\pi\)，那么这个函数的值会通过负数来表示，也就是变成了从 \(x\) 轴正半轴逆时针旋转一个大小为负的角度。这个东西显然也是符合定义的，也比较符合直觉，但是这样会对之后的极角排序带来一定的影响。
极角排序，顾名思义，以极角的大小作为关键字来为向量排序。这个排序的效果是令所有的线与 \(x\) 轴正半轴的角度依次增大。（注意直线的方向性）求半平面交的过程其实比较朴素，就是将直线一条一条地放在平面上，然后围成一个图形。极角排序就可以规整这个“放直线”的顺序，像我们画出凸壳一样逆时针旋转着来放这些直线。
也就是说，我们实际上只需要写一个 sort，然后重定义一下直线排序的关键字，然后里面就是判断直线的极角，小的在前大的在后。
但我们还要考虑一下直线平行（这里的平行指代极角相等而不是斜率相等，还是注意方向性）的情况。我们发现由于是求半平面交，因此如果两条直线平行，那么更右边的那条直线就没用了。因此我们让更右边的直线排在后面，排完序后去掉即可。至于如何判断哪条在右边，可以先不管，看代码的时候依照代码画一下图即可。粗略地讲就是用叉积判一下逆时针还是顺时针。
然后就是“一条一条放直线”的过程。我们在这个过程中实际上干两件事：
- 算一下当前的直线与前一条直线的交点。（由于我们极角排序了，因此最后形成的图形的拐点一定是由我们放的相邻两条直线交成的）
- 删去一定没有贡献的直线。（这里说的没有贡献是指其上的任意一个点都一定不在最后的半平面交内部）
一条线什么时候会完全没用呢？
如图。由于极角排序，因此我们插入直线的顺序是黑到紫到蓝。显然插入了蓝色线以后紫色的线就一定没有任何作用了，删掉了也是一样的。我们发现，还是由于极角排序规定了插入直线的方向性，因此只需要很方便的判定上两条直线的交点是否在当前直线的右边，也就是顺时针方向即可。如果是，那么就可以直接将上一条线去掉。例如图中黑色紫色的交点在蓝色线右边，那紫色线就可以直接去掉了。
还有一种情况是我们已经画了大半了，这个时候插入的线碰到了最前面插入的直线。我们可以用类似于上面的方法判定前面的是否完全没用。可以之后去看代码，不再赘述。
上述操作在实现的时候就是用双端队列即可。
最后的问题是，我们一直在用后插入的线去判断之前插入的线是否没用。最后我们要做的事情就是用之前插入的线去判断后插入的线是否没用。
然后就可以愉快的将所有的拐点合成一个多边形。顺带一提，如果半平面交是一个封闭图形，那么一定是一个凸壳。
旋转卡壳
旋转卡壳求凸包直径的算法，也就是求凸包内最远的两点。
这个算法感性理解很简单，就是对于每一个点“最远”的线段，然后全部取 \(\max\) 即可。这个比较有单调性，随着选的点在转，线段一起转即可。
但是严谨证明实际上需要很多条件，我其实也扯不清楚。这里推荐一个博客：https://www.cnblogs.com/Sinktank/p/18308988 ，说的很清楚，图也很棒。
圆
三点确定一圆
可以用求中垂线的一次函数然后求交点的方法。但是显然误差巨大，同时没有斜率的问题要单独分讨（竖直的线）。
于是我们要抛弃初中方法，拥抱高中数学，掏出圆的方程，然后解出来很大一坨，看起来复杂度是常数，实际慢的离谱。于是我们可以先把算式中的一些重复东西先预处理，重复的浮点数运算是很慢的。
最小圆覆盖
实际上是个超级大暴力，然后分析复杂度期望是线性的。这种看起来很劣的东西是用类似随机化的东西保证的期望复杂度。
其思想就是：
- 以每个点为起点开始枚举，画个以其为圆心半径为 0 的圆。
- 如果某个点不在画的这个圆内，就把圆扩大，直到有三个点已经确定了圆的形态。
显然这种方法的正确性是对的，因为一旦没有了圆其就再找一个更大的圆。如果运气差其将枚举每一种可能的圆。
于是其用什么保证复杂度呢？我们将点集的顺序随机打乱。可以证明复杂度期望会收敛成为线性的。但反正我个人不是很喜欢这种很不符合直觉的算法。可以去看一下代码感受其有多暴力。
code
几乎所有操作都有注释，可以放心食用。同时定义了一些便于书写的函数，如 get(s,t) 就是求 \(s\to t\) 的向量。
同时特别提醒一下叉积是用“^”符号替代。
#include<bits/stdc++.h>
using namespace std;
#define ld double
const ld eps=1e-6,pi=acos(-1.0),inf=2e9+7;
mt19937 rnd(time(0));
int dcmp(ld x){return x<-eps?-1:(x>eps);}
ld Abs(ld x){return dcmp(x)*x;}
namespace bs{//基础的线段向量直线
struct pt{//点或向量
ld x,y;
void rd(){cin>>x>>y;}
void put(){cout<<x<<' '<<y<<'\n';}
};
bool operator < (const pt a,const pt b){return (a.x!=b.x)?a.x<b.x:a.y<b.y;} //以 x 为第一关键字排序，用于求凸包
pt operator * (const pt a,const ld b){return (pt){a.x*b,a.y*b};}//向量的常数运算
pt operator / (const pt a,const ld b){return (pt){a.x/b,a.y/b};}
pt operator + (const pt a,const pt b){return (pt){a.x+b.x,a.y+b.y};}//向量间运算
pt operator - (const pt a,const pt b){return (pt){a.x-b.x,a.y-b.y};}
ld operator * (const pt a,const pt b){return a.x*b.x+a.y*b.y;}//点积
ld operator ^ (const pt a,const pt b){return a.x*b.y-a.y*b.x;}//叉积
bool operator == (const pt a,const pt b){return (dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0);}//相等
ld len(pt a){return sqrt(a.x*a.x+a.y*a.y);}//模长
ld ang(pt a){return atan2(a.y,a.x);} //与 x 轴的夹角
ld dis(pt a,pt b){return len(a-b);} //两点间距离
int side(pt a,pt b){return dcmp(b^a);} //判断向量 a 在 b 的哪边，1 逆，-1 顺，0 共线
pt get(pt a,pt b){return b-a;} //得到两点间的向量，b 为尾
pt rotate(pt a,ld theta){ //向量转动一个特定角度
ld c=cos(theta),s=sin(theta);
return (pt){a.x*c-a.y*s,a.x*s+a.y*c};
}
struct seg{//线段
pt s,t;
};
seg get_seg(pt a,pt b){return (seg){a,b};}
ld len(seg a){return dis(a.s,a.t);} //线段长度
ld dis(seg a,pt b){//点到线段距离
pt s=a.s,t=a.t;
if(dcmp(get(s,b)*get(s,t))<0||dcmp(get(t,b)*get(t,s))<0)return min(dis(s,b),dis(t,b));
return (get(s,b)^get(t,b))/len(a);
}
pt close_pt(seg a,pt b){//线段上最近的点
pt s=a.s,t=a.t;
if(dcmp(get(s,b)*get(s,t))<0||dcmp(get(t,b)*get(t,s))<0)return dis(s,b)<dis(t,b)?s:t;
pt v=get(a.s,a.t);
return (v*(((b-a.s)*v)/(v*v)))+a.s;
}
struct line{//直线
pt s,t; //半平面统一在逆时针方向（左边）
};
line get_line(pt a,pt b){return (line){a,b};}
int side(line a,pt b){return side(b-a.s,a.t-a.s);} //判断点 b 在直线 a 的哪边 注意直线是有方向的
pt cross_line(line a,line b){ //直线交点
pt A=a.s,B=a.t,C=b.s,D=b.t;
return A+(get(A,B)*((get(A,C)^get(A,D))/(get(A,B)^get(C,D))));
}
pt ver_line(line a,pt b){ //一个点向直线做垂线的垂足
pt v=get(a.s,a.t); //直线的方向向量
return (v*(((b-a.s)*v)/(v*v)))+a.s;
}
pt sym_line(line a,pt b){ //一个点关于直线的对称点
return (ver_line(a,b)*2)-b;
}
bool operator < (const line a,const line b){ //极角排序的关键字 注意到这里的排序并非严格的按照极角从小到大的顺序排列，只是相对位置不变（因为极角为负的会在前）
ld a1=ang(get(a.s,a.t)),a2=ang(get(b.s,b.t));//注意在 sort 的关键字中不能出现等号
if(dcmp(a2-a1)!=0)return dcmp(a2-a1)>0; //按极角本身大小排序
return dcmp(get(a.s,b.t)^get(a.s,a.t))>0; //如果 a 在 b 的逆时针方向，那么 b 就没有用就排在前面然后在去重的时候去掉
}
}
using namespace bs;
namespace poly{//多边形
struct polygon{
vector <pt> pts;//原则上逆时针排序
void rd(int n){for(int i=1;i<=n;i++){pt tmp;tmp.rd();pts.push_back(tmp);}}
void put(){cout<<"poly_put\n";for(pt tmp:pts)tmp.put();}
pt operator [](int x){return pts[x];}
int nxt(int x){return x<pts.size()-1?x+1:0;}
ld C(){//多边形周长
ld ans=0;
for(int i=0;i<pts.size();i++)ans+=dis(pts[i],pts[nxt(i)]);
return ans;
}
ld S(){//多边形面积
ld ans=0;
for(int i=0;i<pts.size();i++)ans+=pts[i]^pts[nxt(i)];
return Abs(ans)/2;
}
};
}
using namespace poly;
namespace tb{//凸包
void andrew(polygon &p){ //求点集凸包
sort(p.pts.begin(),p.pts.end());
vector <pt> q(p.pts.size()<<1);
int top=-1,siz=p.pts.size();q[++top]=p[0];
for(int i=1;i<siz;i++){
while(top&&side(get(q[top-1],q[top]),get(q[top-1],p[i]))>=0)top--;
q[++top]=p[i];
}
int now=top;
for(int i=siz-2;i>=0;i--){
while(top>now&&side(get(q[top-1],q[top]),get(q[top-1],p[i]))>=0)top--;
q[++top]=p[i];
}
if(siz>1)top--;
p.pts.clear();for(int i=0;i<=top;i++)p.pts.push_back(q[i]);
}
ld kk(polygon p){ //旋转卡壳求凸包直径
ld ans=0;int beg=0,siz=p.pts.size(),tmp=0;
for(int i=0;i<siz;i++){ //找到对于 0 最远的线段
if(dis(get_seg(p[i],p[p.nxt(i)]),p[0])>tmp)
beg=i,tmp=dis(get_seg(p[i],p[p.nxt(i)]),p[0]);
}
for(int i=0;i<siz;i++){
while(dis(get_seg(p[beg],p[p.nxt(beg)]),p[i])<dis(get_seg(p[p.nxt(beg)],p[p.nxt(p.nxt(beg))]),p[i]))beg=p.nxt(beg);
ans=max({ans,dis(p[beg],p[i]),dis(p[p.nxt(beg)],p[i])});
}
return ans;
}
polygon si(vector <line> Q){ //求半平面交
int siz=Q.size(),len=0;sort(Q.begin(),Q.end());//极角排序
vector <line> q(siz+1),lst(siz+1);vector <pt> p(siz+1);//p 表示交点坐标
for(int i=0;i<siz;i++){ //去重
if(i&&dcmp(ang(get(Q[i].s,Q[i].t))-ang(get(Q[i-1].s,Q[i-1].t)))==0)continue;
q[++len]=Q[i];
}
int l=1,r=0; //双端队列
for(int i=1;i<=len;i++){
while(l<r&&dcmp(get(p[r],q[i].t)^get(p[r],q[i].s))>0)--r;
while(l<r&&dcmp(get(p[l+1],q[i].t)^get(p[l+1],q[i].s))>0)++l;
lst[++r]=q[i];if(r>l)p[r]=cross_line(lst[r],lst[r-1]);
}
while(l<r&&dcmp(get(p[r],lst[l].t)^get(p[r],lst[l].s))>0)--r; //由于前面只有后插入的来弹出前插入的，于是这里要用之前插入的去弹出后插入的
// while(l<r&&dcmp(get(p[l+1],lst[r].t)^get(p[l+1],lst[r].s))>0)++l;
p[r+1]=cross_line(lst[l],lst[r]),++r; //不知道为什么用 r++ 会错，感觉判定有点奇怪（实际上是因为加的时候仍然是 p_r 而不是 p_r+1 被赋值）
polygon P;for(int i=l+1;i<=r;i++)P.pts.push_back(p[i]);
return P;
}
}
using namespace tb;
namespace circle{//圆
struct cir{
pt o;ld r;
bool in(pt a){return dcmp(r-dis(a,o))>=0;}
};
cir get_cir(pt o,ld r){return (cir){o,r};} //直接得到圆
cir get_cir(pt u,pt v,pt w){ //用圆的解析式三点确定一圆，精度相对高
ld A,B,C,D;
ld tmp1=v.x-w.x,tmp2=v.x*w.y,tmp3=w.x*v.y,tmp4=(u.x*u.x+u.y*u.y),tmp5=(v.x*v.x+v.y*v.y),tmp6=(w.x*w.x+w.y*w.y),tmp7=v.y-w.y;
A=u.x*(tmp7)-u.y*(tmp1)+tmp2-tmp3;
B=tmp4*(-tmp7)+tmp5*(u.y-w.y)+tmp6*(v.y-u.y);
C=tmp4*(tmp1)+tmp5*(w.x-u.x)+tmp6*(u.x-v.x);
D=tmp4*(tmp3-tmp2)+tmp5*(u.x*w.y-w.x*u.y)+tmp6*(v.x*u.y-u.x*v.y);
return get_cir((pt){-B/(2.0*A),-C/(2.0*A)},sqrt((B*B+C*C-4*A*D)/(4*A*A)));
}
cir cir_cov(vector <pt> p){ //最小圆覆盖
shuffle(p.begin(),p.end(),rnd);int siz=p.size(); //注意要随机打乱才能保证复杂度
cir c=get_cir((pt){0,0},0);
for(int i=0;i<siz;i++){
if(!c.in(p[i])){
c=get_cir(p[i],0);
for(int j=0;j<i;j++){
if(!c.in(p[j])){
c=get_cir((p[i]+p[j])/2.0,dis(p[i],p[j])/2.0);
for(int k=0;k<j;k++)if(!c.in(p[k]))c=get_cir(p[i],p[j],p[k]);
}
}
}
}
return c;
}
}
using namespace circle;
const int N=2e6+7; //记得对答案特判 -0 的情况避免被创飞
signed main(){
ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
return 0;
}
浙公网安备 33010602011771号