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

题目

给出 $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 = ;

const double eps = 1e-;

struct Point

{

    double x, y, z;

    Point(double x = , double y = , double z = ) : 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 = ;

        for(int i=;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 = ;

    //scanf("%d", &T);

    while(scanf("%d", &n) && n)

    {

        //scanf("%d", &n);

        for(int i=;i<n;i++)

            scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z);

        printf("%.8f\n", SA(Point(,,)));

    }

    return ;

}

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

#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-;

int sgn(double x) {

    if(fabs(x) < eps) return ;

    return x <  ? - : ;

}

struct Point {

    double x, y, z;

    Point(double xp=, double yp=, double zp=): 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) ==  && sgn(y - rhs.y) ==  && sgn(z-rhs.z)==; }

    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 = ;

    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-) * T, (Rand()*2.0-1.0) * T, (Rand()*2.0-) * 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 = ; i < n; ++i) p[i].scan();

    double ans = 1e13;

    rep(i, ) {

        Circle cir = MinCircleAnneal(p, 100000.0, 0.999, 3e-);

        ans = min(ans, cir.r);

    }

    printf("%.10f\n", ans);

    return ;

}

2、三分套三分套三分

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

/*

    三分求球的最小覆盖

*/

#include<iostream>

#include<cstdio>

#include<cstring>

#include<algorithm>

#include<cmath>

using namespace std;

typedef long long ll;

const ll mod=;

const double eps=1e-;

int n;

double x[],y[],z[];

double dis3(double a,double b,double c){

    double ans=;

    for(int i=;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=-;

    double r=;

    double ans=;

    while(r-l>=eps){

        double rmid=(r+l)/;

        double lmid=(l+rmid)/;

        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=-;

    double r=;

    while(r-l>=eps){

        double rmid=(r+l)/;

        double lmid=(l+rmid)/;

        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=;i<=n;i++)scanf("%lf%lf%lf",&x[i],&y[i],&z[i]);

    double l=-;

    double r=;

    while(r-l>=eps){

        double rmid=(r+l)/;

        double lmid=(l+rmid)/;

        if(dis(lmid)<dis(rmid)){

            r=rmid;

        }

        else l=lmid;

    }

    printf("%.8f\n",sqrt(dis(l)));

    return ;

}

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

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

最小球覆盖——模拟退火&&三分套三分套三分的相关教程结束。