【BZOJ1013】[JSOI2008] 球形空间产生器（高斯消元）

点此看题面

大致题意： 给定一个$n$维球体上的$n+1$个点，请你求出这个球体的圆心的位置。

---

列出方程

这一看就是一道解方程题。

我们可以设这个球体的圆心的位置为$(x_1,x_2,..x_n)$，并设每个点到圆心的距离为$dis$。

借助题目中给出的公式，我们可以得到以下方程：

$\{\begin{array}{l}\sqrt{(x_1-a_{1,1}{)}^2+(x_2-a_{1,2}{)}^2+...+(x_n-a_{1,n}{)}^2}=dis\\ \sqrt{(x_1-a_{2,1}{)}^2+(x_2-a_{2,2}{)}^2+...+(x_n-a_{2,n}{)}^2}=dis\\ ......\\ \sqrt{(x_1-a_{n+1,1}{)}^2+(x_2-a_{n+1,2}{)}^2+...+(x_n-a_{n+1,n}{)}^2}=dis\end{array}$

---

方程的转化

原方程看起来十分麻烦，又有平方，又有开方，很难解，因此我们要将它转化一下。

将方程两边同时平方，可得：

$\{\begin{array}{l}(x_1-a_{1,1}{)}^2+(x_2-a_{1,2}{)}^2+...+(x_n-a_{1,n}{)}^2=di{s}^2\\ (x_1-a_{2,1}{)}^2+(x_2-a_{2,2}{)}^2+...+(x_n-a_{2,n}{)}^2=di{s}^2\\ ......\\ (x_1-a_{n+1,1}{)}^2+(x_2-a_{n+1,2}{)}^2+...+(x_n-a_{n+1,n}{)}^2=di{s}^2\end{array}$

但是，这些方程全部都是二次方程，好像非常难做。

因此，我们考虑将每个方程展开：

$\{\begin{array}{l}{x}_1^2-2a_{1,1}{x}_1+a_{1,1}^2+x_2^2-2a_{1,2}{x}_2+a_{1,2}^2+...+x_n^2-2a_{1,n}{x}_n+a_{1,n}^2=di{s}^2\\ x_1^2-2a_{2,1}{x}_1+a_{2,1}^2+x_2^2-2a_{2,2}{x}_2+a_{2,2}^2+...+x_n^2-2a_{2,n}{x}_n+a_{2,n}^2=di{s}^2\\ ......\\ x_1^2-2a_{n+1,1}{x}_1+a_{n+1,1}^2+x_2^2-2a_{n+1,2}{x}_2+a_{n+1,2}^2+...+x_n^2-2a_{n+1,n}{x}_n+a_{n+1,n}^2=di{s}^2\end{array}$

这时候，我们显然可以看出每个方程中左边都有$x_1^2+x_2^2+...+x_n^2$，右边都有$di{s}^2$，不难想到，将第$1\sim n$个方程分别减去第$n+1$个方程，便可以得到一个新的方程组，而且是一次的：

$\{\begin{array}{l}2(a_{n+1,1}-a_{1,1})\cdot x_1+...+2(a_{n+1,n}-a_{1,n})\cdot x_n=a_{n+1,1}^2-a_{1,1}^2+...+a_{n+1,n}^2-a_{1,n}^2\\ 2(a_{n+1,1}-a_{2,1})\cdot x_1+...+2(a_{n+1,n}-a_{2,n})\cdot x_n=a_{n+1,1}^2-a_{2,1}^2+...+a_{n+1,n}^2-a_{2,n}^2\\ ...\\ 2(a_{n+1,1}-a_{n,1})\cdot x_1+...+2(a_{n+1,n}-a_{n,n})\cdot x_n=a_{n+1,1}^2-a_{n,1}^2+...+a_{n+1,n}^2-a_{n,n}^2\end{array}$

由于$a$数组是题目中给出的，我们就可以直接高斯消元了。

$Link$

高斯消元详见博客高斯消元入门

---

代码

#include<bits/stdc++.h>
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))
#define abs(x) ((x)<0?-(x):(x))
#define LL long long
#define ull unsigned long long
#define N 10
using namespace std;
int n;
namespace Gauss
{
	const double eps=1e-10;//eps是一个极小值，防止精度误差
	double a[N+5][N+5],s[N+5];
	inline void swap(double &x,double &y)
	{
		double t=x;x=y,y=t;
	}
	inline void GetData()//将读入的数据转化为方程的系数
	{
		register int i,j;
		for(i=1;i<=n+1;++i) for(j=1;j<=n;++j) scanf("%lf",&a[i][j]);
		for(i=1;i<=n;++i) for(s[i]=0,j=1;j<=n;++j) 
			s[i]+=a[n+1][j]*a[n+1][j]-a[i][j]*a[i][j],a[i][j]=2*(a[n+1][j]-a[i][j]);
	}
	inline void Find_line(int x)//找到一个行数大于等于x且第x个元素系数不为0的方程，将其移至第x行
	{
		register int i=x,j;
		while(i<=n&&fabs(a[i][x])<eps) ++i;
		for(j=1;j<=n;++j) swap(a[x][j],a[i][j]);
	}
	inline void PrintAns()
	{
		register int i,j,k;
		for(i=1;i<=n;++i)
		{
			for(Find_line(i),j=i+1;j<=n;++j)//消去[i+1~n]中每一行第i个元素
			{
				register double delta=-a[j][i]/a[i][i];
				for(s[j]+=s[i]*delta,k=i;k<=n;++k) a[j][k]+=a[i][k]*delta;
			}
		}
		for(i=n;i;--i) for(s[i]/=a[i][i],j=i-1;j;--j) s[j]-=a[j][i]*s[i];//计算出第i个未知数的值，并将第i个元素的值代入第1~i-1行的式子中消去第i个未知数
		for(i=1;i<=n;++i) printf("%.3lf ",s[i]);//输出每一个未知数的值
	}
}
int main()
{
	return scanf("%d",&n),Gauss::GetData(),Gauss::PrintAns(),0;
}