POJ - 2069 Super Star 最小球覆盖裸题（模拟退火）

传送门

---

题目大意

给出三维坐标系的若干个点，在坐标系范围内找到一个点作为圆心，在给出的点中距离它最远的点的距离作为半径作出一个球，该球内包含了其余所有的给定点。

解题思路

本题就是最小球覆盖的裸题

解法一

这题的精度卡的贼死，我调了一个小时模拟退火的参数还是过不了，但是网上普遍流传的一种过法目前看不太懂为什么这么写，参考了$ACdreamer$的博客：

#include <cstdio>
#include <math.h>
#include <iostream>
#include <algorithm>

using namespace std;

const double delta = 0.98;
const double eps = 1e-8;
const double dinf = 1e300;

struct Point {
    double x, y, z;
} p[2020];

int n;

int dcmp(double d) {
    if (fabs(d) < 0) return 0;
    return d > 0 ? 1 : -1;
}

double dis(Point a, Point b) {
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z));
}

double search() {
    Point s = p[0];
    double T = 100.0, ans = dinf;
    while (T > eps) {
        int k = 1;
        for (int i = 1; i <= n; i++) {
            if (dcmp(dis(s, p[i]) - dis(s, p[k])) > 0)
                k = i;
        }
        double d = dis(s, p[k]);
        ans = min(ans, d);
        s.x += (p[k].x - s.x) / d * T;
        s.y += (p[k].y - s.y) / d * T;
        s.z += (p[k].z - s.z) / d * T;
        T *= delta;
    }
    return ans;
}

int main() {
    while (~scanf("%d", &n) && n) {
        for (int i = 1; i <= n; i++) {
            scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);
        }
        printf("%.5lf\n", search());
    }
    return 0;
}

思路二

红书《ACM国际大学生程序设计竞赛——算法与实现》（俞勇，清华大学出版社）上面的模板：

此代码也没有搞懂，但是准确度极高，跑2018南京站那道题也完全OK。只是对于本题$eps$需要调的高一些才能过。

#include <cstdio>
#include <math.h>
#include <iostream>
#include <algorithm>
#include <time.h>

using namespace std;
const double eps = 1e-14;

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

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

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

void minbool(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;
                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 > eps) {
            nouter = 1;
            outer[0] = pt[i];
            minbool(i);
        }
    }
    return sqrt(radius);
}

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

顺便附上我调了很久的模拟退火（WA）：

#include <cstdio>
#include <math.h>
#include <iostream>
#include <algorithm>

using namespace std;

const double T = 45;
const double delta = 0.999;
const double eps = 1e-14;
const double dinf = 1e300;

struct Point {
    double x, y, z;

    Point(double a = 0, double b = 0, double c = 0) : x(a), y(b), z(c) {}
} p[105];

int n;
Point s;

double dis(Point a, Point b) {
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z));
}

int f(Point cur) {
    int ans = 0;
    for (int i = 1; i < n; i++) {
        if (dis(cur, p[i]) > dis(cur, p[ans])) {
            ans = i;
        }
    }
    return ans;
}

double ans;

void SA() {
    double t = T;
    Point tmp = s, nxt;
    while (t > eps) {
        nxt.x = tmp.x + (rand() * 2 - RAND_MAX) * t;
        nxt.y = tmp.y + (rand() * 2 - RAND_MAX) * t;
        nxt.z = tmp.z + (rand() * 2 - RAND_MAX) * t;
        int k = f(nxt);
        double res = dis(nxt, p[k]);
        double del = res - ans;
        if (del < 0) {
            ans = res;
            s = tmp = nxt;
        } else if (exp(-del / t) * RAND_MAX > rand()) {
            tmp = nxt;
        }
        t *= delta;
    }
}

int main() {
    while (~scanf("%d", &n) && n) {
        srand(23333333);
        srand(rand()), srand(rand());
        double sumx = 0, sumy = 0, sumz = 0;
        for (int i = 0; i < n; i++) {
            scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);
            sumx += p[i].x, sumy += p[i].y, sumz += p[i].z;
        }
        s = Point(sumx / n, sumy / n, sumz / n);
        ans = dinf;
        SA();
        printf("%.5lf\n", ans);
    }
    return 0;
}

其他可行但是被卡的解法

三分套三分套三分，做过2018南京那道最小球覆盖应该比较好找到这个解法，那题的通过代码，我改了三分次数，太小会$WA$，稍大一点就$TLE$，还是这题时间卡的太紧了。

TLE代码

//
// Created by Happig on 2020/11/11
//
#include <cstdio>
#include <math.h>
#include <iostream>
#include <algorithm>

using namespace std;
#define fi first
#define se second
#define pb push_back
#define ins insert
#define Vector Point
#define ENDL "\n"
#define lowbit(x) (x&(-x))
#define mkp(x, y) make_pair(x,y)
#define mem(a, x) memset(a,x,sizeof a);
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
typedef pair<double, double> pdd;
const double eps = 1e-8;
const double pi = acos(-1.0);
const int inf = 0x3f3f3f3f;
const double dinf = 1e300;
const ll INF = 1e18;
const int Mod = 1e9 + 7;
const int maxn = 2e5 + 10;

struct Point {
    double x, y, z;

    Point(double a = 0, double b = 0, double c = 0) : x(a), y(b), z(c) {}
} p[105];

const int T = 30;
const double L = 0;
const double R = 100;
double lx, rx, ly, ry, lz, rz;
double ans;
int n;

double dis(Point a, Point b) {
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z));
}

double f(double x, double y, double z) {
    double res = 0;
    Point s(x, y, z);
    for (int i = 0; i < n; i++) {
        res = max(res, dis(s, p[i]));
    }
    return res;
}

double tri_z(double x, double y) {
    lz = L, rz = R;
    double midl, midr;
    for (int i = 1; i <= T; i++) {
        midl = lz + (rz - lz) / 3;
        midr = rz - (rz - lz) / 3;
        if (f(x, y, midl) > f(x, y, midr)) {
            lz = midl;
        } else rz = midr;
    }
    return f(x, y, lz);
}

double tri_y(double x) {
    ly = L, ry = R;
    double midl, midr;
    for (int i = 1; i <= T; i++) {
        midl = ly + (ry - ly) / 3;
        midr = ry - (ry - ly) / 3;
        if (tri_z(x, midl) > tri_z(x, midr)) {
            ly = midl;
        } else ry = midr;
    }
    return f(x, ly, lz);
}

int main() {
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    //ios_base::sync_with_stdio(0),cin.tie(0),cout.tie(0);
    while (~scanf("%d", &n) && n) {
        for (int i = 0; i < n; i++) {
            scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);
        }
        lx = L, rx = R;
        double midl, midr;
        for (int i = 1; i <= T; i++) {
            midl = lx + (rx - lx) / 3;
            midr = rx - (rx - lx) / 3;
            if (tri_y(midl) > tri_y(midr)) {
                lx = midl;
            } else rx = midr;
        }
        printf("%.5lf\n", f(lx, ly, lz));
    }
    return 0;
}