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

http://www.lydsy.com/JudgeOnline/problem.php?id=1013

只要列出方程组就能套高斯来解了。

显然距离相等，所以开不开平方都无所谓。

b表示圆心，可列

sigma((x[i][j]-b[j])^2)=sigma((x[i+1][j]-b[j])^2)

化简得

sigma(2*b[j]*(x[i+1][j]-x[i][j]))=sigma(x[i+1][j]^2-x[i][j]^2)

然后就得到n个等式，而且题目保证有解，就套高斯就行了。

第一次学高斯消元啊，其实就是在一个系数矩阵（叫做啥增广矩阵），每一次将第i行的下边的第i列的系数全部消除，最后得到一个倒三角矩阵，然后回代就是了。

也就是说

1 2 3

2 3 4

的矩阵，当前行在1，我们首先要消掉第二行的第一列，就相当于第一行*2（这个2就是2/1得来），然后第二行减去第一行。

就是消元嘛。。

在这里第i行i列的元素叫做主元素，而我们就是要将所有大于第i行的行将这一列的通过主元素消除。这里需要注意，主元素不要为0，且不要很小，要不然会严重影响精度（就是如果很小的话无法体现分母了）。那么我们在每一次消元时，要找第i列最大的，然后和当前行交换，这样能得到最大的主元素。

高斯中判断无解和无限解的问题还待研究。

#include <cstdio>
#include <cstring>
#include <cmath>
#include <string>
#include <iostream>
#include <algorithm>
#include <queue>
using namespace std;
#define rep(i, n) for(int i=0; i<(n); ++i)
#define for1(i,a,n) for(int i=(a);i<=(n);++i)
#define for2(i,a,n) for(int i=(a);i<(n);++i)
#define for3(i,a,n) for(int i=(a);i>=(n);--i)
#define for4(i,a,n) for(int i=(a);i>(n);--i)
#define CC(i,a) memset(i,a,sizeof(i))
#define read(a) a=getint()
#define print(a) printf("%d", a)
#define dbg(x) cout << (#x) << " = " << (x) << endl
#define printarr2(a, b, c) for1(_, 1, b) { for1(__, 1, c) cout << a[_][__]; cout << endl; }
#define printarr1(a, b) for1(_, 1, b) cout << a[_] << '\t'; cout << endl
inline const int getint() { int r=0, k=1; char c=getchar(); for(; c<'0'||c>'9'; c=getchar()) if(c=='-') k=-1; for(; c>='0'&&c<='9'; c=getchar()) r=r*10+c-'0'; return k*r; }
inline const int max(const int &a, const int &b) { return a>b?a:b; }
inline const int min(const int &a, const int &b) { return a<b?a:b; }const int N=15;
double A[N][N], a[N][N];
int n;void gauss() {for1(i, 1, n-1) {int pos=i;for1(j, i+1, n) if(abs(A[pos][i])<abs(A[j][i])) pos=j;for1(j, 1, n+1) swap(A[i][j], A[pos][j]);for1(j, i+1, n) {double y=A[j][i]/A[i][i];for1(k, i, n+1) A[j][k]-=y*A[i][k];}}for3(i, n, 1) {for1(j, i+1, n) A[i][n+1]-=A[j][n+1]*A[i][j];A[i][n+1]/=A[i][i];}
}
int main() {read(n);for1(i, 1, n+1) for1(j, 1, n) scanf("%lf", &a[i][j]);for1(i, 1, n) {for1(j, 1, n) A[i][j]=2*(a[i+1][j]-a[i][j]);for1(j, 1, n) A[i][n+1]+=a[i+1][j]*a[i+1][j]-a[i][j]*a[i][j];}gauss();for1(i, 1, n-1) printf("%.3lf ", A[i][n+1]); printf("%.3lf\n", A[n][n+1]);return 0;
}

Description

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

Input

第一行是一个整数，n。接下来的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

HINT

数据规模：

对于40%的数据，1<=n<=3

对于100%的数据，1<=n<=10

提示：给出两个定义：

1、球心：到球面上任意一点距离都相等的点。

2、距离：设两个n为空间上的点A, B的坐标为(a1, a2, …, an), (b1, b2, …, bn)，则AB的距离定义为：dist = sqrt( (a1-b1)^2 + (a2-b2)^2 + … + (an-bn)^2 )

Source