【BZOJ1336】Alien（Balkan2002）-最小圆覆盖

测试地址：Alien
做法：本题需要用到最小圆覆盖。
我们知道最小圆一定是点集中某三个点的外接圆，因此我们要找到这三个点。
求最小圆覆盖的算法是这样的：
先枚举一个i$i$，表示现在要求包含前i$i$个点的最小圆。
对于每一个i$i$，如果它不在当前求出的最小圆中，则把当前最小圆修改成点i$i$，然后继续进行下面的操作，否则继续枚举i$i$。
枚举j$j$，表示现在要求包含第i$i$个点，且包含前j$j$个点的最小圆。
对于每一个j$j$，如果它不在当前求出的最小圆中，则把当前最小圆修改成以i,j$i,j$两点连成的线段为直径的圆，然后继续进行下面的操作，否则继续枚举j$j$。
枚举k$k$，表示现在要求包含第i,j$i,j$个点，且包含前k$k$个点的最小圆。
对于每一个k$k$，如果它不在当前求出的最小圆中，则把当前最小圆修改成i,j,k$i,j,k$三点组成的三角形的外接圆，否则继续枚举k$k$。
以上算法看似是O(n3)$O(n^3)$的，但在点排列随机的基础上（注意一开始一定要将点随机排列！），算法实际上是期望O(n)$O(n)$的（因为一个点在外接圆外的概率比较小），于是我们利用一些简单的计算几何技巧就可以解决这一题了。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
int n;
struct point
{
    double x,y;
    int val;
}p[100010];
double nowx,nowy,nowr;

bool cmp(point a,point b)
{
    return a.val<b.val;
}

double dis(double x1,double y1,double x2,double y2)
{
    return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}

bool check(int i)
{
    return dis(p[i].x,p[i].y,nowx,nowy)<=nowr;
}

double multi(point a,point b)
{
    return a.x*b.y-b.x*a.y;
}

point inter(point p,point v,point q,point w)
{
    point ans,u;
    u.x=p.x-q.x,u.y=p.y-q.y;
    double t=multi(w,u)/multi(v,w);
    ans.x=p.x+t*v.x,ans.y=p.y+t*v.y;
    return ans;
}

int main()
{
    srand(19260817);

    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        scanf("%lf%lf",&p[i].x,&p[i].y);
        p[i].val=rand();
    }
    sort(p+1,p+n+1,cmp);

    nowx=p[1].x,nowy=p[1].y,nowr=0.0;
    for(int i=2;i<=n;i++)
        if (!check(i))
        {
            nowx=p[i].x,nowy=p[i].y,nowr=0.0;
            for(int j=1;j<i;j++)
                if (!check(j))
                {
                    nowx=(p[i].x+p[j].x)/2.0;
                    nowy=(p[i].y+p[j].y)/2.0;
                    nowr=dis(p[i].x,p[i].y,p[j].x,p[j].y)/2.0;
                    for(int k=1;k<j;k++)
                        if (!check(k))
                        {
                            point p1,p2,u,v,x;
                            p1.x=(p[i].x+p[j].x)/2.0,p1.y=(p[i].y+p[j].y)/2.0;
                            p2.x=(p[j].x+p[k].x)/2.0,p2.y=(p[j].y+p[k].y)/2.0;
                            u.x=p[i].y-p[j].y,u.y=p[j].x-p[i].x;
                            v.x=p[j].y-p[k].y,v.y=p[k].x-p[j].x;
                            x=inter(p1,u,p2,v);
                            nowx=x.x,nowy=x.y,nowr=dis(nowx,nowy,p[i].x,p[i].y);
                        }
                }
        }

    printf("%.6lf\n%.6lf %.6lf",nowr,nowx,nowy);

    return 0;
}