「JSOI 2008」球形空间产生器

Description

题目链接：BZOJ 1013

有一个球形空间产生器能够在 $n$ 维空间中产生一个坚硬的球体。现在，你被困在了这个 $n$ 维球体中，你只知道球面上 $n+1$ 个点的坐标 $(x_{i,1},x_{i,2},\dots,x_{i,n})$，你需要以最快的速度确定这个 $n$ 维球体的球心坐标，以便于摧毁这个球形空间产生器。答案保留 $3$ 位小数。

数据范围：$1\le n\le 10$，$\vert x_{i,j}\vert\le 20000$

---

Solution

得到 $n+1$ 个等式：

$$\begin{cases} (x_1-p_{1,1})^2+(x_2-p_{1,2})^2+(x_3-p_{1,3})^2+\cdots+(x_n-p_{1,n})^2=R^2 \\ (x_1-p_{2,1})^2+(x_2-p_{2,2})^2+(x_3-p_{2,3})^2+\cdots+(x_n-p_{2,n})^2=R^2 \\ \vdots \\ (x_1-p_{n,1})^2+(x_2-p_{n,2})^2+(x_3-p_{n,3})^2+\cdots+(x_n-p_{n,n})^2=R^2 \\ (x_1-p_{n+1,1})^2+(x_2-p_{n+1,2})^2+(x_3-p_{n+1,3})^2+\cdots+(x_n-p_{n+1,n})^2=R^2 \end{cases}$$

然后用第 $n+1$ 个等式去消去前 $n$ 个等式的 ${x_i}^2$ 项，这样就得到了一个标准的 $n$ 元 $1$ 次方程组，直接高斯消元求出 $n​$ 个未知数。

时间复杂度：$O(n^3)$

---

Code

#include <cstdio>
#include <cmath>
#include <algorithm>

const int N=15;
int n;
double p[N][N],a[N][N],b[N],x[N];

void Gauss(int n) {
    for(int i=1;i<=n;++i) {
        int p=i;
        for(int k=i+1;k<=n;++k) if(fabs(a[k][i])>fabs(a[p][i])) p=k;
        if(i!=p) std::swap(a[i],a[p]),std::swap(b[i],b[p]);
        for(int k=i+1;k<=n;++k) {
            double d=a[k][i]/a[i][i];
            b[k]-=d*b[i];
            for(int j=i;j<=n;++j) a[k][j]-=d*a[i][j];
        }
    }
    for(int i=n;i>=1;--i) {
        for(int j=i+1;j<=n;++j) b[i]-=x[j]*a[i][j];
        x[i]=b[i]/a[i][i];
    }
}
void init() {
    for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) {
        a[i][j]=2.0*(p[n+1][j]-p[i][j]),b[i]+=p[n+1][j]*p[n+1][j]-p[i][j]*p[i][j];
    }
}
int main() {
    scanf("%d",&n);
    for(int i=1;i<=n+1;++i) for(int j=1;j<=n;++j) scanf("%lf",&p[i][j]);
    init();
    Gauss(n);
    for(int i=1;i<=n;++i) printf("%.3lf%c",x[i]," \n"[i==n]);
    return 0;
}