秋游时候的一句话激起了我学这个算法的念头（
大家好，今天我们来学习用高斯消元法解线性方程组
初 一 数 学 既 视 感
废话少说，我们现在开始！
假设你有一个 \(n\) 个方程 \(n\) 个未知数的方程组：
\(\begin{cases}a_{1,1}x_1+a_{1,2}x_2+\dots+a_{1,n}x_n=a_{1,n+1}\\a_{2,1}x_1+a_{2,2}x_2+\dots+a_{2,n}x_n=a_{2,n+1}\\\dots\\a_{n,1}x_1+a_{n,2}x_2+\dots+a_{n,n}x_n=a_{n,n+1}\end{cases}\)
其中 \(a_{i,j}\) 表示每一项的系数。
那么你就可以用矩阵 \(A\) 表示这个方程组：
\(\begin{bmatrix}a_{1,1}\ a_{1,2}\ &\cdots&a_{1,n}\ a_{1,n+1}\\a_{2,1}\ a_{2,2}\ &\cdots&a_{2,n}\ a_{2,n+1}\\\dots\\a_{n,1}\ a_{n,2}\ &\cdots&a_{n,n}\ a_{n,n+1}\end{bmatrix}\)
我们的目标是将这个方程组变为
\(\begin{cases}x_1=t_1\\x_2=t_2\\x_3=t_3\\\dots\\x_n=t_n\end{cases}\)
的形式，也就是说，我们要将矩阵变为最后一列为 \(t_i\)，对角线全是 \(1\) ，其它全是 \(0\) 的形式。
现在问题的关键是，怎样进行转化呢？
在初一的时候，yyq曾经告诉我们，解一次方程组可以用加减消元法和代入消元法
那么问题来了，在我们在实现高斯消元的时候究竟是用加减消元法还是代入消元法呢？
答案是，先用加减消元法消成下三角矩阵，再用代入消元法反代。
不妨举个例子，假设有以下的三元线性方程组
\(\begin{cases}2x+9y-5z=10\\4x+20y+z=24\\x-\dfrac{1}{2}y+3z=8\end{cases}\)
我们先把 \(x\) 消掉，注意，方便起见，我们把 \(1,2,3\) 式中 \(x\) 前面的系数都化为 \(1\)，原方程转化为：
\(\begin{cases}x+\dfrac{9}{2}y-\dfrac{5}{2}z=5\\x+5y+\dfrac{1}{4}z=6\\x-\dfrac{1}{2}y+3z=8\end{cases}\)
再拿 \(2,3\) 式分别与 \(1\) 式子相减：
\(\begin{cases}x+\dfrac{9}{2}y-\dfrac{5}{2}z=5\\\quad\quad\dfrac{1}{2}y+\dfrac{11}{4}z=1\\\quad\ -5y+\dfrac{11}{2}z=3\end{cases}\)
然后我们再消 \(y\)，注意，1 式已经消过了，就不能再动它了，所以我们把 \(2,3\) 式前 \(y\) 的系数都化为 \(1\) 并且相减。
\(\begin{cases}x+\dfrac{9}{2}y-\dfrac{5}{2}z=5\\\quad\quad y+\dfrac{11}{2}z=2\\\quad\quad y-\dfrac{11}{10}z=-\dfrac{3}{5}\end{cases}\)
\(\begin{cases}x+\dfrac{9}{2}y-\dfrac{5}{2}z=5\\\quad\quad y+\dfrac{11}{2}z=2\\\quad\quad\quad-\dfrac{33}{5}z=-\dfrac{13}{5}\end{cases}\)
最后消 \(z\)（其实准确来说只能叫做系数化 \(1\)）
\(\begin{cases}x+\dfrac{9}{2}y-\dfrac{5}{2}z=5\\\quad\quad y+\dfrac{11}{2}z=2\\\quad\quad\quad\quad\quad z=\dfrac{13}{33}\end{cases}\)
至此，第 \(i\) 个式子中未知数 \(i\) 前面的系数都变成了 \(1\)，并且 \(a_{j,i}=0\ (j>i)\)，我们称之为“下三角矩阵”
然后我们回代，把 \(z=\dfrac{13}{33}\) 回代到 2 式可以得到 \(y\) 的值，再把 \(y,z\) 的值代入 1 式可以得到 \(x\) 的值。
最终我们得到方程的解为：
\(\begin{cases}x=\dfrac{889}{132}\\y=-\dfrac{1}{6}\\z=\dfrac{13}{33}\end{cases}\)
上述部分的代码如下：
#include <bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define fz(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
#define ffe(it,v) for(__typeof(v.begin()) it=v.begin();it!=v.end();it++)
#define fill0(a) memset(a,0,sizeof(a))
#define fill1(a) memset(a,-1,sizeof(a))
#define fillbig(a) memset(a,63,sizeof(a))
#define pb push_back
#define mp make_pair
typedef pair<int,int> pii;
typedef long long ll;
const int MAXN=100+5;
const double EPS=1e-3;
int n;double a[MAXN][MAXN],ans[MAXN];
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++)
cin>>a[i][j];
for(int i=1;i<=n;i++){//elimination by adding&minusing
for(int j=i+1;j<=n+1;j++) a[i][j]/=a[i][i];a[i][i]=1;
for(int j=i+1;j<=n;j++){
for(int k=i+1;k<=n+1;k++){
a[j][k]-=a[j][i]*a[i][k];
} a[j][i]=0;
}
}
ans[n]=a[n][n+1];
for(int i=n-1;i;i--){//elimination by substitution
ans[i]=a[i][n+1];
for(int j=i+1;j<=n;j++) ans[i]-=(a[i][j]*ans[j]);
}
for(int i=1;i<=n;i++) printf("%.2lf\n",ans[i]);
return 0;
}
可是，yyq当时还讲过，对于方程组
\(\begin{cases}a_1x+b_1y=c_1\\a_2x+b_2y=c_2\end{cases}\)
当 \(\dfrac{a_1}{a_2}=\dfrac{b_1}{b_2}=\dfrac{c_1}{c_2}\) 时，方程有无数组解，当 \(\dfrac{a_1}{a_2}=\dfrac{b_1}{b_2}\neq\dfrac{c_1}{c_2}\) 时，方程无解。
碰到这种情况怎么办呢？
不妨再举个例子。
\(\begin{cases}x+2y+3z+4w=2\\x+3y+5z+6w=5\\3x+8y+13z+16w=12\\x+4y+7z+10w=11\end{cases}\)
我们把 \(x\) 消完以后得到：
\(\begin{cases}x+2y+3z+4w=2\\y+2z+2w=3\\2y+4z+4w=6\\2y+4z+6w=9\end{cases}\)
如果此时你再把 \(y\) 消掉：
\(\begin{cases}x+2y+3z+4w=2\\y+2z+2w=3\\0=0\\2w=3\end{cases}\)
你会惊奇地发现，\(z\) 没了！
这意味着第三个方程其实是没有用的，实际上只有 \(3\) 个方程 \(4\) 个未知数。
继续把该方程组化为下三角矩阵，可以得到：
\(\begin{bmatrix}1&2&3&4&|&2\\0&1&2&2&|&3\\0&0&0&1&|&\dfrac{3}{2}\\0&0&0&0&|&0\end{bmatrix}\)
所以该方程组有 \(1\) 个自由元，就是 \(z\)。
那怎么判无解的情况呢？
比如说我们把上面的方程组中 “\(3x+8y+13z+16w=12\)” 改为 “\(3x+8y+13z+16w=13\)”。那么你会得到 \(0=1\)，这显然是不可能的，所以方程组就无解了。
更通俗地说，如果化成下三角矩阵后：
- 不存在某一行，它 \(n\) 个未知数前的系数都为 \(0\)，方程有唯一解。
- 有一行 \(n\) 个未知数前的系数都为 \(0\) 以及和都为 \(0\)，方程有无数组解。
- 有一行 \(n\) 个未知数前的系数都为 \(0\)，而和不为 \(0\)，那就无解了。
如果有无穷组解，那么方程自由元的个数就是下三角矩阵中全为 \(0\) 的行的个数。
把上面的代码进行一些魔改就可以过洛谷上的模板 P3389
#include <bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define fz(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
#define ffe(it,v) for(__typeof(v.begin()) it=v.begin();it!=v.end();it++)
#define fill0(a) memset(a,0,sizeof(a))
#define fill1(a) memset(a,-1,sizeof(a))
#define fillbig(a) memset(a,63,sizeof(a))
#define pb push_back
#define mp make_pair
typedef pair<int,int> pii;
typedef long long ll;
const int MAXN=100+5;
const double EPS=1e-3;
int n;double a[MAXN][MAXN],ans[MAXN];
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++)
cin>>a[i][j];
for(int i=1;i<=n;i++){
int t=i;
for(int j=i+1;j<=n;j++){
if(fabs(a[j][i])>fabs(a[t][i])) t=j;
}
if(fabs(a[t][i])<EPS) return puts("No Solution"),0;
for(int j=1;j<=n+1;j++) swap(a[t][j],a[i][j]);
for(int j=i+1;j<=n+1;j++) a[i][j]/=a[i][i];a[i][i]=1;
for(int j=i+1;j<=n;j++){
for(int k=i+1;k<=n+1;k++){
a[j][k]-=a[j][i]*a[i][k];
} a[j][i]=0;
}
}
ans[n]=a[n][n+1];
for(int i=n-1;i;i--){
ans[i]=a[i][n+1];
for(int j=i+1;j<=n;j++) ans[i]-=(a[i][j]*ans[j]);
}
for(int i=1;i<=n;i++) printf("%.2lf\n",ans[i]);
return 0;
}
以上都是 \(n\) 个方程 \(n\) 个未知数的情况，但如果把范围扩大一点，变成 \(m\) 个方程 \(n\) 个未知数的情况怎么办呢？
其实也大差不差，我们先尝试消 \(x_1\)，如果能消（在剩余的方程中存在某个方程，\(x_1\) 前面的系数不为 \(0\)）那就消，否则继续看 \(x_2\)，直到消完 \(x_n\) 或把 \(m\) 个方程用完。
例如
\(\begin{cases}a+b+c+d+e+f=1\\2a+2a+2c+3d+3e+4f=3\\3a+3b+3c+5d+5e+6f=5\end{cases}\)
我们先消 \(a\)。
\(\begin{cases}a+b+c+d+e+f=1\\d+e+2f=1\\2d+2e+3f=2\end{cases}\)
然后发现 \(b\) 和 \(c\) 都消不了，去消 \(d\)。
\(\begin{cases}a+b+c+d+e+f=1\\d+e+2f=1\\-f=0\end{cases}\)
\(e\) 也消不了，消 \(f\)。
\(\begin{cases}a+b+c+d+e+f=1\\d+e+2f=1\\f=0\end{cases}\)
至此，已将原方程组消成下三角矩阵。
异或方程组
异或方程组是形如
\(\begin{cases}a_{1,1}x_1\operatorname{xor}a_{1,2}x_2\operatorname{xor}\dots\operatorname{xor}a_{1,n}x_n=a_{1,n+1}\\a_{2,1}x_1\operatorname{xor}a_{2,2}x_2\operatorname{xor}\dots\operatorname{xor}a_{2,n}x_n=a_{2,n+1}\\\dots\\a_{n,1}x_1\operatorname{xor}a_{n,2}x_2\operatorname{xor}\dots\operatorname{xor}a_{n,n}x_n=a_{n,n+1}\end{cases}\)
的方程组。
其实与普通的方程组大差不差，也是把它消成下三角矩阵然后回代就可以了。
由于不少题目都要用到这个东西，故把它单独拎出来提了一下。
例题
P4035 [JSOI2008]球形空间产生器
模板题，把球心坐标设为 \((x_1,x_2,\dots,x_n)\)，列出 \(n+1\) 个方程，然后相邻两方程相减可以消去 \(x_i^2\) 项。
P3232 [HNOI2013]游走
记得秋游的时候 ycx 跟我说过“你是在随机游走吗”
然后 10 天之后我就做到了这道题。
根据贪心的思想，我们肯定要按每条边经过的期望次数从大到小排序并依次编号 \(1,2,\dots,m\)。
设 \(f_i\) 为经过点 \(i\) 的期望次数，\(deg_i\) 为点 \(i\) 的度数。
如果已知 \(f_i\)，那么边 \((u,v)\) 经过的次数的期望值为 \(\dfrac{f_u}{deg_u}+\dfrac{f_v}{deg_v}\)
接下来我们的任务就是求 \(f_i\)。
稍微想一想即可列出转移方程 \(f_i=\sum\limits_{(j,i)\in E}\dfrac{f_j}{deg_j}\)
也就是你枚举点 \(i\) 是从哪个点来的并乘上从 \(j\) 走到 \(i\) 的概率就是经过点 \(i\) 的期望。
但是点 \(1\) 和点 \(n\) 比较特殊，由于游走是从点 \(1\) 开始的，经过点 \(i\) 的期望次数应当加 \(1\)。而到达点 \(n\) 后就停止游走了，所以 \(f_n=0\)，并且转移的时候也不能从 \(n\) 转移来。
高斯消元就可以求出 \(f_i\)，时间复杂度 \(n^3\)。
832E Vasya and Shifts
设 \(x_i\) 为第 \(i\) 个字符串使用次数（\(x_i \in [0,4]\)）
很显然可以根据结果字符串上每一位的字符列出 \(m\) 个方程出来。
如果方程有解，假设方程中有 \(r\) 个自由元，那么答案就是 \(5^{r}\)（每个自由元的值都可以随便取）；否则答案为 \(0\)。
典型的变成 \(m\) 个方程 \(n\) 个未知数的问题，直接上板子（大雾。
还有一点，那就是对于每组询问都 \(n^3\) 跑一遍会 TLE，但注意到每组询问只改变常数项的值，于是我们可以把 \(q\) 组常数项的值一起放到矩阵中，形成一个 \(m\times(n+q)\) 的矩阵并对其进行高斯消元。
时间复杂度 \((n+q)\times m^2\)