计算几何算法汇总

· · 个人记录

0.前面的话

首先,笔者是一个初三蒟蒻,数学(几何)还特别菜。

所以如果什么地方讲错了或者比较naive,请轻喷。

解析法由于斜率不存在等等麻烦的情况,容易写错,所以本文主要讨论向量法。

(如果您是初中巨佬,还不熟知向量的应用,建议自学一下高中数学的相关内容)

向量法由于三角函数等等原因,跑的比较慢,所以计算几何的出题人一般不会卡常数。

我们实现函数的时候要以代码补偿为主,不要过度追求常数而把自己绕晕。

本文中所有的证明均为方便记忆之用,可能并不严格或者漏了情况。

这篇文章内容多了之后说不定会分块。

1.简单函数

#define db double
const db eps=1e-?; //根据需要调整,一般在 1e-6 ~ 1e-12 之间

众所周知,浮点数运算是有误差的,那么就需要一个符号函数,忽略较小的误差。

int sign(double x)
{return x>eps ? 1 : (x<-eps ? -1 : 0);}

有加减和数乘几种运算。

len是求模长。

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);}
};

封装在结构体里面,|运算是内积,^运算是叉积(小心优先级!)

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就是这两个向量围成的三角形面积。

可以用来算面积,判断方向等等。

更多内容请见必修四

把某个点绕原点逆时针旋转一个弧度,可以使用复数相乘,幅角相加来推导。

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};
}

可以表示直线,射线和线段。

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);}
};
//判断点是否在直线上
//之所以写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。

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 开始弄三个向量,叉左边右边,把中间横着的缩放。

//计算点到直线距离,利用叉积+面积法 
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不在蓝色区域内 ),判断方法是点积。

有效部分约长1.7K.

#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;
}

默写纸:

#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.简单练习

看到最多 10 个点,很明显这是一道爆搜题。

具体的思路就是暴力dfs,每加入一个点就判断有没有和前面已经练出来的线段相交。

注意:暴力dfs的话每种合法方案会被统计 2 次(顺逆时针各一次),最后要除掉。

由于常数大,需要预处理所有 O(n^2) 条线段的两两相交情况。

评测记录

题意 :按照顺时针顺序给出一个凸多边形,让你移动它的顶点使得它把它变成一个非凸多边形。

指定一个实数 D ,然后将每个顶点最多移动 D 个单位距离,问 D 最小是多少。

这是一道基础的贪心题,显然要把凸包变凹,等价于制造三点共线,因为共线之后再往里移动一点点就不是个凸包了。

可以证明,选择凸包上相邻的三个点并令其共线是较优的,我们算出一个点到两边点连成直线的距离,然后除以 2 (因为三个点可以同时移动)取\min即可。

评测记录

容易证明,最优路径的拐点肯定都是在缺口边缘。

由于数据范围极小,根本不用在意复杂度。

我们把缺口之间(包括起点终点)两两连边,然后跑folyd即可。

问题在于如何判断某两个点之间能否直线到达。

把有缺口的墙拆成 3n 条线段,然后暴力判相交即可。

注意把线段在端点处削短一点,否则可能有在端点处算作相交的情况。

(可以不造四个墙壁[滑稽])

评测记录

题意:找4个点形成矩形并且面积最大，问面积。点数\leq 1500。

首先,四个点形成矩形的充要条件是 : 对角线相等且互相平分。

那么我们枚举每两个点,产生的线段按照中点和长度分类,可以通过排序实现。

被分到同一类的对角线,任选两条都可以弄成一个矩形,但是如果暴力枚举的话可能会超时。

把所有候选对角线放在一起大概是这个样子:

我们发现,对于一条候选对角线,肯定是选择端点尽量远的另一条。

这个东西是单峰的,满足某种决策单调性,所以两个指针一起转这个圈圈就好了。

总复杂度 O(n^2logn) ,瓶颈在排序分类。

[todo]

3.算法进阶

模板题:P2742 【模板】二维凸包

这个东西很重要,而且你会在很多貌似和计算几何没有关系的地方见到它。

定义 : 包住平面上某个点集的周长最小的简单多边形,可以证明一定是凸多边形。

感性理解:橡皮筋套钉子。

凸包的求解方法:

首先先按照水平序排序 ( x 相同比 y ,否则比 x )。

可以利用这个凸性来构造,维护一个单调栈,从左往右扫,如果遇见一个点和栈顶的两个点成顺时针夹角,那么就弹掉一个点,最后把新来的点加上。

这样就能保证我们求出了一个斜率递增(或者说,逆时针转)的下凸壳。

反过来再扫一遍则得到上凸壳(同样是逆时针转)。

(注意有些题目有重点,必须去掉)

评测记录

#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://jvruo.com/usr/uploads/2018/02/2986768821.gif (图可能会挂,给个链接)

首先证明一个引理(重要) : 只用考虑斜率恰好与凸包某条边相同的直线。

容易发现,斜率恰好与凸包某条边相同的直线卡住凸包,必然在某一侧同时卡住两个点。

(为了避免卡住多个点的情况,求凸包时不应保留多个共线的点)

否则必然是一边卡住一个点。

那我们可以把两条直线逆时针转一下,转到卡住多一个点为止,那么原来的两个点仍然在直线上,直线的斜率变化得与凸包某条边相同。

对踵点对的个数为 O(n) ,证明如下:

由上面的引理,我们只用考虑斜率恰好与凸包某条边相同的直线,这样的直线有 n 条。

那么直线在对面只会卡住 1 或 2 个点,产生 O(1) 个对踵点对。

所以对踵点对总量为 O(n)*O(1)=O(n) ,证毕。

容易发现,平面最远点对∈对踵点对,采用反证法容易证明。

