高斯消元¶

约 543 个字 77 行代码预计阅读时间 3 分钟

求解线性方程组和异或方程组，求矩阵的逆

用的是高斯约旦消元法，比普通高斯消元法精度更高，代码更简洁

时间复杂度¶

$O(n^3)$

稀疏矩阵可以观察稀疏的地方并优化

线性方程组¶

形式如下

$\begin{cases}a_{11}x_1+a_{12}x_2+...+a_{1n}x_n=b_1\\ a_{21}x_1+a_{22}x_2+...+a_{2n}x_n=b_2\\...\\a_{n1}x_1+a_{n2}x_2+...+a_{nn}x_n=b_n\end{cases}$

写成矩阵形式，AX=B

$A=\begin{vmatrix}a_{11}&a_{12}&...&a_{1n}\\a_{21}&a_{22}&...&a_{2n}\\...\\a_{n1}&a_{n2}&...&a_{nn}\end{vmatrix},X=\begin{vmatrix}x_1\\x_2\\...\\x_n\end{vmatrix},B=\begin{vmatrix}b_1\\b_2\\...\\b_n\end{vmatrix}$

笔算如何求解在”线性代数“这门课里已经讲过了，不再叙述

程序求解就是增广矩阵行初等行变换为行最简形：

$\begin{vmatrix}a_{11}&a_{12}&...&a_{1n}&b_1\\a_{21}&a_{22}&...&a_{2n}&b_2\\...\\a_{n1}&a_{n2}&...&a_{nn}&b_n\end{vmatrix}\rightarrow\begin{vmatrix}1&0&...&0&c_1\\0&1&...&0&c_2\\...\\0&0&...&1&c_n\end{vmatrix}$

这样 $X$ 就等于

$\begin{vmatrix}c_1\\c_2\\...\\c_n\end{vmatrix}$

矩阵求逆¶

逆矩阵在“线性代数”已讲过，不再叙述

对 $n$ 阶矩阵 $A$ ，构造 $n\times 2n$ 的矩阵 $(A,I_n)$

高斯消元化成 $(I_n,A^{-1})$

如果左半部分不能化成单位矩阵，则矩阵 $A$ 不可逆

异或方程组¶

异或运算满足交换律和结合律，所以可以使用类似高斯消元的方法解决方程组

形式如下,且所有系数（即a和b）均为0或1

$\begin{cases}a_{11}x_1\oplus a_{12}x_2\oplus ...\oplus a_{1n}x_n=b_1\\ a_{21}x_1\oplus a_{22}x_2\oplus ...\oplus a_{2n}x_n=b_2\\...\\a_{n1}x_1\oplus a_{n2}x_2\oplus ...\oplus a_{nn}x_n=b_n\end{cases}$

应使用“异或消元”而非“加减消元”，且不需要进行乘除改变系数（因为系数均为0和1）。

注意到异或方程组的增广矩阵是01矩阵，可以用bitset优化

注意处理无解，无穷解情况

-1无解，0无穷解，1唯一解

线性方程组¶

const int N=1e2+1;
const double eps=1e-6;
double a[N][N],ans[N];
int GaussianElimination(int n)
{
    int i,maxx,j,k,cnt=0;
    double temp;
    for(i=0;i<n;++i)
    {
        maxx=cnt;
        for(j=cnt+1;j<n;++j)
        {
            if(fabs(a[j][i])>fabs(a[maxx][i]))
                maxx=j;
        }
        if(fabs(a[maxx][i])<eps)
            continue;
        if(maxx!=cnt)
            swap(a[maxx],a[cnt]);
        for(j=0;j<n;++j)
        {
            if(j==cnt)
                continue;
            temp=a[j][i]/a[cnt][i];
            for(k=i+1;k<=n;++k)
                a[j][k]-=a[cnt][k]*temp;
        }
        ++cnt;
    }
    if(cnt<n)
    {
        for(;cnt<n;++cnt)
        {
            if(fabs(a[cnt][n])>eps)
                return -1;
        }
        return 0;
    }
    for(i=0;i<n;++i)
        ans[i]=a[i][n]/a[i][i];
    return 1;
}

异或方程组¶

const int N=1e2+3;
bitset<N> a[N];
bitset<N> ans;
int GaussElimination(int n, int m)
{
    int i,j,cnt=0;
    for(i=0;i<n;++i)
    {
        for(j=cnt;j<m && !a[j][i];++j){}
        if(j==m)
            continue;
        if(j != cnt)
            swap(a[j], a[cnt]);
        for(j=0;j<m;j++)
        {
            if(cnt==j)
                continue;
            if(a[j][i])
                a[j]^=a[cnt];
        }
        ++cnt;
    }
    if(cnt<n)
    {
        for(;cnt<n;++cnt)
        {
            if(a[cnt][n])
                return -1;
        }
        return 0;
    }
    for (i = 0; i < n; i++)
        ans[i] = a[i][n];
    return 1;
}