[HDU3007]Buried memory（最小圆覆盖）

题目描述

传送门

题解

三个点能确定一个圆，答案一定为某一个三角形的外接圆
有一种非常神奇的复杂度为$O(n)$的随机增量法解决这个问题
算法流程

• 先对点集random_shuffle处理一下
• 假设已经有了前i-1个点的最小圆覆盖$C_{i-1}$，检验点$p_i$是否在$C_{i-1}$内，如果在，则$C_i=C_{i-1}$，否则，$p_i$一定在圆$C_i$上，进行下一步
• 包括点$p_i$在内的前j-1个点(j < <script type="math/tex" id="MathJax-Element-141"><</script>i)的最小圆覆盖为$C_{j-1}$，检验点$p_j$是否在$C_{j-1}$内，如果在，则$C_j=C_{j-1}$，否则，$p_j$一定在圆$C_j$上，进行下一步
• 包括点$p_i$和点$p_j$在内的前k-1个点(k < <script type="math/tex" id="MathJax-Element-150"><</script>j)的最小圆覆盖为$C_{k-1}$，检验点$p_k$是否在$C_{k-1}$内，如果在，则$C_k=C_{k-1}$，否则，$p_k$一定在圆$C_k$上。返回

时间复杂度证明
目测复杂度是$O(n^3)$的，实际上是$O(n)$的

• 首先最后一层循环是$O(j)$的
• 考虑倒数第二层循环，这个点有$j\over 3$的概率不在圆内，这个时候要进行一个复杂度为$O(j)$的循环，有$2j\over 3$的概率不进行这个循环，所以期望是$O(1)$，总体复杂度$O(n)$
• 同理可得，第一层循环的复杂度是$O(n)$的

代码

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define N 505

const double pi=acos(-1.0);
const double eps=1e-9;
int dcmp(double x)
{
    if (x<=eps&&x>=-eps) return 0;
    return (x>0)?1:-1;
}
struct Vector
{
    double x,y;
    Vector(double X=0,double Y=0)
    {
        x=X,y=Y;
    }
};
typedef Vector Point;
Vector operator + (Vector a,Vector b) {return Vector(a.x+b.x,a.y+b.y);}
Vector operator - (Vector a,Vector b) {return Vector(a.x-b.x,a.y-b.y);}
Vector operator * (Vector a,double b) {return Vector(a.x*b,a.y*b);}
Vector operator / (Vector a,double b) {return Vector(a.x/b,a.y/b);}

int n;
double r;
Point p[N],c;

double Dot(Vector a,Vector b)
{
    return a.x*b.x+a.y*b.y;
}
double Cross(Vector a,Vector b)
{
    return a.x*b.y-a.y*b.x;
}
double Len(Vector a)
{
    return sqrt(Dot(a,a));
}
Vector rotate(Vector a,double rad)
{
    return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
Point GLI(Point P,Vector v,Point Q,Vector w)
{
    Vector u=P-Q;
    double t=Cross(w,u)/Cross(v,w);
    return P+v*t;
}
Point center(Point A,Point B,Point C)
{
    Vector v,w;
    v=B-A,w=C-A;
    v=rotate(v,pi/2),w=rotate(w,pi/2);
    return GLI((A+B)/2,v,(A+C)/2,w);
}
void mcc()
{
    random_shuffle(p+1,p+n+1);
    c=p[1],r=0;
    for (int i=2;i<=n;++i)
        if (dcmp(Len(c-p[i])-r)>0)
        {
            c=p[i],r=0;
            for (int j=1;j<i;++j)
                if (dcmp(Len(c-p[j])-r)>0)
                {
                    c=(p[i]+p[j])/2;
                    r=Len(c-p[i]);
                    for (int k=1;k<j;++k)
                        if (dcmp(Len(c-p[k])-r)>0)
                        {
                            c=center(p[i],p[j],p[k]);
                            r=Len(c-p[i]);
                        }
                }
        }
}
int main()
{
    while (~scanf("%d",&n))
    {
        if (!n) break;
        for (int i=1;i<=n;++i)
            scanf("%lf%lf",&p[i].x,&p[i].y);
        mcc();
        printf("%.2lf %.2lf %.2lf\n",c.x,c.y,r);
    }
    return 0;
}