hdu 3932,3007 最小圆覆盖 随机增量法

题目链接：

http://acm.hdu.edu.cn/showproblem.php?pid=3932

http://acm.hdu.edu.cn/showproblem.php?pid=3007

都是最小圆覆盖的题目，学习一下随机增量法。

盗用别人的讲解，慢慢理解。

最小覆盖圆， 增量法
假设圆O是前i-1个点得最小覆盖圆，加入第i个点，如果在圆内或边上则什么也不做。否,新得到的最小覆盖圆肯定经过第i个点。
然后以第i个点为基础（半径为0），重复以上过程依次加入第j个点，若第j个点在圆外，则最小覆盖圆必经过第j个点。
重复以上步骤（因为最多需要三个点来确定这个最小覆盖圆，所以重复三次）。遍历完所有点之后，所得到的圆就是覆盖所有点得最小圆。

证明可以考虑这么做：
最小圆必定是可以通过不断放大半径，直到所有以任意点为圆心，半径为半径的圆存在交点，此时的半径就是最小圆。所以上述定理可以通过这个思想得到。这个做法复杂度是O(n)的，当加入圆的顺序随机时，因为三点定一圆，所以不在圆内概率是3/i,求出期望可得是O(n)。

hdu 3932代码：

/* ***********************************************Author :rabbitCreated Time :2014/4/20 9:39:31File Name :8.cpp************************************************ */#pragma comment(linker, "/STACK:102400000,102400000")#include <stdio.h>#include <iostream>#include <algorithm>#include <sstream>#include <stdlib.h>#include <string.h>#include <limits.h>#include <string>#include <time.h>#include <math.h>#include <queue>#include <stack>#include <set>#include <map>using namespace std;#define INF 0x3f3f3f3f#define eps 1e-10#define pi acos(-1.0)typedef long long ll;int dcmp(double x){    if(fabs(x)<eps)return 0;    return x>0?1:-1;}struct Point{    double x,y;    Point(double _x=0,double _y=0){        x=_x;y=_y;    }};Point operator + (Point a,Point b){    return Point(a.x+b.x,a.y+b.y);}Point operator - (Point a,Point b){    return Point(a.x-b.x,a.y-b.y);}Point operator * (Point a,double p){    return Point(a.x*p,a.y*p);}Point operator / (Point a,double p){    return Point(a.x/p,a.y/p);}bool operator < (const Point &a,const Point &b){    return a.x<b.x||(a.x==b.x&&a.y<b.y);}bool operator == (const Point &a,const Point &b){    return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;}double Dot(Point a,Point b){    return a.x*b.x+a.y*b.y;}double Length(Point a){    return sqrt(Dot(a,a));}double Angle(Point a,Point b){    return acos(Dot(a,b)/Length(a)/Length(b));}double angle(Point a){    return atan2(a.y,a.x);}double Cross(Point a,Point b){    return a.x*b.y-a.y*b.x;}Point vecunit(Point x){    return x/Length(x);}Point Normal(Point x){    return Point(-x.y,x.x);}Point Rotate(Point a,double rad){    return Point(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));}struct Line{    Point p,v;    double ang;    Line(){}    Line(Point P,Point v):p(P),v(v){        ang=atan2(v.y,v.x);    }    bool operator < (const Line &L) const {        return ang<L.ang;    }    Point point(double a){        return p+(v*a);    }};bool SegmentIntersection(Point a1,Point a2,Point b1,Point b2){    double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1,b2-a1),           c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-b1);    return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;}struct Circle{    Point c;    double r;    Circle(){}    Circle(Point c,double r):c(c),r(r){}    Point point(double a){        return Point(c.x+cos(a)*r,c.y+sin(a)*r);    }};Point p[100010];Circle CircumscribedCircle(Point p1,Point p2,Point p3){    double Bx=p2.x-p1.x,By=p2.y-p1.y;    double Cx=p3.x-p1.x,Cy=p3.y-p1.y;    double D=2*(Bx*Cy-By*Cx);    double cx=(Cy*(Bx*Bx+By*By)-By*(Cx*Cx+Cy*Cy))/D+p1.x;    double cy=(Bx*(Cx*Cx+Cy*Cy)-Cx*(Bx*Bx+By*By))/D+p1.y;    Point p=Point(cx,cy);    return Circle(p,Length(p1-p));}Circle Min_cover_circle(int n){    Circle ans;ans.c=p[0];ans.r=0;    random_shuffle(p,p+n);    for(int i=1;i<n;i++)        if(Length(ans.c-p[i])+eps>ans.r){            ans.c=p[i];ans.r=0;   //第一个点            for(int j=0;j<i;j++)                if(Length(ans.c-p[j])+eps>ans.r){//第二个点                    ans.c=(p[i]+p[j])/2;ans.r=Length(ans.c-p[j]);                    for(int k=0;k<j;k++)                        if(Length(ans.c-p[k])+eps>ans.r)//第三个点                            ans=CircumscribedCircle(p[i],p[j],p[k]);                }        }    return ans;}int main(){      int n;double X,Y;    while(~scanf("%lf%lf%d",&X,&Y,&n)){          for(int i=0;i<n;i++)              scanf("%lf%lf",&p[i].x,&p[i].y);          Circle ans=Min_cover_circle(n);        printf("(%.1f,%.1f).\n%.1f\n",ans.c.x,ans.c.y,ans.r);      }      return 0;  }  

