BZOJ - 1013 高斯消元

n维空间中给出n+1个球面上的点求圆心坐标(x0,x1,...xn-1)
任选其中一个点坐标如第一个点(a0,b0...z0)
(x0-a0)^2+(x1-b0)2+...=r^2
对于剩下的n个点都与上面的式子作差,把高次方的消去得到线性方程组
2(a1-a0)x0+2(b1-b0)x1+2(c1-c0)x2+...=dis1-dis0
...
共n行n列,
然后构造Ax=b丢进去高斯消元就行了
//其实我只是来试试板子(来自挑战程序设计竞赛)
//用vector来构造矩阵有点不适应

#include<bits/stdc++.h>
#define rep(i,j,k) for(int i=j;i<=k;i++)
using namespace std;
const int maxn = 1e5+11;
const double eps = 1e-8;
typedef vector<double> vec;
typedef vector<vec> mat;
vec gauss(const mat& A,const vec& b){
	int n=A.size();
	mat B(n,vec(n+1));//n element with vec(n+1)
	rep(i,0,n-1)rep(j,0,n-1)B[i][j]=A[i][j];
	rep(i,0,n-1) B[i][n]=b[i];//增广 
	rep(i,0,n-1){//未知数系数绝对值最大的移至第i行 
		int now=i;
		rep(j,i,n-1){
			if(abs(B[j][i])>abs(B[now][i])){
				now=j;
			}
		}
		swap(B[i],B[now]);
		if(abs(B[i][i])<eps)return vec();
		rep(j,i+1,n) B[i][j]/=B[i][i];
		rep(j,0,n-1) if(i!=j){
			rep(k,i+1,n){
				B[j][k]-=B[j][i]*B[i][k];
			}
		} 
	}
	vec x(n);
	rep(i,0,n-1) x[i]=B[i][n];
	return x; 
}
int main(){
	int n;
	double a0[22];
	while(scanf("%d",&n)!=EOF){
		mat A(n,vec(n));//n行n列 
		vec b(n);//n行1列 
		rep(i,0,n-1) scanf("%lf",&a0[i]);
		double dis0=0;
		rep(i,0,n-1) dis0+=a0[i]*a0[i];
		rep(i,0,n-1){
			double dis2=0;
			rep(j,0,n-1){
				scanf("%lf",&A[i][j]);
				dis2+=A[i][j]*A[i][j];
				A[i][j]=2*(A[i][j]-a0[j]);
			}
			b[i]=dis2-dis0;
		}
		vec x=gauss(A,b);
		for(int i=0;i<x.size();i++){
			if(i==x.size()-1) printf("%.3lf\n",x[i]);
			else printf("%.3lf ",x[i]);
		}
	}
	return 0;
}