bzoj1013/luogu4035 球形空间生成器 (高斯消元)

根据各点到圆心的距离相等，可以列出有N个等号的方程

假设圆心坐标是(x,y,z,...)的话，$x^2,y^2,z^2$等是可以消掉的

于是整理一下，就变成了N元1次方程组，有N个方程、而且保证是相容的

高斯消元的话，就是拿着第一式去把剩下的第一项都消了，再拿第二式把剩下的第二项都消了，...到最后方程组呈一个阶梯状，然后从最后一点点往回带就能解出来

复杂度$O(n^3)$

#include<bits/stdc++.h>
#define pa pair<int,int>
#define CLR(a,x) memset(a,x,sizeof(a))
using namespace std;
typedef long long ll;
const int maxn=15;

inline ll rd(){
    ll x=0;char c=getchar();int neg=1;
    while(c<'0'||c>'9'){if(c=='-') neg=-1;c=getchar();}
    while(c>='0'&&c<='9') x=x*10+c-'0',c=getchar();
    return x*neg;
}

int N;
double pos[maxn][maxn],a[maxn][maxn],b[maxn];

int main(){
    //freopen("","r",stdin);
    int i,j,k;
    N=rd();
    for(i=1;i<=N+1;i++){
        for(j=1;j<=N;j++){
            scanf("%lf",&pos[i][j]);
        }
        if(i>1){
            for(j=1;j<=N;j++){
                a[i-1][j]=2*(pos[i][j]-pos[i-1][j]);
                a[i-1][N+1]+=pos[i][j]*pos[i][j]-pos[i-1][j]*pos[i-1][j];
            }
        }
    }
    for(i=1;i<=N;i++){
        for(j=i+1;j<=N;j++){
            double p=a[j][i]/a[i][i];
            for(k=i;k<=N+1;k++){
                a[j][k]-=p*a[i][k];
            }
        }
    }
    for(i=N;i;i--){
        double r=a[i][N+1];
        for(j=N;j>i;j--){
            r-=b[j]*a[i][j];
        }
        b[i]=r/a[i][i];
    }
    for(i=1;i<=N;i++)
        printf("%.3lf ",b[i]);
    return 0;
}