[Voronoi图 || 模拟退火 特判 || 圆的离散化] Ural 1520 Empire strikes back

论文：顾研《浅谈随机化思想在几何问题中的应用》

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

与POJ 1379 Run Away差别不大

以下摘录特判

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<ctime>
using namespace std;

const int N=1005;
const double eps=1e-5;
const double PI=acos(-1.0);

inline int dcmp(double a,double b){
  if (fabs(a-b)<eps) return 0;
  if (a>b) return 1; return -1;
}

inline double sqr(double x){ return x*x; }

inline double ran(double x){
  return x*(rand()%10000)/10000;
}

inline double rnd(){
  return 1.0*(rand()%10000)/10000;
}

struct Point{
  double x,y;
  Point(double x=0,double y=0):x(x),y(y) { }
  void read(){ scanf("%lf%lf",&x,&y); }
  double len(){ return sqrt(x*x+y*y); }
  bool operator < (const Point &B) const { return x==B.x?y<B.y:x<B.x; }
  friend Point operator - (Point A,Point B){ return Point(A.x-B.x,A.y-B.y); }
  friend double Dis(Point A,Point B){ return (A-B).len(); }
}P[N];


double R; int n;
double Ans;

inline double F(Point p){
  double ret=1e130;
  for (int i=1;i<=n;i++)
    ret=min(ret,Dis(p,P[i]));
  return ret;
}

const double Eps=1e-8;
const double beta=.9; 
const int Test=15;

inline double SA(){
  Point sp,np; double theta;
  double ans,ret;
  theta=ran(2*PI); 
  sp=Point(sin(theta)*R,cos(theta)*R); ans=F(sp);
  for (double T=R;T>Eps;T*=beta){
    Point pre=sp;
    for (int j=1;j<=30;j++){
      theta=ran(2*PI);
      np=Point(pre.x+sin(theta)*T,pre.y+cos(theta)*T);
      if (np.len()>R) continue;
      ret=F(np);
      if (ret>ans)
	sp=np,ans=ret;
    }	
  }
  return ans;
}

double best; 
Point tp;

int main(){
  Point tem;
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  srand(time(0));
  scanf("%d%lf",&n,&R);
  for (int i=1;i<=n;i++) P[i].read();
  Ans=0;
  for (int i=1;i<=n;i++)
    if (P[i].len()>eps){
      tem.x=P[i].x/P[i].len()*R;
      tem.y=P[i].y/P[i].len()*R;
      Ans=max(Ans,F(tem));
      tem.x=-tem.x; tem.y=-tem.y;
      Ans=max(Ans,F(tem));
    }
  for (int i=1;i<=n;i++)
    for (int j=1;j<=n;j++)
      if (!(!dcmp(P[i].x,P[j].x) && !dcmp(P[i].y,P[j].y))){
	if (!dcmp(P[i].y,P[j].y)){
	  tem.x=(P[i].x+P[j].x)/2;
	  tem.y=sqrt(R*R-tem.x*tem.x);
	  Ans=max(Ans,F(tem));
	  tem.y=-tem.y;
	  Ans=max(Ans,F(tem));
	}else{
	  double k=(P[i].x-P[j].x)/(P[j].y-P[i].y);
	  double b=(P[i].y+P[j].y)/2-k*(P[i].x+P[j].x)/2;
	  tem.x=(-2*k*b+sqrt(4*k*k*b*b-4*(k*k+1)*(b*b-R*R)))/(2*k*k+2);
	  tem.y=k*tem.x+b;
	  Ans=max(Ans,F(tem));
	  tem.x=(-2*k*b-sqrt(4*k*k*b*b-4*(k*k+1)*(b*b-R*R)))/(2*k*k+2);
	  tem.y=k*tem.x+b;
	  Ans=max(Ans,F(tem));	  
	}
      }
  for (int t=1;t<=Test;t++)
    Ans=max(Ans,SA());
  printf("%.5lf\n",Ans);
  return 0;
}

UPD 另一篇论文《圆的离散化》中讲了另一个方法

