[FROM WOJ]#3004 [JSOI2008]球形空间产生器sphere

#3004 [JSOI2008]球形空间产生器sphere

题面
火星人不能忍受地球人对他们的歧视，终于发明了一种非常强大的武器： 球 形空间产生器。球形空间产生器能产生一个 N 维球体屏障，而且这个屏障是坚 不可摧的，被困在球体内的地球人就被切断了与外界的联系。 Js08 现在就被困在 了屏障中，情况十分危急，必须尽快找出并摧毁球形空间产生器。 Js08 经过摸索 和碰壁，给出了球体上 N+1 个点的坐标，希望你能够帮 Js08 找出球形空间产生 器的位置——它总是位于球形空间的球心。
球体定义为到空间中定点距离等于定长的点的集合。此外，空间中两点的距 离定义为：
dist = sqrt( (a1-b1)^2 + (a2-b2)^2 + … + (an-bn)^2 )

输入
输入文件第一行为整数 N，代表了空间的维度。
接下来 N+1 行，每行 N 个实数代表一个坐标。输入数据精确到小数点后 6 位。且其绝对值都不超过20000。
输入数据保证输出结果唯一。

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

样例输入
4
24.720172 14.560868 7.592580 25.639614
24.275281 8.136164 21.498302 -12.570364
33.948090 10.669582 4.705742 -2.648426
-22.717309 3.519597 16.222675 2.425559
-9.593546 18.349904 -15.757415 10.414070

样例输出
5.360 5.443 6.078 4.669

提示
数据规模：
对于40%的数据，1<=n<=3
对于100%的数据，1<=n<=10
提示：给出两个定义：
1、 球心：到球面上任意一点距离都相等的点。
2、 距离：设两个n为空间上的点A, B的坐标为$(a1,a2,\dots ,an),(b1,b2,\dots ,bn)$，则AB的距离定义为：$dist=sqrt((a1-b1{)}^2+(a2-b2{)}^2+\dots +(an-bn{)}^2)$

SOL
我们可以把每个点的坐标看成包含n个信息的集合，用{Rn}（化竞勿喷）表示n维球体半径坐标
对于两个点x，y，存在：
$$\sum _{i=1}^{i<=n}(Xi-Ri{)}^2=\sum _{i=1}^{i<=n}(Yi-Ri{)}^2$$
化简之后可以得到：
$$\sum _{i=1}^{i<=n}(Xi-Yi{)}^2=\sum _{i=1}^{i<=n}2(Xi-Yi)\ast Ri$$
左边是一个常数，于是就可以转成n元一次方程组
然后高斯消元。

代码：

#include<bits/stdc++.h>
#define pc putchar
#define N 15
using namespace std;
int n;
double in[N][N],a[N][N],ans[N],zero=1e-6;
inline void clean(int k){
	double register p;
	int register x=k,mx=fabs(a[k][k]);
	for(int register i=k+1;i<=n;i++)if(fabs(a[i][k])>mx)mx=fabs(a[i][k]),x=i;//注意用fabs 
	if(x!=k)for(int register i=k;i<=n+1;i++)swap(a[x][i],a[k][i]);
	if(fabs(a[k][k])<=zero)return;
	for(int register i=k+1;i<=n;i++){
		p=a[i][k]/a[k][k];
		for(int register j=k;j<=n+1;j++)a[i][j]-=p*a[k][j];
	}
}
inline void solve(){
	bool flag=1;
	for(int register i=n;i>0;i--){
		for(int register j=n;j>i;j--)a[i][n+1]-=a[i][j]*ans[j];
		ans[i]=a[i][n+1]/a[i][i];
	}
	for(int register i=1;i<=n;i++)printf("%.3lf ",ans[i]);
}
inline void gauss(){
	for(int register i=1;i<=n;i++)clean(i);
	solve();
}
signed main(){
	scanf("%d",&n);
	for(int register i=1;i<=n+1;i++)for(int register j=1;j<=n;j++)scanf("%lf",&in[i][j]);
	for(int register i=1;i<=n;i++){
		double register sum=0.000000;
		for(int register j=1;j<=n;j++)a[i][j]=2.000000*(in[i][j]-in[i+1][j]),sum+=in[i][j]*in[i][j]-in[i+1][j]*in[i+1][j];
		a[i][n+1]=sum;
	}
	gauss();
	return 0;
}