洛谷 P4035 [JSOI2008]球形空间产生器

洛谷 P4035 [JSOI2008]球形空间产生器

思路

高斯消元

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

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

\[\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})\)

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

\[\sum\limits_{j = 1}^{n}(a_{i_1,j} - x_j)^2 = \sum\limits_{j = 1}^{n}(a_{i_2,j} - x_j)^2 \]

展开式子得

\[\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)) \]

进一步化简

\[\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)) \]

\[\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\)个线性方程组，则最后的矩阵为

\[\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;
}