题目大意
最小球覆盖裸题
模拟退火
这题和 P O J 2069 POJ2069 POJ2069最大的区别是,这题的精度卡不紧,使用一般的模拟退火就能写过
#include <cstdio>
#include <math.h>
#include <iostream>
#include <algorithm>
using namespace std;
const double T = 100000;
const double delta = 0.998;
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 SA() {
double t = T, ans = dinf;
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;
}
return ans;
}
int main() {
while (~scanf("%d", &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);
printf("%.8lf\n", SA());
}
return 0;
}
这份代码是红书上的模板,也可以过 P O J 2069 POJ2069 POJ2069
#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;
int n, 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 < n; i++) {
if (dist(res, pt[i]) - radius > eps) {
nouter = 1;
outer[0] = pt[i];
minbool(i);
}
}
return sqrt(radius);
}
void solve() {
scanf("%d", &n);
for (int i = 0; i < n; i++) {
scanf("%lf%lf%lf", &pt[i].x, &pt[i].y, &pt[i].z);
}
printf("%.8lf\n", smallest_bool());
}
int main() {
solve();
return 0;
}
三分套三分套三分
这题和二维上的费马点经典问题很类似,费马点可以三分套三分,类似的此题可以三分套三分套三分,结论需要记住:
当确定 y , z y,z y,z时,函数图像随 x x x是一个单峰函数;当确定 x , z x,z x,z时,函数图像随 y y y是一个单峰函数;当确定 x , y x,y x,y时,函数图像随 z z z是一个单峰函数
因为数据范围和精度都没有很大,那么三分五十次即可(三分一百次 2.8 s 2.8s 2.8s险过)
//
// Created by Happig on 2020/11/11
//
#include <bits/stdc++.h>
#include <unordered_map>
#include <unordered_set>
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 = 50;
const double L = -100000;
const double R = 100000;
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);
scanf("%d", &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("%.8lf\n", f(lx, ly, lz));
return 0;
}