BZOJ2961:共点圆 圆的反演+cdq分治

思路:

(现在我还没有AC这道题,先记录一下思路吧T^T)

昨天膜拜了xhr犇的论文,get了一系列新姿势.

总的来说,保证修改相互独立的数据结构问题可以用如下三种方法来水:

1.对时间分治

题目要求:修改独立,允许离线

我们利用对时间分治,将问题转化为先给出若干修改,然后进行若干询问的形式.这样我们就可以先对修改进行预处理,随后回答每个询问,这样就避免了动态修改.

2.二进制分组

题目要求:修改独立,强制在线

我们按照顺序进行操作,同时维护已经完成的修改的二进制分组,当在后面加入一个修改的时候,我们暴力重构为新的分组;当回答询问时,我们在之前的\(O(\log n)\)的分组中分别累加答案.这样,我们就通过一个\(\log\)的代价避免了动态修改.

3.整体二分

这个与这道题无关就先不说了.

 

那么来说这道题.

首先有两种基本思路:

[1]利用圆的反演转化为动态半平面交问题

我们以圆心为反演中心,合适长度为半径进行反演,那么每个圆都被反演成一条不过原点的直线(在圆内区域形成一个半平面).我们对每个询问点也进行反演,不难发现若询问点在所有圆的并集之内,则反演点应该在所有直线的半平面交之内.

那么现在问题转化为:支持两个操作,加入一个半平面,询问一个点在不在之前插入的半平面的半平面交之内.

{1}貌似可以拿平衡树维护,\(O(n\log n)\),我表示跪了.

{2}按照时间分治,在一次分治内,我们处理出前一半操作中插入的半平面的半平面交,随后去判断后一半操作中的询问点是否在这个半平面交中,若不在,则表示这个询问点不在所有的圆中.于是这样我们可以做到\(O(n\log ^2n)\).

(至于如何快速判断一个点在不在一个凸包中,我们可以任取凸包内的一个点,将凸包上的点按照关于这个点的极角序排序,对于询问点,我们二分得到询问点在凸包上的哪两个点之间,随后利用\(O(1)\)判断询问点在不在三角形中.这样我们就做到了\(O(\log n)\)询问.)

{3}利用二进制分组,维护\(O(\log n)\)个分组的半平面交,容易发现重构半平面交的复杂度为\(O(n\log ^2n)\),对于每一个询问的复杂度也为\(\log ^2n\),因此总复杂度为\(O(n\log ^2n)\).

[2]通过等式变换转化为另一个问题

由于题目保证圆通过原点,我们考虑点\((x_0,y_0)\)在圆心为\((ox,oy)\)的圆内部或边界上的条件:

\[(x_0-ox)^2+(y_0-oy)^2\leq{ox^2+oy^2}\]

\[x_0^2+y_0^2\leq{2x_0{ox}+2y_0{oy}}\]

这样相当于是每一个询问是给定一个半平面,然后问之前出现过的所有圆形成的点是否都在半平面之内.并且要支持加点.

{1}利用平衡树做动态凸包,\(O(n\log n)\),同样是给跪了.

{2}按时间分治,对于每一次分治,我们预处理出前一半操作中的点形成的凸包,然后去判断后面一半操作中的询问半平面是否完全包含凸包.

我们没有必要作凸包,只需\(O(n)\)维护上下凸壳即可.对于后面的询问,我们可以直接二分,也可以预先将所有询问排序随后线性扫描.

这样就可以做到\(O(n\log ^2n)\)或是\(O(n\log n)\).

{3}二进制分组,我们将之前的\(n\)次加点操作分为\(O(\log n)\)组,然后对于每一次询问在每一个维护的结构中二分.这个方法支持强制在线,不过时间复杂度为\(O(n\log ^2n)\).

 

我目前写的做法是[1]{2},不过貌似精度死得很惨.反演半径居然会影响答案.

现在用[2]{2}的方法AC了.时间复杂度\(O(n\log n)\).

#include<cstdio>
#include<cstring>
#include<cctype>
#include<iostream>
#include<algorithm>
using namespace std;
 
typedef double f2;
inline f2 sqr(const f2&x){return x*x;}
 
#define N 500010
struct Case{
    int qte,lab;f2 x,y;
    bool operator<(const Case&B)const{
        if(qte!=B.qte)return qte<B.qte;
        if(qte==0)return x<B.x||(x==B.x&&y<B.y);
        if(qte==1)return(-x/y)<(-B.x/B.y);
    }
}q[N],ql[N],qr[N],up[N],down[N];
int cntl,cntr,topup,topdown;
 
struct Point{
    f2 x,y;
    Point(){}
    Point(f2 _x,f2 _y):x(_x),y(_y){}
    Point operator-(const Point&B)const{return Point(B.x-x,B.y-y);}
};
f2 cross(const Point&A,const Point&B){return A.x*B.y-A.y*B.x;}
struct Line{
    f2 A,B,C;
};
inline Line Getline(const Case&c){
    Line l;
    l.A=2*c.x,l.B=2*c.y,l.C=sqr(c.x)+sqr(c.y);
    return l;
}
inline bool onleft(const Line&l,const Point&p){
    return l.A*p.x+l.B*p.y>=l.C;
}
bool OK[N];
 
