[note] 高斯消元
高斯消元是一种在 \(O(n^3)\) 的时间复杂度下解线性方程组的算法，方程组形如：
\[\left\{
\begin{aligned}
&a_{1,1}x_{1}+a_{1,2}x_{2}+\dots+a_{1,n}x_{n}=b_1\\
&a_{2,1}x_{1}+a_{2,2}x_{2}+\dots+a_{2,n}x_{n}=b_2\\
&\dots\\
&a_{n,1}x_{1}+a_{n,2}x_{2}+\dots+a_{n,n}x_{n}=b_n
\end{aligned}
\right.
\]
算法思路
1. 高斯消元
发现线性方程组的系数和等号右侧的常数是解方程的关键，因此我们可以把系数和常数拎出来，形成一个矩阵。
比如：
\[\left\{
\begin{aligned}
3x_1+7x_2+8x_3&=14\\
6x_1-x_2-2x_3&=7\\
-12x_1+x_2-x_3&=-10\\
\end{aligned}
\right.
\]
可以改写为：
\[\left(
\begin{array}{ccc|c}
3& 7& 8& 14\\
6& -1& -2& 7\\
-12& 1& -1& -10
\end{array}
\right)
\]
这样的矩阵被称为增广矩阵。
回忆解多元一次方程组的过程，我们发现我们可以对方程组进行如下操作：
- 交换任意两个方程的位置；
- 将一个方程加上另一个方程；
- 将整个方程乘上一个数。
对应到增广矩阵中，就是：
- 交换任意两行；
- 将某一行的数分别加上另一行对应的数；
- 将某一行的所有数乘上同一个数。
高斯消元的目的，就是通过上面的三种操作，将增广矩阵变为形式如下的上三角矩阵，得到 \(x_n\) 的值，然后依次代入得到整个方程的解：
\[\left(
\begin{array}{ccc|c}
\times& \times& \times& \times\\
\times& \times& \times& \times\\
\times& \times& \times& \times
\end{array}
\right)
\Rightarrow
\left(
\begin{array}{ccc|c}
\times& \times& \times& \times\\
0& \times& \times& \times\\
0& 0& \times& \times
\end{array}
\right)
\Rightarrow
\left(
\begin{array}{ccc|c}
\times& 0& 0& \times\\
0& \times& 0& \times\\
0& 0& \times& \times
\end{array}
\right)
\]
考虑每次对一个未知数进行消元（即一列），设当前未知数为 \(x_i\)（第 \(i\) 列），我们按照以下步骤消元：
1. 无解判断
若这一列从第 \(i\) 行往下都是 \(0\)，说明这个未知数无解或为任意解，判无解。
2. 选主元
选择这一列从第 \(i\) 行往下系数绝对值最大的那一个作为主元，将这一行与第 \(i\) 行交换
比如：
\[\
\left(
\begin{array}{ccc|c}
3& 7& 8& 14\\
6& -1& -2& 7\\
-12& 1& -1& -10
\end{array}
\right)
\Rightarrow
\left(
\begin{array}{ccc|c}
\underline{3}& 7& 8& 14\\
\underline{6}& -1& -2& 7\\
\boxed{-12}& \underline{1}& \underline{-1}& \underline{-10}
\end{array}
\right)
\Rightarrow
\left(
\begin{array}{ccc|c}
-12& 1& -1& -10\\
6& -1& -2& 7\\
3& 7& 8& 14
\end{array}
\right)
\]
这一操作的目的是减少误差（原因见下文）。
3. 消元
对于第 \(j\) 行（\(j>i\)），为了使 \(a_{j,i}\) 消为 \(0\) ，我们需要将第 \(j\) 行减去第 \(i\) 行（主元行）的 \(\frac{a_{j,i}}{a_{i,i}}\) 倍，这样 \(a_{j,i}-a_{i,i}\times\frac{a_{j,i}}{a_{i,i}}=0\)。
注意，消第 \(i\) 个未知数时，只需要将第 \(i\) 列第 \(i\) 行下面的数消为 \(0\) 即可。
\[\begin{aligned}
&\left(
\begin{array}{ccc|c}
-12& 1& -1& -10\\
6& -1& -2& 7\\
3& 7& 8& 14
\end{array}
\right)\\
\Rightarrow
&\left(
\begin{array}{ccc|c}
\color{red}{-12}& \color{red}{1}& \color{red}{-1}& \color{red}{-10}\\
6-(\textcolor{red}{-12})\times\textcolor{blue}{\frac{6}{-12}}& -1-\textcolor{red}{1}\times\textcolor{blue}{\frac{6}{-12}}& -2-(\textcolor{red}{-1})\times\textcolor{blue}{\frac{6}{-12}}& 7-(\textcolor{red}{-10})\times\textcolor{blue}{\frac{6}{-12}}\\
3-(\textcolor{red}{-12})\times\textcolor{green}{\frac{3}{-12}}& 7-\textcolor{red}{1}\times\textcolor{green}{\frac{3}{-12}}& 8-(\textcolor{red}{-1})\times\textcolor{green}{\frac{3}{-12}}& 14-(\textcolor{red}{-10})\times\textcolor{green}{\frac{3}{-12}}
\end{array}
\right)\\
\Rightarrow
&\left(
\begin{array}{ccc|c}
-12& 1& -1& -10\\
\underline{0}& -0.5& -2.5& 2\\
\underline{0}& 7.25& 7.75& 11.5
\end{array}
\right)
\end{aligned}
\]
从左往右对每个未知数（每一列）进行以上三步操作，得到上三角矩阵：
\[\left(
\begin{array}{ccc|c}
-12& 1& -1& -10\\
0& 7.25& 7.75& 11.5\\
0& 0& -\frac{114}{58}& \frac{81}{29}
\end{array}
\right)
\]
4. 得到答案
将上三角矩阵写成方程式，有
\[\left\{
\begin{aligned}
-12x_1+x_2-x_3&=-10\\
7.25x_2+7.75x_3&=11.5\\
-\frac{114}{58}x_3&=\frac{81}{29}
\end{aligned}
\right.
\]
发现最后一个未知数的解我们已经知道了，将解依次代回就可以得到方程组的解了。
\[\left\{
\begin{aligned}
x_1&=\frac{-10-x_2+x_3}{-12}=1.210526\\
x_2&=\frac{11.5-7.75x_3}{7.25}=3.105263\\
x_3&=-1.421053\\
\end{aligned}
\right.
\]
5.误差问题
考虑这样一组数据：
\[\left\{
\begin{aligned}
0.3 \times 10^{-11}x_1+x_2&=0.7\\
x_1+x_2&=0.9
\end{aligned}
\right.
\]
如果以第一行为主元进行消元，则有：
\[\left(
\begin{array}{cc|c}
0.3\times 10^{-11}& 1& 0.7\\
0& 1-\frac{10}{3}\times10^{11}\times1& 0.9-\frac{10}{3}\times10^{11}\times0.7
\end{array}
\right)
\]
会得到：
\[\left\{
\begin{aligned}
x_1&=0.000000\\
x_2&=0.700000
\end{aligned}
\right.
\]
但如果以第二行为主元消元，则有：
\[\left(
\begin{array}{cc|c}
1& 1& 0.9\\
0& 1-0.3\times10^{-11}\times1& 0.7-0.3\times10^{-11}\times0.9
\end{array}
\right)
\]
会得到：
\[\left\{
\begin{aligned}
x_1&=0.200000\\
x_2&=0.700000
\end{aligned}
\right.
\]
明显更接近正解。
所以，在高斯消元中，应选择系数绝对值最大的一行作为主元行，使 \(\frac{a_{j,i}}{a_{i,i}}\) 中的 \(a_{i,i}\) 绝对值更大，原分数值更小，减少计算时的误差。
时间复杂度 \(O(n^3)\)。
模版代码
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
const double eps=1e-7;
int n;
double a[105][105],ans[105];
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)
{
for(int j=1;j<=n+1;j++)
{
scanf("%lf",&a[i][j]);
}
}
for(int i=1;i<=n;i++)//枚举消去的未知数
{
//找主元行
int maxn=i;
for(int j=i+1;j<=n;j++)
{
if(fabs(a[j][i])>fabs(a[maxn][i])) maxn=j;
}
//判无解
if(fabs(a[maxn][i])<eps)
{
printf("No Solution\n");
return 0;
}
//将主元行交换到第 i 行
for(int j=1;j<=n+1;j++)
{
swap(a[i][j],a[maxn][j]);
}
//消元
for(int j=i+1;j<=n;j++)
{
for(int k=n+1;k>=i;k--)//注意倒序，否则 a[j][i] 的值会被提前修改
{
a[j][k]-=a[i][k]*(a[j][i]/a[i][i]);
}
}
}
//代回，得到答案
for(int i=n;i>=1;i--)
{
ans[i]=a[i][n+1];
for(int j=n;j>i;j--) ans[i]-=a[i][j]*ans[j];
ans[i]/=a[i][i];
}
for(int i=1;i<=n;i++)
{
printf("%.2lf\n",ans[i]);
}
return 0;
}
2. 高斯-约旦消元法
别看名字很高大上，实际上就是高斯消元法的一个小优化。
刚刚的消元过程中，我们的增广矩阵的变化为：
\[\left(
\begin{array}{ccc|c}
\times& \times& \times& \times\\
\times& \times& \times& \times\\
\times& \times& \times& \times
\end{array}
\right)
\Rightarrow
\left(
\begin{array}{ccc|c}
\times& \times& \times& \times\\
0& \times& \times& \times\\
0& 0& \times& \times
\end{array}
\right)
\Rightarrow
\left(
\begin{array}{ccc|c}
\times& 0& 0& \times\\
0& \times& 0& \times\\
0& 0& \times& \times
\end{array}
\right)
\]
而高斯-约旦消元法则是跳过上三角矩阵，直接得到答案。
\[\left(
\begin{array}{ccc|c}
\times& \times& \times& \times\\
\times& \times& \times& \times\\
\times& \times& \times& \times
\end{array}
\right)
\Rightarrow
\left(
\begin{array}{ccc|c}
1& 0& 0& \times\\
0& 1& 0& \times\\
0& 0& 1& \times\\
\end{array}
\right)
\]
具体来讲，对于第 \(i\) 个未知数，在选完主元行后，先将主元行整体除以 \(a_{i,i}\) 使 \(a_{i,i}\)（即对角线）等于 \(1\)，然后将每一行都减去主元行的 \(a_{j,i}\) 倍即可。
时间复杂度 \(O(n^3)\)。
模版代码
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
const double eps=1e-7;
int n;
double a[105][105],ans[105];
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)
{
for(int j=1;j<=n+1;j++)
{
scanf("%lf",&a[i][j]);
}
}
for(int i=1;i<=n;i++)//枚举消去的未知数
{
//找主元行
int maxn=i;
for(int j=i+1;j<=n;j++)
{
if(fabs(a[j][i])>fabs(a[maxn][i])) maxn=j;
}
//判无解
if(fabs(a[maxn][i])<eps)
{
printf("No Solution\n");
return 0;
}
//将主元行交换到第 i 行
for(int j=1;j<=n+1;j++)
{
swap(a[i][j],a[maxn][j]);
}
//将 a[i][i] 消为 1
for(int j=n+1;j>=i;j--) a[i][j]/=a[i][i];
//消元
for(int j=1;j<=n;j++)//注意！！从 1 开始，对每一行进行消元
{
if(j!=i)
for(int k=n+1;k>=i;k--)//注意倒序，否则 a[j][i] 的值会被提前修改
{
a[j][k]-=a[i][k]*a[j][i];
}
}
}
//直接得到答案
for(int i=1;i<=n;i++)
{
printf("%.2lf\n",a[i][n+1]);
}
return 0;
}
3.异或方程组
异或方程组是形式如下的方程组，其中的 \(a\)，\(b\)，\(x\) 均为 \(0\) 或 \(1\):
\[\left\{
\begin{aligned}
&a_{1,1}x_1 \operatorname{xor}
a_{1,2}x_2
\operatorname{xor}
\dots
\operatorname{xor}
a_{1,n}x_n=b_1\\
&a_{2,1}x_1 \operatorname{xor}
a_{2,2}x_2
\operatorname{xor}
\dots
\operatorname{xor}
a_{2,n}x_n=b_2\\
&\dots\\
&a_{n,1}x_1 \operatorname{xor}
a_{n,2}x_2
\operatorname{xor}
\dots
\operatorname{xor}
a_{n,n}x_n=b_n
\end{aligned}
\right.
\]
回忆上文，之所以线性方程组的增广矩阵可以高斯消元，是因为它可以：
- 交换任意两个方程的位置；
- 将一个方程加上另一个方程；
- 将整个方程乘上一个数。
推广到异或方程组，由于异或运算有交换律和结合律，所以它可以：
- 交换任意两个方程的位置；
- 将一个方程异或上一个方程；
- 将整个方程异或上一个数。
因此，异或方程组可以通过高斯消元解决。
使用 std::bitset 进行优化，复杂度为 \(O(\frac{n^3}{\omega})\)。
模版代码
int n;
bitset<N> a[N];
vector<bool> gauss()
{
for(int i=1;i<=n;i++)
{
int maxn=i;
for(int j=i+1;j<=n;j++)
{
if(a[j][i]>a[maxn][i]) maxn=j;
}
if(a[maxn][i]==0) return vector<bool>(0);
if(i!=maxn) swap(a[i],a[maxn])
for(int j=1;j<=n;j++)
{
if(j!=i&&a[j][i]) a[j]^=a[i];
}
}
vector<bool> ans(n+1);
for(int i=1;i<=n;i++) ans[i]=a[i][n+1];
return ans;
}
浙公网安备 33010602011771号