题目
给出 $N(1 leq N leq 100)$ 个点的坐标 $x_i,y_i,z_i$（$-100000 leq x_i,y_i,z_i leq 100000$），求包围全部点的最小的球。
分析
方法一：模拟退火
模拟退火是 解决最小球覆盖的经典方法，效果也非常好。
随机得到球的中心，如果更小的半径或设定的概率，则转移。（详细解释见链接）
//这个代码严格说不是模拟退火
有一个事实：最小球的球心，它不然是一个确定的点，就是距它最远的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