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丟進去高斯消元就行了
//其實我只是來試試板子(來自挑戰程序設計競賽)
//用vector來構造矩陣有點不適應

#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;
}

BZOJ - 1013 高斯消元