1336: [Balkan2002]Alien最小圆覆盖
Time Limit: 1 Sec Memory Limit: 162 MBSec Special Judge
Submit: 1473 Solved: 648
[Submit][Status][Discuss]
Description
Input
先给出点的个数N,2<=N<=100000,再给出坐标Xi,Yi.(-10000.0<=xi,yi<=10000.0)
Output
Sample Input
6
8.0 9.0
4.0 7.5
1.0 2.0
5.1 8.7
9.0 2.0
4.5 1.0
8.0 9.0
4.0 7.5
1.0 2.0
5.1 8.7
9.0 2.0
4.5 1.0
Sample Output
5.00
5.00 5.00
5.00 5.00
随机增量法求最小圆覆盖。
三重循环。
令ci为前i个点的覆盖圆,新加入一个点i+1时,若其在圆内,跳过,若其在圆外,修改圆心使i+1在圆c(i+1)上。
检查之前的点,令ci为前i个点的覆盖圆,且点j在圆周上,若第i+1个点无法被圆覆盖,修改圆心使点i+1和点j都在圆周上。
检查之前的点,令ci为前i个点的覆盖圆,且点j和点k在圆周上,若第i+1个点无法被圆覆盖,修改圆心使点i+1和点j、点k都在圆周上
这算法倒是还能理解,但是求外心的几何算法表示看不懂。这个技能还是等高二再解锁吧。
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
const double eps=1e-;
const int mxn=;
int n;
struct point{
double x,y;
friend point operator +(const point a,const point b){
return (point){a.x+b.x,a.y+b.y};
}
friend point operator -(const point a,const point b){
return (point){a.x-b.x,a.y-b.y};
}
friend point operator /(const point a,double b){
return (point){a.x/b,a.y/b};
}
}p[mxn];
inline double dis(point a,point b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
} point center(point a,point b,point c){//返回三角形外心
double a1,a2,b1,b2,c1,c2;
point ans;
a1=*(b.x-a.x);b1=*(b.y-a.y);c1=(b.x*b.x)-(a.x*a.x)+(b.y*b.y)-(a.y*a.y);
//c1=(a1*a1+b1*b1)/2
a2=*(c.x-a.x);b2=*(c.y-a.y);c2=(c.x*c.x)-(a.x*a.x)+(c.y*c.y)-(a.y*a.y);
//c2=(a2*a2+b2*b2)/2
if(fabs(a1)<eps){
ans.y=c1/b1;
ans.x=(c2-ans.y*b2)/a2;
}
else if(fabs(b1)<eps){
ans.x=c1/a1;
ans.y=(c2-ans.x*a2)/b2;
}
else{
ans.x=(c2*b1-c1*b2)/(a2*b1-a1*b2);
ans.y=(c2*a1-c1*a2)/(b2*a1-b1*a2);
}
return ans;
}
int main(){
scanf("%d",&n);
int i,j,k;
for(i=;i<=n;i++){
scanf("%lf%lf",&p[i].x,&p[i].y);
}
random_shuffle(p+,p+n+);
point t=p[];
double r=0.0;
for(i=;i<=n;i++)//
if(dis(t,p[i])>r+eps){
t=(p[i]+p[])/;//默认圆心,等待增量
r=dis(p[i],t);//半径
for(j=;j<i;j++)//
if(dis(t,p[j])>r+eps){//若有点在圆外,更新圆心
t=(p[i]+p[j])/;
r=dis(t,p[i]);
for(k=;k<j;k++){//最多三点确定一圆
if(dis(p[k],t)>r+eps){
t=center(p[i],p[j],p[k]);
r=dis(p[i],t);
}
}
}
}
printf("%.10lf\n%.10lf %.10lf",r,t.x,t.y);
return ;
}