BZOJ1013 球形空间产生器sphere （高斯消元）

题目大意

有一个球形空间产生器能够在n维空间中产生一个坚硬的球体。现在，你被困在了这个n维球体中，你只知道球面上n+1个点的坐标，你需要以最快的速度确定这个n维球体的球心坐标，以便于摧毁这个球形空间产生器。

---

题解

以二维为例，设圆心坐标为 $x0$ 和 $y0$，若有两点的坐标( $x1$ , $y1$)、( $x2$ , $y2$)，则有以下关系：

$$(x1-x0)^2+(y1-y0)^2=(x2-x0)^2+(y2-y0)^2$$
展开得：
$$x1^2+x0^2-2x1x0+y1^2+y0^2-2y1y0=x2^2+x0^2-2x2x0+y2^2+y0^2-2y2y0$$
整理得：
$$x1^2-2x1x0+y1^2-2y1y0=x2^2-2x2x0+y2^2-2y2y0$$
移项：
$$2(x2-x1)x0+2(y2-y1)y0=x2^2-x1^2+y2^2-y1^2$$

因为 $x1$ $y1$ 和 $x2$ $y2$ 已知，也就是说都是常数项，那么就可以将n+1个n维点转化为n+1个n次方程，用高斯消元直接得出答案。

---

代码

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
using namespace std;

#define sqr(k) ((k)*(k))
const int maxn=11;
int n,m;
double a[maxn][maxn],f[maxn][maxn];
double ans[maxn];

inline void Gauss() {
    for(int i=1;i<=n;++i) {
        int cer=i;
        for(int j=cer+1;j<=n;++j) if(fabs(f[cer][i])<fabs(f[j][i])) cer=j;
        for(int j=i;j<=m;++j)
            swap(f[cer][j],f[i][j]);
        for(int j=i+1;j<=n;++j) {
            double mul=f[j][i]/f[i][i];
            for(int k=i;k<=m;++k)
                f[j][k]=f[j][k]-f[i][k]*mul;
        }
    }
    for(int i=n;i>=1;--i) {
        double cur=f[i][m];
        for(int j=n;j>i;--j)
            cur-=ans[j]*f[i][j];
        ans[i]=cur/f[i][i];
    }
    return ;
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("input.txt","r",stdin);
    freopen("output.txt","w",stdout);
#endif

    scanf("%d",&n);
    m=n+1;
    for(int i=0;i<=n;++i)
        for(int j=1;j<=n;++j)
            scanf("%lf",&a[i][j]);

    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j) {
            f[i][j]=(a[i][j]-a[i-1][j])*2;
            f[i][m]+=sqr(a[i][j])-sqr(a[i-1][j]);
        }
    Gauss();

    for(int i=1;i<n;++i) printf("%.3lf ",ans[i]);
    printf("%.3f\n",ans[n]);

    return 0;
}