题目描述

有一个球形空间产生器能够在 $n$ 维空间中产生一个坚硬的球体。

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

注意： 数据保证有唯一解。

输入格式

第一行是一个整数 $n$。

接下来的 $n+1$ 行，每行有 $n$ 个实数，表示球面上一点的 $n$ 维坐标。

每一个实数精确到小数点后 $6$ 位，且其绝对值都不超过 $20000$。

输出格式

有且只有一行，依次给出球心的 $n$ 维坐标（$n$ 个实数），两个实数之间用一个空格隔开。

每个实数精确到小数点后 $3$ 位。

数据范围

$1\le n\le 10$

输入样例

2
0.0 0.0
-1.0 1.0
1.0 0.0

输出样例：

0.500 1.500

关于这道题中的球体，有两个显而易见的性质：

1. 球心是到球面上任意一点距离都相等的点。
2. 两个 $n$ 维空间上的点的距离可表示为 $\sqrt{(a_1-b_1)^2+(a_2-b_2)^2+\cdots+(a_n-b_n)^2}$。

根据这两个性质，就可以列出 $n+1$ 个 $n$ 元 $2$ 次方程，
其中球心表示为 $\lbrace x_1,x_2,\cdots,x_n \rbrace$，半径表示为 $R$，球面上的点 $i$ 表示为 $\lbrace a_{i,1},a_{i,2},\cdots,a_{i,n} \rbrace$，$1 \le i \le n+1$。

$\sqrt{\sum_{j=1}^n {(a_{i,j}-x_j)}^2}=R$
$\sum_{j=1}^n {(a_{i,j}-x_j)}^2=R^2$

现在还不能使用高斯消元，高斯消元只能解线性方程组，再来转化。
使用平方差公式 $(a-b)^2=a^2+b^2-2ab$ 转化得：
$\sum_{j=1}^n ({a_{i,j}}^2+{x_j}^2-2a_{i,j}x_j)=R^2$

都写出来：

$ \begin{cases} \sum_{j=1}^n ({a_{1,j}}^2+{x_j}^2-2a_{1,j}x_j)=R^2\\\ \sum_{j=1}^n ({a_{2,j}}^2+{x_j}^2-2a_{2,j}x_j)=R^2\\\ \cdots\\\ \sum_{j=1}^n ({a_{n+1,j}}^2+{x_j}^2-2a_{n+1,j}x_j)=R^2\\\ \end{cases} $

再把相邻的两个式子相减：

$ \begin{cases} \sum_{j=1}^n ({a_{2,j}}^2+{x_j}^2-2a_{2,j}x_j)-\sum_{j=1}^n ({a_{1,j}}^2+{x_j}^2-2a_{1,j}x_j)=0\\\ \sum_{j=1}^n ({a_{3,j}}^2+{x_j}^2-2a_{3,j}x_j)-\sum_{j=1}^n ({a_{2,j}}^2+{x_j}^2-2a_{2,j}x_j)=0\\\ \cdots\\\ \sum_{j=1}^n ({a_{n+1,j}}^2+{x_j}^2-2a_{n+1,j}x_j)-\sum_{j=1}^{n} ({a_{n,j}}^2+{x_j}^2-2a_{n,j}x_j)=0\\\ \end{cases} $

化简：

$ \begin{cases} \sum_{j=1}^n ({a_{2,j}}^2-{a_{1,j}}^2+2x_j(a_{1,j}-a_{2,j}))=0\\\ \sum_{j=1}^n ({a_{3,j}}^2-{a_{2,j}}^2+2x_j(a_{2,j}-a_{3,j}))=0\\\ \cdots\\\ \sum_{j=1}^n ({a_{n+1,j}}^2-{a_{n,j}}^2+2x_j(a_{n,j}-a_{n+1,j}))=0\\\ \end{cases} $

把变量和常数分开：

$ \begin{cases} \sum_{j=1}^n ({a_{2,j}}^2-{a_{1,j}}^2)=2x_j(a_{2,j}-a_{1,j})\\\ \sum_{j=1}^n ({a_{3,j}}^2-{a_{2,j}}^2)=2x_j(a_{3,j}-a_{2,j})\\\ \cdots\\\ \sum_{j=1}^n ({a_{n+1,j}}^2-{a_{n,j}}^2)=2x_j(a_{n+1,j}-a_{n,j})\\\ \end{cases} $

现在我们得到了一个线性方程组，可以使用高斯消元来解。
套进板子里就好了。

代码

代码里有注释

#include <cmath>
#include <iomanip>
#include <iostream>

using namespace std;

const double EPS = 1e-8;
int N;
double A[19][19];  //输入数组
double B[19][19];  //高斯消元数组

void In () {
    cin >> N;

    for (int i = 0; i < N + 1; i++) {
        for (int m = 0; m < N; m++) {
            cin >> A[i][m];
        }
    }
}

void Gauss () {  //高斯消元
    int r = 0;
    int t;

    for (int c = 0; c < N; c++) {
        t = r;

        for (int i = r; i < N; i++) {  //找在该列绝对值最大数所在的行
            if (fabs(B[i][c]) > fabs(B[t][c])) {
                t = i;
            }
        }

        if (fabs(B[t][c]) < EPS) {  //B[t][c]是否为0
            continue;
        }

        for (int i = c; i < N + 1; i++) {
            swap(B[t][i], B[r][i]);  //与最上面一行交换
        }

        for (int i = N; i > c - 1; i--) {
            B[r][i] /= B[r][c];  //把B[r][c]改成1
        }

        for (int i = r + 1; i < N; i++) {
            if (fabs(B[i][c]) > EPS) {
                for (int m = N; m > c - 1; m--) {
                    B[i][m] -= B[r][m] * B[i][c];  //把该列其它数改成0
                }
            }
        }

        r++;
    }

    for (int i = N - 1; i > -1; i--) {
        for (int m = i + 1; m < N; m++) {
            B[i][N] -= B[i][m] * B[m][N];  //计算答案
        }
    }
}

void Solve () {
    for (int i = 0; i < N; i++) {
        for (int m = 0; m < N; m++) {
            B[i][m] = 2 * (A[i + 1][m] - A[i][m]);
            B[i][N] += A[i + 1][m] * A[i + 1][m] - A[i][m] * A[i][m];
        }
    }

    Gauss();
}

void Out () {
    for (int i = 0; i < N; i++) {
        cout << fixed << setprecision(3) << B[i][N] << " ";
    }

    cout << endl;
}

int main () {

    In();
    Solve();
    Out();

    return 0;
}