题意:
三维空间给出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;
}