[BZOJ 1013][JSOI 2008] 球形空间产生器

题目大意

给定 \(n\) 维空间上一个球面上的 \(n + 1\) 个点,求其圆心。数据保证有解。

\(1 \leqslant n \leqslant 10\)

\(|坐标| \leqslant 20,000\)

题目链接

BZOJ 1013

题解

考虑两个点 \(P\)\(Q\) ,记圆心为 \(O\) ,有 \[ \begin{align} \sum_{i = 1}^{n} (P_i - O_i)^2 &= \sum_{i = 1}^{n} (Q_i - O_i)^2 \\ \sum_{i = 1}^{n} P_i^2 - 2 P_i O_i &= \sum_{i = 1}^{n} Q_i^2 - 2 Q_i O_i \\ \sum_{i = 1}^{n} 2(P_i - Q_i)O_i &= \sum_{i = 1}^{n} P_i^2 - Q_i^2 \end{align} \] 发现二次项被约掉,于是可以高斯消元解方程。

代码

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
#include <cstdio>
#include <cmath>
#include <algorithm>
const int MAXN = 10;
const double EPS = 1e-7;
double a[MAXN][MAXN + 1];
bool GaussJordan(int n) {
for (int i = 0; i < n; i++) {
int max = i;
for (int j = i + 1; j < n; j++) if (fabs(a[j][i]) > fabs(a[max][i])) max = j;
if (fabs(a[max][i]) < EPS) return false;
if (max != i) for (int j = i; j <= n; j++) std::swap(a[i][j], a[max][j]);
for (int j = 0; j < n; j++) if (i != j) {
for (int k = n; k >= i; k--) a[j][k] -= a[i][k] / a[i][i] * a[j][i];
}
}
return true;
}
int main() {
int n;
scanf("%d", &n);
for (int i = 0; i <= n; i++) for (int j = 0; j < n; j++) {
double x;
scanf("%lf", &x);
if (i != n) a[i][j] += 2 * x, a[i][n] += x * x;
if (i != 0) a[i - 1][j] -= 2 * x, a[i - 1][n] -= x * x;
}
if (!GaussJordan(n)) return puts("-1"), 0;
for (int i = 0; i < n; i++) printf("%.3lf%c", a[i][n] / a[i][i], i == n - 1 ? '\n' : ' ');
return 0;
}