更好的阅读体验：https://www.cnblogs.com/To-Carpe-Diem/p/19525292

大意

给你一个 $n$ 维空间，给我 $n + 1$ 个点，你要确定这个球的中心点的坐标。

思路

原始方程（$n + 1$ 个）对于每一个点 $i$，到球心 $(x_1, \dots, x_n)$ 的距离平方相等：$\sum_{j=1}^n (p_{i,j} - x_j)^2 = R^2$

降次消元（作差）用相邻两个方程 $(i+1)$ 和 $(i)$ 相减，消去二次项 $x_j^2$ 和 $R^2$：$$\sum_{j=1}^n (p_{i+1,j} - x_j)^2 - \sum_{j=1}^n (p_{i,j} - x_j)^2 = 0$$展开并整理：$$\sum_{j=1}^n (p_{i+1,j}^2 - 2p_{i+1,j}x_j) - \sum_{j=1}^n (p_{i,j}^2 - 2p_{i,j}x_j) = 0$$

标准线性方程将未知数 $x_j$ 移至左侧，常数移至右侧：$$\sum_{j=1}^n 2(p_{i+1,j} - p_{i,j})x_j = \sum_{j=1}^n (p_{i+1,j}^2 - p_{i,j}^2)$$

增广矩阵 $[A|B]$

$$\left[\begin{array}{ccc|c}2(p_{2,1}-p_{1,1}) & \dots & 2(p_{2,n}-p_{1,n})  & \sum (p_{2,j}^2 - p_{1,j}^2) \\ 2(p_{3,1}-p_{2,1}) & \dots & 2(p_{3,n}-p_{2,n})  & \sum (p_{3,j}^2 - p_{2,j}^2) \\ \vdots & \ddots & \vdots  & \vdots \\ 2(p_{n+1,1}-p_{n,1}) & \dots & 2(p_{n+1,n}-p_{n,n})  & \sum (p_{n+1,j}^2 - p_{n,j}^2) \end{array}\right]$$

代码

#include<iostream>
#include<cmath>
using namespace std;

const int MAXN = 15;
const double eps = 1e-6;

int n;
double a[MAXN][MAXN];
double p[MAXN][MAXN];

bool Jornh(){
    for(int i = 1;i <= n;i ++){
        int piv = i;
        for(int k = i + 1;k <= n;k ++){
            if(fabs(a[k][i]) > fabs(a[piv][i])){
                piv = k;
            }
        }
        for(int j = 1;j <= n + 1;j ++){
            swap(a[i][j], a[piv][j]);
        }
        if(fabs(a[i][i]) < eps){
            return false;
        }
        for(int k = 1;k <= n;k ++){
            if(k == i) continue;
            double fac = a[k][i] / a[i][i];
            for(int j = i;j <= n + 1;j ++){
                a[k][j] -= fac * a[i][j];
            }
        }
    }
    for(int i = 1;i <= n;i ++){
        a[i][n + 1] /= a[i][i];
    }
    return true;
}

int main(){

    cin >> n;
    for(int i = 1;i <= n + 1;i ++){
        for(int j = 1;j <= n;j ++){
            cin >> p[i][j];
        }
    }

    for(int i = 1;i <= n;i ++){
        double res = 0;
        for(int j = 1;j <= n;j ++){
            a[i][j] = 2.0 * (p[i + 1][j] - p[i][j]);
            res += (p[i + 1][j] * p[i + 1][j] - p[i][j] * p[i][j]);
        }
        a[i][n + 1] = res;
    }

    if(Jornh()){
        for(int i = 1;i <= n;i ++){
            printf("%.3lf ", a[i][n + 1]);
        }
    }
    return 0;
}