• [loj2586]选圆圈


    下面先给出比较简单的KD树的做法——

    根据圆心建一棵KD树,然后模拟题目的过程,考虑搜索一个圆

    剪枝:如果当前圆[与包含该子树内所有圆的最小矩形]都不相交就退出

    然而这样的理论复杂度是$o(n^2)$,所以会被出题人卡了

    但是如果将坐标系旋转45度,即对于$(x,y)$,变为$((x-y)/\sqrt{2},(x+y)/\sqrt{2})$,就不会被卡了

    然后因为可能相切,精度要求高,eps大约要取$10^{-3}$

     1 #include<bits/stdc++.h>
     2 using namespace std;
     3 #define N 3000005
     4 #define k (l+r>>1)
     5 #define s (sqrt(2))
     6 #define eps (1e-3)
     7 #define sqr(x) (1LL*(x)*(x))
     8 int n,t,id[N],ans[N],ls[N],rs[N];
     9 struct P{
    10     double x,y;
    11 }mn[N],mx[N];
    12 struct ji{
    13     int r,id;
    14     P a;
    15     bool operator < (const ji &p){
    16         if (t)return a.x<p.a.x;
    17         return a.y<p.a.y;
    18     }
    19 }a[N];
    20 bool cmp(int x,int y){
    21     return (a[x].r>a[y].r)||(a[x].r==a[y].r)&&(a[x].id<a[y].id);
    22 }
    23 P rotate(P x){
    24     return P{(x.x-x.y)/s,(x.x+x.y)/s};
    25 }
    26 double dis(P x,P y){
    27     return sqr(x.x-y.x)+sqr(x.y-y.y);
    28 }
    29 double calc(int t,P x){
    30     P y;
    31     if (x.x<mn[t].x)y.x=mn[t].x;
    32     else
    33         if (x.x<=mx[t].x)y.x=x.x;
    34         else y.x=mx[t].x;
    35     if (x.y<mn[t].y)y.y=mn[t].y;
    36     else
    37         if (x.y<=mx[t].y)y.y=x.y;
    38         else y.y=mx[t].y;
    39     return dis(x,y);
    40 }
    41 void up(int l,int r){
    42     mn[k].x=min(mn[ls[k]].x,mn[rs[k]].x);
    43     mn[k].y=min(mn[ls[k]].y,mn[rs[k]].y);
    44     mx[k].x=max(mx[ls[k]].x,mx[rs[k]].x);
    45     mx[k].y=max(mx[ls[k]].y,mx[rs[k]].y);
    46     if (ans[a[k].id])return;
    47     mn[k].x=min(mn[k].x,a[k].a.x-a[k].r);
    48     mn[k].y=min(mn[k].y,a[k].a.y-a[k].r);
    49     mx[k].x=max(mx[k].x,a[k].a.x+a[k].r);
    50     mx[k].y=max(mx[k].y,a[k].a.y+a[k].r);
    51 }
    52 int build(int l,int r,int p){
    53     if (l>r)return 0;
    54     t=p;
    55     nth_element(a+l,a+k,a+r+1);
    56     ls[k]=build(l,k-1,p^1);
    57     rs[k]=build(k+1,r,p^1);
    58     up(l,r);
    59     return k;
    60 } 
    61 void query(int l,int r,int x){
    62     if (l>r)return;
    63     if (calc(k,a[x].a)>sqr(a[x].r)+eps)return;
    64     if ((!ans[a[k].id])&&(dis(a[x].a,a[k].a)<=sqr(a[x].r+a[k].r)+eps))ans[a[k].id]=a[x].id;
    65     query(l,k-1,x);
    66     query(k+1,r,x);
    67     up(l,r);
    68 }
    69 int main(){
    70     scanf("%d",&n);
    71     for(int i=1;i<=n;i++){
    72         scanf("%lf%lf%d",&a[i].a.x,&a[i].a.y,&a[i].r);
    73         a[i].a=rotate(a[i].a);
    74         a[i].id=id[i]=i;
    75     }
    76     mn[0].x=mn[0].y=2e9+1;
    77     mx[0].x=mx[0].y=-2e9-1;
    78     build(1,n,0);
    79     sort(id+1,id+n+1,cmp);
    80     for(int i=1;i<=n;i++)
    81         if (!ans[a[id[i]].id])query(1,n,id[i]);
    82     for(int i=1;i<=n;i++)printf("%d ",ans[i]);
    83 }
    View Code

     当然,上面的做法并不是正解,下面给出更巧妙的做法——

    对于圆$c_{i}$,设圆心为$(x_{i},y_{i})$,半径为$r_{i}$,显然$c_{i}$和$c_{j}$相交当且仅当$(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}\le (r_{i}+r_{j})^{2}$

    考虑Subtask4,即所有圆半径都相同的情况(假设半径都为$r$)

    建立新坐标系,原坐标系中$(x,y)$对应新坐标系中$(\lfloor\frac{x}{r}\rfloor,\lfloor\frac{y}{r}\rfloor)$

    记$(x'_{i},y'_{i})$为$(x_{i},y_{i})$在新坐标系中对应的点,即$(\lfloor\frac{x_{i}}{r}\rfloor,\lfloor\frac{y_{i}}{r}\rfloor)$

    此时,考虑两个圆$c_{i}$和$c_{j}$的相交与$(x'_{i},y'_{i})$的关系:

    1.若$|x'_{i}-x'_{j}|>2$或$|y'_{i}-y'_{j}|>2$,则该维在原坐标系中相差大于$2r$,显然这两个圆不交

    2.若$(x'_{i},y'_{i})=(x'_{j},y'_{j})$,则两维坐标在原坐标系中相差都小于等于$r$,显然这两个圆相交

    当搜索$c_{i}$时,根据第1个性质,其仅需要遍历所有$|x'_{i}-x'_{j}|,|y'_{i}-y'_{j}|\le 2$的圆$c_{j}$,不妨先枚举$(x'_{j},y'_{j})$,再枚举所有对应该位置的$j$,这个将所有点按照$(x'_{i},y'_{i})$这个二元组排序即可做到

    (删除可以用类似指针的方式维护,但实际上重复访问已经被删除的圆不影响复杂度)

    分析此时的复杂度,不妨考虑每一个$j$会被搜索到的次数,注意到搜索$c_{i}$后,根据第2个性质,$(x'_{i},y'_{i})$上所有点都会被选,那么之后就不会从$(x'_{i},y'_{i})$搜到$j$了,因此总复杂度为$o(n)$

    另外由于前面的排序,总复杂度为$o(n\log n)$,可以通过

    考虑原题,实际上若$r\ge \max_{i=1}^{n}r_{i}$,上面的第1个性质仍然成立,但第2个性质并不一定成立

    但是,有一个类似的性质:若$(x'_{i},y'_{i})=(x'_{j},y'_{j})$且$r_{i},r_{j}>\frac{r}{2}$,则这两个圆相交

    初始令$r=\max_{i=1}^{n}r_{i}$,若当前$\max_{i=1}^{n}r_{i}\le \frac{r}{2}$(不包括删去的圆),则将$r$变为$\frac{r}{2}$并重新计算$(x'_{i},y'_{i})$

    来分析此时的复杂度:对于每一个$r$,同样可以证明每一个$j$不会重复从同一个$(x'_{i},y'_{i})$访问$r$(若两次访问,由于是同一个$r$,这两个圆半径都大于$\frac{r}{2}$,根据上面的性质应该被删除了),复杂度仍为$o(n\log n)$

    现在,复杂度的瓶颈已经是排序了,暴力sort复杂度为$o(n\log^{2}n)$,我们将$r$的初值改为$2^{30}$,那么每一次都可以在上一次基础上拆分,即单次排序变为$o(n)$,总排序复杂度为$o(n\log n)$

    最终,总复杂度为$o(n\log n)$,但常数很大,反正被卡了QAQ

     1 #include<bits/stdc++.h>
     2 using namespace std;
     3 #define N 300005 
     4 #define ll long long
     5 int n,R,ans[N];
     6 struct Circle{
     7     int x,y,r,id;
     8     int xx()const{
     9         if (x>=0)return x/R;
    10         return (x-R+1)/R;
    11     }
    12     int yy()const{
    13         if (y>=0)return y/R;
    14         return (y-R+1)/R;
    15     }
    16     bool operator < (const Circle &k)const{
    17         return (xx()<k.xx())||(xx()==k.xx())&&(yy()<k.yy());
    18     }
    19 }a[N],b[N];
    20 vector<Circle>v[2];
    21 ll sqr(int k){
    22     return 1LL*k*k;
    23 }
    24 bool cmp1(Circle x,Circle y){
    25     return (x.xx()<y.xx())||(x.xx()==y.xx())&&(x.y<y.y);
    26 }
    27 bool cmp2(Circle x,Circle y){
    28     return (x.r>y.r)||(x.r==y.r)&&(x.id<y.id);
    29 }
    30 bool check(Circle x,Circle y){
    31     return sqr(x.x-y.x)+sqr(x.y-y.y)<=sqr(x.r+y.r);
    32 }
    33 void find(int x,int y,Circle k){
    34     Circle o;
    35     o.x=min(abs(x),(1<<30)/R)*R;
    36     if (x<0)o.x=-o.x;
    37     o.y=min(abs(y),(1<<30)/R)*R;
    38     if (y<0)o.y=-o.y;
    39     int l=lower_bound(a+1,a+n+1,o)-a;
    40     int r=upper_bound(a+1,a+n+1,o)-a-1;
    41     for(int i=l;i<=r;i++)
    42         if ((!ans[a[i].id])&&(check(k,a[i])))ans[a[i].id]=k.id;
    43 }
    44 int main(){
    45     scanf("%d",&n);
    46     for(int i=1;i<=n;i++){
    47         scanf("%d%d%d",&a[i].x,&a[i].y,&a[i].r);
    48         a[i].id=i;
    49     }
    50     R=(1<<30);
    51     sort(a+1,a+n+1,cmp1);
    52     memcpy(b,a,sizeof(b));
    53     sort(b+1,b+n+1,cmp2);
    54     for(int i=1;i<=n;i++)
    55         if (!ans[b[i].id]){
    56             while (b[i].r<=R/2){
    57                 for(int j=1,lst=1;j<=n;j++)
    58                     if ((j==n)||(a[j].xx()<a[j+1].xx())){
    59                         int x=a[lst].xx();
    60                         v[0].clear(),v[1].clear();
    61                         R/=2;
    62                         for(int k=lst;k<=j;k++)v[a[k].xx()-x*2].push_back(a[k]);
    63                         for(int k=0;k<v[0].size();k++)a[lst++]=v[0][k];
    64                         for(int k=0;k<v[1].size();k++)a[lst++]=v[1][k];
    65                         R*=2;
    66                     }
    67                 R/=2;
    68             }
    69             for(int j=-2;j<=2;j++)
    70                 for(int k=-2;k<=2;k++)find(b[i].xx()+j,b[i].yy()+k,b[i]);
    71         }
    72     for(int i=1;i<=n;i++)printf("%d ",ans[i]);
    73 }
    View Code
  • 相关阅读:
    oracle函数 exp(y)
    oracle函数 power(x,y)
    oracle函数 floor(x)
    oracle函数 ceil(x)
    oracle函数 ABS(x)
    简明Python3教程(A Byte of Python 3)
    C#实现窗口最小化到系统托盘
    简明Python3教程 4.安装
    ubuntu
    Javascript 笔记与总结(2-6)var
  • 原文地址:https://www.cnblogs.com/PYWBKTDA/p/11536022.html
Copyright © 2020-2023  润新知