n 维空间中给出 n+1 个球面上的点求圆心坐标 (x0,x1,…xn-1)

任选其中一个点坐标如第一个点 (a0,b0…z0)

(x0-a0)^2+(x1-b0)^2+…=r^2

对于剩下的 n 个点都与上面的式子作差，把高次方的消去得到线性方程组

2(a1-a0)x0+2(b1-b0)x1+2(c1-c0)x2+…=dis1-dis0

…

共 n 行 n 列，

然后构造 Ax=b 丢进去高斯消元就行了

//其实我只是来试试板子 (来自挑战程序设计竞赛)

#include<bits/stdc++.h>
#define rep(i,j,k) for(int i=j;i<=k;i++)
using namespace std;
const int maxn = 1e5+11;
const double eps = 1e-8;
typedef vector<double> vec;
typedef vector<vec> mat;
vec gauss(const mat& A,const vec& b){
    int n=A.size();
    mat B(n,vec(n+1));//n element with vec(n+1)
    rep(i,0,n-1)rep(j,0,n-1)B[i][j]=A[i][j];
    rep(i,0,n-1) B[i][n]=b[i];//增广 
    rep(i,0,n-1){//未知数系数绝对值最大的移至第 i 行 
        int now=i;
        rep(j,i,n-1){
            if(abs(B[j][i])>abs(B[now][i])){
                now=j;
            }
        }
        swap(B[i],B[now]);
        if(abs(B[i][i])<eps)return vec();
        rep(j,i+1,n) B[i][j]/=B[i][i];
        rep(j,0,n-1) if(i!=j){
            rep(k,i+1,n){
                B[j][k]-=B[j][i]*B[i][k];
            }
        } 
    }
    vec x(n);
    rep(i,0,n-1) x[i]=B[i][n];
    return x; 
}
int main(){
    int n;
    double a0[22];
    while(scanf("%d",&n)!=EOF){
        mat A(n,vec(n)); // n 行 n 列 
        vec b(n); // n 行 1 列 
        rep(i,0,n-1) scanf("%lf",&a0[i]);
        double dis0=0;
        rep(i,0,n-1) dis0+=a0[i]*a0[i];
        rep(i,0,n-1){
            double dis2=0;
            rep(j,0,n-1){
                scanf("%lf",&A[i][j]);
                dis2+=A[i][j]*A[i][j];
                A[i][j]=2*(A[i][j]-a0[j]);
            }
            b[i]=dis2-dis0;
        }
        vec x=gauss(A,b);
        for(int i=0;i<x.size();i++){
            if(i==x.size()-1) printf("%.3lf\n",x[i]);
            else printf("%.3lf ",x[i]);
        }
    }
    return 0;
}