Gym 101981D 2018南京

最新推荐文章于 2021-07-24 10:58:01 发布

二分抄代码

最新推荐文章于 2021-07-24 10:58:01 发布

阅读量622

点赞数

分类专栏：计算几何三维几何

本文链接：https://blog.csdn.net/liufengwei1/article/details/89303612

版权

计算几何同时被 2 个专栏收录

48 篇文章

订阅专栏

三维几何

1 篇文章

订阅专栏

本文介绍了一种求解三维空间中点集最小球覆盖的算法实现，通过递归方式找到能够覆盖所有点的最小球体，适用于点数不多于100的情况。算法首先初始化球体，然后通过检查每个点是否在当前球体内来逐步调整球体的大小和位置。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

最小球覆盖。。。小红书上有O(n)的板子，但这题n只有100..抄了就过了

似乎三分套三分套三分也能过？

#include<bits/stdc++.h>
#define maxl 110
#define eps 1e-8
struct point
{
  double x,y,z;
  point(double a=0,double b=0,double c=0)
  {
    x=a;y=b;z=c;
  }
};

int npoint,nouter;
point pt[maxl],outer[4],res;
double radius,tmp,ans;

inline double dist(point p1,point p2)
{
  double dx=p1.x-p2.x,dy=p1.y-p2.y,dz=p1.z-p2.z;
  return (dx*dx+dy*dy+dz*dz);
}

inline double dot(point p1,point p2)
{
  return p1.x*p2.x+p1.y*p2.y+p1.z*p2.z;
}

inline void ball()
{
  point 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.x=(outer[0].x+outer[1].x)/2;
      res.y=(outer[0].y+outer[1].y)/2;
      res.z=(outer[0].z+outer[1].z)/2;
      radius=dist(res,outer[0]);
      break;
    case 3:
      for(int i=0;i<2;i++)
	{
	  q[i].x=outer[i+1].x-outer[0].x;
	  q[i].y=outer[i+1].y-outer[0].y;
	  q[i].z=outer[i+1].z-outer[0].z;
	}
      for(int i=0;i<2;i++)
	for(int j=0;j<2;j++)
	  m[i][j]=dot(q[i],q[j])*2;
      for(int 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.x=outer[0].x+q[0].x*L[0]+q[1].x*L[1];
      res.y=outer[0].y+q[0].y*L[0]+q[1].y*L[1];
      res.z=outer[0].z+q[0].z*L[0]+q[1].z*L[1];
      radius=dist(res,outer[0]);
      break;
    case 4:
      for(int i=0;i<3;i++)
	{
	  q[i].x=outer[i+1].x-outer[0].x;
	  q[i].y=outer[i+1].y-outer[0].y;
	  q[i].z=outer[i+1].z-outer[0].z;
	  sol[i]=dot(q[i],q[i]);
	}
      for(int i=0;i<3;i++)
	for(int 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(int j=0;j<3;j++)
	{
	  for(int 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(int i=0;i<3;i++)
	    m[i][j]=dot(q[i],q[j])*2;
	}
      res=outer[0];
      for(int i=0;i<3;i++)
	{
	  res.x+=q[i].x*L[i];
	  res.y+=q[i].y*L[i];
	  res.z+=q[i].z*L[i];
	}
      radius=dist(res,outer[0]);
    }
}

inline void minball(int n)
{
  ball();
  if(nouter<4)
    for(int i=0;i<n;i++)
      if(dist(res,pt[i])-radius>eps)
	{
	  outer[nouter]=pt[i];
	  ++nouter;
	  minball(i);
	  --nouter;
	  if(i>0)
	    {
	      point Tt=pt[i];
	      memmove(&pt[1],&pt[0],sizeof(point)*i);
	      pt[0]=Tt;
	    }
	}
}

inline double smallest_ball()
{
  radius=-1;
  for(int i=0;i<npoint;i++)
    if(dist(res,pt[i])-radius>eps)
      {
	nouter=1;
	outer[0]=pt[i];
	minball(i);
      }
  return sqrt(radius);
}

inline void prework()
{
  for(int i=0;i<npoint;i++)
    scanf("%lf%lf%lf",&pt[i].x,&pt[i].y,&pt[i].z);
}

inline void mainwork()
{
  ans=smallest_ball();
}

inline void print()
{
  printf("%.5f\n",ans);
}

int main()
{
  while(~scanf("%d",&npoint))
    {
      prework();
      mainwork();
      print();
    }
  return 0;
}