2018南京Gym - 101981D - Country Meow 【最小球覆盖】

https://codeforces.com/gym/101981/attachments

红书模板。

#include "bits/stdc++.h"

using namespace std;
const double esp = 1e-8;

int npoint, nouter;

struct Tpoint {
    double x, y, z;
};

Tpoint pt[200000], outer[4], res;
double radius, tmp;

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 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]) < esp)
                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) < esp) 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]);
    }
}

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

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

int main() {
    scanf("%d", &npoint);
    for (int i = 0; i < npoint; i++) {
        scanf("%lf%lf%lf", &pt[i].x, &pt[i].y, &pt[i].z);
    }
    double ans = smallest_bool();
    printf("%.12f\n", ans);
    return 0;
}