[圆的离散化 几何] 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;
}