hdu 3007代码：

/* ***********************************************Author :rabbitCreated Time :2014/4/20 9:39:31File Name :8.cpp************************************************ */#pragma comment(linker, "/STACK:102400000,102400000")#include <stdio.h>#include <iostream>#include <algorithm>#include <sstream>#include <stdlib.h>#include <string.h>#include <limits.h>#include <string>#include <time.h>#include <math.h>#include <queue>#include <stack>#include <set>#include <map>using namespace std;#define INF 0x3f3f3f3f#define eps 1e-10#define pi acos(-1.0)typedef long long ll;int dcmp(double x){    if(fabs(x)<eps)return 0;    return x>0?1:-1;}struct Point{    double x,y;    Point(double _x=0,double _y=0){        x=_x;y=_y;    }};Point operator + (Point a,Point b){    return Point(a.x+b.x,a.y+b.y);}Point operator - (Point a,Point b){    return Point(a.x-b.x,a.y-b.y);}Point operator * (Point a,double p){    return Point(a.x*p,a.y*p);}Point operator / (Point a,double p){    return Point(a.x/p,a.y/p);}bool operator < (const Point &a,const Point &b){    return a.x<b.x||(a.x==b.x&&a.y<b.y);}bool operator == (const Point &a,const Point &b){    return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;}double Dot(Point a,Point b){    return a.x*b.x+a.y*b.y;}double Length(Point a){    return sqrt(Dot(a,a));}double Angle(Point a,Point b){    return acos(Dot(a,b)/Length(a)/Length(b));}double angle(Point a){    return atan2(a.y,a.x);}double Cross(Point a,Point b){    return a.x*b.y-a.y*b.x;}Point vecunit(Point x){    return x/Length(x);}Point Normal(Point x){    return Point(-x.y,x.x);}Point Rotate(Point a,double rad){    return Point(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));}struct Line{    Point p,v;    double ang;    Line(){}    Line(Point P,Point v):p(P),v(v){        ang=atan2(v.y,v.x);    }    bool operator < (const Line &L) const {        return ang<L.ang;    }    Point point(double a){        return p+(v*a);    }};bool SegmentIntersection(Point a1,Point a2,Point b1,Point b2){    double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1,b2-a1),           c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-b1);    return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;}struct Circle{    Point c;    double r;    Circle(){}    Circle(Point c,double r):c(c),r(r){}    Point point(double a){        return Point(c.x+cos(a)*r,c.y+sin(a)*r);    }};Point p[100010];Circle CircumscribedCircle(Point p1,Point p2,Point p3){    double Bx=p2.x-p1.x,By=p2.y-p1.y;    double Cx=p3.x-p1.x,Cy=p3.y-p1.y;    double D=2*(Bx*Cy-By*Cx);    double cx=(Cy*(Bx*Bx+By*By)-By*(Cx*Cx+Cy*Cy))/D+p1.x;    double cy=(Bx*(Cx*Cx+Cy*Cy)-Cx*(Bx*Bx+By*By))/D+p1.y;    Point p=Point(cx,cy);    return Circle(p,Length(p1-p));}Circle Min_cover_circle(int n){    Circle ans;ans.c=p[0];ans.r=0;    random_shuffle(p,p+n);    for(int i=1;i<n;i++)        if(Length(ans.c-p[i])+eps>ans.r){            ans.c=p[i];ans.r=0;            for(int j=0;j<i;j++)                if(Length(ans.c-p[j])+eps>ans.r){                    ans.c=(p[i]+p[j])/2;ans.r=Length(ans.c-p[j]);                    for(int k=0;k<j;k++)                        if(Length(ans.c-p[k])+eps>ans.r)                            ans=CircumscribedCircle(p[i],p[j],p[k]);                }        }    return ans;}int main(){      int n;    while(~scanf("%d",&n)&&n){          for(int i=0;i<n;i++)              scanf("%lf%lf",&p[i].x,&p[i].y);          Circle ans=Min_cover_circle(n);        printf("%.2lf %.2lf %.2lf\n",ans.c.x,ans.c.y,ans.r);      }      return 0;  }