P4035 [JSOI2008]球形空间产生器

题意描述

洛谷

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

数据范围： (nleq 100)。

solution

高斯消元裸题。

设球心的第 (i) 维坐标为 (x_i), 然后根据球心的性质可得： (dis(a[i],x) = dis(a[i+1],x)) 。

展开可得：

((a_{i,1}-x_1)^2 + (a_{i,2}-x_2)^2 + .... + (a_{i,n}-x_n)^2 = (a_{{i+1},1}-x_1)^2 + (a_{i+1,2}-x_2)^2 + .... + (a_{i+1,n}-x_n)^2)

化简在合并同类项可得：

$2(a_{i+1,1}-a_{i,1})x_1 + 2(a_{i+1,2}-a_{i,2})x_2 + .... + 2(a_{i+1,n}-a_{i,n})x_n = a_{i,1}^2 - a_{i+1,1}^2 + a_{i,2}^2 - a_{i+1,2}^2 + ....+a_{i,n}^2 - a_{i+1,n}^2 $

预处理一个 (sum[i]) 表示 (displaystylesum_{j=1}^{n} a_{i,j}^2), 上面的式子可以变为：

$2(a_{i+1,1}-a_{i,1})x_1 + 2(a_{i+1,2}-a_{i,2})x_2 + .... + 2(a_{i+1,n}-a_{i,n})x_n = sum[i] - sum[i+1] $

然后我们就可以得到 (n) 个等式，把 (x_i) 当成未知数，用高斯消元求解即可。

时间复杂度： (O(n^3))

Code

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath> 
using namespace std;
const int N = 2010;
const double eps = 1e-8;
int n;
double sum[N],p[N][N],a[N][N];
inline int read()
{
    int s = 0, w = 1; char ch = getchar();
    while(ch < '0' || ch > '9'){if(ch == '-') w = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9'){s = s * 10 + ch - '0'; ch = getchar();}
    return s * w;
}
void Gauss()//高斯消元
{
	for(int i = 1; i <= n; i++)
	{
		int pai = i;
		while(pai <= n && fabs(a[pai][i]) < eps) pai++;
		for(int j = 1; j <= n+1; j++) swap(a[i][j],a[pai][j]);
		double rate = a[i][i];
		for(int j = 1; j <= n+1; j++) a[i][j] /= rate;
		for(int j = 1; j <= n; j++)
		{
			if(j == i) continue;
			rate = a[j][i] / a[i][i];
			for(int k = i; k <= n+1; k++) a[j][k] -= a[i][k] * rate;
		}
	}
}
int main()
{
	scanf("%d",&n);
	for(int i = 1; i <= n+1; i++)
	{
		for(int j = 1; j <= n; j++)
		{
			scanf("%lf",&p[i][j]); 
			sum[i] += p[i][j] * p[i][j];
		}
	}
	for(int i = 1; i <= n; i++)
	{
		for(int j = 1; j <= n; j++)
		{
			a[i][j] = 2 * (p[i+1][j]-p[i][j]);
		}
		a[i][n+1] = sum[i+1]-sum[i];//构建n个等式
	}
	Gauss();
	for(int i = 1; i <= n; i++) printf("%.3lf ",a[i][n+1]);
	printf("
");
	return 0;
}