[BZOJ 1013][JSOI 2008] 球形空间产生器sphere 题解(高斯消元)

[BZOJ 1013][JSOI 2008] 球形空间产生器sphere

Description

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

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

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

Solution

1.考虑构造对于不同的店同样结构的方程。

因为是一个球,所以每个点到球心的距离都相等,

我们设这个半径为R,球心坐标为O(x1,x2,....,xn);

那么对于每一个点P(ai1,ai2,...,ain):我们易得

sqrt ( ( ai1 - x1 ) ^ 2 + ( ai2 - x2 ) ^ 2 + ... + ( ain - xn ) ^2 ) = R;

2.考虑构造方程组

将上式两侧平方再展开,得

-2 ( ai1 * x1 + ai2 * x2 + ... + ain * xn )
+( ai1 ^ 2 + ai2 ^ 2 + ... + ain ^ 2 )
+( x1 ^ 2 + x2 ^ 2 + ... + xn ^ 2 )
= R ^ 2;

这时我们看数据,给出n+1个点,n个点就可以构造该方程组,那多给的一个点是用来干什么的呢(设选第一个点为这个点)?

没错!消掉重复出现的部分( x1 ^ 2 + x2 ^ 2 + ... + xn ^ 2 ) 和R ^ 2,

即令其他所有的点的方程都减掉多的一个点的方程,整理得到其他方程格式为:

2 ( ai1 * x1 + ai2 * x2 + ... + ain * xn )
= ( ai1 ^ 2 + ai2 ^ 2 + ... + ain ^ 2 )
- ( a11 ^ 2 + a12 ^ 2 + ... + a1n ^ 2 )
- 2 ( a11 * x1 + a12 * x2 + ... + a1n * xn ) ;

右侧是常数,左侧展开就是一个愉快的高斯消元方程组,解方程组即可。

这里使用的高斯消元法是一种比较毒瘤的方法,详解参考我的随笔:http://www.cnblogs.com/COLIN-LIGHTNING/p/8981923.html

Code

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;

const int max_n=15;
double a[max_n][max_n+1],v[max_n],del[max_n],x;
int n,w[max_n]; 

inline double sqr(double x){return x*x;} 

void init(){
	for(int i=1;i<=n;++i)scanf("%lf",&del[i]);
	for(int i=1;i<=n;++i){
		for(int j=1;j<=n;++j){
			scanf("%lf",&x);
			a[i][n+1]+=sqr(x)-sqr(del[j]);
			a[i][j]=2*(x-del[j]);
		}
	}
}

void gauss(){
	double eps=1e-6;
	for(int i=1;i<=n;++i){//enumerate the equation;
		int p=0;	      //Record the position of the largest number; 
		double mx=0; 	  //Recording the largest number;
		for(int j=1;j<=n;++j)
			if(fabs(a[i][j])-eps>mx){
				mx=fabs(a[i][j]);p=j;//fabs() returns the absolute value of float; 
			}
		w[i]=p;
		for(int j=1;j<=n;++j)
			if(i!=j){		//other equations
				double t=a[j][p]/a[i][p];
				for(int k=1;k<=n+1;++k)//n+1 is important
					a[j][k]-=a[i][k]*t;
			}
	}
	for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]];
	for(int i=1;i<=n;++i) printf("%.3lf ",v[i]);
}

int main(){
	scanf("%d",&n);
	init();
	gauss();
	return 0; 
}
posted @ 2018-05-02 21:11  COLINGAO  阅读(349)  评论(0编辑  收藏  举报