BZOJ 1013 | 一份写了一堆注释的高斯消元题解

题意

给出\(n\)维直角坐标系中\(n + 1\)个点的坐标，它们都在一个\(n\)维球面上，求球心坐标。

题解

设球面上某两个点坐标为\((a_1, a_2, ... a_n)\)和\((b_1, b_2, ... b_n)\)，则可以列出方程:

\[(x_1 - a_1)^2 + (x_2 - a_2)^2 + ... + (x_n - a_n)^2 = (x_1 - b_1)^2 + (x_2 - b_2)^2 + ... + (x_n - b_n)^2 \]

括号打开化简得

\[2*(a_1 - b_1)x_1 + 2*(a_2 - b_2)x_2 + ... + 2*(a_n - b_n)x_n = a_1^2 - b_1^2 + a_2^2 - b_2^2 + ... + a_n^2 - b_n^2 \]

这样可得n个方程，然后高斯消元即可。

#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#define space putchar(' ')
#define enter putchar('\n')
using namespace std;
typedef long long ll;
template <class T>
void read(T &x){
    char c;
    bool op = 0;
    while(c = getchar(), c < '0' || c > '9')
	if(c == '-') op = 1;
    x = c - '0';
    while(c = getchar(), c >= '0' && c <= '9')
	x = x * 10 + c - '0';
    if(op) x = -x;
}
template <class T>
void write(T x){
    if(x < 0) putchar('-'), x = -x;
    if(x >= 10) write(x / 10);
    putchar('0' + x % 10);
}

const int N = 20;
int n;
double g[N][N], f[N][N];

int main(){
    /*
      下面以样例输入为例，模拟一下高斯消元的过程。
      样例输入：
      2
      0.0 0.0
      -1.0 1.0
      1.0 0.0
    */
    scanf("%d", &n);
    for(int i = 1; i <= n + 1; i++)
	for(int j = 1; j <= n; j++)
	    scanf("%lf", &g[i][j]);
    //首先，我们构造一个矩阵来表示n个方程。
    for(int i = 1; i <= n; i++)
	for(int j = 1; j <= n; j++){
	    f[i][j] = 2 * (g[i][j] - g[i + 1][j]);
	    f[i][n + 1] += g[i][j] * g[i][j] - g[i + 1][j] * g[i + 1][j];
	}
    /*
      上一步中，我们构造出了n个方程，分别是：
      2 * x[1] - 2 * x[2] = -2
      -4 * x[1] + 2 * x[2] = 1
      写成矩阵的形式就是：
      2 -2 -2
      -4 2 1
      这个矩阵存在f[][]数组中。
    */
    for(int i = 1; i <= n; i++){
	//这次循环，我们以x[i]为主元。
	int l = i;
	for(int j = i + 1; j <= n; j++)
	    if(fabs(f[j][i]) > fabs(f[l][i])) l = j;
	if(l != i)
	    for(int j = i; j <= n + 1; j++)
		swap(f[i][j], f[l][j]);
	//上一步中，我们找出了主元系数绝对值最大的一个方程，把它换到第i行（据说这么做精度能高一些）
	/*
	  程序第一次执行到这里的时候，矩阵变成了
	  -4 2 1
	  2 -2 -2
	 */
	for(int j = n + 1; j >= i; j--) //因为循环内要用到当前的f[i][i]，所以f[i][i]最后修改
	    f[i][j] /= f[i][i]; //将主元系数变为1，方程中其他项系数也等比例扩大/缩小。
	for(int j = i + 1; j <= n; j++)
	    for(int k = n + 1; k >= i; k--) //因为循环内要用到当前的f[j][i]，所以f[j][i]最后修改
		f[j][k] -= f[i][k] * f[j][i];
	/*
	  将第i个方程的每项系数扩大/缩小，使得主元系数和第j个方程主元系数相同
	  然后第j个方程 -= 第i个方程，这样第j个方程就消去了当前主元。
	*/
	/*
	  程序第一次进行到这一步的时候，矩阵变成了
	  1 -0.5 -0.25
	  0 -1 -1.5
	  第二次则变成了
	  1 -0.5 -0.25
	  0 1 1.5
	 */
    }
    for(int i = n; i; i--) //循环到i的时候，f[i][n + 1]表示的已经是最后的解x[i]
	for(int j = 1; j < i; j++)
	    f[j][n + 1] -= f[j][i] * f[i][n + 1]; //将x[i]带入到第j个方程中
    for(int i = 1; i <= n; i++)
	printf("%.3lf%c", f[i][n + 1], " \n"[i == n]);

    return 0;
}