[圆的离散化 几何] POJ 1688 Dolphin Pool

调了一天半 就是过不去 心累 先挖个坑吧

论文：高逸涵《与圆有关的离散化方法》

#include<cstdio>#include<cstdlib>#include<algorithm>#include<cmath>#include<cstring>#define cl(x) memset(x,0,sizeof(x))using namespace std;typedef long double ld;typedef pair<ld,ld> abcd;const ld eps=1e-11;inline int dcmp(ld a,ld b){  if (fabs(a-b)<eps) return 0;  //if (a==b) return 0;  if (a<b) return -1; return 1;}inline ld sqr(ld x){ return x*x; }inline bool Pd(abcd A,abcd B){  return !(dcmp(A.second,B.first)<0 || dcmp(B.second,A.first)<0);}const int N=55;namespace Tset{  int fat[N*N*N],rank[N*N*N],ncnt;  inline int New(){    ++ncnt; fat[ncnt]=ncnt; rank[ncnt]=0; return ncnt;  }  inline int Fat(int u){    return u==fat[u]?u:fat[u]=Fat(fat[u]);  }  inline bool Merge(int x,int y){    x=Fat(x); y=Fat(y); if (x==y) return 0;    if (rank[x]>rank[y]) swap(x,y);    if (rank[x]==rank[y]) rank[y]++;    fat[x]=y; return 1;  }  inline void clear(){    ncnt=0; //cl(fat); cl(rank);  }}struct Point{    ld x,y;    Point(ld x=0,ld y=0):x(x),y(y) { }    void read(){      int ix,iy; scanf("%d%d",&ix,&iy); x=ix; y=iy;   }  };struct Circle{  Point c; ld r;  void read(){    int ir; c.read(); scanf("%d",&ir); r=ir;  }  bool cross(ld h,ld &x1,ld &x2){    ld l=fabs(h-c.y);    if (dcmp(l,r)>0) return 0;    ld tem=sqrt(r*r-l*l);    x1=c.x-tem; x2=c.x+tem; return 1;  }  friend bool Cross(Circle C1,Circle C2,Point &p1,Point &p2){      ld r1=C1.r,r2=C2.r;      ld l2=sqr(C1.c.x-C2.c.x)+sqr(C1.c.y-C2.c.y);      if (dcmp(l2,sqr(r1+r2))>0 || dcmp(l2,sqr(r1-r2))<0) return 0;      ld l=sqrt(l2);      ld cs0=(C2.c.x-C1.c.x)/l,sn0=(C2.c.y-C1.c.y)/l;      ld cs1=(r1*r1+l2-r2*r2)*0.5/r1/l,sn1=sqrt(1.0-cs1*cs1);      p1=Point(C1.c.x+r1*(cs0*cs1-sn0*sn1),C1.c.y+r1*(cs0*sn1+cs1*sn0));      p2=Point(C1.c.x+r1*(cs0*cs1+sn0*sn1),C1.c.y+r1*(-cs0*sn1+cs1*sn0));      return 1;    }}C[N];int n,Ans;ld sx[N*N*N]; int icnt;int last[N*N],cur[N*N];abcd down[N*N],up[N*N];int dcnt,ucnt;int tag[N*N*N];inline void Calc(ld h,abcd *lst,int &cnt){  ld t1,t2;  abcd tem[N*N],ret[N*N]; int tcnt=0,rcnt=0;  for (int i=1;i<=n;i++)    if (C[i].cross(h,t1,t2))      tem[++tcnt]=abcd(t1,t2);  sort(tem+1,tem+tcnt+1);  ld l=tem[1].first,r=tem[1].second;  for (int i=2;i<=tcnt;i++)    if (dcmp(tem[i].first,r)>0)      ret[++rcnt]=abcd(l,r),l=tem[i].first,r=tem[i].second;    else      r=max(r,tem[i].second);  if (tcnt) ret[++rcnt]=abcd(l,r);  cnt=rcnt+1;  lst[1].first=-1e10; lst[rcnt+1].second=1e10;  for (int i=1;i<=rcnt;i++)    lst[i].second=ret[i].first,lst[i+1].first=ret[i].second;}  int main(){  const ld Eps=1e-9;  using namespace Tset;  int Q; Point t1,t2;  freopen("t.in","r",stdin);  freopen("t.out","w",stdout);  scanf("%d",&Q);  while (Q--){    scanf("%d",&n);    icnt=0; Ans=0;    for (int i=1;i<=n;i++) C[i].read(),sx[++icnt]=C[i].c.y-C[i].r,sx[++icnt]=C[i].c.y+C[i].r;     for (int i=1;i<=n;i++) for (int j=i+1;j<=n;j++) if (Cross(C[i],C[j],t1,t2)) sx[++icnt]=t1.y,sx[++icnt]=t2.y;    sort(sx+1,sx+icnt+1);    int tem=1;    for (int i=2;i<=icnt;i++) if (sx[i]-sx[i-1]>Eps) sx[++tem]=sx[i];    icnt=tem;    dcnt=0; last[++dcnt]=New();    for (int i=1;i<=icnt;i++){      Calc(sx[i]-Eps,down,dcnt);      Calc(sx[i]+Eps,up,ucnt);      for (int j=1;j<=ucnt;j++){cur[j]=New();for (int k=1;k<=dcnt;k++)  if (Pd(down[k],up[j])){    Merge(cur[j],last[k]);    tag[Fat(last[k])]=i;  }      }      for (int k=1;k<=dcnt;k++)if (tag[Fat(last[k])]!=i)  Ans++;      for (int j=1;j<=ucnt;j++)last[j]=cur[j];    }    printf("%d\n",Ans);    clear(); cl(tag); cl(sx); cl(C);  }  return 0;}