洛谷P4035 [JSOI2008]球形空间产生器（高斯消元）
球啊球 @xzz_233 qaq
高斯消元模板题，关键在于将已知条件转化为方程组。
可以发现题目要求的未知量有\(n\)个，题目却给了我们\(n+1\)个点的坐标，这其中必有玄机。
由高中数学知识可以知道，三点定圆（二维），四点定球（三维）······以此类推，应该是\(n+1\)个点才能确定一个\(n\)维空间下的球。
那么隐藏的另一个关键未知量在哪里呢？
想想圆的标准方程\((x-x_0)^2+(y-y_0)^2=r^2\)，除了圆心坐标，半径不也对这个圆起到决定性作用么？
接下来，额外设一个未知量——球的半径\(r\)，开始试着对条件式进行变换。
对于\(n+1\)个点，它们与球心的距离是定值\(r\)，那么我们可以得到形式如下的\(n+1\)个方程（a为球面点坐标，x为球心坐标）
\[(a_1-x_1)^2+(a_2-x_2)^2+···+(a_n-x_n)^2=r^2
\]
显然我们要把已知量和未知量分开，于是展开，移项
\[a_1^2-2a_1x_1+x_1^2+a_2^2-2a_2x_2+x_2^2+···+a_n^2-2a_nx_n+x_n^2=r^2
\]
\[2a_1x_1+2a_2x_2+···+2a_nx_n+r^2-x_1^2-x_2^2-···-x_n^2=a_1^2+a_2^2+···+a_n^2
\]
发现\(r^2-x_1^2-x_2^2-···-x_n^2\)与\(a\)无关，所以考虑换元，设\(t=r^2-x_1^2-x_2^2-···-x_n^2\)（实际上我们并不用求\(r\)）
终于，我们可以看到一个关于\(x_1,x_2,···,x_n,t\)的\((n+1)\)元方程组了，上高斯消元
具体实现看代码
#include<cmath>
#include<cstdio>
#define R register
#define init for(i=1;i<=n;++i)ne[i-1]=pr[i+1]=i
const int N=19;
int p[N],pr[N],ne[N];
double a[N][N];
int main(){
R int n,i,j,k,x;
R double mx,d;
scanf("%d",&n);++n;
init;//链表初始化，为了实现交换行
for(i=1;i<=n;++i){
for(j=1;j<n;++j){
scanf("%lf",&d);//处理系数
a[i][j]=d*2;a[i][n+1]+=d*d;
}
a[i][n]=1;//t的系数为1
}
for(k=1;k<=n;++k){
mx=0;//蒟蒻没有交换主元，而是交换行
//这样做防掉精度的效果可能不如交换主元
for(i=ne[0];i;i=ne[i])
if(mx<fabs(a[i][k]))
mx=fabs(a[i][k]),x=i;
d=a[p[k]=x][k];//选择当前a最大的一行
pr[ne[pr[x]]=ne[x]]=pr[x];
for(j=1;j<=n+1;++j)
a[x][j]/=d;
for(i=ne[0];i;i=ne[i])
for(d=a[i][j=k];j<=n+1;++j)
a[i][j]-=d*a[x][j];
}//高斯消元
init;
for(k=n;k;--k){
d=a[x=p[k]][n+1];
pr[ne[pr[x]]=ne[x]]=pr[x];
for(i=ne[0];i;i=ne[i])
a[i][n+1]-=d*a[i][k];
}//回代
for(k=1;k<n;++k)
printf("%.3f ",a[p[k]][n+1]);
puts("");
return 0;
}
