计算几何选讲
[SHOI2012]信用卡凸包
通过平移可以将最终的图形分成一个以信用卡四个角的圆心为顶点的凸包以及一个圆，分别计算即可。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
int n,top,m;
double ans=0,A,B,R,l,phi;
struct AS{
double x,y;
}a[100005],p[100005];
inline int read(){
int ret=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')f=-f;ch=getchar();}
while(ch<='9'&&ch>='0')ret=ret*10+ch-'0',ch=getchar();
return ret*f;
}
double check(AS A1,AS A2,AS B1,AS B2){
return (A2.x-A1.x)*(B2.y-B1.y)-(B2.x-B1.x)*(A2.y-A1.y);
}
double d(AS A,AS B){
return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y));
}
bool cmp(AS x,AS y){
double tmp=check(a[1],x,a[1],y);
if(tmp>0)return 1;
if(tmp==0&&(d(a[1],x)<d(a[1],y)))return 1;
return 0;
}
void add(double x,double y){
n++;
a[n].x=x;a[n].y=y;
if(a[n].y<a[1].y)swap(a[1],a[n]);
}
int main(){
scanf("%d%lf%lf%lf",&m,&A,&B,&R);A-=2*R;B-=2*R;
l=sqrt(A*A+B*B)/2;
phi=atan(A/B);
for(int i=1;i<=m;i++){
double x,y,theta;scanf("%lf%lf%lf",&x,&y,&theta);
double dx=cos(theta+phi)*l,dy=sin(theta+phi)*l;
add(x+dx,y+dy);add(x-dx,y-dy);
dx=cos(theta-phi)*l;dy=sin(theta-phi)*l;
add(x+dx,y+dy);add(x-dx,y-dy);
}
sort(a+2,a+1+n,cmp);
p[top=1]=a[1];
for(int i=2;i<=n;i++){
while(top>1&&check(p[top-1],p[top],p[top],a[i])<=0)top--;
top++;p[top]=a[i];
}
p[++top]=p[1];
for(int i=1;i<top;i++){
ans+=d(p[i],p[i+1]);
}
printf("%.2lf\n",ans+3.141592653589793*2*R);
return 0;
}
[NOI2005] 月下柠檬树
可以通过 \(\text{Simpson}\) 公式逼近
点击查看代码
#include<bits/stdc++.h>
using namespace std;
#define EPS 1e-7
double alpha;
int n;
struct circle{
double x, r;
}p[1000];
struct tan_line{
double k, b, left, right;
}q[1000];
inline double Gougu(double a,double b){
return sqrt(a*a-b*b);
}
inline void get_tan(int x,int y){
if (fabs(p[x].r-p[y].r)<EPS){
q[x].left=p[x].x;q[x].right=p[y].x;
q[x].k=0;q[x].b=p[x].r;
return;
}
double dx=p[y].x-p[x].x,dr=fabs(p[x].r-p[y].r);
double ly,ry;
if(p[x].r>p[y].r){
q[x].left=p[x].x+p[x].r*dr/dx;
q[x].right=p[y].x+(q[x].left-p[x].x)*p[y].r/p[x].r;
ly=Gougu(p[x].r,q[x].left-p[x].x);ry=Gougu(p[y].r,q[x].right-p[y].x);
q[x].k=(ly-ry)/(q[x].left-q[x].right);
q[x].b=ly-q[x].left*q[x].k;
}
else {
q[x].right=p[y].x-p[y].r*dr/dx;
q[x].left=p[x].x-(p[y].x-q[x].right)*p[x].r/p[y].r;
ly=Gougu(p[x].r,q[x].left-p[x].x);ry=Gougu(p[y].r,q[x].right-p[y].x);
q[x].k=(ly-ry)/(q[x].left-q[x].right);
q[x].b=ly-q[x].left*q[x].k;
}
}
inline double F(double x){
double ans=0.0;
for (int i=1;i<=n;i++){
if (x<p[i].x+p[i].r&&x>p[i].x-p[i].r)ans=(ans>Gougu(p[i].r,x-p[i].x)?ans:Gougu(p[i].r,x-p[i].x));
}
for(int i=1;i<=n;i++){
if(x>=q[i].left&&x<=q[i].right)ans=(ans>q[i].k*x+q[i].b?ans:q[i].k*x+q[i].b);
}
return ans;
}
inline double simpson(double a,double b){
double c=a+(b-a)/2.0;
return (b-a)*(F(a)+4.0*F(c)+F(b))/6.0;
}
inline double asr(double a,double b,double ans){
double c=a+(b-a)/2.0;
double left = simpson(a,c),right=simpson(c,b);
if(fabs(left+right-ans)<EPS)return left+right;
else return asr(a,c,left)+asr(c,b,right);
}
int main(){
scanf("%d%lf",&n,&alpha);
alpha=1.0/tan(alpha);
scanf("%lf", &p[1].x);
p[1].x*=alpha;
for(int i=2;i<=n+1;i++){
scanf("%lf",&p[i].x);
p[i].x*=alpha;p[i].x+=p[i-1].x;
}
for(int i=1;i<=n;i++)scanf("%lf",&p[i].r);
++n;
p[n].r=0.0;
for(int i=1;i<=n-1;i++)get_tan(i,i+1);
double ll=p[1].x-p[1].r,rr=p[n].x;
for(int i=1;i<=n;i++){
rr=rr>(p[i].x+p[i].r)?rr:(p[i].x+p[i].r);
ll=ll<(p[i].x-p[i].r)?ll:(p[i].x-p[i].r);
}
printf("%.2lf\n",2.0*asr(ll,rr,simpson(ll,rr)));
return 0;
}
「BZOJ3707」圈地
可以发现枚举完两点后最优点就是最接近该直线的点。
如果将该直线看作 \(y\) 轴用两边 \(x\) 坐标的绝对值最小的点更新即可。
于是可以将所有斜率排序后不断旋转 \(y\) 轴并维护当前按 \(x\) 从小到大的序列。
发现对于直线 \((x,y)\)，两点 \(x\) 坐标的相对关系仅在 \(y\) 轴斜率超过直线斜率后才会变化。
因此每处理完一条直线，将两点的顺序交换即可。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
int n;
struct node{
double x,y;
node(double _x=0,double _y=0){
x=_x;y=_y;
}
inline node operator -(const node &b)const{
return node(x-b.x,y-b.y);
}
inline double operator *(const node &b)const{
return x*b.y-y*b.x;
}
inline bool operator <(const node &b)const{
if(x!=b.x)return x<b.x;
return y<b.y;
}
}a[1005];
struct line{
int x,y;double k;
line(int _x,int _y,double _k){
x=_x;y=_y;k=_k;
}
inline bool operator <(const line &b)const{
return k>b.k;
}
};
inline double S(int i,int j,int k){
return fabs((a[j]-a[i])*(a[k]-a[i]));
}
int loc[1005],id[1005];
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)scanf("%lf%lf",&a[i].x,&a[i].y);
sort(a+1,a+n+1);
vector<line> vec;
for(int i=1;i<=n;i++){
for(int j=i+1;j<=n;j++){
double k=a[i].x==a[j].x?1e18:(a[i].y-a[j].y)/(a[i].x-a[j].x);
vec.push_back(line(i,j,k));
}
}
sort(vec.begin(),vec.end());
for(int i=1;i<=n;i++)loc[i]=i;
for(int i=1;i<=n;i++)id[i]=i;
double ans=1e18;
for(auto it:vec){
int x=loc[it.x],y=loc[it.y];
if(x>y)swap(x,y);
if(x!=1)ans=min(ans,S(id[x-1],it.x,it.y));
if(y!=n)ans=min(ans,S(id[y+1],it.x,it.y));
swap(id[x],id[y]);
swap(loc[it.x],loc[it.y]);
}
printf("%.2lf\n",ans/2);
return 0;
}
[HNOI2007]最小矩形覆盖
旋转卡壳。
点击查看代码
[Balkan2002]Alien最小圆覆盖
我们知道最小圆一定是点集中某三个点的外接圆，因此我们要找到这三个点。
求最小圆覆盖的算法是这样的：
先枚举一个 \(i\)，表示现在要求包含前 \(i\) 个点的最小圆。
对于每一个 \(i\)，如果它不在当前求出的最小圆中，则把当前最小圆修改成点 \(i\)，然后继续进行下面的操作，否则继续枚举 \(i\)。
枚举 \(j\)，表示现在要求包含第 \(i\) 个点，且包含前 \(j\) 个点的最小圆。
对于每一个 \(j\)，如果它不在当前求出的最小圆中，则把当前最小圆修改成以 \(i,j\) 两点连成的线段为直径的圆，然后继续进行下面的操作，否则继续枚举 \(j\)。
枚举 \(k\)，表示现在要求包含第 \(i,j\) 个点，且包含前 \(k\) 个点的最小圆。
对于每一个 \(k\)，如果它不在当前求出的最小圆中，则把当前最小圆修改成 \(i,j,k\) 三点组成的三角形的外接圆，否则继续枚举 \(k\)。
以上算法看似是 \(O(n^3)\) 的，但在点排列随机的基础上，算法实际上是期望 \(O(n)\) 的（因为一个点在外接圆外的概率比较小）。
点击查看代码
#include <bits/stdc++.h>
using namespace std;
int n;
struct node{
double x,y;
node(double _x=0,double _y=0){
x=_x;y=_y;
}
inline node operator -(const node &b)const{
return node(x-b.x,y-b.y);
}
inline node operator +(const node &b)const{
return node(x+b.x,y+b.y);
}
inline node operator /(double v){
return node(x/v,y/v);
}
inline double operator *(const node &b)const{
return x*b.y-y*b.x;
}
}a[100010];
double nowx,nowy,nowr;
inline double dis(double x1,double y1,double x2,double y2){
return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}
inline bool check(int i){
return dis(a[i].x,a[i].y,nowx,nowy)<=nowr;
}
inline node inter(node p,node v,node q,node w){
node u=p-q;
double t=(w*u)/(v*w);
return node(p.x+t*v.x,p.y+t*v.y);
}
int main(){
srand(time(0));
scanf("%d",&n);
for(int i=1;i<=n;i++)scanf("%lf%lf",&a[i].x,&a[i].y);
random_shuffle(a+1,a+n+1);
nowx=a[1].x,nowy=a[1].y,nowr=0.0;
for(int i=2;i<=n;i++){
if(!check(i)){
nowx=a[i].x,nowy=a[i].y,nowr=0.0;
for(int j=1;j<i;j++){
if(!check(j)){
nowx=(a[i].x+a[j].x)/2.0;
nowy=(a[i].y+a[j].y)/2.0;
nowr=dis(a[i].x,a[i].y,a[j].x,a[j].y)/2.0;
for(int k=1;k<j;k++){
if(!check(k)){
node p1=(a[i]+a[j])/2,p2=(a[j]+a[k])/2,u=a[i]-a[j],v=a[j]-a[k],x;
swap(u.x,u.y);u.y=-u.y;
swap(v.x,v.y);v.y=-v.y;x=inter(p1,u,p2,v);
nowx=x.x,nowy=x.y,nowr=dis(nowx,nowy,a[i].x,a[i].y);
}
}
}
}
}
}
printf("%.10lf\n%.10lf %.10lf",nowr,nowx,nowy);
return 0;
}
[HNOI2012]射箭
设抛物线方程为 \(y=ax^2+bx(a<0,b>0)\)，我们想要求出一组 \(a,b\) 使得它尽可能满足更多的要求。这个显然可以二分答案。
如何 \(\text{check}\) 当前的 \(\text{mid}\) 是否合法呢？每一个限制条件形如 \(y_{i_1}\le ax_i^2+bx_i\le y_{i_2}\) ，也就是 \(\dfrac{y_{i_1}}{x_i}\le ax_i+b\le \dfrac{y_{i_2}}{x_i}\) 。把 \(a,b\) 看成自变量，实际上每个不等式就是一个半平面，我们需要求出半平面交。
如果不是在线插入的话，可以做到 \(O(n\log n)\) 。
我们用有向直线表示半平面，以下默认半平面在有向直线的左侧。
对有向直线按方向向量的极角排序，维护一个双端队列，存储当前构成半平面的直线以及相邻两直线的交点。
每次加入一条有向直线，如果队首/队尾的交点在直线右侧（用叉积判）则弹掉队首/队尾的直线。
为什么这样是对的呢？因为加入直线的单调性，所以要被弹出的直线一定在队首或队尾。
需要注意的细节：
-
加入直线时，先弹队尾，再弹队首。
-
最后还要检查队尾交点是否在队首直线的右侧，如果是也要弹掉。
-
特判平行直线，在右侧的要弹掉。
-
如果题目给出的半平面不一定有限制边界，则应该手动加入一个 \(\text{inf}\) 边界。
点击查看代码
[NOI2011] 智能车比赛
很容易就能够看出，通过两点之间线段最短说明路程一定是折线或线段，如果是折线的话拐点一定是位于矩形的顶点，与此同时拐点位于重合的边上，那么每个矩形产生两个点，这一张联通图最多只有 \((2\times n+2)\) 个点。
令 \(dp[i]\) 表示起点到点 \(i\) 的最短距离，写出状态转移式：\(dp[i]=min(dp[i],dp[j]+len(i,j))\)（线段 \(ij\) 合法）。
在起点固定的情况下，我们的直线行走路径的斜率会被限制在一个区间里，那么直接在枚举下一个点的过程中（或固定终点，枚举上一个点），不断更新合法斜率区间，来判断当前这个转移是否合法。复杂度就是 \(O(n^2)\) 。
点击查看代码
#include<bits/stdc++.h>
#define ERO (1e7)
#define ESP (1e-12)
#define XX(x) ((x)*(x))
#define fabs(x) ((x)>0?(x):-(x))
#define ET0(x) (fabs(x)<=ESP)
using namespace std;
struct Point{
double x,y;
Point operator + (Point p){
p.x+=x;p.y+=y;return p;
}
Point operator - (Point p){
p.x=-p.x+x;p.y=-p.y+y;
return p;
}
inline Point operator *(double lambda){
return {x*lambda,y*lambda};
}
inline double operator *(Point p){
return x*p.x+y*p.y;
}
inline double operator %(Point p){
return x*p.y-y*p.x;
}
inline bool operator ==(Point p){
return ET0(x-p.x)&&ET0(y-p.y);
}
};
struct Line{
Point p1,p2;
};
inline double Dist(Point p1,Point p2){
return sqrt(XX(p1.x-p2.x)+XX(p1.y-p2.y));
}
Point Intersection(Line l1,Line l2){
Point re;
double c1x=l1.p2.x-l1.p1.x,c2x=l2.p2.x-l2.p1.x,c1y=l1.p2.y-l1.p1.y,c2y=l2.p2.y-l2.p1.y;
re.x=(c1y*c2x*l1.p1.x-c1x*c2y*l2.p1.x+c1x*c2x*(l2.p1.y-l1.p1.y))/(c1y*c2x-c1x*c2y);
re.y=(c1x*c2y*l1.p1.y-c1y*c2x*l2.p1.y+c1y*c2y*(l2.p1.x-l1.p1.x))/(c1x*c2y-c1y*c2x);
return re;
}
Line Door[2018];
Point path[2018],Rect[2018][2];
double sx,sy,tx,ty,speed,lastmindis=0,mindis=1e9;
int n;
int main(){
srand(time(NULL));
scanf("%d",&n);
scanf("%lf%lf%lf%lf",&Rect[1][0].x,&Rect[1][0].y,&Rect[1][1].x,&Rect[1][1].y);
for(int i=2;i<=n;i++){
scanf("%lf%lf%lf%lf",&Rect[i][0].x,&Rect[i][0].y,&Rect[i][1].x,&Rect[i][1].y);
Door[i-1].p1={Rect[i][0].x,max(Rect[i-1][0].y,Rect[i][0].y)},
Door[i-1].p2={Rect[i][0].x,min(Rect[i-1][1].y,Rect[i][1].y)};
}
scanf("%lf%lf%lf%lf%lf",&sx,&sy,&tx,&ty,&speed);
if(sx>tx)swap(sx,tx),swap(sy,ty);
int Pos1=1,Pos2=n-1;
while(Pos1<=n-1&&Door[Pos1].p1.x<=sx)++Pos1;
if(Door[Pos1-1].p1.x==sx&&(Door[Pos1-1].p1.y>sy||Door[Pos1-1].p2.y<sy)&&Rect[Pos1-2][0].y<=sy&&sy<=Rect[Pos1-2][1].y)--Pos1;
while(Pos2>=1&&Door[Pos2].p1.x>=tx)--Pos2;
if(Door[Pos2+1].p1.x==tx&&(Door[Pos2+1].p1.y>ty||Door[Pos2+1].p2.y<ty)&&Rect[Pos2+2][0].y<=ty&&ty<=Rect[Pos2+2][1].y)++Pos2;
path[Pos1-1]={sx,sy};
path[Pos2+1]={tx,ty};
for(int i=Pos1;i<=Pos2;i++){
path[i]=(Door[i].p1+Door[i].p2)*0.5;
}
while(lastmindis!=mindis){
lastmindis=mindis;
mindis=0;
for(int i=Pos1;i<=Pos2;i++){
Line l;l.p1=path[i-1];l.p2=path[i+1];
Point p=Intersection(l,Door[i]);
if(p.y>=Door[i].p1.y&&p.y<=Door[i].p2.y) path[i]=p;
else {
if((path[i]-path[i-1])%(path[i+1]-path[i])>0)path[i]=Door[i].p2;
else path[i]=Door[i].p1;
}
mindis+=Dist(path[i],path[i-1]);
}
mindis+=Dist(path[Pos2],path[Pos2+1]);
}
printf("%.8lf\n",mindis/speed);
return 0;
}
[ZJOI2008]无序运动
首先不考虑翻转，对于两个等长的序列，如果任意两条相邻边的边长比以及夹角都相等，那么就匹配。
为了避免实数运算，边长比可以上下平方，然后约分。夹角可以用叉积和点积的最简比值来表示，注意上下符号都要保留。
然后将这些信息离散化，转化成数字串匹配问题，建出用 \(\text{map}\) 表存边的 \(\text{AC}\) 自动机后即可解决。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
int n,m,k,tot=0,fail[200005],last[200005],cnt[200005],ans[200005],tp[200005];
struct point{
int x,y;
point(int _x=0,int _y=0){
x=_x;y=_y;
}
inline point operator -(const point &rhs){
return point(x-rhs.x,y-rhs.y);
}
inline int operator *(const point &rhs){
return x*rhs.y-y*rhs.x;
}
inline int operator ^(const point &rhs){
return x*rhs.x+y*rhs.y;
}
}po[200005];
struct node{
int a,b,c,d;
bool operator <(const node &rhs)const{
if(a!=rhs.a)return a<rhs.a;
if(b!=rhs.b)return b<rhs.b;
if(c!=rhs.c)return c<rhs.c;
return d<rhs.d;
}
}s[200005];
int gcd(int x,int y){
return y?gcd(y,x%y):x;
}
inline node change(point a,point b,point c){
node an;
an.a=(b.x-a.x)*(b.x-a.x)+(b.y-a.y)*(b.y-a.y);
an.b=(c.x-b.x)*(c.x-b.x)+(c.y-b.y)*(c.y-b.y);
an.c=(b-a)*(c-b),an.d=(b-a)^(c-b);
int t=gcd(an.a,an.b);
if(t!=0)an.a/=t,an.b/=t;
t=gcd(abs(an.c),abs(an.d));
if(t!=0)an.c/=t,an.d/=t;
return an;
}
map<node,int> nex[200005];
vector<int> ed[200005];
inline void ins(int len,int id){
int now=0;
for(int i=1;i<=len;i++){
auto it=nex[now].find(s[i]);
if(it==nex[now].end())now=nex[now][s[i]]=++tot;
else now=it->second;
}
ed[now].push_back(id);
}
inline void build(){
queue<int> que;
for(auto it:nex[0])que.push(it.second);
while(!que.empty()){
int now=que.front();
que.pop();
for(auto it:nex[now]){
node a=it.first;
int b=it.second,c=fail[now];
while(c&&nex[c].find(a)==nex[c].end())c=fail[c];
if(nex[c].find(a)!=nex[c].end())c=nex[c][a];
fail[b]=c;
last[b]=ed[c].empty()?last[c]:c;
que.push(b);
}
}
}
inline void solve(){
int now=0;
for(int i=1;i<=n-2;i++){
int x=now;
while(x&&nex[x].find(s[i])==nex[x].end())x=fail[x];
if(nex[x].find(s[i])!=nex[x].end())x=nex[x][s[i]];
now=x;
for(;x;x=last[x])cnt[x]++;
}
}
int main(){
scanf("%d%d",&n,&m);
for(int i=1;i<=m;i++){
int falg=1;
scanf("%d",&k);
for(int j=1;j<=k;j++)scanf("%d%d",&po[j].x,&po[j].y);
if(k<=2){
ans[i]=n-1;
continue;
}
for(int j=2;j<k;j++){
s[j-1]=change(po[j-1],po[j],po[j+1]);
if(s[j-1].c)falg=0;
}
if(falg)tp[i]=1;
ins(k-2,i);
}
build();
for(int i=1;i<=n;i++)scanf("%d %d",&po[i].x,&po[i].y);
for(int i=2;i<n;i++)s[i-1]=change(po[i-1],po[i],po[i+1]);
solve();
for(int i=1;i<=n;i++)po[i].x=-po[i].x;
for(int i=2;i<n;i++)s[i-1]=change(po[i-1],po[i],po[i+1]);
solve();
for(int i=1;i<=tot;i++){
for(int j=0;j<ed[i].size();j++)ans[ed[i][j]]+=cnt[i];
}
for(int i=1;i<=m;i++)printf("%d\n",ans[i]/(tp[i]+1));
return 0;
}
「雅礼集训 2017 Day10」拍苍蝇
首先可以用一个矩形去套这个多边形，那么我们只要枚举这个矩形的左下角就可以枚举完所有多边形的位置了
我们先对每一个 \(x\) 坐标开一个 \(\text{bitset}\)，表示这个 \(x\) 坐标里哪些 \(y\) 坐标处有苍蝇。然后再处理出矩形中哪些位置会被覆盖，这个同样可以枚举 \(x\) 坐标，然后对于所有线段，如果它穿过这个 \(x\) 坐标，就用一个\(\text{stack}\) 存起来，然后把所有 \(\text{stack}\) 里的 \(\text{sort}\) 一下，乱搞就好了。
点击查看代码
