BZOJ - 1013. 球形空间产生器

#题面

#题目描述

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

#输入格式

第一行是一个整数 $n$。接下来的 $n + 1$ 行，每行有 $n$ 个实数，表示球面上一点的 $n$ 维坐标。每一个实数精确到小数点后 $6$ 位，且其绝对值都不超过 $20000$。

#输出格式

有且只有一行，依次给出球心的 $n$ 维坐标（ $n$ 个实数），两个实数之间用一个空格隔开。每个实数精确到小数点后 $3$ 位。数据保证有解。你的答案必须和标准输出一模一样才能够得分。

#输入输出样例

样例输入 #1

2
0.0 0.0
-1.0 1.0
1.0 0.0

样例输出 #1

0.500 1.500

#数据范围与约定

对于 $100\%$ 的数据，$1 \leq n \leq 10$。

#提示

给出两个定义：

1. 球心：到球面上任意一点距离都相等的点。
2. 距离：设两个 n 维空间上的点 $A, B$ 的坐标为 $(a_1, a_2, \cdots , a_n), (b_1, b_2, \cdots , b_n)$，则 $AB$ 之间的距离定义为：$\mathit{dist} = \sqrt{ \sum_{i = 1}^{n} (a_i - b_i)^2 }$。

#思路

设球的半径为 $R$，球心坐标为 $(x_1, x2, \ldots, x_n)$。那么根据球体的性质，显然有：

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

那么展开来写可以得到一个方程组的形式：

$$\begin{cases} \begin{aligned} (a_{1, 1} - x_1)^2 + (a_{1, 2} - x_2)^2 + \cdots + (a_{1, n} - x_n)^2 &= R^2 \\ (a_{2, 1} - x_1)^2 + (a_{2, 2} - x_2)^2 + \cdots + (a_{2, n} - x_n)^2 &= R^2 \\ \vdots \\ (a_{n + 1, 1} - x_1)^2 + (a_{n + 1, 2} - x_2)^2 + \cdots + (a_{n + 1, n} - x_n)^2 &= R^2 \\ \end{aligned} \end{cases}$$

但这个方程组是由 $n + 1$ 个 $n$ 元二次方程构成的，不是线性方程组。本着「出题人绝对不会多给 1 Byte 有用数据」的原则，我们可以察觉到 $n + 1$ 中的异样之处。通过将相邻两个方程作差, 得到的结果是：

$$\begin{cases} \begin{aligned} 2(a_{2, 1} - a_{1, 1}) x_1 + 2(a_{2, 2} - a_{1, 2}) x_2 + \cdots + 2(a_{2, n} - a_{1, n}) x_n &= a_{2, 1}^2 + a_{2, 2}^2 + \cdots + a_{2, n}^2 - a_{1, 1}^2 - a_{1, 2}^2 - \cdots - a_{1, n}^2 \\ 2(a_{3, 1} - a_{2, 1}) x_1 + 2(a_{3, 2} - a_{2, 2}) x_2 + \cdots + 2(a_{3, n} - a_{2, n}) x_n &= a_{3, 1}^2 + a_{3, 2}^2 + \cdots + a_{3, n}^2 - a_{2, 1}^2 - a_{2, 2}^2 - \cdots - a_{2, n}^2 \\ \vdots \\ 2(a_{n + 1, 1} - a_{n, 1}) x_1 + 2(a_{n + 1, 2} - a_{n, 2}) x_2 + \cdots + 2(a_{n + 1, n} - a_{n, n}) x_n &= a_{n + 1, 1}^2 + a_{n + 1, 2}^2 + \cdots + a_{n + 1, n}^2 - a_{n, 1}^2 - a_{n, 2}^2 - \cdots - a_{n, n}^2 \\ \end{aligned} \end{cases}$$

可以转化成增广矩阵的形式：

$$\left[ \begin{array}{cccc|c} 2(a_{2, 1} - a_{1, 1}) & 2(a_{2, 2} - a_{1, 2}) & \cdots & 2(a_{2, n} - a_{1, n}) & \sum_{j = 1}^{n} (a_{2, j}^2 - a_{1, j}^2) \\ 2(a_{3, 1} - a_{2, 1}) & 2(a_{3, 2} - a_{2, 2}) & \cdots & 2(a_{3, n} - a_{2, n}) & \sum_{j = 1}^{n} (a_{3, j}^2 - a_{2, j}^2) \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 2(a_{n + 1, 1} - a_{n, 1}) & 2(a_{n + 1, 2} - a_{n, 2}) & \cdots & 2(a_{n + 1, n} - a_{n, n}) & \sum_{j = 1}^{n} (a_{n + 1, j}^2 - a_{n, j}^2) \\ \end{array} \right]$$

对该矩阵进行高斯消元即可。

#代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

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

using std::cin;
using std::cout;
const char endl = '\n';

const int N = 15;
const double eps = 1e-6;

int n;
double a[N][N], b[N][N];

int main() {
    std::ios::sync_with_stdio(false);

    cin >> n;

    for (int i = 1; i <= n + 1; i++) {
        for (int j = 1; j <= n; j++) {
            cin >> a[i][j];
        }
    }

    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            b[i][j] = 2.0 * (a[i][j] - a[i + 1][j]);
            b[i][n + 1] += std::pow(a[i][j], 2) - std::pow(a[i + 1][j], 2);
        }
    }

    for (int i = 1; i <= n; i++) {
        for (int j = i; j <= n; j++) {
            if (std::abs(b[j][i]) > eps) {
                std::swap(b[i], b[j]);
            }
        }

        for (int j = 1; j <= n; j++) {
            if (i == j) continue;
            double x = b[j][i] / b[i][i];
            for (int k = i; k <= n + 1; k++) {
                b[j][k] -= b[i][k] * x;
            }
        }
    }

    for (int i = 1; i <= n; i++) {
        cout << std::fixed << std::setprecision(3) << b[i][n + 1] / b[i][i] << ' ';
    }
    cout << endl;

    return 0;
}