那我们求出所有的对踵点对,然后对距离取\max即可。

问题在于如何求对踵点对?

(以下内容可以结合DP的决策单调性优化总结 / 常用优化手段 / 单峰性-单指针跳跃-1D来理解)

我们逆时针枚举边,然后看看对面有什么被卡住了。

我们发现,由于凸包的凸性,对踵点对一定是单调逆时针变化的。

假设我们已经求出了上一次的对踵点对,我们把边逆时针转动时,另一边的对踵点也会逆时针转动。

而且,由于凸包的凸性,逆时针方向的点到该直线的距离一定是个单峰函数,而且峰顶就是对踵点对。

转动点的时候,我们可以计算这个点到该边的距离,如果小于逆时针方向的下一个点,就跳到下一个点,如果比相邻的两个点都牛逼,那么这就是对踵点。

图解:

可能同时卡住对面的两个点,这时,需要特判产生4双对踵点对。

我们考虑一个现象 : 假如把所有的点轻微扰动,可以认为卡住对面的两个点这种情况不存在。

而轻微扰动对大多数东西的影响不大,所以我们在大多数题目里面,可以不考虑卡住对面的两个点这种情况。

比如这道题叫我们求的是最远点对,那么可以不特判。

还要注意凸包只有两个点的情况。

#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;
}

评测记录

本文只讨论凸包的闵可夫斯基和。如图： ( \rm From : WinXP )

不难发现两个凸壳 A,B 的闵可夫斯基和也是凸壳，而且其顶点一定是 A,B 各取一点相加而得。

结论 : 将两个凸包上的边按照极角序顺次连接即可得到答案。(具体见代码)

如图： ( \rm From : xzyxzy )

注意可能会出现三点共线。

题目 : P4557 [JSOI2018]战争

题意 : 给出两个凸包 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) 中，求闵可夫斯基和即可。

写法中,分别求上下凸壳,分别判定。

#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]凸多边形 /【模板】半平面交

容易转化成下列问题:

给出若干条有向直线,每条直线的右边算作一个半平面,询问所有半平面的交集面积大小,保证交集有界。

也容易证明交集一定是个凸包,感性理解就是 : 削苹果不停刀一定削不出坑来。

斜率优化有两种写法,一种是基于凸包,另一种就是基于半平面交的,我个人认为后者更易于理解,这也说明两者有很多相似之处。

先讲左凸壳 : 按照幅角排序(可以atan2或者手写斜率),然后依次加入,维护一个栈。

假如上一个交点不在新加入的半平面内(叉积判断即可),那么可以把上一条线弹出,如图:

这样逆时针转一圈就得到了凸壳,最后要判一下首尾两条线,如图:

还是把交点在班平面外的线都弹掉即可。

如果保证答案有界,最后至少剩下三条线,否则无解。

某些题目里需要自行加边界。

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)了。

#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;
}

评测记录

4.进阶题

如果 L=0 ,那么直接求凸包就好了。

否则,我们可以把点变成半径为 L 的柱子,然后再围起来。

不难发现,有弧度的部分合起来是一个半径为 L 的正圆。

那么剩下的直线部分,平移之后能够拼成凸包,答案就是凸包周长+2\pi L。

评测记录

和上一题类似,也是求出凸包周长再加上2\pi r,注意要去重(ch\acute ong)点。

评测记录

题意:给你平面上n个点,点有点权(可正可负),要求你选择三个点形成劣角。

然后将会得到劣角内部所有点点权和的收益,注意,选择的三个点不计入收益。

最大化这个东西。n\leq2000

枚举这个角的顶点,然后移动坐标系再极角排序,这样子能选择的东西就变成序列了,使用单队维护一下便可。

这是一道结合了网络流与计算几何的题目。

首先建立网络流模型 :

二分时间限制,巫妖连上汇点,容量为 : 攻击次数=(时间/冷却)+1

精灵连上能够攻击他的巫妖,容量为 1 ,源点连精灵,容量为 1 。

如果满流则方案可行。

剩下的问题就是判断某个巫妖能否攻击某个精灵。

场地里面有唯一的障碍物 : 树木。可以抽象成圆形。

我们判断一下巫妖和精灵连成的线段后没有和某个圆相交就好了(注意距离还要小于R)。

直接求线段到圆的距离,如果小于等于半径则认为有交。

没有实体在树里面,所以求端点距离可以省了。

复杂度O(n^3+\text{网络流}log?)

评测记录

题意:给出 n 个点,选出四个使得面积最大。

容易证明,这四个点一定都在凸包上。

(如果凸包上只有三个点,那么直接枚举第四个点即可,正确性显然)

然后我们固定一点, O(n) 枚举一条对角线(对点)

剩下的两个点就在两侧,他们对面积的贡献是与对角线连成的三角形的面积,由于该三角形底一定,那么要最大化三角形的高。

我们发现,随着对角线的转动,这两个点肯定不会往回转,因为往回转会使得高变小。

此外,高在逆时针方向上还是个单峰函数。

那么四个点一起转圈圈就好了,复杂度 O(n^2)。

还有基于对踵点的 O(n\log n) 做法:

注意到最优解中对角线一定是一对对踵点,所以旋转卡壳先求出所有候选对角线。

然后让候选对角线转逆时针圈圈,另外的两个点类似旋转卡壳一样转,即可线性扫描完毕。

评测记录 (O(n\log n))

能够感知,最优解一定有一条边与凸包上的一条边重合。

我们就依次枚举这条边,然后在其余三个方向上求出最远点,这也是具有单调性的,可以四个点一起转。

复杂度 O(n\log n)。

[Todo]

毒瘤凸包题,不想多说。

评测记录