poj2069/2018南京（最小球覆盖）

上交NB！

POJ2069 题目链接：http://poj.org/problem?id=2069

2018南京（gym101981 D） 题目链接：https://codeforces.com/gym/101981

POJ2069 题意：比较裸，给你n个恒星在三维空间的位置，让你找到包含所有恒星的最小球体（每个恒星看成一个点），输出球体的半径。

gym101981 D 题意：给你n个城市在三维空间的位置，现在要建立一个总部，使得总部与n个城市之间的距离尽可能的小，输出总部与最远城市之间的距离，那么这个题就可以转换成最小球覆盖问题。

上交NB！上交的代码两个题都可以过，时间复杂度为，网上POJ2069题解的代码不能过gym的，gym题解的代码不能过POJ2069的，但是上交的代码两个题都可以过！但是需要在精度问题上修改一下。

gym上说的误差不超过1e-3，我们eps开到1e-3,1e-4都可以过。

POJ上说误差不超过0.00001，但是eps得开到1e-11才可以过，？.jpg

代码：

# define _CRT_SECURE_NO_WARNINGS
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <string.h>
#include <vector>
using namespace std;
const double eps = 1e-11;
//poj2069 eps = 1e-11 0.00001
//gym 
//O（n）
struct Tpoint {
  double x, y, z;
  Tpoint() {}
  Tpoint(double x, double y, double z) :x(x), y(y), z(z) {}
};

inline double dist(Tpoint p1, Tpoint 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(Tpoint p1, Tpoint p2) {
  return p1.x * p2.x + p1.y * p2.y + p1.z * p2.z;
}
void ball(Tpoint outer[], Tpoint& res, double& radious, int& nouter) {
  Tpoint q[3];
  double m[3][3], sol[3], L[3], det;
  int i, j;
  res.x = res.y = res.z = radious = 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;
    radious = 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];
    radious = 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];
    }
    radious = dist(res, outer[0]);
  }
}
void minball(vector<Tpoint>pt, Tpoint outer[], int n, Tpoint& res, double& radious, int& nouter) {
  ball(outer, res, radious, nouter);
  if (nouter < 4)
    for (int i = 0; i < n; ++i)
      if (dist(res, pt[i]) - radious > eps) {
        outer[nouter] = pt[i];
        ++nouter;
        minball(pt, outer, i, res, radious, nouter);
        --nouter;
        if (i > 0) {
          Tpoint Tt = pt[i];
          memmove(&pt[1], &pt[0], sizeof(Tpoint) * i);
          pt[0] = Tt;
        }
      }
}

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

signed main() {
  int npoint;//点数
  Tpoint res;//球心坐标
  double radious;//球的半径
  int nouter = 0;
  while (~scanf("%d", &npoint) && npoint) {
    vector<Tpoint>pt;
    Tpoint outer[4];
    for (int i = 0; i < npoint; i++) {
      double x, y, z;
      scanf("%lf%lf%lf", &x, &y, &z);
      pt.push_back(Tpoint(x, y, z));
    }
    printf("%.5f\n", smallest_ball(pt, outer, npoint, res, radious, nouter));
  }
  return 0;
}