inline void Solve(int l,int r){
    if(l==r)return;
    int mid=(l+r)>>1;
    register int i,j;
    for(cntl=cntr=0,i=l;i<=r;++i)if(q[i].lab<=mid)ql[++cntl]=q[i];else qr[++cntr]=q[i];
    int id=l;for(i=1;i<=cntl;++i)q[id++]=ql[i];for(i=1;i<=cntr;++i)q[id++]=qr[i];
    int cnt=0;for(i=l;i<=mid;++i)cnt+=(q[i].qte==0);
    if(cnt){
        for(topup=topdown=0,i=l;i<=mid&&q[i].qte==0;++i){
            while(topup>1&&(q[i].y-up[topup].y)*(up[topup].x-up[topup-1].x)>=(up[topup].y-up[topup-1].y)*(q[i].x-up[topup].x))--topup;
            up[++topup]=q[i];
            while(topdown>1&&(q[i].y-down[topdown].y)*(down[topdown].x-down[topdown-1].x)<=(down[topdown].y-down[topdown-1].y)*(q[i].x-down[topdown].x))--topdown;
            down[++topdown]=q[i];
        }
        int ins;
        for(ins=mid+1;ins<=r;++ins)if(q[ins].qte==1)break;
        if(ins<=r){
            i=ins;
            for(j=1;j<=topdown;++j){
                for(;i<=r&&(j==topdown||-q[i].x/q[i].y*(down[j+1].x-down[j].x)<=(down[j+1].y-down[j].y));++i)
                    if(!onleft(Getline(q[i]),Point(down[j].x,down[j].y)))OK[q[i].lab]=0;   
            }
            i=r;
            for(j=1;j<=topup;++j){
                for(;i>=ins&&(j==topup||-q[i].x/q[i].y*(up[j+1].x-up[j].x)>=up[j+1].y-up[j].y);--i)
                    if(!onleft(Getline(q[i]),Point(up[j].x,up[j].y)))OK[q[i].lab]=0;
            }
        }
    }
    Solve(l,mid),Solve(mid+1,r);
}
 
int main(){
    int n;scanf("%d",&n);
    register int i;for(i=1;i<=n;++i)scanf("%d%lf%lf",&q[i].qte,&q[i].x,&q[i].y),q[i].lab=i;
    for(i=1;i<=n;++i)if(q[i].qte==1)OK[i]=1;
    for(i=1;i<=n&&q[i].qte==1;++i)OK[i]=0;
    sort(q+1,q+n+1);
    Solve(1,n);
    for(i=1;i<=n;++i)if(q[i].qte==1)puts(OK[i]?"Yes":"No");
    return 0;
}

(另外此题真的非常卡精度T^T)

Codeforces77E Martian Food 计算几何+圆的反演

思路:

设两个比较小的圆分别是A,B.

我们将最大圆以及其中一个圆A的交点取出,并以这个点为反演中心,以合适的长度为反演半径,进行反演.

容易发现大圆和A反演之后都成了一条直线.而圆B由于与这两个圆都有唯一交点,那么反演之后也一定都有唯一交点,因此只能是夹在两条直线中间.

不妨以反演中心为原点,三个圆的直径为x轴(方向任意)建立平面直角坐标系,显然B反演之后的圆纵坐标为0.

那么如果我们再放上别的圆(假设一直放在上方),由于它与这三个圆都有唯一交点,所以依然是被卡在两条直线中间,而且在刚才反演的那个圆上方.

这些圆的大小相同,因此我们可以\(O(1)\)得到第\(k\)个圆的圆心半径.

接下来只要再反演回来就好了.

#include<cstdio>
#include<cstring>
#include<cctype>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;

typedef double f2;

struct Point{
    f2 x,y;
    Point(){}
    Point(f2 _x,f2 _y):x(_x),y(_y){}
}P1,P2;
f2 sqr(const f2&x){return x*x;}
f2 dis(const Point&A,const Point&B){return sqrt(sqr(A.x-B.x)+sqr(A.y-B.y));}

inline void GetPointOnDiameter(f2 ox,f2 oy,f2 r,f2 k,f2 b){
    f2 deltax=r/sqrt(k*k+1);
    P1=Point(ox-deltax,k*(ox-deltax)+b);
    P2=Point(ox+deltax,k*(ox+deltax)+b);
}

Point Invert(f2 ox,f2 oy,f2 r,Point A){
    Point o(ox,oy);f2 d=dis(o,A),_d=r*r/d;
    f2 x=ox+(A.x-ox)*(_d/d);
    f2 y=oy+(A.y-oy)*(_d/d);
    return Point(x,y);
}

int main(){
    int Testcase;cin>>Testcase;
    f2 R,r,lx,rx,invr,ox,oy,totalr;int k;
    while(Testcase--){
        scanf("%lf%lf%d",&R,&r,&k);
        invr=r/2;
        rx=invr*invr/(2*r);
        lx=invr*invr/(2*R);
        totalr=(rx-lx)/2,ox=(lx+rx)/2,oy=2*k*totalr;
        GetPointOnDiameter(ox,oy,totalr,oy/ox,0);
        printf("%.10lf\n",dis(Invert(0,0,invr,P1),Invert(0,0,invr,P2))/2);
    }return 0;
}