【圆的反演变换+cdq分治】共点圆

http://www.tsinsen.com/ViewGProblem.page?gpid=A1381

因为所有的圆都过原点，因此可以以原点为反演中心，将所有的圆反演为一条直线，然后根据反演的性质，可以发现圆的内部会被反演到一个半平面上去，因此每次询问就相当于询问点是否在半平面的交之中，这个可以动态插半平面询问，也可以cdq分治来离线回答；询问的时候二分点在凸包的哪块三角形中，然后叉积一下就可以判断了。

我采用了第二种方法，然后开始了艰苦卓绝地卡常数的历程，历经各种常数优化，终于把最后两个点都卡在了3s整...

#include <cstdio>#include <cstdlib>#include <cstring>#include <iostream>#include <algorithm>#include <cmath>#include <ctime>#define sqr(x) ((x)*(x))const double pi=acos(-1.0);const double eps=1e-8;const double oo=1e15,LR=100000.0;using namespace std;struct point{    double x,y,z,_,d;    point () {}    point (double X,double Y) {x=X,y=Y,z=1;}    point operator -(const point &p) const {        point e;        e.x=x-p.x,e.y=y-p.y,e.z=1;        return e;    }    point operator +(const point &p) const {        point e;        e.x=x+p.x,e.y=y+p.y,e.z=1;        return e;        }        point operator *(const double a) const {            point e;            e.x=x*a,e.y=y*a,e.z=1;            return e;        }        void print() {            printf("X=%f Y=%f Z=%f\n",(double)x,(double)y,(double)z);        }        double dist() {return (sqr(x)+sqr(y));}}p[500050],b[500050];int ans[500050],a[500050],n,st[500050],u[500050],c[500050];inline double cr(point e,point r) {return e.x*r.y-e.y*r.x;}inline double dot(point &e,point &r) {return e.x*r.x+e.y*r.y+e.z*r.z;}inline void cross(point &p,point &q,point &e){    e.x=p.y*q.z-p.z*q.y;    e.y=p.z*q.x-p.x*q.z;    e.z=p.x*q.y-p.y*q.x;}inline void ori(point &a){    a._=atan2(a.y,a.x);    a.d=sqr(a.z)/(sqr(a.x)+sqr(a.y));}inline bool cmp(int i,int j){    double tmp=p[i]._-p[j]._;    if (fabs(tmp)>pi/2) return tmp<-eps;    tmp=cr(p[i],p[j]);    if (fabs(tmp)>eps) return tmp>eps;    return p[i].d<p[j].d;}inline bool check(point L,point T,point I){    point p;    cross(L,T,p);    if (dot(p,I)>-eps) return 1;    return 0;}int search(point b[],int tot,point &P,point &Q,point &U){        if (cr(P,U)<-eps || cr(Q,U)>eps) return 0;        int l=2,r=tot-1;point pp;         for (;l<=r;) {            int mid=((l+r)>>1);            pp=b[mid]-b[1];            if (cr(pp,U)>eps) l=mid+1;            else r=mid-1;        }        return l-1;}point OB;bool cmp2(point i,point j){    double tmp=i._-j._;        if (fabs(tmp)>pi/2) return tmp<-eps;    tmp=cr(i-OB,j-OB);    if (fabs(tmp)>eps) return tmp>eps;    return i.d<j.d;}void modify(point b[],int tot,int a[],int sum) {    for (int i=2;i<=tot;i++)        if (b[i].y<b[1].y-eps || (fabs(b[i].y-b[1].y)<eps && b[i].x<b[1].x-eps)) swap(b[1],b[i]);    for (int i=2;i<=tot;i++)         b[i]._=atan2(b[i].y-b[1].y,b[i].x-b[1].x);    OB=b[1];    sort(b+2,b+tot+1,cmp2);    b[++tot]=b[1];    point P=b[2]-b[1],Q=b[tot-1]-b[1];    for (int i=1;i<=sum;i++) {    point U=p[a[i]]-b[1];        int k=search(b,tot,P,Q,U);        if (!k) ans[a[i]]=1;        point e,r;        e=b[k+1]-b[k],r=p[a[i]]-b[k];        if (cr(e,r)<-eps) ans[a[i]]=1;    }}void inter(int u[],int tot,int a[],int sum){    for (int i=n+1;i<=n+4;i++) u[++tot]=i;    sort(u+1,u+tot+1,cmp);    int h,r;    st[h=r=1]=u[1];     for (int i=2;i<=tot;i++) {        if (fabs(cr(p[u[i]],p[u[i-1]]))<eps) continue;        for (;(h<r) && (!check(p[st[r-1]],p[st[r]],p[u[i]]));r--) ;        for (;(h<r) && (!check(p[st[h]],p[st[h+1]],p[u[i]]));h++) ;        st[++r]=u[i];    }    for (;(h<r) && (!check(p[st[r-1]],p[st[r]],p[st[h]]));r--) ;    st[r+1]=st[h];    int cnt=0;    for (int i=h;i<=r;i++) {        point e;        cross(p[st[i]],p[st[i+1]],e);        e.x/=e.z,e.y/=e.z,e.z=1;        b[++cnt]=e;    }    modify(b,cnt,a,sum);}void calc(int l,int r){    if (l==r) return ;    int mid=(l+r)>>1;    calc(l,mid);    calc(mid+1,r);    int tot=0,sum=0;    for (int i=l;i<=mid;i++)        if (!c[i]) u[++tot]=i;    for (int i=mid+1;i<=r;i++)        if (c[i] && !ans[i]) a[++sum]=i;    if (sum && tot) inter(u,tot,a,sum);}void origin(){    int tot=n;    ++tot,p[tot].x=1,p[tot].y=0,p[tot].z=oo;    ++tot,p[tot].x=0,p[tot].y=1,p[tot].z=oo;    ++tot,p[tot].x=-1,p[tot].y=0,p[tot].z=oo;    ++tot,p[tot].x=0,p[tot].y=-1,p[tot].z=oo;    for (int i=1;i<=tot;i++) ori(p[i]);    }void point_reverse(point &p){        double len=LR/p.dist();        p.x*=len,p.y*=len;}void circle_reverse(double x,double y,point &P){        point p,q,O(x,y);        p.x=x+x,p.y=y+y,p.z=1;        q.x=x-y,q.y=y+x,q.z=1;        point_reverse(p),point_reverse(q);        cross(p,q,P);        point_reverse(O);        if (dot(P,O)<-eps) {            P.x=-P.x,P.y=-P.y,P.z=-P.z;        }}int main(){    scanf("%d",&n);    for (int i=1;i<=n;i++) {        double x,y;        scanf("%d%lf%lf",&c[i],&x,&y);        if (!c[i]) circle_reverse(x,y,p[i]);        else {            p[i].x=x,p[i].y=y;            point_reverse(p[i]);        }    }    origin();    calc(1,n);    for (int i=1;i<=n;i++)        if (c[i]) ans[i]=1;        else break;    for (int i=1;i<=n;i++)        if (c[i]) {            if (ans[i]) printf("No\n");            else printf("Yes\n");        }    return 0;}

