_bzoj1013 [JSOI2008]球形空间产生器sphere【高斯消元】

传送门：http://www.lydsy.com/JudgeOnline/problem.php?id=1013

保存高斯消元模版。

ps，这一题的英文名字是ヨスガノソラ的开发商~^_^

#include <cstdio>
#include <cstring>
#include <cmath>

const int maxn = 15;

int n;
double a0[maxn], s0, t, sqr_s, geass[maxn][maxn], tem[maxn];

int main(void) {
	//freopen("in.txt", "r", stdin);
	scanf("%d", &n);
	for (int j = 1; j <= n; ++j) {
		scanf("%lf", a0 + j);
		s0 += a0[j] * a0[j];
	}
	for (int i = 1; i <= n; ++i) {
		sqr_s = 0;
		for (int j = 1; j <= n; ++j) {
			scanf("%lf", &t);
			geass[i][j] = (t - a0[j]) * 2;
			sqr_s += t * t;
		}
		geass[i][n + 1] = sqr_s - s0;
	}
	
	int p;
	double xs;
	for (int i = 1; i <= n; ++i) {
		p = i;
		for (int j = i + 1; j <= n; ++j) {
			if (fabs(geass[j][i]) > fabs(geass[p][i])) {
				p = j;
			}
		}
		if (p != i) {
			memcpy(tem, geass[p], sizeof tem);
			memcpy(geass[p], geass[i], sizeof tem);
			memcpy(geass[i], tem, sizeof tem);
		}
		
		for (int j = i + 1; j <= n; ++j) {
			xs = geass[j][i] / geass[i][i];
			for (int k = i; k <= n + 1; ++k) {
				geass[j][k] -= geass[i][k] * xs;
			}
		}
	}
	
	for (int i = n; i; --i) {
		geass[i][n + 1] /= geass[i][i];
		geass[i][i] = 1.0;
		for (int j = i - 1; j; --j) {
			geass[j][n + 1] -= geass[j][i] * geass[i][n + 1];
			geass[j][i] = 0.0;
		}
	}
	
	for (int i = 1; i < n; ++i) {
		printf("%.3f ", geass[i][n + 1]);
	}
	printf("%.3f\n", geass[n][n + 1]);
	return 0;
}