题目就是已知 $n$ 维空间中球面上 $n + 1$ 个点 $X$ 的坐标 $(x_1,x_2,…,x_n)$ 求球心 $O(v_1, v_2,…,v_n)$ 的坐标。
$n$ 维空间中欧氏距离是 $\sqrt{\sum_\limits{i=1}^n(v_i-x_i)^2}$
根据球面上的点 $x$ 到球心的距离相等,题目给了 $n+1$ 个这样的式子:
$$
用 x_{i,j} 表示向量 X_i 的第 j 个分量
$$
$$
\sum_\limits{i=1}^n(v_{1,i}-x_{1,i})^2 = R^2
$$
$$
\sum_\limits{i=1}^n(v_{2,i}-x_{2,i})^2 = R^2
$$
$$…$$
$$\sum_\limits{i=1}^n(v_{n+1,i}-x_{n+1,i})^2 \ = R^2$$
显然不是一个线性方程组,相邻两式相减,消去平方项得到:
$$
\sum_\limits{i=1}^n(2v_{1,i}-x_{1,i}-x_{2,i})(x_{2,i}-x_{1,i}) = 0
$$
$$
\sum_\limits{i=1}^n(2v_{2,i}-x_{2,i}-x_{3,i})(x_{3,i}-x_{2,i}) = 0
$$
$$
…
$$
$$
\sum_\limits{i=1}^n(2v_{n,i}-x_{n,i}-x_{n+1,i})(x_{n+1,i}-x_{n,i}) = 0
$$
把 $v$ 拿出来,系数缩起来,令 $b_i = x_{i+1,j} - x_{i,j}$
$$
\sum_\limits{i=1}^n2v_{1,i}(x_{2,i}-x_{1,i}) - \sum_\limits{i=1}^n(x_{1,i}+x_{2,i})(x_{2,i}-x_{1,i}) = 0
$$
$$
\sum_\limits{i=1}^n2b_i\times v_{1,i} - \sum_\limits{i=1}^n(x_{1,i}+x_{2,i})\times b_i = 0
$$
得到
$$
\sum_\limits{i=1}^n2b_i\times v_{1,i} = \sum_\limits{i=1}^n(x_{1,i}+x_{2,i})\times b_i
$$
$$
…
$$
$$
\sum_\limits{i=1}^n2b_i\times v_{n,i} = \sum_\limits{i=1}^n(x_{n,i}+x_{n+1,i})\times b_i \\
$$
竖着看,是一个由 $b^T=[b_1,b_2,…,b_n]、v_i^T=[v_{i,1},v_{i,2},…,v_{i,n}]、x_i^T=[x_{i,1},x_{i,2},…,x_{i,n}]$ 组成的线性方程组,跑一边高斯消元即可
#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#include "../template/debug.h"
#else
#define debug(...) 42
#endif
const double eps = 1e-9;
bool IsZero(double v) {
return abs(v) < 1e-9;
}
enum GAUSS_MODE {
DEGREE, ABS
};
int solution = 0;
template <typename T>
void GaussianElimination(vector<vector<T>>& a, int limit, GAUSS_MODE mode = DEGREE) {
if (a.empty() || a[0].empty()) {
return;
}
int h = static_cast<int>(a.size());
int w = static_cast<int>(a[0].size());
for (int i = 0; i < h; i++) {
assert(w == static_cast<int>(a[i].size()));
}
assert(limit <= w);
vector<int> deg(h);
for (int i = 0; i < h; i++) {
for (int j = 0; j < w; j++) {
deg[i] += !IsZero(a[i][j]);
}
}
int r = 0;
for (int c = 0; c < limit; c++) {
int id = -1;
for (int i = r; i < h; i++) {
if (!IsZero(a[i][c]) && (id == -1 || (mode == DEGREE && deg[i] < deg[id]) || (mode == ABS && abs(a[id][c]) < abs(a[i][c])))) {
id = i;
}
}
if (id == -1) {
continue;
}
if (id > r) {
swap(a[r], a[id]);
swap(deg[r], deg[id]);
for (int j = c; j < w; j++) {
a[id][j] = -a[id][j];
}
}
vector<int> nonzero;
for (int j = c; j < w; j++) {
if (!IsZero(a[r][j])) {
nonzero.push_back(j);
}
}
T inv_a = 1 / a[r][c];
for (int i = r + 1; i < h; i++) {
if (IsZero(a[i][c])) {
continue;
}
T coeff = -a[i][c] * inv_a;
for (int j : nonzero) {
if (!IsZero(a[i][j])) deg[i]--;
a[i][j] += coeff * a[r][j];
if (!IsZero(a[i][j])) deg[i]++;
}
}
++r;
}
for (r = h - 1; r >= 0; r--) {
for (int c = 0; c < limit; c++) {
if (!IsZero(a[r][c])) {
T inv_a = 1 / a[r][c];
for (int i = r - 1; i >= 0; i--) {
if (IsZero(a[i][c])) {
continue;
}
T coeff = -a[i][c] * inv_a;
for (int j = c; j < w; j++) {
a[i][j] += coeff * a[r][j];
}
}
break;
}
}
}
}
template <typename T>
vector<T> SolveLinearSystem(vector<vector<T>> a, const vector<T>& b, int w) {
int h = static_cast<int>(a.size());
assert(h == static_cast<int>(b.size()));
if (h > 0) {
assert(w == static_cast<int>(a[0].size()));
}
for (int i = 0; i < h; i++) {
a[i].push_back(b[i]);
}
GaussianElimination(a, w);
vector<T> x(w, 0);
for (int i = 0; i < h; i++) {
for (int j = 0; j < w; j++) {
if (!IsZero(a[i][j])) {
x[j] = a[i][w] / a[i][j];
break;
}
}
}
return x;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
int n;
cin >> n;
vector<vector<double>> mat(n + 1, vector<double>(n, 0));
for (int i = 0; i < n + 1; ++i) {
for (int j = 0; j < n; ++j) {
cin >> mat[i][j];
}
}
vector<vector<double>> coeff(n, vector<double>(n, 0));
vector<double> b(n);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
coeff[i][j] = (mat[i][j] - mat[i + 1][j]) * 2;
b[i] += mat[i][j] * mat[i][j] - mat[i + 1][j] * mat[i + 1][j];
}
}
auto res = SolveLinearSystem(coeff, b, n);
for (auto x : res) {
cout << fixed << setprecision(3) << x << ' ';
}
return 0;
}