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

题目传送门

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

---

前置技能

下面简要说一说有关高斯消元和与其相关的线性代数最基本的知识。

假如要求解线性方程组$Ax=b$，其中的$A$表示系数矩阵（一般是方阵），$x$表示未知数向量，$b$表示常数向量（默认是列向量）。

很明显，在$A$矩阵可逆时$x=A^{-1}b$，假如能够求出$A$矩阵的逆矩阵，做一个矩乘就解出来了。

有这样一个定理：
矩阵$A$可逆的充要条件是$A\rightarrow E$，其中那个箭头表示初等行变换（由于我打不出那个符号所以用箭头代替），$E$是单位矩阵（只有对角线为1）。

初等行变换就是那三条基本操作：①换行②乘k③+其他行乘k。

那我们要求$x$，就要先算出$A^{-1}$。其实还有更简单的方法：将$b$放在$A$的右侧，构造出新的矩阵$(A,b)$，我们称其为增广矩阵。有$(A,b)\rightarrow(E,x)$，于是直接对增广矩阵变换就行了。小学的二元一次方程组怎么解，多元一次就怎么解。

而高斯消元就是通过对增广矩阵进行初等行变换解方程组的过程。

高斯消元有一种写法是消成对角线矩阵。另外一种消出上三角矩阵然后回带。第二种会快一点。我写的是简单的第一种，又称为高斯—约旦消元。

在实现上，为了减小误差，对于每一个未知数应该找系数绝对值最大的那行将其保留。这叫做列主元高消。

如何判断$n$元方程组的解的情况呢？假设$A$是方阵吧。

①方程有唯一解$R(A)=R(A,b)=n$;

②方程无解$R(A)<R(A,b)$;

③方程有不定解$R(A)=R(A,b)<n$.（有自由变元，可求通解）

$R(A)$表示$A$的秩，就是行变换成行阶梯型矩阵后矩阵的非零行数。

更简单的实现其实只要$A[i][i]=0$就不是唯一解了。

---

本题思路

每两个点的坐标搞出一条方程，将二次项约掉，剩下一个线性方程。将n条方程的系数和常数放进矩阵的对应位置直接高消即可。

时间复杂度，也是高消的复杂度：$O(n^3)$。

---

代码

#include <bits/stdc++.h>
#define maxn 15

using namespace std;

int n;
double pos[maxn][maxn], A[maxn][maxn];

void Gauss(){
    for(int i = 1; i <= n; i++){
        int x = i;
        for(int j = i+1; j <= n; j++)
            if(fabs(A[j][i]) > fabs(A[x][i]))  x = i;

        swap(A[x], A[i]);

        for(int j = n+1; j >= i; j--)  A[i][j] /= A[i][i];

        for(int j = 1; j <= n; j++){
            if(i == j)  continue;
            for(int k = n+1; k >= i; k--)  A[j][k] -= A[j][i] * A[i][k];
        }
    }
}

int main(){

    scanf("%d", &n);

    for(int i = 1; i <= n+1; i++)
        for(int j = 1; j <= n; j++)
            scanf("%lf", &pos[i][j]);   

    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++){
            A[i][j] = 2 * (pos[i+1][j] - pos[i][j]);
            A[i][n+1] += pos[i+1][j] * pos[i+1][j] - pos[i][j] * pos[i][j];
        }

    Gauss();

    for(int i = 1; i <= n; i++)
        printf("%.3lf ", A[i][n+1]);

    return 0;
} 

---