[LOJ6360]复燃「恋之埋火」

Description

古明地恋(koishi)和小石子(koishi)是好朋友。 ​ 旧地狱的空中散布着许多颗小石子。恋恋想找出一个位置，使得这个位置离最远的小石子的距离尽可能小。

需要注意的是，这里的空间可能是高维空间。

“在这幻想乡里，可不能被常理所束缚！”——东风谷早苗

对于100%的数据，n≤20000,m≤5,0≤n\leq 20000,m\leq 5,0\leqn≤20000,m≤5,0≤所有坐标≤104\leq 10^4≤10​4​​

Solution

恋恋世界第一！
毒瘤的出题zjtjz题出的瘤毒
最小圆覆盖，不过是高维的
然而随机增量法的分析还是正确的，问题变成了给定m维空间中的k个点，在这k个点所在的k-1维平面上找一个点使得其到这k个点的距离相等。
设这k个点是x1.x2..xk，不妨令x1是原点，其余的k-1个向量的线性组合就是其所在平面，对于每一个向量其系数是未知数，那么就有k-1个未知数，并且有k-1个方程(到x1和xi的距离相等)，那么就可以高斯消元了。

Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
using namespace std;

typedef double db;

const int N=2e4+5;

int n,m;
db eps=1e-7;

struct P{
    db a[5];
    friend P operator + (P a,P b) {
        P c;
        fo(i,0,m-1) c.a[i]=a.a[i]+b.a[i];
        return c;
    }
    friend P operator - (P a,P b) {
        P c;
        fo(i,0,m-1) c.a[i]=a.a[i]-b.a[i];
        return c;
    }
    friend P operator * (P a,db b) {
        P c;
        fo(i,0,m-1) c.a[i]=a.a[i]*b;
        return c;
    }
}p[N],c[N];

struct C{P c;db r;}an;

db dot(P a,P b) {db res=0;fo(i,0,m-1) res+=a.a[i]*b.a[i];return res;}

db dis(P a,P b) {return sqrt(dot(a-b,a-b));}

db a[10][10];

void gauss(int n) {
    fo(i,1,n) {
        fo(j,i+1,n) 
            if (fabs(a[j][i])>fabs(a[i][i])) 
                fo(k,i,n+1) swap(a[i][k],a[j][k]);
        if (fabs(a[i][i])<eps) continue;
        fo(j,i+1,n) {
            db w=a[j][i]/a[i][i];
            fo(k,i,n+1) a[j][k]-=a[i][k]*w;
        }
    }
    fd(i,n,1) {
        if (fabs(a[i][i])<eps) continue;
        fo(j,1,i-1) {
            db w=a[j][i]/a[i][i];
            fo(k,i,n+1) a[j][k]-=a[i][k]*w;
        }
    }
    fo(i,1,n) 
        if (fabs(a[i][i])<eps) a[i][n+1]=0;
        else a[i][n+1]/=a[i][i];
}

void get_circle(int n) {
    if (n==0) {an.c=p[1],an.r=0;return;}
    if (n==1) {an.c=c[1],an.r=0;return;}
    if (n==2) {
        fo(i,0,m-1) an.c.a[i]=(c[1].a[i]+c[2].a[i])/2;
        an.r=dis(an.c,c[1]);
        return;
    }
    fo(i,2,n) {
        P A=c[i]-c[1];
        fo(j,2,n) {
            P B=c[j]-c[1];
            a[i-1][j-1]=2*dot(A,B);
        }
        a[i-1][n]=dot(A,A);
    }
    gauss(n-1);
    an.c=c[1];
    fo(i,2,n) an.c=an.c+(c[i]-c[1])*a[i-1][n];
    an.r=dis(an.c,c[1]);
}

void solve(int n,int cnt) {
    get_circle(cnt);
    fo(i,1,n)
        if (dis(p[i],an.c)>an.r+eps) {
            c[cnt+1]=p[i];
            solve(i-1,cnt+1);
        }
}

int main() {
    scanf("%d%d",&n,&m);
    fo(i,1,n) fo(j,0,m-1) scanf("%lf",&p[i].a[j]);
    srand(20180817);
    random_shuffle(p+1,p+n+1);
    solve(n,0);
    fo(i,0,m-1) printf("%.6lf ",an.c.a[i]);
    return 0;
}