最小球覆盖——模拟退火&&三分套三分套三分

题目

给出 $N(1 leq N leq 100)$ 个点的坐标 $x_i,y_i,z_i$（$-100000 leq x_i,y_i,z_i leq 100000$），求包围全部点的最小的球。

2018南京区域赛D题

分析

方法一：模拟退火

模拟退火是 解决最小球覆盖的经典方法，效果也非常好。

随机得到球的中心，如果更小的半径或设定的概率，则转移。（详细解释见链接）

//这个代码严格说不是模拟退火

有一个事实：最小球的球心，它不然是一个确定的点，就是距它最远的4个点且等距

于是，我们任选一个点作为初始球心，再在点集中找到距它最远的点，我们让球心靠近最远的点，同时使步长逐渐变小，不断重复此过程，就可以让球心到达正确的位置

#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <queue>
#include <stack>
#include <map>
#include <cmath>
#include <set>
#define LL long long
using namespace std;
const int MAXN = 50;
const double eps = 1e-6;
struct Point
{
    double x, y, z;
    Point(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) { }
};
Point operator - (Point A, Point B)
{
    return Point(A.x - B.x, A.y - B.y, A.z - B.z);
}
double dis(Point A, Point B)
{
    double x = A.x - B.x;
    double y = A.y - B.y;
    double z = A.z - B.z;
    return sqrt(x * x + y * y + z * z);
}
int n;
Point p[MAXN];
double SA(Point start)
{
    double delta = 100.0;
    double ans = 1e20;
    while(delta > eps)
    {
        int d = 0;
        for(int i=0;i<n;i++)
        {
            if(dis(p[i], start) > dis(p[d], start))
                d = i;
        }
        double r = dis(start, p[d]);
        ans = min(ans, r);
        start.x += (p[d].x - start.x) / r * delta;
        start.y += (p[d].y - start.y) / r * delta;
        start.z += (p[d].z - start.z) / r * delta;
        delta *= 0.98;
    }
    return ans;
}
int main()
{
    int T, kcase = 1;
    //scanf("%d", &T);
    while(scanf("%d", &n) && n)
    {
        //scanf("%d", &n);
        for(int i=0;i<n;i++)
            scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);
        printf("%.8f
", SA(Point(0,0,0)));
    }
    return 0;
}

 这有一个严格按模拟退火的定义写的（精度较低，但更好理解

#include<bits/stdc++.h>
#define rep(i, n) for(int i = 1, i##_end_ = (n); i <= i##_end_; ++i)
using namespace std;
typedef pair<int, int> pii;
typedef long long ll;


const double eps = 1e-12;
int sgn(double x) {
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}
struct Point {
    double x, y, z;
    Point(double xp=0, double yp=0, double zp=0): x(xp), y(yp), z(zp) { }
    Point operator + (const Point& rhs) const { return Point(x+rhs.x, y+rhs.y, z+rhs.z); }
    Point operator - (const Point& rhs) const { return Point(x-rhs.x, y-rhs.y, z-rhs.z); }
    Point operator * (const double& k) const { return Point(x*k, y*k, z*k); }
    Point operator / (const double& k) const { return Point(x/k, y/k, z/k); }
    bool operator < (const Point& rhs) const { return x < rhs.x || (x==rhs.x && y<rhs.y) || (x==rhs.x&&y==rhs.y&&z<rhs.z); }
    bool operator == (const Point& rhs) const {return sgn(x - rhs.x) == 0 && sgn(y - rhs.y) == 0 && sgn(z-rhs.z)==0; }
    void scan() { scanf("%lf%lf%lf", &x, &y, &z); }
};
typedef Point Vector;

double dot(Vector x, Vector y) { return x.x*y.x + x.y*y.y + x.z*y.z; }
double length(Vector x) { return sqrt(dot(x, x)); }
double dist2(Point A, Point B) { return dot(A - B, A - B); }
// Circle
struct Circle {
    Point o;
    double r;
    Circle(Point O, double R): o(O), r(R) { }
};


mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());

double Eval(const vector<Point>& pt, Point o) {
    double res = 0;
    for(auto g : pt) res = max(res, dist2(g, o));
    return res;
}
uniform_real_distribution<double> rgen(0.0, 1.0);
double Rand(){ return rgen(rng); }

Circle MinCircleAnneal(const vector<Point>& pt, double T, double dec, double ed) {
    Point pcur, pbest, pnew;

    int sz = pt.size();
    for(auto g : pt) pcur = pcur + g;
    pbest = pcur = pcur / sz;

    double vcur = Eval(pt, pcur), vnew, vbest = vcur;

    while(T > ed) {
        pnew = pcur + Point((Rand()*2.0-1) * T, (Rand()*2.0-1.0) * T, (Rand()*2.0-1) * T);
        vnew = Eval(pt, pnew);
        if(vnew <= vbest) vbest = vcur = vnew, pbest = pcur = pnew;
        if(vnew <= vcur || Rand() < exp(-(vnew-vcur)/T))
            vcur = vnew, pcur = pnew;
        T *= dec;
    }

    return Circle(pbest, sqrt(vbest));
}

int n;
int main() {
    scanf("%d", &n);
    vector<Point> p(n);
    for(int i = 0; i < n; ++i) p[i].scan();
    double ans = 1e13;
    rep(i, 40) {
        Circle cir = MinCircleAnneal(p, 100000.0, 0.999, 3e-7);
        ans = min(ans, cir.r);
    }
    printf("%.10f
", ans);

    return 0;
}

2、三分套三分套三分

代码也很好理解，直接看代码吧

/*
    三分求球的最小覆盖
*/
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long ll;
const ll mod=9999973;
const double eps=1e-4;
int n;
double x[105],y[105],z[105];
double dis3(double a,double b,double c){
    double ans=0;
    for(int i=1;i<=n;i++){
        ans=max(ans,(x[i]-a)*(x[i]-a)+(y[i]-b)*(y[i]-b)+(z[i]-c)*(z[i]-c));
    }
    return ans;
}
double dis2(double a,double b){
    double l=-100000;
    double r=100000;
    double ans=0;
    while(r-l>=eps){
        double rmid=(r+l)/2;
        double lmid=(l+rmid)/2;
        if(dis3(a,b,lmid)<dis3(a,b,rmid)){
            r=rmid;
        }
        else l=lmid;
    }
    return dis3(a,b,l);
}
double dis(double a){
    double l=-100000;
    double r=100000;
    while(r-l>=eps){
        double rmid=(r+l)/2;
        double lmid=(l+rmid)/2;
        if(dis2(a,lmid)<dis2(a,rmid)){
            r=rmid;
        }
        else l=lmid;
    }
    return dis2(a,l);
}
int main(){
   // int n;
    scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%lf%lf%lf",&x[i],&y[i],&z[i]);
    double l=-100000;
    double r=100000;
    while(r-l>=eps){
        double rmid=(r+l)/2;
        double lmid=(l+rmid)/2;
        if(dis(lmid)<dis(rmid)){
            r=rmid;
        }
        else l=lmid;
    }
    printf("%.8f
",sqrt(dis(l)));
    return 0;
}

不管用上面哪种方法，都需要结合题目的精度要求，调节参数，达到精度和时间的平衡。

1. http://www.voidcn.com/article/p-ffokglab-uy.html

2. https://www.cnblogs.com/MekakuCityActor/p/10613934.html

3. https://www.cnblogs.com/heisenberg-/p/6827790.html