[JSOI2008]球形空间产生器

题目描述

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

输入输出格式

输入格式：

第一行是一个整数 n (1<=N=10)。接下来的 n+1 行，每行有 n 个实数，表示球面上一点的 n 维坐标。每一个实数精确到小数点后 6位，且其绝对值都不超过 20000。

输出格式：

有且只有一行，依次给出球心的 n 维坐标（ n 个实数），两个实数之间用一个空格隔开。每个实数精确到小数点后 位。数据保证有解。你的答案必须和标准输出一模一样才能够得分。

输入输出样例

输入样例#1：

2
0.0 0.0
-1.0 1.0
1.0 0.0

输出样例#1：

0.500 1.500

说明

提示：给出两个定义：

---

高维球长啥样啊 ==

只能拿三维的来想><

我们知道球上的任意一点到球心的距离都等于半径

然后我们又有给定提示的第二个公式

我们就可以列方程辣

我们设给的n+1个点为\(A,B,C\)

设球心为\(x\)

显然我们可以列一个方程(以第一个点为例)

\(\sum_{i=1}^{n}{(A[i]-X[i])^2}=r^2\)

然后我们化简一下

\(\sum_{i=1}^{n}{A[i]^2} - 2*\sum_{i=1}^{n}{A[i]*X[i]}+\sum_{i=1}^{n}{X[i]^2}=r^2\)

然后我们发现这是个二次的方程,而且还不知道r

那高斯消元个锤子

然后我们发现一个问题

由n+1个点组成的方程组第1项是常数,每个方程的第3项都相同

所以我们可以作差

我们还是以第一项和第二项为例

这样式子就成了\(\sum_{i=1}^{n}{2(A[i]-B[i])X[i]}=\sum_{i=1}^{n}{A[i]^2-B[i]^2}\)

然后就化成了n个方程,并且只有n个未知数,并且最高次项是1次

就可以高斯消元辣

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
const int M = 15 ;
const double eps = 1e-7 ;
using namespace std ;
int n ;
double B[M][M] , p[M][M] , sum[M] ;
inline void Gauss() {
	for(int i = 1 ; i <= n ; i ++) {
		int Tab = i ;
		for(int j = i ; j <= n ; j ++)
		  if(fabs(fabs(B[j][i]) - fabs(B[Tab][i])) <= eps)
		    Tab = j ;
		for(int j = 1 ; j <= n + 1 ; j ++) swap(B[i][j] , B[Tab][j]) ;
	    for(int j = i + 1 ; j <= n + 1 ; j ++) B[i][j] /= B[i][i] ;
	    for(int j = 1 ; j <= n ; j ++) 
	      if(i != j)
	        for(int k = i + 1 ; k <= n + 1 ; k ++)
	          B[j][k] -= B[j][i] * B[i][k] ;
	}
	for(int i = 1 ; i <= n ; i ++) printf("%.3lf ",B[i][n + 1]) ;
}
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 ++)
			B[i][j] = ( p[i][j] - p[i + 1][j] ) * 2 ;
	    B[i][n + 1] = sum[i] - sum[i + 1] ;
	}
	Gauss() ; 
	return 0 ;
}