洛谷 P4035 [JSOI2008]球形空间产生器-优快云博客

本文链接：https://blog.youkuaiyun.com/Loceaner/article/details/105831499

本文介绍了一种使用高斯消元法解决在n维球形空间中找到一个点，该点到给定的n+1个点距离相等的问题。通过将几何问题转化为线性代数问题，最终形成一个线性方程组，并利用高斯消元法求解未知点坐标。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

思路

高斯消元

题意：在 $n$ 维的球形空间中给定 $n + 1$ 个点，求到所有 $n + 1$ 个点的距离相等的点的坐标

由题意易知我们要求出在 $n$ 维空间中的一个点 $x_1,x_2,x_3,…x_n)$ ，满足：

$∀i∈[1,n+1],∑j=1n(ai,j−xj)2=R\forall i\in[1,n+ 1],\sum\limits_{j = 1}^{n}(a_{i,j} - x_j)^2 = R$

其中 $R$ 是一个常数，第 $i$ 个点的坐标为 $a_{i,1},a_{i,2},a_{i,3},…a_{i,n})$

假设有 $i1,i2∈[1,n+1],i1≠i2i_1,i2\in[1,n+1],i1≠i2$ ，由 $∑j=1n(ai1,j−xj)2=R\sum\limits_{j = 1}^{n}(a_{i_1,j} - x_j)^2 = R$ 和 $∑j=1n(ai2,j−xj)2=R\sum\limits_{j = 1}^{n}(a_{i_2,j} - x_j)^2 = R$ 得

$∑j=1n(ai1,j−xj)2=∑j=1n(ai2,j−xj)2\sum\limits_{j = 1}^{n}(a_{i_1,j} - x_j)^2 = \sum\limits_{j = 1}^{n}(a_{i_2,j} - x_j)^2$

展开式子得

$∑j=1n(ai1,j2+xj2−2∗(ai1,j∗xj))=∑j=1n(ai2,j2+xj2−2∗(ai2,j∗xj))\sum\limits_{j = 1}^{n}(a_{i_1,j}^2 + x_j^2 - 2 *(a_{i_1,j}* x_j))=\sum\limits_{j = 1}^{n}(a_{i_2,j}^2 + x_j^2 - 2 *(a_{i_2,j}* x_j))$

进一步化简

$∑j=1n(ai1,j2−2∗(ai1,j∗xj))=∑j=1n(ai2,j2−2∗(ai2,j∗xj))\sum\limits_{j = 1}^{n}(a_{i_1,j}^2 - 2 *(a_{i_1,j}* x_j))=\sum\limits_{j = 1}^{n}(a_{i_2,j}^2 - 2 *(a_{i_2,j}* x_j))$

$∑j=1n2∗(ai1,j−ai2,j)xj=∑j=1n(ai1,j2−ai2,j2)\sum\limits_{j = 1}^{n} 2 *(a_{i_1,j} - a_{i_2,j})x_j= \sum\limits_{j = 1}^{n}(a_{i_1,j}^2 - a_{i_2,j}^2)$

这样就转化成了一个线性方程组，由此可以用每个 $i$ 和 $i + 1$ 两两组合，得出 $n$ 个线性方程组，则最后的矩阵为

$…2(an,n−an+1,n)∑j=1n(an,j2−an+1,j2)]\begin{bmatrix}2(a_{1,1}-a_{2,1})\ 2(a_{1,2}-a_{2,2})\ … 2(a_{1,n}-a_{2,n})\sum\limits_{j=1}^{n}(a^2_{1,j}-a^2_{2,j})\\2(a_{2,1}-a_{3,1})\ 2(a_{2,2}-a_{3,2})\ … 2(a_{2,n}-a_{3,n})\sum\limits_{j=1}^{n}(a^2_{2,j}-a^2_{3,j})\\2(a_{3,1}-a_{4,1})\ 2(a_{3,2}-a_{4,2})\ … 2(a_{3,n}-a_{4,n})\sum\limits_{j=1}^{n}(a^2_{3,j}-a^2_{4,j})\\…\\2(a_{n,1}-a_{n+1,1})\ 2(a_{n,2}-a_{n+1,2})\ … 2(a_{n,n}-a_{n+1,n})\sum\limits_{j=1}^{n}(a^2_{n,j}-a^2_{n+1,j})\end{bmatrix}$

对此矩阵进行高斯消元求解即可，由于保证有解，所以直接做就好了

代码

/*
Author:loceaner
*/
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define eps 1e-8
using namespace std;

const int A = 22;
const int B = 1e6 + 11;
const int mod = 1e9 + 7;
const int inf = 0x3f3f3f3f;

inline int read() {
	char c = getchar(); int x = 0, f = 1;
	for( ; !isdigit(c); c = getchar()) if(c == '-') f = -1;
	for( ; isdigit(c); c = getchar()) x = x * 10 + (c ^ 48);
	return x * f;
}


int n;
double a[A][A], G[A][A], b[A];
//a是输入的点坐标，G是矩阵序数数组，b是方程右边的常数

int main() {
	n = read();
	for (int i = 1; i <= n + 1; i++) 
		for (int j = 1; j <= n; j++) scanf("%lf", &a[i][j]);
	for (int i = 1; i <= n; i++) 
		for (int j = 1; j <= n; j++) {
			G[i][j] = 2 * (a[i][j] - a[i + 1][j]);
			b[i] = b[i] + a[i][j] * a[i][j] - a[i + 1][j] * a[i + 1][j];
		} 
	for (int i = 1; i <= n; i++) {
		for (int j = i; j <= n; j++) {
			if (fabs(G[i][j]) > eps) {
				for (int k = 1; k <= n; k++) swap(G[i][k], G[j][k]);
				swap(b[i], b[j]);
			}
		}
		for (int j = 1; j <= n; j++) {
			if (i == j) continue;
			double tmp = G[j][i] / G[i][i];
			for (int k = i; k <= n; k++) G[j][k] -= G[i][k] * tmp;
			b[j] -= b[i] * tmp;
		}
	}
	for (int i = 1; i <= n; i++) printf("%.3lf ", b[i] / G[i][i]);
	return 0;
}