BZOJ 1336: [Balkan2002]Alien最小圆覆盖 随机增量法

1336: [Balkan2002]Alien最小圆覆盖

Time Limit: 1 Sec  Memory Limit: 162 MBSec  Special Judge
Submit: 1806  Solved: 805
[Submit][Status][Discuss]

Description

给出N个点，让你画一个最小的包含所有点的圆。

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

Sample Output

5.00
5.00 5.00

神奇·随机增量法

注意 1337 要求输出三位小数1336 要求精度极高。。。

求三角形外心是抄的。。。要是考再YY吧。。。

#include<cmath>
#include<ctime>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<complex>
#include<iostream>
#include<algorithm>
#include<iomanip>
#include<vector>
#include<string>
#include<bitset>
#include<map>
#include<set>
using namespace std;

typedef double db;

inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch<='9'&&ch>='0'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
	return x*f;
}
void print(int x)
{if(x<0)putchar('-'),x=-x;if(x>=10)print(x/10);putchar(x%10+'0');}

const int N=100100;
const db eps=1e-10;

struct P{db x,y;}p[N];

inline db dis(const P &x,const P &y){return sqrt((x.x-y.x)*(x.x-y.x)+(x.y-y.y)*(x.y-y.y));}

P geto(P a,P b,P c)
{
	P res;db a1,a2,b1,b2,c1,c2;
	a1=2*(b.x-a.x);b1=2*(b.y-a.y);c1=b.x*b.x-a.x*a.x+b.y*b.y-a.y*a.y;
	a2=2*(c.x-a.x);b2=2*(c.y-a.y);c2=c.x*c.x-a.x*a.x+c.y*c.y-a.y*a.y;
	if(abs(a1)<eps)res.y=c1/b1,res.x=(c2-res.y*b2)/a2;
	else if(fabs(b1)<eps)res.x=c1/a1,res.y=(c2-res.x*a2)/b2;
	else res.x=(c2*b1-c1*b2)/(a2*b1-a1*b2),res.y=(c2*a1-c1*a2)/(b2*a1-b1*a2);
	return res;
}

int main()
{
	int n=read();
	register int i,j,k;
	for(i=1;i<=n;++i)scanf("%lf%lf",&p[i].x,&p[i].y);
	random_shuffle(p+1,p+1+n);
	P o=p[1];db r=0;
	for(i=1;i<=n;++i)
	{
		if(dis(o,p[i])<r||abs(dis(o,p[i])-r)<eps)continue;
		o.x=(p[i].x+p[1].x)/2;o.y=(p[i].y+p[1].y)/2;r=dis(p[i],p[1])/2;
		for(j=2;j<i;++j)
		{
			if(dis(o,p[j])<r||abs(dis(o,p[j])-r)<eps)continue;
			o.x=(p[i].x+p[j].x)/2;o.y=(p[i].y+p[j].y)/2;r=dis(p[i],p[j])/2;
			for(k=1;k<j;++k)
			{
				if(dis(o,p[k])<r||abs(dis(o,p[k])-r)<eps)continue;
				o=geto(p[i],p[j],p[k]);
				r=dis(o,p[i]);
			}
		}
	}
	printf("%.10lf\n%.10lf %.10lf\n",r,o.x,o.y);
	return 0;
}
/*
4
1 0
0 1
0 -1
-1 0

1.000

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

5.00
5.00 5.00
*/