updata：把交点存一下会快一点

#include <cstdio>#include <cstdlib>#include <cstring>#include <iostream>#include <algorithm>#include <cmath>#include <ctime>#define sqr(x) ((x)*(x))const double pi=acos(-1.0);const double eps=1e-8;const double oo=1e15,LR=100000.0;using namespace std;struct point{    double x,y,z,_,d;    point () {}    point (double X,double Y) {x=X,y=Y,z=1;}    point operator -(const point &p) const {        point e;        e.x=x-p.x,e.y=y-p.y,e.z=1;        return e;    }    point operator +(const point &p) const {        point e;        e.x=x+p.x,e.y=y+p.y,e.z=1;        return e;        }        point operator *(const double a) const {            point e;            e.x=x*a,e.y=y*a,e.z=1;            return e;        }        void print() {            printf("X=%f Y=%f Z=%f\n",(double)x,(double)y,(double)z);        }        double dist() {return (sqr(x)+sqr(y));}}p[500050],b[500050],St[500050];int ans[500050],a[500050],n,st[500050],u[500050],c[500050];inline double cr(point e,point r) {return e.x*r.y-e.y*r.x;}inline double dot(point &e,point &r) {return e.x*r.x+e.y*r.y+e.z*r.z;}inline void cross(point &p,point &q,point &e){    e.x=p.y*q.z-p.z*q.y;    e.y=p.z*q.x-p.x*q.z;    e.z=p.x*q.y-p.y*q.x;}inline void ori(point &a){    a._=atan2(a.y,a.x);    a.d=sqr(a.z)/(sqr(a.x)+sqr(a.y));}inline bool cmp(int i,int j){    double tmp=p[i]._-p[j]._;    if (fabs(tmp)>pi/2) return tmp<-eps;    tmp=cr(p[i],p[j]);    if (fabs(tmp)>eps) return tmp>eps;    return p[i].d<p[j].d;}inline bool check(point &p,point I){    if (dot(p,I)>eps) return 1;    return 0;}inline bool equ(point &p) {return fabs(p.x)<eps && fabs(p.y)<eps && fabs(p.z)<eps;}int search(point b[],int tot,point &P,point &Q,point &U){        if (cr(P,U)<-eps || cr(Q,U)>eps) return 0;        int l=2,r=tot-1;point pp;         for (;l<=r;) {            int mid=((l+r)>>1);            pp=b[mid]-b[1];            if (cr(pp,U)>eps) l=mid+1;            else r=mid-1;        }        return l-1;}point OB;bool cmp2(point i,point j){    double tmp=i._-j._;        if (fabs(tmp)>pi/2) return tmp<-eps;    tmp=cr(i-OB,j-OB);    if (fabs(tmp)>eps) return tmp>eps;    return i.d<j.d;}void modify(point b[],int tot,int a[],int sum) {    for (int i=2;i<=tot;i++)        if (b[i].y<b[1].y-eps || (fabs(b[i].y-b[1].y)<eps && b[i].x<b[1].x-eps)) swap(b[1],b[i]);    for (int i=2;i<=tot;i++)         b[i]._=atan2(b[i].y-b[1].y,b[i].x-b[1].x);    OB=b[1];    sort(b+2,b+tot+1,cmp2);    b[++tot]=b[1];    point P=b[2]-b[1],Q=b[tot-1]-b[1];    for (int i=1;i<=sum;i++) {    point U=p[a[i]]-b[1];        int k=search(b,tot,P,Q,U);        if (!k) ans[a[i]]=1;        point e,r;        e=b[k+1]-b[k],r=p[a[i]]-b[k];        if (cr(e,r)<-eps) ans[a[i]]=1;    }}void inter(int u[],int tot,int a[],int sum){    for (int i=n+1;i<=n+4;i++) u[++tot]=i;    sort(u+1,u+tot+1,cmp);    int h,r;   st[h=r=1]=u[1];    for (int i=2;i<=tot;i++) {       if (fabs(cr(p[u[i]],p[u[i-1]]))<eps) continue;       for (;(h<r) && (!check(St[r],p[u[i]]));r--) ;       for (;(h<r) && (!check(St[h+1],p[u[i]]));h++) ;       st[++r]=u[i];       cross(p[st[r-1]],p[st[r]],St[r]);       if (equ(St[r])) return ;   }   for (;(h<r) && (!check(St[r],p[st[h]]));r--) ;   for (;(h<r) && (!check(St[h+1],p[st[r]]));h++) ;   if (r-h+1<=2) return ;    st[r+1]=st[h];    int cnt=0;    for (int i=h;i<=r;i++) {        point e;        cross(p[st[i]],p[st[i+1]],e);        e.x/=e.z,e.y/=e.z,e.z=1;        b[++cnt]=e;    }    modify(b,cnt,a,sum);}void calc(int l,int r){    if (l==r) return ;    int mid=(l+r)>>1;    calc(l,mid);    calc(mid+1,r);    int tot=0,sum=0;    for (int i=l;i<=mid;i++)        if (!c[i]) u[++tot]=i;    for (int i=mid+1;i<=r;i++)        if (c[i] && !ans[i]) a[++sum]=i;    if (sum && tot) inter(u,tot,a,sum);}void origin(){    int tot=n;    ++tot,p[tot].x=1,p[tot].y=0,p[tot].z=oo;    ++tot,p[tot].x=0,p[tot].y=1,p[tot].z=oo;    ++tot,p[tot].x=-1,p[tot].y=0,p[tot].z=oo;    ++tot,p[tot].x=0,p[tot].y=-1,p[tot].z=oo;    for (int i=1;i<=tot;i++) ori(p[i]);    }void point_reverse(point &p){        double len=LR/p.dist();        p.x*=len,p.y*=len;}void circle_reverse(double x,double y,point &P){        point p,q,O(x,y);        p.x=x+x,p.y=y+y,p.z=1;        q.x=x-y,q.y=y+x,q.z=1;        point_reverse(p),point_reverse(q);        cross(p,q,P);        point_reverse(O);        if (dot(P,O)<-eps) {            P.x=-P.x,P.y=-P.y,P.z=-P.z;        }}int main(){    scanf("%d",&n);    for (int i=1;i<=n;i++) {        double x,y;        scanf("%d%lf%lf",&c[i],&x,&y);        if (!c[i]) circle_reverse(x,y,p[i]);        else {            p[i].x=x,p[i].y=y;            point_reverse(p[i]);        }    }    origin();    calc(1,n);    for (int i=1;i<=n;i++)        if (c[i]) ans[i]=1;        else break;    for (int i=1;i<=n;i++)        if (c[i]) {            if (ans[i]) printf("No\n");            else printf("Yes\n");        }    return 0;}