高斯消元 学习笔记

模板及讲解

高斯消元是什么

求解多元一次方程组的一个消元方法,主要借助矩阵来完成。

例题

Luogu 2455

我们发现给定的方程组解只和系数有关,所以将系数和常数提取出来构造矩阵:

$$
\begin{cases}
2x+y+z=1\\
6x+2y+z=-1\\
-2x+2y+z=7
\end{cases}
$$
则构造矩阵:
$$
\left[
\begin{array}{ccc|c}
2 & 1 & 1&1 \\
6 & 2 & 1&-1\\
-2&2&1&7
\end{array}
\right]
$$

高斯消元的操作有(初等行变换):
1、交换两行。
2、将一行用非零数乘。
3、将一行的若干倍加到另一行上。

也就是差不多是小学学的加减消元带入消元的操作。

高斯消元目标:
将原矩阵转化为一个阶梯形矩阵,即
$$
\left[
\begin{array}{ccc|c}
2 & 1 & 1&1 \\
0 & -1& -2&-4\\
0&0&-4&-4
\end{array}
\right]
$$
将后面的行带入上面的行,可得一个简化阶梯形矩阵
$$
\left[
\begin{array}{ccc|c}
1 & 0 & 0&-1 \\
0 & 1& 0&2\\
0&0&1&1
\end{array}
\right]
$$
直接可以读出每个元素的解。

高斯消元一般步骤:
逐步处理每一个$x_i$,将$i$列除了$i$行的其他系数都变为0,即我们可以先保证$i$行$i$列不为0,然后以这一行为基准将其他行$i$位置的系数都变为0。用初等行变换的交换可以满足第一个条件,用初等行变换将一行的若干倍加到另一行上可以做到第二个条件。做完之后只有对角线和常数有值,每一行的常数除以对角线的值即为第$i$个元素的解。

无解情况:系数都为0,常数不为0
无穷解情况:系数都为0,常数也为0

两个概念:
1、主元:确定唯一的未知量
2、自由元:不确定值未知量,对应情况为该行全为0,有无穷解

注意在判这两个之前要先判有没有不是对角线上的系数有值,如果有则这一行不进行检查。

#include<cstdio> 
#include<cstring>
#include<algorithm>
#include<iostream>
#include<vector>
#include<queue>
#include<cmath>
#define ms(i, j) memset(i, j, sizeof i)
#define LL long long
#define db double
#define fir first
#define sec second
#define mp make_pair
using namespace std;

namespace flyinthesky {

    const int MAXN = 105;
    const db eps = 1e-12; 

    int n;
    db c[MAXN][MAXN], b[MAXN];
    // c, b 为消元矩阵 

    void clean() {
        ms(b, 0);
    }
    int solve() {
        clean();
        scanf("%d", &n);
        for (int i = 1; i <= n; ++i) {
            for (int j = 1; j <= n; ++j) scanf("%lf", &c[i][j]);
            scanf("%lf", &b[i]);
        }
        for (int i = 1; i <= n; ++i) { // 处理 xi[i] 的系数 
            for (int j = i; j <= n; ++j) { // 从后面没处理的行找 xi[i] 系数不为 0 的行交换 
                if (fabs(c[j][i]) > eps) {
                    swap(b[i], b[j]);
                    for (int k = 1; k <= n; ++k) swap(c[i][k], c[j][k]);
                }
            }
            if (fabs(c[i][i]) < eps) continue ; // 注意 
            for (int j = 1; j <= n; ++j) { // 消去其他方程 xi[i] 的系数 
                if (i == j) continue ;
                db rat = c[j][i] / c[i][i];
                for (int k = 1; k <= n; ++k) c[j][k] -= rat * c[i][k];
                b[j] -= rat * b[i];
            }
        }
        int bj = 1;
        for (int i = 1; i <= n; ++i) {
            int fl = 0;
            for(int j = 1; j <= n; ++j){
                if (i != j && fabs(c[i][j]) > eps) {fl = 1; break ;}
            }
            if (fl) continue ; // 必须判 
            if (fabs(c[i][i]) < eps && fabs(b[i]) > eps) bj = -1;
            if (fabs(c[i][i]) < eps && fabs(b[i]) < eps && bj != -1) bj = 0; 
        }
        if (bj != 1) return printf("%d\n", bj), 0; 
        for (int i = 1; i <= n; ++i) {
            if (fabs(b[i] / c[i][i]) < eps) printf("x%d=0\n", i); else printf("x%d=%.2lf\n", i, b[i] / c[i][i]);
        }
        return 0;
    }
}
int main() {
    flyinthesky::solve();
    return 0;
}

