若$R=0$,那麼顯然答案爲離原點最遠的點到原點的距離。算法
不然若全部點都在原點,那麼顯然答案爲$R$。dom
不然考慮二分答案$mid$,檢查$mid$是否可行。ui
那麼每一個點根據對應圓交,能夠覆蓋圓上的一部分,每一個可行方案均可以經過平移使得恰好卡住某個交點。spa
枚舉每一個交點,算出圓上$n$個位置的座標,而後匈牙利算法判斷是否存在完美匹配,時間複雜度$O(n^4\log w)$,不能承受。blog
注意到這個圖是個稠密圖,因此能夠用bitset對匈牙利進行加速,作到$O(\frac{n^3}{32})$每次匹配。it
另外一方面,能夠先枚舉一個點$x$,而後再二分答案$mid$,判斷是否有可行方案使得$x$恰好匹配$x$和圓的交點。io
在這裏,顯然只須要在以前答案$ans$的基礎之上往下二分,若是$ans-eps$不可行那麼就沒有繼續二分的必要。class
即:設$f[x]$表示$x$獲得的最優解,若$f[x]$不是$f[1,x]$的最小值,那麼就沒有繼續二分的必要。基礎
考慮將讀入的$n$個點隨機打亂,那麼$f[x]$是$f[1,x]$的最小值的機率爲$\frac{1}{x}$,一共只有指望$O(\log n)$個$x$有二分的必要。im
檢查次數驟降爲$O(n+\log n\log w)$,時間複雜度$O(\frac{(n+\log n\log w)n^3}{32})$。
注意要特判$mid$太小或者過大致使$x$與圓沒有交點的狀況。
#include<cstdio> #include<algorithm> #include<cstdlib> #include<cmath> using namespace std; typedef unsigned int U; const int N=205,M=7; const double eps=1e-9,pi=acos(-1.0); int n,m,R,i,j;double lim,ans,rot[N][2],b[N][2]; int f[N];U v[M],g[N][M]; struct P{int x,y;}a[N]; inline int sgn(double x){ if(x>eps)return 1; if(x<-eps)return -1; return 0; } bool find(int x){ for(int i=0;i<=m;i++){ U t=v[i]&g[x][i]; while(t){ int j=i<<5|__builtin_ctz(t); v[i]^=1U<<(j&31); if(f[j]<0||find(f[j]))return f[j]=x,1; t-=t&-t; } } return 0; } inline bool check(int A,int B,double C){ double d=sqrt(a[A].x*a[A].x+a[A].y*a[A].y); double l=(a[A].x*a[A].x+a[A].y*a[A].y+C*C-R*R)/(2*d); double h=sqrt(max(C*C-l*l,0.0)); double bx=-a[A].x/d,by=-a[A].y/d; double px=a[A].x+bx*l,py=a[A].y+by*l; by=-by; swap(bx,by); bx*=h,by*=h; if(B==0)px+=bx,py+=by;else px-=bx,py-=by; int i,j; b[0][0]=px,b[0][1]=py; for(i=1;i<n;i++){ b[i][0]=px*rot[i][1]-py*rot[i][0]; b[i][1]=px*rot[i][0]+py*rot[i][1]; } C*=C; for(i=0;i<n;i++){ for(j=0;j<=m;j++)g[i][j]=0; for(j=0;j<n;j++)if(sgn((a[i].x-b[j][0])*(a[i].x-b[j][0])+(a[i].y-b[j][1])*(a[i].y-b[j][1])-C)<=0)g[i][j>>5]|=1U<<(j&31); } for(i=0;i<n;i++)f[i]=-1; for(i=0;i<n;i++){ for(j=0;j<=m;j++)v[j]=~0U; if(!find(i))return 0; } return 1; } int main(){ scanf("%d%d",&n,&R); m=(n-1)>>5; for(i=0;i<n;i++)scanf("%d%d",&a[i].x,&a[i].y); if(!R){ int ans=0; for(i=0;i<n;i++)ans=max(ans,a[i].x*a[i].x+a[i].y*a[i].y); double ret=sqrt(ans); return printf("%.15f",ret),0; } random_shuffle(a,a+n); for(i=1;i<n;i++){ double o=pi*2*i/n; rot[i][0]=sin(o),rot[i][1]=cos(o); } for(i=0;i<n;i++){ int t=a[i].x*a[i].x+a[i].y*a[i].y; double val; if(t<=R)val=R-sqrt(t);else val=sqrt(t)-R; lim=max(lim,val); ans=max(ans,sqrt(t)+R); } for(i=0;i<n;i++)if(a[i].x||a[i].y)for(j=0;j<2;j++){ double l=lim,r=max(min(sqrt(a[i].x*a[i].x+a[i].y*a[i].y)+R,ans-eps),lim); if(!check(i,j,r))continue; while(l+eps<r){ double mid=(l+r)/2; if(check(i,j,mid))r=ans=mid;else l=mid; } } return printf("%.15f",ans),0; }