P4035 [JSOI2008]球形空间产生器

我们可以根据球面上任意一点到球心的距离相等来列方程组

设球心为\((x_1,x_2,x_3……x_n)\)，题目中给出的第\(i\)个点的坐标为\((a_{i.1},a_{i,2}……a_{i,n})\)

可以得出\(\sum\limits_{j=0}^n(a_{i,j}-x_j)^2=c\)

这个方程组是由\(n\)个\(n\)元\(2\)次方程所得,不是线性方程组

我们可以通过相邻的两个方程做差，将它变成线性方程组

这里可以自己用二维的为例算一下^_

最后得到的矩阵为

\(\begin{bmatrix}2(a_{1,1}-a_{2,1})\ 2(a_{1,2}-a_{2,2})\ …… \sum\limits_{j=1}^{n}(a^2_{1,j}-a^2_{2,j})\\2(a_{2,1}-a_{3,1})\ 2(a_{2,2}-a_{3,2})\ …… \sum\limits_{j=1}^{n}(a^2_{2,j}-a^2_{3,j})\\2(a_{3,1}-a_{4,1})\ 2(a_{3,2}-a_{4,2})\ …… \sum\limits_{j=1}^{n}(a^2_{3,j}-a^2_{4,j})\\……\\2(a_{n,1}-a_{n+1,1})\ 2(a_{n,2}-a_{n+1,2})\ …… \sum\limits_{j=1}^{n}(a^2_{n,j}-a^2_{n+1,j})\end{bmatrix}\)

题目中保证了方程组有唯一解，所有直接对上述矩阵进行高斯消元即可

关于初等行交换的一些问题

• 初等行交换可以防止当前\(i\)位置出现\(0\)详见这个讨论
• 有\(double\)进行高斯消元时，初等行交换可以减小误差

进行高斯消元时对两行消
需要除以其中一行的第一个数字
如果这个数字接近\(0\)
那得到的数字就会很大
\(double\)精度会出问题
如果把还没有消的行中第一个数字比较大的换上来
除的时候就可以避免精度问题

以上内容来自学长wxyww有个\(nb\)的学长是真滴爽，什么时候我也能像学长这么\(nb\)

#include<bits/stdc++.h>
using namespace std;
const int N=25;
double a[N][N],b[N],c[N][N];
int n; 
int main()
{
	freopen("sphere9.in","r",stdin);
//	freopen(".out","w",stdout);
	cin>>n;
	for(int i=1;i<=n+1;++i)
		for(int j=1;j<=n;++j)
			cin>>a[i][j];
	for(int i=1;i<=n;++i)
	{
		for(int j=1;j<=n;++j)
		{
			c[i][j]=2*(a[i][j]-a[i+1][j]);
			b[i]+=a[i][j]*a[i][j]-a[i+1][j]*a[i+1][j];
		}
	}
	for(int i=1;i<=n;i++)
	{
		for(int j=i;j<=n;++j)
		{
			if(fabs(c[j][i])>1e-8)
			{
				for(int k=1;k<=n;++k) swap(c[i][k],c[j][k]);
				swap(b[i],b[j]);
			}
		}
		for(int j=1;j<=n;++j)
		{
			if(i==j) continue;
			double rate=c[j][i]/c[i][i];
			for(int k=i;k<=n;++k) c[j][k]-=c[i][k]*rate;
			b[j]-=b[i]*rate;
		}
	}
	for(int i=1;i<=n;++i) printf("%.13lf \n",b[i]/c[i][i]);
	return 0;
}