题意:给你100个三维空间里的点,让你求一个点,使得他到所有点距离最大的值最小,也就是让你找一个最小的球覆盖掉这n个点
题解:红书模板题,这题也因为数据小,精度也不高,所以也可以用随机算法,模拟退火之类的
1 #include<bits/stdc++.h> 2 using namespace std; 3 int npoint,nouter; 4 struct Tpoint 5 { 6 double x,y,z; 7 }; 8 Tpoint pt[200000],outer[4],res; 9 #define eps 1e-9 10 double radius,tmp; 11 inline double dist(Tpoint p1,Tpoint p2) 12 { 13 double dx=p1.x-p2.x,dy=p1.y-p2.y,dz=p1.z-p2.z; 14 return (dx*dx+dy*dy+dz*dz); 15 } 16 inline double dot(Tpoint p1,Tpoint p2) 17 { 18 return p1.x*p2.x+p1.y*p2.y+p1.z*p2.z; 19 } 20 void ball() 21 { 22 Tpoint q[3]; 23 double m[3][3],sol[3],L[3],det; 24 int i,j; 25 res.x=res.y=res.z=radius=0; 26 switch (nouter) 27 { 28 case 1:res=outer[0];break; 29 case 2: 30 res.x=(outer[0].x+outer[1].x)/2; 31 res.y=(outer[0].y+outer[1].y)/2; 32 res.z=(outer[0].z+outer[1].z)/2; 33 radius=dist(res,outer[0]); 34 break; 35 case 3: 36 for (int i=0;i<2;i++) 37 { 38 q[i].x=outer[i+1].x-outer[0].x; 39 q[i].y=outer[i+1].y-outer[0].y; 40 q[i].z=outer[i+1].z-outer[0].z; 41 } 42 for (int i=0;i<2;i++) 43 for (int j=0;j<2;j++) 44 m[i][j]=dot(q[i],q[j])*2; 45 for (int i=0;i<2;i++) sol[i]=dot(q[i],q[i]); 46 if (fabs(det=m[0][0]*m[1][1]-m[0][1]*m[1][0])<eps) return; 47 L[0]=(sol[0] * m[1][1] - sol[1] * m[0][1])/det; 48 L[1]=(sol[1] * m[0][0] - sol[0] * m[1][0])/det; 49 res.x=outer[0].x+q[0].x*L[0]+q[1].x*L[1]; 50 res.y=outer[0].y+q[0].y*L[0]+q[1].y*L[1]; 51 res.z=outer[0].z+q[0].z*L[0]+q[1].z*L[1]; 52 radius=dist(res,outer[0]); 53 break; 54 case 4: 55 for (int i=0;i<3;i++) 56 { 57 q[i].x=outer[i+1].x-outer[0].x; 58 q[i].y=outer[i+1].y-outer[0].y; 59 q[i].z=outer[i+1].z-outer[0].z; 60 sol[i]=dot(q[i],q[i]); 61 } 62 for (int i=0;i<3;i++) 63 for (int j=0;j<3;j++) m[i][j]=dot(q[i],q[j])*2; 64 det=m[0][0]*m[1][1]*m[2][2] 65 +m[0][1]*m[1][2]*m[2][0] 66 +m[0][2]*m[2][1]*m[1][0] 67 -m[0][1]*m[1][0]*m[2][2] 68 -m[0][2]*m[1][1]*m[2][0] 69 -m[0][0]*m[1][2]*m[2][1]; 70 if (fabs(det)<eps) return ; 71 for (int j=0;j<3;j++) 72 { 73 for (int i=0;i<3;i++) 74 m[i][j]=sol[i]; 75 L[j]=( m[0][0]*m[1][1]*m[2][2] 76 +m[0][1]*m[1][2]*m[2][0] 77 +m[0][2]*m[2][1]*m[1][0] 78 -m[0][2]*m[1][1]*m[2][0] 79 -m[0][1]*m[1][0]*m[2][2] 80 -m[0][0]*m[1][2]*m[2][1] 81 )/det; 82 for (int i=0;i<3;i++) 83 m[i][j]=dot(q[i],q[j])*2; 84 } 85 res=outer[0]; 86 for (int i=0;i<3;i++) 87 { 88 res.x+=q[i].x*L[i]; 89 res.y+=q[i].y*L[i]; 90 res.z+=q[i].z*L[i]; 91 } 92 radius=dist(res,outer[0]); 93 } 94 } 95 void minball(int n) 96 { 97 ball(); 98 if (nouter<4) 99 { 100 for (int i=0;i<n;i++) 101 if (dist(res,pt[i])-radius>eps) 102 { 103 outer[nouter]=pt[i]; 104 ++nouter; 105 minball(i); 106 --nouter; 107 if (i>0) 108 { 109 Tpoint Tt=pt[i]; 110 memmove(&pt[1],&pt[0],sizeof(Tpoint)*i); 111 pt[0]=Tt; 112 } 113 } 114 } 115 } 116 int main() 117 { 118 radius=-1; 119 scanf("%d",&npoint); 120 for (int i=0;i<npoint;i++) 121 { 122 scanf("%lf%lf%lf",&pt[i].x,&pt[i].y,&pt[i].z); 123 } 124 for (int i=0;i<npoint;i++) 125 if (dist(res,pt[i])-radius>eps) 126 { 127 nouter=1; 128 outer[0]=pt[i]; 129 minball(i); 130 } 131 printf("%.10lf ",sqrt(radius)); 132 }