最小圆覆盖 随机增量算法

　最小圆覆盖。神奇的随机算法。当点以随机的顺序加入时期望复杂度是线性的。

　　------------------------------------------------------------------------------------

　　algorithm:

　　A、令Ci表示为前i个点的最小覆盖圆。当加入新点pi时如果pi不在Ci-1里那么pi必定在Ci的边界上。

　　B、再从新考虑这样一个问题，Ci为前i个点最小覆盖圆且p在Ci的的边界上！同理加入新点pi时如果p

　　i不在Ci-1里那么pi必定在Ci的边界上。这时我们就包含了两个点在这个最小圆的边界上。

　　C、再从新考虑这样一个问题，Ci为前i个点最小覆盖圆且有两个确定点再边界上！此时先让

　　O(N)的方法能够判定出最小圆。

　　------------------------------------------------------------------------------------

　　analysis:

　　现在来分析为什么是线性的。

　　C是线性的这是显然的。

　　B<-C的过程中。考虑pi 他在园内的概率为 (i-1)/i 。在圆外的概率为 1/i 所以加入pi的期望复杂度为：(1-i)/i*O(1) +(1/i)*O(i) {前者在园内那么不进入C，只用了O(1)。后者进入C用了O(i)的时间}这样分析出来，复杂度实际上仍旧

　　是线性的。

　　A<-B的过程中。考虑方法相同，这样A<-B仍旧是线性。于是难以置信的最小圆覆盖的复杂度变成了线性的。

　　-------------------------------------------------------------------------------------

　　下面的程序没有先将点随机化，因为数据通常也是随机的= =！

　　1

　　2 #include<iostream>

　　3 #include<cstdio>

　　4 #include<cmath>

　　5 using namespace std;

　　6 struct node{

　　7 double x,y;

　　8 };

　　9 int n;

　　10 node p[200000];

　　11 double r;

　　12 node O;

　　13 double dist(node a,node b)

　　14 {

　　15 return sqrt( (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y) );

　　16 }

　　17 void calc(double a,double b,double c,double d,double e,double f) //给出两条直线ax+by+c=0,dx+ey+f=0 求交点

　　18 { //注意到三角形里两条中垂线不可能平行，所以不会产生除0错误

　　19 O.y=(c*d-f*a)/(b*d-e*a);

　　20 O.x=(c*e-f*b)/(a*e-b*d);

　　21 }

　　22 int main()

　　23 {

　　24 freopen("HYOJ1337.in","r",stdin);

　　25 freopen("HYOJ1337.out","w",stdout);

　　26 scanf("%d",&n);

　　27 for (int i=1;i<=n;++i) scanf("%lf%lf",&p[i].x,&p[i].y);

　　28 O=p[1];r=0; //初始C1

　　29

　　30 for (int i=2;i<=n;++i) //A

　　31 if (dist(O,p[i])>r+1e-6)

　　32 {

　　33 O=p[i];r=0;

　　34 for (int j=1;j<=i-1;++j) //B

　　35 if (dist(O,p[j])>r+1e-6)

　　36 {

　　37 O.x=(p[i].x+p[j].x)/2;O.y=(p[i].y+p[j].y)/2;r=dist(O,p[j]);

　　38 for (int k=1;k<=j-1;++k) //C

　　39 if (dist(O,p[k])>r+1e-6)

　　40 {

　　41 calc(p[j].x-p[i].x,p[j].y-p[i].y,(p[j].x*p[j].x+p[j].y*p[j].y-p[i].x*p[i].x-p[i].y*p[i].y)/2,

　　42 p[k].x-p[i].x,p[k].y-p[i].y,(p[k].x*p[k].x+p[k].y*p[k].y-p[i].x*p[i].x-p[i].y*p[i].y)/2);

　　43 r=dist(O,p[k]);

　　44 }

　　45 }

　　46 }

　　47 printf("%.3lf/n",r);

　　48 return 0;

　　49 }