计算几何算法汇总

# 0.前面的话 首先,笔者是一个初三蒟蒻,数学(几何)还特别菜。 所以如果什么地方讲错了或者比较naive,请轻喷。 解析法由于斜率不存在等等麻烦的情况,容易写错,所以本文主要讨论**向量法**。 (如果您是初中巨佬,还不熟知向量的应用,建议自学一下高中数学的相关内容) 向量法由于三角函数等等原因,跑的比较慢,所以计算几何的出题人一般不会卡常数。 我们实现函数的时候要以**代码补偿**为主,不要过度追求常数而把自己绕晕。 本文中所有的证明均为**方便记忆**之用,可能并不严格或者漏了情况。 这篇文章内容多了之后说不定会分块。 # 1.简单函数 - ### 约定/宏 ```cpp #define db double const db eps=1e-?; //根据需要调整,一般在 1e-6 ~ 1e-12 之间 ``` - ### 符号函数`sign` 众所周知,浮点数运算是有**误差**的,那么就需要一个符号函数,忽略较小的误差。 ```cpp int sign(double x) {return x>eps ? 1 : (x<-eps ? -1 : 0);} ``` - ## 点(向量)结构体`Point` 有加减和数乘几种运算。 len是求模长。 ```cpp struct Point { db x,y; inline void print() {printf("(%.3lf,%.3lf)\n",x,y);} Point operator + (const Point &B) const {return (Point){x+B.x,y+B.y};} Point operator - (const Point &B) const {return (Point){x-B.x,y-B.y};} Point operator * (const double &B) const {return (Point){x*B,y*B};} Point operator / (const double &B) const {return (Point){x/B,y/B};} double len() {return sqr(x*x+y*y);} }; ``` - ## 叉积,内积 封装在结构体里面,`|`运算是内积,`^`运算是叉积(小心优先级!) ```cpp double operator | (const Point &B) const {return x*B.x+y*B.y;} double operator ^ (const Point &B) const {return x*B.y-y*B.x;} ``` 向量内积: $\vec a·\vec b=|a||b|cosθ=a.x*b.x+a.y*b.y$ 根据 $|a||b|cosθ$ 这个式子,向量内积的数值就等于 $a$ 在 $b$ 上的投影与 $b$ 的模长积。 由于包含投影,可以用来判断垂直 : 内积为$0$即垂直。 假如两个向量的方向较为相同,即夹角小于 $90$ 度,那么点积的结果 $>0$ ; 假如两个向量的方向较为不同,即夹角大于 $90$ 度,那么点积的结果 $<0$ ; 两个向量同时旋转相同角度,其内积结果不变(显然)。 向量叉积: $\vec a\times \vec b=a.x*b.y-a.y*b.x$ 几何意义是 $\vec a,\vec b$ 两个向量张成的平行四边形(有向)面积,假如**逆时针则符号为正**,否则符号为负。 取绝对值除以2就是这两个向量围成的三角形面积。  可以用来算面积,判断方向等等。 ~~更多内容请见必修四~~ - ## 旋转 把某个点绕原点逆时针旋转一个弧度,可以使用复数相乘,幅角相加来推导。 ```cpp const double Pi=acos(-1); Point rotate(const Point &A,double len) { double s1=sin(len),s2=cos(len); return (Point){A.x*s2-A.y*s1,A.x*s1+A.y*s2}; } ``` - ## 线`Line` 可以表示直线,射线和线段。 ```cpp db dis(Point A,Point B) {return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y));} struct Line { Point a,b; double len(){return dis(a,b);} }; ``` - ## 点在直线/线段上 ```cpp //判断点是否在直线上 //之所以写l0是因为直线有0个端点 bool onl0(Point A,Line L) {return sign((A-L.a)^(A-L.b))==0;} //先判断是否在直线上,然后在判断是否在四至矩形内 //写l2是因为线段有2个端点 bool onl2(Point A,Line L) { return sign((A-L.a)^(A-L.b))==0&& sign(A.x-min(L.a.x,L.b.x))>=0&& sign(A.y-min(L.a.y,L.b.y))>=0&& sign(A.x-max(L.a.x,L.b.x))<=0&& sign(A.y-max(L.a.y,L.b.y))<=0; } ``` `sign((A-L.a)^(A-L.b))==0` 意思就是点 $A$ 是否在 $L$ 所在的直线上,具体原理 : 众所周知,点 $C-D$ 会形成向量 $\overrightarrow{DC}$ 那么这个式子的意思就是:判断 $\overrightarrow{AL_a},\overrightarrow{AL_b}$ 这两个向量叉积是否为$0$。 也就相当于判断 $a,L_a,L_b$ 三点围成的三角形面积。不难发现 : 三点共线 $\Longleftrightarrow$ 面积为$0$。 如图:  此后的四个条件是判断$A$是否在如右图矩阵中,注意可以视情况不取等号(一般要取)。 - ## 求两直线交点(常用) 调用该函数前请先排除两直线**重合**,**平行**等情况,否则会出现`nan`。 ```cpp Point inter(Line L1,Line L2) { double ls=(L1.b-L1.a)^(L2.a-L1.a); double rs=(L1.b-L1.a)^(L2.b-L1.a); return L1.a+(L1.b-L1.a)*ls/(ls-rs); } ``` 这东西的原理如下图:  下称 `L1.a→A1, L1.b→B1, L2.a→A2, L2.b→B2`。 算出 $\overrightarrow{A_1B_1}$ 和 $\overrightarrow{A_1A_2}$ 的叉积,由于是逆时针转,得到的数值是 $S\triangle A_1B_1A_2$ 的 $2$ 倍,称作 $ls$。 算出 $\overrightarrow{A_1B_1}$ 和 $\overrightarrow{A_1B_2}$ 的叉积,由于是顺时针转,得到的数值是 $S\triangle A_1B_1B_2$ 的 $-2$ 倍,称作 $rs$。 然后观察这两个三角形,假如把红线当做底,那么他们是等高的,底边比就相当于面积比。 我们把 $\overrightarrow{A_2B_2}$ 这个向量放缩到 $\Large\frac{ls}{ls-rs}$ 倍(注意符号),那么加上 $A_2$ 刚好够到那个点。 以上只是一种情况,经过严谨的讨论其他的情况下也都是正确的(使用向量一般不需要特判)。 记忆方法 : 从 $A1$ 开始弄三个向量,叉左边右边,把中间横着的缩放。 - ## 点到线段/直线的最短距离 ```cpp //计算点到直线距离,利用叉积+面积法 double disl0(Point A,Line L) {return abs((L.a-A)^(L.b-A))/dis(L.a,L.b);} //计算点到线段距离,利用点积判方向 double disl2(Point A,Line L) { if(sign(A-L.a|L.b-L.a)<0||sign(A-L.b|L.a-L.b)<0) return min(dis(A,L.a),dis(A,L.b)); return disl0(A,L); } ``` 如果是直线,求出垂线段长度即可。 方法 : 求出直线上两点与该点围成的三角形的面积,然后除以三角形的底。  线段的话,如果垂线段不落在线段上,那就不算数。 我们就要判断一个点有没有资格使用垂线段。  由上图可知,只有蓝色区域内的点有资格使用垂线段,外面的点直接对到两端点距离取min就好了 要判断深蓝色的两个夹角,( 有一个大于90度$\Longleftrightarrow$不在蓝色区域内 ),判断方法是点积。 - ## 完整模板Code 有效部分约长1.7K. ```cpp #include<algorithm> #include<cstdio> #include<cmath> #define db double using namespace std; const db Pi=acos(-1),eps=1e-8; int sign(db x) {return x<eps ? (-eps<x ? 0 : -1 ) : 1;} struct Point{ db x,y; void print() {printf("(%.3lf,%.3lf)\n",x,y);} Point operator + (const Point &B) const {return (Point){x+B.x,y+B.y};} Point operator - (const Point &B) const {return (Point){x-B.x,y-B.y};} Point operator * (const db &B) const {return (Point){x*B,y*B};} Point operator / (const db &B) const {return (Point){x/B,y/B};} db operator | (const Point &B) const {return x*B.x+y*B.y;} db operator ^ (const Point &B) const {return x*B.y-y*B.x;} }; Point rotate(const Point &A,db len){ db lx=cos(len),ly=sin(len); return (Point){A.x*lx-A.y*ly,A.y*lx+A.x*ly}; } db dis(Point A,Point B) {return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y));} struct Line{ Point a,b; db len(){return dis(a,b);} }; //判断点是否在直线上 //之所以写l0是因为直线有0个端点 bool onl0(Point A,Line L) {return sign((L.a-A)^(L.b-A))==0;} //先判断是否在直线上,然后在判断是否在四至矩形内 //写l2是因为线段有2个端点 bool onl2(Point A,Line L) { return onl0(A,L) &&sign(A.x-min(L.a.x,L.b.x))>=0 &&sign(A.y-min(L.a.y,L.b.y))>=0 &&sign(A.y-max(L.a.y,L.b.y))<=0 &&sign(A.y-max(L.a.y,L.b.y))<=0; } //算出左右面积,然后根据比例决定 Point inter(Line L1,Line L2) { db ls=(L1.b-L1.a)^(L2.a-L1.a), rs=(L1.b-L1.a)^(L2.b-L1.a); return (L2.b-L2.a)*(ls/(ls-rs))+L2.a; } //计算点到直线距离,利用叉积+面积法 db disl0(Point A,Line L) {return abs((L.a-A)^(L.b-A))/L.len();} //计算点到线段距离,利用点积判方向 db disl2(Point A,Line L) { if (sign((A-L.a)|(L.b-L.a))<=0||sign((A-L.b)|(L.a-L.b))<=0) return min(dis(A,L.a),dis(A,L.b)); return disl0(A,L); } void test1() { Point a=(Point){1,2},b=(Point){4,5}; printf("%.3lf ",dis(a,b)); a=rotate(a,Pi/4);b=rotate(b,Pi/2); printf("%.3lf\n",dis(a,b)); //4.243 4.686 puts(""); } void test2() { Point a=(Point){1,2},b=(Point){4,5},c=(Point){3,3},d=(Point){1.5,2.5}; Line l=(Line){(Point){0,1},(Point){3,4}}; printf("%d %d %d %d %d\n",onl2(a,l),onl2(b,l),onl2(c,l),onl0(b,l),onl2(d,l)); //1 0 0 1 1 puts(""); } void test3() { Line l1=(Line){(Point){0,1},(Point){3,4}}; Line l2=(Line){(Point){-2,1},(Point){-1,3}}; inter(l1,l2).print(); //(-4.000,-3.000) l1.b=l1.b*2; inter(l1,l2).print(); //(-4.800,-4.600) puts(""); } void test4() { Point a=(Point){0,2},b={-2,2},c={4,5}; Line l=(Line){(Point){0,1},(Point){3,4}}; printf("%.3lf %.3lf %.3lf\n",disl0(a,l),disl0(b,l),disl0(c,l)); //0.707 2.121 0.000 printf("%.3lf %.3lf %.3lf\n",disl2(a,l),disl2(b,l),disl2(c,l)); //0.707 2.236 1.414 puts(""); } int main() { test1(); test2(); test3(); test4(); return 0; } ``` 默写纸: ```cpp #include<algorithm> #include<cstdio> #include<cmath> #define db double using namespace std; const db Pi=acos(-1),eps=1e-8; int sign(double x) {} struct Point { db x,y; inline void print() {printf("(%.3lf,%.3lf)\n",x,y);} }; Point rotate(const Point &A,double len) { } double dis(Point A,Point B) {} struct Line { }; bool onl0(Point A,Line L) {} bool onl2(Point A,Line L) { } Point inter(Line L1,Line L2) { } double disl0(Point A,Line L) {} double disl2(Point A,Line L) { } void test1() { Point a=(Point){1,2},b=(Point){4,5}; printf("%.3lf ",dis(a,b)); a=rotate(a,Pi/4);b=rotate(b,Pi/2); printf("%.3lf\n",dis(a,b)); //4.243 4.686 puts(""); } void test2() { Point a=(Point){1,2},b=(Point){4,5},c=(Point){3,3},d=(Point){1.5,2.5}; Line l=(Line){(Point){0,1},(Point){3,4}}; printf("%d %d %d %d %d\n",onl2(a,l),onl2(b,l),onl2(c,l),onl0(b,l),onl2(d,l)); //1 0 0 1 1 puts(""); } void test3() { Line l1=(Line){(Point){0,1},(Point){3,4}}; Line l2=(Line){(Point){-2,1},(Point){-1,3}}; inter(l1,l2).print(); //(-4.000,-3.000) l1.b=l1.b*2; inter(l1,l2).print(); //(-4.800,-4.600) puts(""); } void test4() { Point a=(Point){0,2},b={-2,2},c={4,5}; Line l=(Line){(Point){0,1},(Point){3,4}}; printf("%.3lf %.3lf %.3lf\n",disl0(a,l),disl0(b,l),disl0(c,l)); //0.707 2.121 0.000 printf("%.3lf %.3lf %.3lf\n",disl2(a,l),disl2(b,l),disl2(c,l)); //0.707 2.236 1.414 puts(""); } int main() { test1(); test2(); test3(); test4(); return 0; } ``` # 2.简单练习 - [P1153 点和线](https://www.luogu.com.cn/problem/P1153) 看到最多 $10$ 个点,很明显这是一道爆搜题。 具体的思路就是暴力dfs,每加入一个点就判断有没有和前面已经练出来的线段相交。 注意:暴力dfs的话每种合法方案会被统计 $2$ 次(顺逆时针各一次),最后要除掉。 由于常数大,需要预处理所有 $O(n^2)$ 条线段的两两相交情况。 [评测记录](https://www.luogu.com.cn/record/22155996) - [P3744 李彬的几何](https://www.luogu.com.cn/problem/P3744) **题意** :按照顺时针顺序给出一个凸多边形,让你移动它的顶点使得它把它变成一个非凸多边形。 指定一个实数 $D$ ,然后将每个顶点最多移动 $D$ 个单位距离,问 $D$ 最小是多少。 这是一道基础的贪心题,显然要把凸包变凹,等价于制造三点共线,因为共线之后再往里移动一点点就不是个凸包了。 可以证明,选择凸包上相邻的三个点并令其共线是较优的,我们算出一个点到两边点连成直线的距离,然后除以 $2$ (因为三个点可以同时移动)取$\min$即可。 [评测记录](https://www.luogu.com.cn/record/22195192) - [P1354 房间最短路问题](https://www.luogu.com.cn/problem/P1354) 容易证明,最优路径的拐点肯定都是在缺口边缘。 由于数据范围极小,根本不用在意复杂度。 我们把缺口之间(包括起点终点)两两连边,然后跑folyd即可。 问题在于如何判断某两个点之间能否直线到达。 把有缺口的墙拆成 $3n$ 条线段,然后暴力判相交即可。 注意把线段在端点处削短一点,否则可能有在端点处算作相交的情况。 (可以不造四个墙壁[滑稽]) [评测记录](https://www.luogu.com.cn/record/22214359) - [P3217 [HNOI2011]数矩形](https://www.luogu.com.cn/problem/P3217) 题意:找4个点形成**矩形**并且面积最大，问面积。点数$\leq 1500$。 首先,四个点形成矩形的充要条件是 : 对角线相等且互相平分。 那么我们枚举每两个点,产生的线段按照中点和长度分类,可以通过排序实现。 被分到同一类的对角线,任选两条都可以弄成一个矩形,但是如果暴力枚举的话可能会超时。 把所有候选对角线放在一起大概是这个样子:  我们发现,对于一条候选对角线,肯定是选择端点尽量远的另一条。 这个东西是单峰的,满足某种决策单调性,所以两个指针一起转这个圈圈就好了。 总复杂度 $O(n^2logn)$ ,瓶颈在排序分类。 [todo] ------------ # 3.算法进阶 - ## 凸包 模板题:[P2742 【模板】二维凸包](https://www.luogu.com.cn/problem/P2742) 这个东西很重要,而且你会在很多貌似和计算几何没有关系的地方见到它。 **定义** : 包住平面上某个点集的**周长最小**的简单多边形,可以证明一定是凸多边形。 **感性理解:橡皮筋套钉子**。 凸包的求解方法: 首先先按照水平序排序 ( $x$ 相同比 $y$ ,否则比 $x$ )。 可以利用这个凸性来构造,维护一个单调栈,从左往右扫,如果遇见一个点和栈顶的两个点成顺时针夹角,那么就弹掉一个点,最后把新来的点加上。 这样就能保证我们求出了一个斜率递增(或者说,逆时针转)的下凸壳。 反过来再扫一遍则得到上凸壳(同样是逆时针转)。 (注意有些题目有重点,必须去掉) [评测记录](https://www.luogu.com.cn/record/22207317) ```cpp #include<algorithm> #include<cstdio> #include<cmath> #define eps 1e-7 using namespace std; int sign(double x) {return x>eps ? 1 : (x<-eps ? -1 : 0);} struct Point { double x,y; Point operator - (const Point &B) const {return (Point){x-B.x,y-B.y};} Point operator * (const double &B) const {return (Point){x*B,y*B};} double operator ^ (const Point &B) const {return x*B.y-y*B.x;} }p[10500]; bool cmp(Point A,Point B) {return sign(A.x-B.x)==0 ? A.y<B.y : A.x<B.x;} double dis(Point A,Point B) {return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y));} bool anticlock(Point A,Point B,Point C) {return sign((A-B)^(C-B))>=0;} //不写等于号可以抛弃共线的点,在这道题里面没有影响 int n,stk[10500],tot; double ans; int main() { scanf("%d",&n); for (int i=1;i<=n;i++) scanf("%lf%lf",&p[i].x,&p[i].y); sort(p+1,p+n+1,cmp); stk[1]=1;stk[tot=2]=2; for (int i=3;i<=n;i++){ while(tot>=2&&!anticlock(p[i],p[stk[tot]],p[stk[tot-1]])) tot--; stk[++tot]=i; } for (int i=1;i<tot;i++) ans+=dis(p[stk[i]],p[stk[i+1]]); stk[1]=n;stk[tot=2]=n-1; for (int i=n-2;i;i--){ while(tot>=2&&!anticlock(p[i],p[stk[tot]],p[stk[tot-1]])) tot--; stk[++tot]=i; } for (int i=1;i<tot;i++) ans+=dis(p[stk[i]],p[stk[i+1]]); printf("%.2lf",ans); return 0; } ``` - ## 旋转卡壳 经典问题 : **平面最远点对** 题目 : [P1452 Beauty Contest](https://www.luogu.com.cn/problem/P1452) 们明显,平面最远点对一定都在凸包上,所以先求出凸包。 定义凸包上的**对踵点对** : 用两条平行直线卡着凸包转,着两条直线一定会卡住至少两个点,这两个点称为对踵点对。  <https://jvruo.com/usr/uploads/2018/02/2986768821.gif> (图可能会挂,给个链接) 首先证明一个引理(**重要**) : 只用考虑斜率恰好与凸包某条边相同的直线。 容易发现,斜率恰好与凸包某条边相同的直线卡住凸包,必然在某一侧同时卡住两个点。 (为了避免卡住多个点的情况,求凸包时不应保留多个共线的点) 否则必然是一边卡住一个点。 那我们可以把两条直线逆时针转一下,转到卡住多一个点为止,那么原来的两个点仍然在直线上,直线的斜率变化得与凸包某条边相同。 对踵点对的个数为 $O(n)$ ,证明如下: 由上面的引理,我们只用考虑斜率恰好与凸包某条边相同的直线,这样的直线有 $n$ 条。 那么直线在对面只会卡住 $1$ 或 $2$ 个点,产生 $O(1)$ 个对踵点对。 所以对踵点对总量为 $O(n)*O(1)=O(n)$ ,证毕。 容易发现,平面最远点对∈对踵点对,采用反证法容易证明。 那我们求出所有的对踵点对,然后对距离取$\max$即可。 问题在于如何求对踵点对? (以下内容可以结合[DP的决策单调性优化总结 **/** 常用优化手段 **/** 单峰性-单指针跳跃-1D](https://www.luogu.com.cn/blog/command-block/dp-di-jue-ce-dan-diao-xing-you-hua-zong-jie)来理解) 我们逆时针枚举边,然后看看对面有什么被卡住了。 我们发现,由于凸包的凸性,对踵点对一定是单调逆时针变化的。 假设我们已经求出了上一次的对踵点对,我们把边逆时针转动时,另一边的对踵点也会逆时针转动。 而且,由于凸包的凸性,逆时针方向的点到该直线的距离一定是个单峰函数,而且峰顶就是对踵点对。 转动点的时候,我们可以计算这个点到该边的距离,如果小于逆时针方向的下一个点,就跳到下一个点,如果比相邻的两个点都牛逼,那么这就是对踵点。 图解:  可能同时卡住对面的两个点,这时,需要特判产生4双对踵点对。 我们考虑一个现象 : 假如把所有的点轻微扰动,可以认为卡住对面的两个点这种情况不存在。 而轻微扰动对大多数东西的影响不大,所以我们在大多数题目里面,可以不考虑卡住对面的两个点这种情况。 比如这道题叫我们求的是最远点对,那么可以不特判。 还要注意凸包只有两个点的情况。 ```cpp #include<algorithm> #include<cstdio> #include<cmath> #define eps 1e-7 #define MaxN 50500 using namespace std; int sign(double x) {return x>eps ? 1 : (x<-eps ? -1 : 0);} struct Point { double x,y; Point operator - (const Point &B) const {return (Point){x-B.x,y-B.y};} Point operator * (const double &B) const {return (Point){x*B,y*B};} double operator ^ (const Point &B) const {return x*B.y-y*B.x;} }p[MaxN],t[MaxN];int tot; bool cmp(Point A,Point B) {return sign(A.x-B.x)==0 ? A.y<B.y : A.x<B.x;} double dis(Point A,Point B) {return (A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y);} bool anticlock(Point A,Point B,Point C) {return sign((A-B)^(C-B))>0;} double tris(Point A,Point B,Point C) {return abs((B-A)^(C-A));} int n,stk[MaxN],top; void hull() { sort(p+1,p+n+1,cmp); stk[1]=1;stk[top=2]=2; for (int i=3;i<=n;i++){ while(top>=2&&!anticlock(p[i],p[stk[top]],p[stk[top-1]])) top--; stk[++top]=i; } for (int i=1;i<=top;i++)t[++tot]=p[stk[i]]; stk[1]=n;stk[top=2]=n-1; for (int i=n-2;i;i--){ while(top>=2&&!anticlock(p[i],p[stk[top]],p[stk[top-1]])) top--; stk[++top]=i; } for (int i=2;i<top;i++)t[++tot]=p[stk[i]]; } double ans; void rothull() { if (tot==2){ ans=dis(t[1],t[2]); return ; }t[0]=t[tot]; int pos=2; for (int i=1;i<=tot;i++){ while(tris(t[pos%tot+1],t[i],t[i-1])>eps+tris(t[pos],t[i],t[i-1])) pos=pos%tot+1; ans=max(ans,dis(t[pos],t[i])); ans=max(ans,dis(t[pos],t[i-1])); } } int main() { scanf("%d",&n); for (int i=1;i<=n;i++) scanf("%lf%lf",&p[i].x,&p[i].y); hull(); rothull(); printf("%.0lf",ans); return 0; } ``` [评测记录](https://www.luogu.com.cn/record/22387228) - ## 闵可夫斯基和 **定义** ： 两个区域 $A,B$ 的闵可夫斯基和定义为 $\{a+b|a\in A,b\in B\}$. 本文只讨论凸包的闵可夫斯基和。如图： ( $\rm From : WinXP$ )  不难发现两个凸壳 $A,B$ 的闵可夫斯基和也是凸壳，而且其顶点一定是 $A,B$ 各取一点相加而得。 结论 : 将两个凸包上的边按照极角序顺次连接即可得到答案。(具体见代码) 如图： ( $\rm From : xzyxzy$ )  注意可能会出现**三点共线**。 **题目** : [P4557 [JSOI2018]战争](https://www.luogu.com.cn/problem/P4557) **题意** : 给出两个凸包 $A,B$ 和 $q$ 个向量 $\vec x_i$ ，求各个 $A+\vec x_i$ 和 $B$ 有没有交。 令 $a\in A,b\in B$ 则问题变为判定是否存在 $a+x=b$。 移项得到 $x=b-a$ ,又变成判定 $x$ 是否在 $B+(-A)$ 中，求闵可夫斯基和即可。 写法中,分别求上下凸壳,分别判定。 ```cpp #include<algorithm> #include<cstdio> #include<cmath> #define ll long long #define MaxN 100500 using namespace std; struct Point { ll x,y; Point operator + (const Point &B) const {return (Point){x+B.x,y+B.y};} Point operator - (const Point &B) const {return (Point){x-B.x,y-B.y};} Point operator * (const ll &B) const {return (Point){x*B,y*B};} ll operator ^ (const Point &B) const {return x*B.y-y*B.x;} }p1[MaxN],p2[MaxN],t0[MaxN<<1],t1[MaxN<<1]; bool cmpP(const Point &A,const Point &B) {return A.x==B.x ? A.y<B.y : A.x<B.x;} bool cmpK(const Point &A,const Point &B) {return A.y*B.x<B.y*A.x;} bool antic(Point A,Point B,Point C) {return ((A-B)^(C-B))>0;} void hull(Point *p,int n,Point *t0,int &m0,Point *t1,int &m1) { static int stk[MaxN],tot; sort(p+1,p+n+1,cmpP); stk[1]=1;stk[tot=2]=2; for (int i=3;i<=n;i++){ while(tot>=2&&!antic(p[i],p[stk[tot]],p[stk[tot-1]])) tot--; stk[++tot]=i; } for (int i=1;i<tot;i++) t0[++m0]=p[stk[i+1]]-p[stk[i]]; stk[1]=n;stk[tot=2]=n-1; for (int i=n-2;i;i--){ while(tot>=2&&!antic(p[i],p[stk[tot]],p[stk[tot-1]])) tot--; stk[++tot]=i; } for (int i=tot-1;i;i--) t1[++m1]=p[stk[i]]-p[stk[i+1]]; } int n1,n2,m0,m1,q; int main() { scanf("%d%d%d",&n1,&n2,&q); for (int i=1;i<=n1;i++) scanf("%lld%lld",&p1[i].x,&p1[i].y); hull(p1,n1,t0,m0,t1,m1); for (int i=1;i<=n2;i++){ scanf("%lld%lld",&p2[i].x,&p2[i].y); p2[i].x*=-1;p2[i].y*=-1; }hull(p2,n2,t0,m0,t1,m1); sort(t0+1,t0+m0+1,cmpK); t0[0]=p1[1]+p2[1]; for (int i=1;i<=m0;i++) t0[i]=t0[i]+t0[i-1]; int l0=0,r0=m0; while(t0[l0].x==t0[l0+1].x)l0++; while(t0[r0].x==t0[r0-1].x)r0--; sort(t1+1,t1+m1+1,cmpK); reverse(t1+1,t1+m1+1); t1[0]=p1[1]+p2[1]; for (int i=1;i<=m1;i++) t1[i]=t1[i]+t1[i-1]; int l1=0,r1=m1; while(t1[l1].x==t1[l1+1].x)l1++; while(t1[r1].x==t1[r1-1].x)r1--; Point d; for (int i=1;i<=q;i++){ scanf("%lld%lld",&d.x,&d.y); if (d.x<t0[0].x||t0[m0].x<d.x) {puts("0");continue;} int l=l0,r=r0,mid; while(l<r){ mid=(l+r+1)>>1; if (t0[mid].x<=d.x)l=mid; else r=mid-1; }if ((l<r0&&((t0[l+1]-t0[l])^(d-t0[l]))<0) ||l==r0&&d.y<t0[r0].y) {puts("0");continue;} l=l1;r=r1; while(l<r){ mid=(l+r+1)>>1; if (t1[mid].x<=d.x)l=mid; else r=mid-1; }if ((l<r1&&((t1[l+1]-t1[l])^(d-t1[l]))>0) ||l==r1&&d.y>t1[r1].y) {puts("0");continue;} puts("1"); }return 0; } ``` - ## 半平面交 题目 : [P4196 [CQOI2006]凸多边形 /【模板】半平面交](https://www.luogu.com.cn/problem/P4196) 容易转化成下列问题: 给出若干条**有向**直线,每条直线的**右边**算作一个半平面,询问所有半平面的交集面积大小,保证交集有界。 也容易证明交集一定是个凸包,感性理解就是 : 削苹果不停刀一定削不出坑来。 斜率优化有两种写法,一种是基于凸包,另一种就是基于半平面交的,我个人认为后者更易于理解,这也说明两者有很多相似之处。 先讲左凸壳 : 按照幅角排序(可以`atan2`或者手写斜率),然后依次加入,维护一个栈。 假如上一个交点不在新加入的半平面内(叉积判断即可),那么可以把上一条线弹出,如图:  这样逆时针转一圈就得到了凸壳,最后要判一下首尾两条线,如图:  还是把交点在班平面外的线都弹掉即可。 如果保证答案有界,最后至少剩下三条线,否则无解。 某些题目里需要自行加边界。 ```cpp ``` - ## 平面最近点对分治 [P1429 平面最近点对（加强版）](https://www.luogu.com.cn/problem/P1429) 应用范围好像不广的样子,不过鉴于其思想精妙还是特别讲讲。 首先,我们按照水平序排序。 我们现在有 $n$ 个点,然后根据 $x$ 坐标开始分治(被分割线穿过的点随便丢到任意一边)。 求出了左边和右边的最近点对之后,我们更新目前的全局 $ans$。 这里如果枚举左右点对合并,复杂度肯定假了。 显然,$x$ 坐标离分割线如果距离小于 $ans$ 才有可能更优,我们可以只枚举这些点。 但是左右的这些点集仍然有可能是 $O(n)$ 大小的,直接枚举仍然有可能TLE. 上面我们只限制了 $x$ 坐标,没有限制 $y$ 坐标,导致复杂度退化。 我们把离分割线近的那些点拿出来按照 $y$ 坐标排序,然后限制两个点之间的 $y$ 坐标差不超过 $ans$ ,这样子左边的点一定对应右边的一个区间(注意要一直动态更新 $ans$ ,否则复杂度是假的) 可以证明每个点对应点不超过 $9$ 个: 我们的限制框出了一个 $3ans*3ans$ 的矩形,我们往里面放置直径为 $ans$ 的圆,需要保证圆与圆之间不相交,最多能放置几个。容易想象最多 $9$ 个 (反正是$O(1)$个)。 如果每次对 $y$ 排序使用`std::sort`,复杂度是$O(n\log^2n)$,但是由于数据随机跑得很快。 严谨的方法是按照 $y$ 归并排序,然后合并的时候再取出 $x$ 满足限制的那些点,这样子复杂度就是$O(n\log n)$了。 ```cpp #include<algorithm> #include<cstdio> #include<cmath> #define MaxN 200500 #define eps 1e-6 using namespace std; struct Point {double x,y;}p[MaxN],t[MaxN]; bool cmp(Point A,Point B) {return abs(A.x-B.x)<eps ? A.y<B.y : A.x<B.x;} double dis(Point A,Point B) {return (A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y);} double ans; void solve(int l,int r) { if (l==r)return ; int mid=(l+r)>>1; double ml=p[mid].x; solve(l,mid);solve(mid+1,r); int pl=l,pr=mid+1,tot=l-1; while(tot<r){ if (pl<=mid&&(p[pl].y<p[pr].y||pr>r))t[++tot]=p[pl++]; else t[++tot]=p[pr++]; }for (int i=l;i<=r;i++) p[i]=t[i]; tot=0; for (int i=l;i<=r;i++) if ((p[i].x-ml)*(p[i].x-ml)+eps<ans) t[++tot]=p[i]; for (int i=1;i<=tot;i++) for (int j=i+1;j<=tot;j++) if ((t[j].y-t[i].y)*(t[j].y-t[i].y)+eps<ans) ans=min(ans,dis(t[i],t[j])); else break; } int n; int main() { scanf("%d",&n); for (int i=1;i<=n;i++) scanf("%lf%lf",&p[i].x,&p[i].y); sort(p+1,p+n+1,cmp); ans=dis(p[1],p[2]); solve(1,n); printf("%.4lf",sqrt(ans)); return 0; } ``` [评测记录](https://www.luogu.com.cn/record/22222511) # 4.进阶题 - [P2116 城墙](https://www.luogu.com.cn/problem/P2116) 如果 $L=0$ ,那么直接求凸包就好了。 否则,我们可以把点变成半径为 $L$ 的柱子,然后再围起来。 不难发现,有弧度的部分合起来是一个半径为 $L$ 的正圆。 那么剩下的直线部分,平移之后能够拼成凸包,答案就是凸包周长$+2\pi L$。 [评测记录](https://www.luogu.com.cn/record/22210824) - [P3829 [SHOI2012]信用卡凸包](https://www.luogu.com.cn/problem/P3829) 和上一题类似,也是求出凸包周长再加上$2\pi r$,注意要去重($ch\acute ong$)点。 [评测记录](https://www.luogu.com.cn/record/22224474) - SMOJ2919 神仙题.jpg 题意:给你平面上n个点,点有点权(可正可负),要求你选择三个点形成**劣角**。 然后将会得到劣角内部所有点点权和的收益,注意,选择的三个点不计入收益。 最大化这个东西。$n\leq2000$ 枚举这个角的顶点,然后移动坐标系再**极角排序**,这样子能选择的东西就变成序列了,使用单队维护一下便可。 - [P4048 [JSOI2010]冷冻波](https://www.luogu.com.cn/problem/P4048) 这是一道结合了网络流与计算几何的题目。 首先建立网络流模型 : 二分时间限制,巫妖连上汇点,容量为 : 攻击次数=(时间/冷却)+1 精灵连上能够攻击他的巫妖,容量为 $1$ ,源点连精灵,容量为 $1$ 。 如果满流则方案可行。 剩下的问题就是判断某个巫妖能否攻击某个精灵。 场地里面有唯一的障碍物 : 树木。可以抽象成圆形。 我们判断一下巫妖和精灵连成的线段后没有和某个圆相交就好了(注意距离还要小于R)。 直接求线段到圆的距离,如果小于等于半径则认为有交。 没有实体在树里面,所以求端点距离可以省了。 复杂度$O(n^3+\text{网络流}log?)$ [评测记录](https://www.luogu.com.cn/record/22210824) - [P4166 [SCOI2007]最大土地面积](https://www.luogu.com.cn/problem/P4166) 题意:给出 $n$ 个点,选出四个使得面积最大。 容易证明,这四个点一定都在凸包上。 (如果凸包上只有三个点,那么直接枚举第四个点即可,正确性显然) 然后我们固定一点, $O(n)$ 枚举一条对角线(对点) 剩下的两个点就在两侧,他们对面积的贡献是与对角线连成的三角形的面积,由于该三角形底一定,那么要最大化三角形的高。 我们发现,随着对角线的转动,这两个点肯定不会往回转,因为往回转会使得高变小。 此外,高在逆时针方向上还是个单峰函数。 那么四个点一起转圈圈就好了,复杂度 $O(n^2)$。 还有基于对踵点的 $O(n\log n)$ 做法: 注意到最优解中对角线一定是一对对踵点,所以旋转卡壳先求出所有候选对角线。 然后让候选对角线转逆时针圈圈,另外的两个点类似旋转卡壳一样转,即可线性扫描完毕。 [评测记录](https://www.luogu.com.cn/record/22388367) ($O(n\log n)$) - [P3187 [HNOI2007]最小矩形覆盖](https://www.luogu.com.cn/problem/P3187) 能够感知,最优解一定有一条边与凸包上的一条边重合。 我们就依次枚举这条边,然后在其余三个方向上求出最远点,这也是具有单调性的,可以四个点一起转。 复杂度 $O(n\log n)$。 [**Todo**] - [P3680 [CERC2016]凸轮廓线 Convex Contour](https://www.luogu.com.cn/problem/P3680) 毒瘤凸包题,不想多说。 [评测记录](https://www.luogu.com.cn/record/22392904)