异或方程组

高斯消元可以用来解线性的异或方程组,方法与前者类似。

例题:poj 1830

设$x_i$为$i$开关是否(0/1)操作,$a_{i,j}$为$j$开关操作后是否(1/0)会影响$i$开关

则可以列出线性异或方程组

$$
\begin{cases}
a_{1,1}x_1 xor a_{1,2}x_2 xor … xor a_{1,n}x_n = st_1 xor ed_1 \\
a_{2,1}x_1 xor a_{2,2}x_2 xor … xor a_{2,n}x_n = st_2 xor ed_2 \\
\vdots \\
a_{n,1}x_1 xor a_{n,2}x_2 xor … xor a_{n,n}x_n = st_n xor ed_n
\end{cases}
$$

然后对于异或实际上我们也可以进行消元得到简化阶梯形矩阵,并且不需要乘除,因为只有 0 和 1。
可以用 bitset 加速运算且代码更简洁高效。

无解即为形如$0 xor 0 xor … xor 0 =1$的,显然不成立。

#include<cstdio> 
#include<cstring>
#include<algorithm>
#include<iostream>
#include<vector>
#include<queue>
#include<bitset>
#define ms(i, j) memset(i, j, sizeof i)
#define LL long long
#define db double
#define fir first
#define sec second
#define mp make_pair
using namespace std;

namespace flyinthesky { 

    const int MAXN = 35;

    int n, a[MAXN];
    bitset<MAXN> bs[MAXN];

    int ksm(int a, int b) {
        int bs = a, ans = 1;
        while (b) {
            if (b & 1) ans *= bs;
            bs *= bs;
            b >>= 1;
        }
        return ans;
    }

    void clean() {
    }
    int solve() {
        clean();
        cin >> n;
        for (int i = 1; i <= n; ++i) scanf("%d", &a[i]);
        for (int x, i = 1; i <= n; ++i) scanf("%d", &x), a[i] ^= x;
        for (int i = 1; i <= n; ++i) bs[i].reset(), bs[i][0] = a[i], bs[i][i] = 1;
        int x, y;
        while (scanf("%d%d", &x, &y) == 2 && x && y) bs[y][x] = 1;
        int ans = 1, whw = 0;
        for (int i = 1; i <= n; ++i) { // 每列
             if (bs[i][i] == 0) { // 找一个非零的
                 for (int j = i + 1; j <= n; ++j) if (bs[j][i]) {swap(bs[j], bs[i]); break;}
             }
             if (bs[i][0] == 1 && bs[i].count() == 1) {ans = 0; break ;} // 无解
             if (bs[i][i] == 0) ++whw; // 自由元
             for (int j = i + 1; j <= n; ++j) if (bs[j][i]) bs[j] ^= bs[i];
        }
        if (ans == 0) printf("Oh,it's impossible~!!\n"); else printf("%d\n", ksm(2, whw));
        return 0; 
    }
}
int main() {
    int k; cin >> k;
    while (k--) flyinthesky::solve();
    return 0;
}

求解带环期望DP

1、CF 24D:列与列之间不满足无后效性 (注意此时高斯消元做成阶梯矩阵即可,不需要简化,简化直接回带)

线性基

定义

