前言
题目传送门 Luogu P4035 [JSOI2008] 球形空间产生器
这回真是模板题了。
题意简述
有一个 \(n(n\le 10)\) 维球体，给定该球面上的 \(n+1\) 个点的坐标，求球心的坐标（\(n\) 个实数）。答案保留三位小数。
思路
不难看出，本题肯定是要求解一个线性方程组的。问题在于怎么列出这个方程组。我们从二维的情况开始入手分析，设 \(O(p_1, p_2), A(a_1, a_2), B(b_1,b_2), C(c_1,c_2)\) ，于是我们有：
\[\begin{align*}
&OA = OB = OC\\
&(p_1-a_1)^2 + (p_2-a_2)^2\\
=&(p_1-b_1)^2 + (p_2-b_2)^2\\
=&(p_1-c_1)^2 + (p_2-c_2)^2
\end{align*}
\]
进行一些化简，我们得到这个连等式：
\[\begin{align*}
&2p_1a_1-a_1^2+2p_2a_2-a_2^2\\
=&2p_1b_1-b_1^2+2p_2b_2-b_2^2\\
=&2p_1c_1-c_1^2+2p_2c_2-c_2^2
\end{align*}
\]
任取两个，移项，我们得到三个关于 \(p_1, p_2\) 的方程。任取两个联立，我们得到以下方程组：（为了方便，我取第一、二和第二、三个联立）
\[\begin{align*}
\begin{cases}
2(a_1-b_1)p_1 + 2(a_2-b_2)p_2 = a_1^2 + a_2^2 - b_1^2 - b_2^2\\
2(b_1-c_1)p_1 + 2(b_2-c_2)p_2 = b_1^2 + b_2^2 - c_1^2 - c_2^2
\end{cases}
\end{align*}
\]
接下来，我们把所有的系数对应填入增广矩阵，得到如下矩阵：
\[\left[
\begin{array}{cc|c}
2(a_1-b_1) & 2(a_2-b_2) & a_1^2 + a_2^2 - b_1^2 - b_2^2\\
2(b_1-c_1) & 2(b_2-c_2) & b_1^2 + b_2^2 - c_1^2 - c_2^2\\
\end{array}
\right]
\]
观察规律，我们也能轻松地把它扩展到 \(n\) 维：
\[\left[
\begin{array}{ccc|c}
2(a_{1,1}-a_{2,1}) & \cdots & 2(a_{1,n}-a_{2,n}) & \sum_{i=1}^n a_{1,i}^2 - a_{2,i}^2\\
\vdots & \ddots & \vdots & \vdots\\
2(a_{n,1}-a_{n+1,1}) & \cdots & 2(a_{n,n}-a_{n+1,n}) & \sum_{i=1}^n a_{n,i}^2 - a_{n+1,i}^2\\
\end{array}
\right]
\]
最后，我们把这个矩阵丢给高斯-约旦消元即可。
由于题目保证有解，我们也不需要担心特殊情况。
代码
代码如下。
#include<iostream>
#include<cmath>
#define fastio ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
#define MIKU 0
using namespace std;const double eps = 1e-7;int n;
double a[15][15], p[2][15];void Gauss() {for(int i=1; i<=n; i++) {int r = i;for(int j=i+1; j<=n; j++) if(fabs(a[j][i]) > fabs(a[r][i])) r = j;if(r != i) swap(a[i], a[r]);double div = a[i][i];for(int j=i; j<=n+1; j++) a[i][j] /= div;for(int j=1; j<=n; j++) {if(j == i) continue;div = a[j][i];for(int k=i; k<=n+1; k++) a[j][k] -= a[i][k] * div;}}
}int main() {fastio;cin>>n;for(int i=1; i<=n; i++) cin>>p[1][i];for(int i=2; i<=n+1; i++) {for(int j=1; j<=n; j++) cin>>p[2][j];double t = 0;for(int j=1; j<=n; j++) {a[i-1][j] = 2 * (p[1][j] - p[2][j]);t += p[1][j] * p[1][j] - p[2][j] * p[2][j];}a[i-1][n+1] = t;swap(p[1], p[2]);}Gauss();for(int i=1; i<=n; i++) printf("%.3lf ", a[i][n+1]);return MIKU;
}