n维球上的任意一点到球心距离相等，故设球心坐标为(x1,x2,...,xn)

则有公式∑(a_i,j-x_j)²=C

设法消去平方项x²，考虑相邻两项相减

可将公式化为Σ2(a_i,j-a_i+1,j)x_j=Σ(a_i,j²-a_i+1.j²)

此时就可以构造线性方程组了

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<queue>
using namespace std;
double eps=1e-8;
double a[1010][1010],b[1010][1010];
int n;
int flag=0;


bool pd(int k,int x,int y){
    if (fabs(fabs(a[x][k])-fabs(a[y][k]))>eps) return fabs(a[x][k])>fabs(a[y][k]); 
    for (int i=k+1;i<=n;i++){
        if (fabs(fabs(a[x][i])-fabs(a[y][i]))>eps) return fabs(a[x][i])<fabs(a[y][i]); 
    }
    return 0;
}
void print()
{
   int i,j;
   for(i=1;i<=n;i++)
    {
        for(j=1;j<=n+1;j++)
        {
            printf("%f ,",a[i][j]);
        }
        printf("\n");
    }
}

void swaprow(int now,int maxrow)
{
    double tmp;
    for(int i=now;i<=n+1;i++)
    {
        tmp=a[now][i];
        a[now][i]=a[maxrow][i];
        a[maxrow][i]=tmp;
    }
    return;
}

void Gauss()
{
    double tmp;
    int maxrow;
    for(int j=1;j<=n;j++)
    {
        maxrow=j;
        for(int i=j;i<=n;i++)
        {
            if(pd(j,i,maxrow))maxrow=i;
            maxrow=i;
        }//选择主元
        if(maxrow!=j)
        {
            swaprow(j,maxrow);
        }
        for(int i=j+1;i<=n;i++)
        {
            tmp=a[i][j]/a[j][j];
            for(int k=j;k<=n+1;k++)
            {
                a[i][k]-=a[j][k]*tmp;
            }
        }
        //print();
        //cout<<endl;
    }
    for(int i=n;i>=1;i--)
    {
        for(int j=i+1;j<=n;j++)
            a[i][n+1]-=a[i][j]*a[j][n+1];
        if(a[i][n+1]==0&&a[i][i]==0){
            flag=1;
        }
        if(a[i][n+1]!=0&&a[i][i]==0){
            flag=2;
            return;
        }
        a[i][n+1]/=a[i][i];
    }
    return ;
}

int main()
{
    cin>>n;
    for(int i=1;i<=n+1;i++)
    {
        for(int j=1;j<=n;j++)
        {
            cin>>b[i][j];
        }
    }
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
        {
            a[i][j]=2*(b[i][j]-b[i+1][j]);
            a[i][n+1]+=(b[i][j]*b[i][j]-b[i+1][j]*b[i+1][j]);
        }
    }
    Gauss();
    for(int i=1;i<=n;i++)
    {
        if(fabs(a[i][n+1])<eps)a[i][n+1]=0;
        printf("%.3lf ",a[i][n+1]);
    }
    return 0;
}