Gym - 101981D The 2018 ICPC Asia Nanjing Regional Contest D.Country Meow 最小球覆盖
题意:给你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\n",sqrt(radius)); 132 }
Anderyi!
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· .NET开发智能桌面机器人:用.NET IoT库编写驱动控制两个屏幕
· 用纯.NET开发并制作一个智能桌面机器人:从.NET IoT入门开始
· 一个超经典 WinForm,WPF 卡死问题的终极反思
· ASP.NET Core - 日志记录系统(二)
· .NET 依赖注入中的 Captive Dependency
· 在外漂泊的这几年总结和感悟,展望未来
· 博客园 & 1Panel 联合终身会员上线
· 支付宝事故这事儿,凭什么又是程序员背锅?有没有可能是这样的...
· https证书一键自动续期,帮你解放90天限制
· 在 ASP.NET Core WebAPI如何实现版本控制?