2018-2019 ACM-ICPC, Asia Nanjing Regional Contest D Country Meow（最小覆盖球）

题目链接

题意：

三维空间给出n个点，在空间中求一个点到这些点的最长距离最短，输出距离

思路： 板子题

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const double eps=1e-10;
struct Tpoint{
    double x,y,z;

};
Tpoint operator +(Tpoint a,Tpoint b){
    Tpoint c;
    c.x=a.x+b.x;
    c.y=a.y+b.y;
    c.z=a.z+b.z;
    return c;
}
Tpoint operator -(Tpoint a,Tpoint b){
    Tpoint c;
    c.x=a.x-b.x;
    c.y=a.y-b.y;
    c.z=a.z-b.z;
    return c;
}
Tpoint operator /(Tpoint a,double k){
    Tpoint c;
    c.x=a.x/k;
    c.y=a.y/k;
    c.z=a.z/k;
    return c;
}
Tpoint operator *(Tpoint a,double k){
    Tpoint c;
    c.x=a.x*k;
    c.y=a.y*k;
    c.z=a.z*k;
    return c;
}
int npoint,nouter;
Tpoint pt[200000],outer[4],res;
double radius,tmp;

double dist2(Tpoint a,Tpoint b){
    return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z);
}
double dot(Tpoint a,Tpoint b){
    return a.x*b.x+a.y*b.y+a.z*b.z;
}
void ball(){
    Tpoint q[3];
    double m[3][3],sol[3],l[3],det;
    int i,j;
    res.x=res.y=res.z=radius=0;
    switch(nouter){
    case 1:
        res=outer[0];break;
    case 2:
        res=(outer[0]+outer[1])/2;
        radius=dist2(res,outer[0]);
        break;
    case 3:
        for(i=0;i<2;i++) q[i]=outer[i+1]-outer[0];
        for(i=0;i<2;i++)
            for(j=0;j<2;j++)
                m[i][j]=dot(q[i],q[j])*2;
        for(i=0;i<2;i++) sol[i]=dot(q[i],q[i]);
        if(fabs(det=m[0][0]*m[1][1]-m[0][1]*m[1][0])<eps) return;
        l[0]=(sol[0]*m[1][1]-sol[1]*m[0][1])/det;
        l[1]=(sol[1]*m[0][0]-sol[0]*m[1][0])/det;
        res=outer[0]+q[0]*l[0]+q[1]*l[1];
        radius=dist2(res,outer[0]);
        break;
    case 4:
        for(i=0;i<3;i++) q[i]=outer[i+1]-outer[0],sol[i]=dot(q[i],q[i]);
        for(i=0;i<3;i++)
            for(j=0;j<3;j++)
                m[i][j]=dot(q[i],q[j])*2;
        det=m[0][0]*m[1][1]*m[2][2]
        +m[0][1]*m[1][2]*m[2][0]
        +m[0][2]*m[2][1]*m[1][0]
        -m[0][2]*m[1][1]*m[2][0]
        -m[0][1]*m[1][0]*m[2][2]
        -m[0][0]*m[1][2]*m[2][1];
        if(fabs(det)<eps) return;
        for(j=0;j<3;j++){
            for(i=0;i<3;i++) m[i][j]=sol[i];
            l[j]=(m[0][0]*m[1][1]*m[2][2]
                  +m[0][1]*m[1][2]*m[2][0]
                  +m[0][2]*m[2][1]*m[1][0]
                  -m[0][2]*m[1][1]*m[2][0]
                  -m[0][1]*m[1][0]*m[2][2]
                  -m[0][0]*m[1][2]*m[2][1])/det;
            for(i=0;i<3;i++) m[i][j]=dot(q[i],q[j])*2;
        }
        res=outer[0];
        for(i=0;i<3;i++) res=res+q[i]*l[i];
        radius=dist2(res,outer[0]);
    }
}
void minball(int n){
    ball();
    if(nouter<4)
        for(int i=0;i<n;i++)
            if(dist2(res,pt[i])-radius>eps){
                outer[nouter++]=pt[i];
                minball(i);
                --nouter;
                if(i>0){
                    Tpoint Tt=pt[i];
                    memmove(&pt[1],&pt[0],sizeof(Tpoint)*i);
                    pt[0]=Tt;
                }
            }
}
int main()
{
    cin>>npoint;
    for(int i=0;i<npoint;i++) cin>>pt[i].x>>pt[i].y>>pt[i].z;
    random_shuffle(pt,pt+npoint);
    radius=-1;
    for(int i=0;i<npoint;i++)
        if(dist2(res,pt[i])-radius>eps)
            nouter=1,outer[0]=pt[i],minball(i);
    printf("%.5lf\n",sqrt(radius));
    return 0;
}