APIO2018 T2,一道神奇的 $\texttt{k-D Tree}$ 搜索剪枝题,或者一道毒瘤的 平衡树 + CDQ 分治 题。
题目链接:Luogu P4631/BZOJ 5465/LibreOJ 2586/APIO2018 T2。
题目
题目描述
在平面上,有 $n$ 个圆,记为 $c_1,c_2,\cdots,c_n$。我们尝试对这些圆运行这个算法:
- 找到这些圆中半径最大的。如果有多个半径最大的圆,选择编号最小的。记为 $c_i$。
- 删除 $c_i$ 及与其有交集的所有圆。两个圆有交集当且仅当平面上存在一个点,这个点同时在这两个圆的圆周上或圆内。(原文直译:如果平面上存在一个点被这两个圆所包含,我们称这两个圆有交集。一个点被一个圆包含,当且仅当它位于圆内或圆周上。)
- 重复上面两个步骤直到所有的圆都被删除。
当 $c_i$ 被删除时,若循环中第一步选择的圆是 $c_j$ ,我们说 $c_i$ 被 $c_j$ 删除。对于每个圆,求出它是被哪一个圆删除的。
数据范围
变量 | 数据范围 |
---|---|
$n$ | $[1,3\times 10^5]$ |
$x,y$ | $[-10^9,10^9]$ |
$r$ | $[1,10^9]$ |
时空限制
题目编号 | 时间限制 | 空间限制 |
---|---|---|
Luogu P4631 | $3\text{s}$ | $1000\text{MiB}$ |
BZOJ 5465 | $?\text{s}$ | $1024\text{MiB}$ |
LibreOJ 2586 | $3\text{s}$ | $1024\text{MiB}$ |
APIO2018 T2 | $3\text{s}$ | $1024\text{MiB}$ |
题解
本题目有多种玄学方法。
方法一
思路
暴力
考虑暴力,肯定要先以半径为第一关键字降序排列,编号为第二关键字升序排列。然后对于每个圆,求出它能删除那些圆,并记录答案,后面只要对没有答案的进行查找即可。最坏的时间复杂度是 $\Theta(n^2)$。
下面考虑优化这个过程。
必要条件的转化
我们不难发现,如果圆 $i$ 与圆 $j$ 相交(或者说接触),那么从它们的 $x$ 轴投影和 $y$ 轴投影可以看出:区间 $[x _ i-r _ i,x _ i+r _ i]$ 与区间 $[x _ j-r _ j,x _ j+r _ j]$ 有交,$y$ 轴同理。
所以,坐标对应的投影区间有交是两圆相交的必要条件。
所以,我们可以把圆的外部抽象出一个外接正方形,如果一个圆与这个外接正方形没有公共部分,必然也与这个圆没有公共部分,而正方形和圆可以用 $\texttt{k-D Tree}$ 来储存。
用 $\texttt{k-D Tree}$ 可以通过洛谷上的数据,但是 $\texttt{k-D Tree}$ 毕竟是暴力数据结构,LibreOJ 和 UOJ 上好像有卡 $\texttt{k-D Tree}$ 的数据,不过旋转一下点就解决了。
代码
#include<bits/stdc++.h> using namespace std; #define reg register #define eps 1e-6 #define getchar()(p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++) static char buf[100000],*p1=buf,*p2=buf; inline int read(void){ reg bool f=false; reg char ch=getchar(); reg int res=0; while(ch<'0'||'9'<ch)f|=(ch=='-'),ch=getchar(); while('0'<=ch&&ch<='9')res=10*res+ch-'0',ch=getchar(); return f?-res:res; } long double sina,cosa; //旋转角度的 sin 值和 cos 值 inline long double sqr(reg long double x){ //求平方和 return x*x; } const int MAXN=300000+5; int n,ans[MAXN]; namespace kD_Tree{ //k-D Tree #define mid ( ( (l) + (r) ) >> 1 ) struct Point{ long double x[2],r; int id; inline void Read(reg int i){ static double X,Y; id=i; X=read(),Y=read(); x[0]=X*cosa-Y*sina; //旋转 x[1]=X*sina+Y*cosa; //旋转 r=read(); //读入半径 return; } inline bool operator<(const Point& a)const{ return r==a.r?id<a.id:r>a.r; } }; int root,tot; int lson[MAXN],rson[MAXN]; long double Min[MAXN][2],Max[MAXN][2]; Point p[MAXN],val[MAXN]; int WD; inline bool cmp(const Point& a,const Point& b){ return a.x[WD]<b.x[WD]; } inline void pushup(reg int k){ //维护信息 reg long double x=val[k].x[0],y=val[k].x[1],R=val[k].r; Min[k][0]=x-R;Max[k][0]=x+R; Min[k][1]=y-R;Max[k][1]=y+R; reg int l=lson[k],r=rson[k]; if(l){ Min[k][0]=min(Min[k][0],Min[l][0]); Max[k][0]=max(Max[k][0],Max[l][0]); Min[k][1]=min(Min[k][1],Min[l][1]); Max[k][1]=max(Max[k][1],Max[l][1]); } if(r){ Min[k][0]=min(Min[k][0],Min[r][0]); Max[k][0]=max(Max[k][0],Max[r][0]); Min[k][1]=min(Min[k][1],Min[r][1]); Max[k][1]=max(Max[k][1],Max[r][1]); } return; } inline int Build(reg int l,reg int r,reg int wd){ //建树 if(l>r) return 0; reg int k=++tot; WD=wd; nth_element(p+l,p+mid,p+r+1,cmp); val[k]=p[mid]; lson[k]=Build(l,mid-1,wd^1); rson[k]=Build(mid+1,r,wd^1); pushup(k); return k; } inline bool check(reg int k,const Point& a){ //有交集 return sqr(a.r+val[k].r)-(sqr(a.x[0]-val[k].x[0])+sqr(a.x[1]-val[k].x[1]))>=-eps; } inline bool far(reg int k,reg const Point& a){ //完全脱离 return (a.x[0]+a.r<Min[k][0])||(a.x[0]-a.r>Max[k][0])||(a.x[1]+a.r<Min[k][1])||(a.x[1]-a.r>Max[k][1]); } inline void Query(reg int k,reg const Point& a){ if(far(k,a)) return; if(!ans[val[k].id]&&check(k,a)) ans[val[k].id]=a.id; if(lson[k]) Query(lson[k],a); if(rson[k]) Query(rson[k],a); return; } } using namespace kD_Tree; int main(void){ n=read(); srand(time(0)); sina=sin((long double)rand()),cosa=sqrt(1-sina*sina); for(reg int i=1;i<=n;++i) p[i].Read(i); root=Build(1,n,0); sort(p+1,p+n+1); for(reg int i=1;i<=n;++i) if(!ans[p[i].id]) Query(root,p[i]); for(reg int i=1;i<=n;++i) printf("%d%c",ans[i],i==n?'\n':' '); return 0; }
方法二
思路
防止被卡的方法
考虑到坐标对应的投影区间有交是两圆相交的必要条件。不难发现,对于我们要查询的一个圆,它的半径必然小于之前的所有圆(我们以半径为降序排过序了)。
而且,前面的圆能够不被删去,它们必定两两不相交,且半径大于圆 $c_i$。
那么我们可以对前面的圆维护扫描线,每个圆和当前的直线 $x=x_0$ 相交两次,可以用括号表示。而且由于这些圆两两不相交,括号相对次序不会变。
由于前面的圆半径都比它大,因此和它有交的圆必然经过 $x=x_i-r_i$ 或 $x=x_i+r_i$ 或 $y=y_i−r_i$ 或 $y=y_i+r_i$。所以我们可以在做扫描线时,查询这四个位置的平衡树上,当前圆的前驱后继。
应对多重询问
有很多个询问,那就加上一个 CDQ 分治 咯。
时间复杂度为 $\Theta(n\log_2^2n)$。
代码
#include<bits/stdc++.h> using namespace std; #define reg register #define INF 0X7F7F7F7F typedef long long ll; #define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++) static char buf[100000],*p1=buf,*p2=buf; inline int read(void){ reg bool f=false; reg char ch=getchar(); reg int res=0; while(ch<'0'||'9'<ch)f|=(ch=='-'),ch=getchar(); while('0'<=ch&&ch<='9')res=10*res+ch-'0',ch=getchar(); return f?-res:res; } const int MAXN=300010; struct Circle{ ll x,y,r; int id; inline void Read(reg int i){ x=read(),y=read(),r=read(); id=i; return; } }; struct Event{ ll t; int op,v; inline Event(reg ll t=0,reg int op=0,reg int v=0):t(t),op(op),v(v){ return; } }; inline bool cmp(const Circle& a,const Circle& b){ return (a.r==b.r)?(a.id<b.id):(a.r>b.r); } inline bool cmp2(const Event& a,const Event& b){ return a.t<b.t; } int n; Circle a[MAXN]; int ans[MAXN]; int final[MAXN]; int b[MAXN]; Event E[2*MAXN]; set<pair<int,int>/**/> S; inline ll sqr(reg ll x){ return x*x; } inline bool inter(reg int x,reg int y){ return sqr(a[x].x-a[y].x)+sqr(a[x].y-a[y].y)<=sqr(a[x].r+a[y].r); } inline void Solve(reg int l,reg int r){ if(l==r){ if(ans[l]==INF){ ans[l]=l; b[l]=1; } return; } reg int mid=(l+r)>>1; Solve(l,mid); reg int cnt=0; for(reg int i=l;i<=mid;++i) if(b[i]){ E[++cnt]=Event(3*(a[i].x-a[i].r)-2,1,i); E[++cnt]=Event(3*(a[i].x+a[i].r),2,i); } for(reg int i=mid+1;i<=r;++i){ E[++cnt]=Event(3*(a[i].x-a[i].r)-1,3,i); E[++cnt]=Event(3*(a[i].x+a[i].r)-1,3,i); } sort(E+1,E+cnt+1,cmp2); for(reg int i=1;i<=cnt;++i) if(E[i].op==1) S.insert(make_pair(a[E[i].v].y,E[i].v)); else if(E[i].op==2) S.erase(make_pair(a[E[i].v].y,E[i].v)); else{ auto it=S.lower_bound(make_pair(a[E[i].v].y,0)); if(it!=S.end()){ int x=it->second; if(inter(x,E[i].v)) ans[E[i].v]=min(ans[E[i].v],x); } if(it!=S.begin()){ --it; int x=it->second; if(inter(x,E[i].v)) ans[E[i].v]=min(ans[E[i].v],x); } } cnt=0; for(reg int i=l;i<=mid;++i) if(b[i]){ E[++cnt]=Event(3*(a[i].y-a[i].r)-2,1,i); E[++cnt]=Event(3*(a[i].y+a[i].r),2,i); } for(reg int i=mid+1;i<=r;++i){ E[++cnt]=Event(3*(a[i].y-a[i].r)-1,3,i); E[++cnt]=Event(3*(a[i].y+a[i].r)-1,3,i); } sort(E+1,E+cnt+1,cmp2); for(reg int i=1;i<=cnt;++i) if(E[i].op==1) S.insert(make_pair(a[E[i].v].x,E[i].v)); else if(E[i].op==2) S.erase(make_pair(a[E[i].v].x,E[i].v)); else{ auto it=S.lower_bound(make_pair(a[E[i].v].x,0)); if(it!=S.end()){ int x=it->second; if(inter(x,E[i].v)) ans[E[i].v]=min(ans[E[i].v],x); } if(it!=S.begin()){ --it; int x=it->second; if(inter(x,E[i].v)) ans[E[i].v]=min(ans[E[i].v],x); } } Solve(mid+1,r); return; } int main(void){ n=read(); ll Minx=INF,Miny=INF; for(reg int i=1;i<=n;++i){ a[i].Read(i); Minx=min(Minx,a[i].x),Miny=min(Miny,a[i].y); } for(reg int i=1;i<=n;++i){ a[i].x-=Minx-1; a[i].y-=Miny-1; } sort(a+1,a+n+1,cmp); memset(ans,0X7F,sizeof(ans)); Solve(1,n); for(reg int i=1;i<=n;++i) final[a[i].id]=a[ans[i]].id; for(reg int i=1;i<=n;++i) printf("%d%c",final[i],i==n?'\n':' '); return 0; }
近期评论