题目链接
洛谷
LOJ
题目分析
在 n 维空间中找到一个点,使其到所有点的距离相等。
最容易想到的是模拟退火,在空间中找点,计算和每个点的距离,然后求方差,方差小于先前最小值则更新答案。
还可以将欧氏距离的方程转化成 n 元一次方程组,进行高斯消元。
首先看看二维下的情况,对于空间内一点 \(O(a_0,b_0)\),该点到 \((a_1,b_1)\) \((a_2,b_2)\) \((a_3,b_3)\) 三个点距离相等,则有
\[
\begin{split}
(a_1-a_0)^2 + (b_1-b_0)^2 = r^2\\
(a_2-a_0)^2 + (b_2-b_0)^2 = r^2\\
(a_3-a_0)^2 + (b_3-b_0)^2 = r^2
\end{split}
\]
分别相减得到
\[
\begin{split}
2 \times (a_1 - a_2) \times a_0 + 2 \times (b_1 - b_2) \times b_0 = (a_1^2-a_2^2) + (b_1^2 - b_2^2) \\
2 \times (a_2 - a_3) \times a_0 + 2 \times (b_2 - b_3) \times b_0 = (a_2^2-a_3^2) + (b_2^2 - b_3^2)
\end{split}
\]
这就是我们需要的方程组了。
同理可以推广到高维,接着使用高斯消元即可解决。
代码实现
高斯消元:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
| #include <bits/stdc++.h>
using namespace std;
int n; double ipt[20][20]; double m[20][20];
int main() { cin >> n; for (int i = 1; i <= n; i++) { cin >> ipt[1][i]; }
for (int i = 2; i <= n + 1; i++) { for (int j = 1; j <= n; j++) { cin >> ipt[i][j]; m[i - 1][j] = 2 * (ipt[i - 1][j] - ipt[i][j]); m[i - 1][n + 1] += ipt[i - 1][j] * ipt[i - 1][j] - ipt[i][j] * ipt[i][j]; } }
for (int i = 1; i <= n; i++) { int maxn = i; for (int j = 1; j <= n; j++) { if (fabs(m[j][i]) > fabs(m[maxn][i])) maxn = j; } for (int j = 1; j <= n + 1; j++) { swap(m[i][j], m[maxn][j]); } for (int j = n + 1; j >= i; j--) { m[i][j] = m[i][j] / m[i][i]; } for (int j = 1; j <= n; j++) { if (i == j) continue; double tmp = m[j][i] / m[i][i]; for (int k = 1; k <= n + 1; k++) { m[j][k] -= m[i][k] * tmp; } } } for (int i = 1; i <= n; i++) { printf("%.3lf ", m[i][n + 1]); } return 0; } }
|
模拟退火:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
| #include <bits/stdc++.h>
using namespace std;
int n; double m[20][20]; double cur[20]; double ans[20]; double dis[20]; double tmp[20]; double min_var = 1e9;
double calc(double arr[]) { double sum = 0; for (int i = 1; i <= n + 1; i++) { dis[i] = 0; for (int j = 1; j <= n; j++) { dis[i] += (m[i][j] - arr[j]) * (m[i][j] - arr[j]); } dis[i] = 666 * dis[i]; sum += dis[i]; }
double aver = sum / (n + 1); double var = 0; for (int i = 1; i <= n + 1; i++) { var += (aver - dis[i]) * (aver - dis[i]); }
if (var < min_var) { min_var = var; for (int i = 1; i <= n; i++) { ans[i] = cur[i]; } }
return var; }
void anneal() { double t = 10000; unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); mt19937 rand_num(seed); uniform_real_distribution<double> range(-1, 1); while (t > 1e-6) { for (int i = 1; i <= n; i++) { tmp[i] = cur[i] + range(rand_num) * t; } double d1 = calc(tmp); double d2 = calc(cur); double delta = d1 - d2; if (exp(-delta / t) > fabs(range(rand_num))) { for (int i = 1; i <= n; i++) { cur[i] = tmp[i]; } } t *= 0.9999; } }
int main() { cin >> n; for (int i = 1; i <= n + 1; i++) { for (int j = 1; j <= n; j++) { cin >> m[i][j]; } } anneal(); for (int i = 1; i <= n; i++) printf("%.3lf ", ans[i]); return 0; }
|