bzoj 1013: [JSOI2008]球形空间产生器sphere（高斯消元）

1013: [JSOI2008]球形空间产生器sphere

Time Limit: 1 Sec   Memory Limit: 162 MB
Submit: 5598   Solved: 2943
[ Submit][ Status][ Discuss]

Description

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

Input

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

Output

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

Sample Input

2

0.0 0.0

-1.0 1.0

1.0 0.0

Sample Output

0.500 1.500

设球的半径为r，圆心坐标为(x1, x2, x3, …, xn)

那么可以列出n+1个方程

(a1-x1)^2+(a2-x2)^2+...+(an-xn)^2==r^2

(b1-x1)^2+(b2-x2)^2+...+(bn-xn)^2==r^2

……

从第二项开始，每项都减去第一项，即可得到一个n元方程组

2*(a1-b1)x1+2*(a2-b2)*x2+...+2*(an-bn)*xn==a1^2-b1^2+a2^2-b2^2+...+an^2-bn^2

2*(a1-c1)x1+2*(a2-c2)*x2+...+2*(an-cn)*xn==a1^2-c1^2+a2^2-c2^2+...+an^2-bc^2

……

之后高斯消元求解 

#include<stdio.h>
#include<math.h>
#include<algorithm>
using namespace std;
double Jz[15][15], c[15];
void Guess(int n, int m);
int main(void)
{
	int n, i, j;
	double a[15], b[15], sum;
	while(scanf("%d", &n)!=EOF)
	{
		for(i=1;i<=n;i++)
			scanf("%lf", &a[i]);
		for(i=1;i<=n;i++)
		{
			sum = 0;
			for(j=1;j<=n;j++)
			{
				scanf("%lf", &b[j]);
				Jz[i][j] = 2*(a[j]-b[j]);
				sum += a[j]*a[j]-b[j]*b[j];
			}
			c[i] = sum;
		}
		Guess(n, n);
		printf("%.3f", c[1]);
		for(i=2;i<=n;i++)
			printf(" %.3f", c[i]);
		printf("\n");
	}
	return 0;
}

void Guess(int n, int m)
{
	int i, j, p, q, r;
	p = q = 1;
	while(p<=n && q<=m)
	{
		r = p;
		for(i=p+1;i<=n;i++)
		{
			if(fabs(Jz[i][q])>fabs(Jz[r][q]))
				r = i;
		}
		for(i=q;i<=m;i++)
			swap(Jz[r][i], Jz[p][i]);
		swap(c[r], c[p]);
		for(i=q+1;i<=m;i++)
			Jz[p][i] /= Jz[p][q];
		c[p] /= Jz[p][q];
		for(i=1;i<=n;i++)
		{
			if(i!=p && Jz[i][q])
			{
				for(j=q+1;j<=m;j++)
					Jz[i][j] -= Jz[p][j]*Jz[i][q];
				c[i] -= c[p]*Jz[i][q];
				Jz[i][q] = 0;
			}
		}
		q++, p++;
	}
}