一下内容部分来自 线性基学习笔记 | Menci’s Blog

异或和:设 $S$ 为一无符号整数集(下文中除非特殊说明,集合均指无符号整数集),定义其异或和为$S_1 xor S_2 xor \cdots xor S_{|S|}$。

线性相关:集合中存在一个元素$S_j$,使得集合中其他元素可以通过异或和的方式得到$S_j$,则$S$线性相关。若不存在,则$S$线性无关。

线性基:称集合 $B$ 是集合 $S$ 的线性基,当且仅当:

  • $B$ 线性无关
  • $B$ 中元素可以通过异或和方式得到$S$中的所有值

线性基性质

  • $B$ 是极小的满足线性基性质的集合,它的任何真子集都不可能是线性基
  • $S$ 中的任意元素都可以唯一表示为 $B$ 中若干个元素异或起来的结果。

构造

构造出来的线性基相当于于所有插入值高斯消元后的简化阶梯矩阵。

设维护的数的二进制长度为$L$,线性基数组为$a_i(0 \leq i \leq L)$
为了保证线性无关,$a_i$的第$i$位若为$1$,则其他线性基上第$i$位都为$0$,并且$a_i$高于$i$位上的值都必须为$0$(否则会影响后面线性基的”其他线性基上第$i$位都为$0$”)。

所以构造方法即为(插入数$t$):
1、将$t$倒序枚举二进制位$i$,然后
2、若$t$的$i$位为0,则continue
否则:

  • 如果$a_i$不为$0$,则用$a_i$消去$t$的$i$位上的二进制$1$,然后continue
  • 如果$a_i$为$0$,那么$t$就能插进$i$位。但是要先用$[1, i)$的线性基消去$t$的相应位置的$1$,以保证”其他线性基上第$i$位都为$0$”,然后再用得到的线性基消去$[i+1, L]$线性基的$i$位上的$1$,因为一路上我们都将高位消去,所以这里不用担心影响其他位的情况。

例题:Luogu 3812

给定$n$个整数(数字可能重复),求在这些数中选取任意个,使得他们的异或和最大。

对于构造的线性基,$a_i$的第$i$位若是$1$,则将答案异或$a_i$,这样一定能得到更大的值,因为这个位置的$1$之后不会再被异或消去,而因为二进制的性质,贪心选取更高位数的$1$,所以将所有线性基异或起来即为答案。

#include<cstdio> 
#include<cstring>
#include<algorithm>
#include<iostream>
#include<vector>
#include<set>
#include<cmath>
#define ms(i, j) memset(i, j, sizeof i)
#define LL long long
#define db double
#define fir first
#define sec second
#define mp make_pair
using namespace std;

namespace flyinthesky {

    const LL MAXLEN = 51;

    LL n, t, a[55];

    void clean() {
        ms(a, 0);
    }
    int solve() {

        clean();

        scanf("%lld", &n);

        for (LL o = 1; o <= n; ++o) {
            scanf("%lld", &t);
            for (LL i = MAXLEN; i >= 0; --i) { // 倒序枚举二进制
                if (!((t >> i) & 1)) continue ; // t 当前位不是 1
                if (a[i]) t ^= a[i]; // 用 a[i] 消去 t 的 i 位上的 1
                  else { 
                    for (LL j = 1; j < i; ++j) if ((t >> j) & 1) t ^= a[j]; // 用 a[j] 消去 t 低于 i 位上的 1
                    for (LL j = i + 1; j <= MAXLEN; ++j) if ((a[j] >> i) & 1) a[j] ^= t; // 用 t 消去 a[j] 的 i 位上的 1
                    a[i] = t; // 插入线性基
                    break ;
                }
            }
        }
        LL ans = 0;
        for (LL i = MAXLEN; i >= 0; --i) ans ^= a[i];

        printf("%lld\n", ans);

        return 0;
    }
}
int main() {
    flyinthesky::solve();
    return 0;
}
------ 本文结束 ------