#include<cstdio>  
#include<cstdlib>  
#include<algorithm>  
#include<cmath>
#include<cstring>
#define cl(x) memset(x,0,sizeof(x))
#define L first  
#define R second  
using namespace std;  
typedef double ld;  
typedef pair<ld,ld> abcd;  
  
const int N=305;  
  
inline ld sqr(ld x){ return x*x; }

const ld PI=acos(-1.0);  
const ld eps=1e-10;

inline int dcmp(ld a,ld b){  
  if (fabs(a-b)<eps) return 0;  
  return a<b?-1:1;  
}  
  
inline ld K(ld x,ld y){  
  x+=eps;  
  if (x<0) return atan(y/x)+PI;  
  if (y<0) return atan(y/x)+PI*2;  
  return atan(y/x);   
}  

struct Circle{  
  ld x,y,r;  
  friend ld dist(Circle A,Circle B){  
    return sqrt(sqr(A.x-B.x)+sqr(A.y-B.y));  
  }  
  void read(){  
    double _x,_y; scanf("%lf%lf",&_x,&_y); x=_x,y=_y;
  }  
  bool Cover(Circle B){  
    return dcmp(B.r+dist(*this,B),r)<=0;  
  }  
}C[N];  

int n; ld Ra;  
int cover[N];  
abcd S[N*N]; int cnt;  

inline bool check(ld mid){
  ld ix,iy,dd,tmp;
  for (int i=1;i<=n;i++) C[i].r=mid;
  for (int i=1;i<=n;i++){  
    cnt=0; 
    for (int j=1;j<=n && !cover[i];j++){  
      if (i==j || cover[j]) continue;  
      if (C[j].Cover(C[i])) { cover[i]=1; continue; }  
      if (C[i].Cover(C[j])) continue;  
      if (dcmp(dd=dist(C[i],C[j]),C[i].r+C[j].r)<0){  
	++cnt;  
	tmp=K(C[j].x-C[i].x,C[j].y-C[i].y);    
	ix=(sqr(C[i].r)-sqr(C[j].r)+sqr(dd))/2/dd;    
	iy=sqrt(sqr(C[i].r)-ix*ix);  
	S[cnt].L=tmp-K(ix,iy);    
	S[cnt].R=tmp+K(ix,iy);    
	if (S[cnt].L<0) S[cnt].L+=2*PI,S[cnt].R+=2*PI;
      }
    }
    if (!C[0].Cover(C[i])){
      int j=0;
      ++cnt;  dd=dist(C[i],C[j]);
      tmp=K(C[j].x-C[i].x,C[j].y-C[i].y);    
      ix=(sqr(C[i].r)-sqr(C[j].r)+sqr(dd))/2/dd;    
      iy=sqrt(sqr(C[i].r)-ix*ix);  
      S[cnt].L=tmp-K(ix,iy);    
      S[cnt].R=tmp+K(ix,iy);    
      if (S[cnt].L<0) S[cnt].L+=2*PI,S[cnt].R+=2*PI;
      S[cnt].R-=2*PI; swap(S[cnt].L,S[cnt].R);
      if (S[cnt].L<0) S[cnt].L+=2*PI,S[cnt].R+=2*PI;
    }
    if (cover[i]) continue;  
    sort(S+1,S+cnt+1);   
    int r=0;  
    for (int j=1;j<=cnt;j++)   
      if (r==0 || dcmp(S[j].L,S[r].R)>0)  
	S[++r]=S[j];  
      else if (dcmp(S[r].R,S[j].R)<0)  
	S[r].R=S[j].R;
    int flag=1;
    for (int i=1;i<=r && flag;i++)
      if (dcmp(S[i].R-S[i].L,2*PI)>=0)
	flag=0;
    if (flag) return 0;
  }
  cl(cover);
  return 1;
}

int main(){  
  freopen("t.in","r",stdin);  
  freopen("t.out","w",stdout);  
  scanf("%d%lf",&n,&Ra);  
  C[0].r=Ra;
  for (int i=1;i<=n;i++) C[i].read();
  ld iL=0,iR=2*Ra,MID;
  while (iR-iL>1e-6)
    if (check(MID=(iL+iR)/2))
      iR=MID;
    else
      iL=MID;
  printf("%.5lf\n",(double)iR);
  return 0;
}