Luogu P4035 [JSOI2008] 球形空间产生器 题解

前言

题目传送门 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;
}