高斯消元

1.高斯消元解线性方程组。这个模板只处理n方程n个未知数的，分别处理无解，无穷多解，以及唯一解的情况。结果存在最后一列中（一开始的常数项），

  • 无解返回0
  • 唯一解返回1
  • 无穷多解返回2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
int gauss(vector<vector<db>>& a, int n) {
    int r = 1;                      // 当前行
    for (int c = 1; c <= n; c++) {  // 消元进行到第c列
        // 1.找到c列的最大行t
        int t = r;
        for (int i = r; i <= n; i++)
            if (fabs(a[i][c]) > fabs(a[t][c]))
                t = i;
        if (fabs(a[t][c]) < eps)
            continue;  // c列已全为0

        // 2.把最大行换到上面
        for (int i = c; i <= n + 1; i++) swap(a[t][i], a[r][i]);

        // 3.把当前行r的第一个数，变成1
        for (int i = n + 1; i >= c; i--) a[r][i] /= a[r][c];

        // 4.把当前列c下面的所有数，全部消成0
        for (int i = r + 1; i <= n; i++)
            if (fabs(a[i][c]) > eps)
                for (int j = n + 1; j >= c; j--)
                    a[i][j] -= a[i][c] * a[r][j];
        r++;  // 从下一行开始消元下一列
    }
    if (r <= n) {  // 说明已经提前变成梯形矩阵
        for (int i = r; i <= n; i++) {
            if (fabs(a[i][n + 1]) > eps)
                return 0;
        }  // 左边=0，右边≠0,无解
        return 2;  // 0==0，无穷多解
    }
    // 5.唯一解，从下往上回代，得到方程的解
    for (int i = n; i >= 1; i--)
        for (int j = i + 1; j <= n; j++)
            a[i][n + 1] -= a[i][j] * a[j][n + 1];
    return 1;
}

void solve() {
    int n;
    cin >> n;
    vector b(n + 1, vector<db>(n + 2));
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n + 1; j++)
            cin >> b[i][j];
    int t = gauss(b, n);
    if (t == 0) {
        cout << "No solution " << endl;
    } else if (t == 2) {
        cout << "Infinite group solutions" << endl;
    } else {
        for (int i = 1; i <= n; i++) {
            baoliu(b[i][n + 1], 2);
            cout << endl;
        }
    }
}

异或线性方程组，基本结构和上面一样，在系数变换上面方便许多

解释：异或线性方程组就是常数项和各变量系数都是0/1，并且各变量取值也是0/1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
 gauss(vector<vector<int>>& a, int n, vector<int>& solution) {
    vector<int> freevar;           // 记录自由变量对应的列
    vector<int> pivot(n + 1, -1);  // 记录每一行的主元所在的列
    //solution.assign(n + 1, 0);
    int r = 1;
    for (int c = 1; c <= n; c++) {
        int t = r;
        // 找到当前列中的主元
        for (int i = r; i <= n; i++) {
            if (a[i][c]) {
                t = i;
                break;
            }
        }
        if (!a[t][c]) {
            freevar.push_back(c);
            continue;  // 当前列没有主元，继续到下一列
        }
        pivot[r] = c;  // 第 r 行的主元在 c 列
        if (t != r) {  // 交换行，将主元行放在第 r 行
            for (int i = c; i <= n + 1; i++)
                swap(a[r][i], a[t][i]);
        }
        // 消去主元下方的所有行
        for (int i = r + 1; i <= n; i++) {
            if (a[i][c])
                for (int j = n + 1; j >= c; j--) a[i][j] ^= a[r][j];
        }
        r++;
    }

    // 检查是否有解
    for (int i = r; i <= n; i++) {
        if (a[i][n + 1])
            return 0;  // 无解
    }
    //int tot = 0;
    int rksz = r - 1;  // 这是系数矩阵的秩
    // 自由变量根据题目要求情况去赋值
    for (auto i : freevar) solution[i] =0;
    //------------------------------
    for (int i = rksz; i >= 1; i--) {
        int sum = a[i][n + 1];
        for (int j = 1; j <= n; j++) {
            if (j == pivot[i]) {
                continue;
            }  // 如果不是主元所在的列
            sum ^= (a[i][j] * solution[j]);  // 右边已经求出来的，左边自由变量遗留
        }
        solution[pivot[i]] = sum;  // 求解对应的主元变量
    }
    assert(rksz <= n);
    if (rksz < n)
        return 2;  // 无穷多解
    return 1;      // 唯一解
}
// int t = gauss(b, n, sol);

解决n个方程m个方程AX=BAX=B问题的模板：处理了无穷多解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
db b[N][N];  // 增广矩阵
int id[N];   // id[j]=i：表示第j列的答案最终在第i行被计算
int rid[N];  // 第i行可以算出第j列的主元
void print(int n, int m) {
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= m + 1; j++) {
            cerr << setiosflags(ios::left) << setw(5) << b[i][j];
        }
        cerr << endl;
    }
    cerr << "---------------------" << endl;
}
int rksz;
int freenum;
vector<db> ans;
int gauss(db a[][N], int n, int m) {  // n个方程，m个未知数。默认多余自由变量为0，记录映射关系
    int r = 1;                        // 当前行
    for (int c = 1; c <= m; c++) {    // 消元进行到第c列
        // 1.找到c列的最大行t
        int t = r;
        for (int i = r; i <= n; i++)
            if (fabs(a[i][c]) > fabs(a[t][c]))
                t = i;
        if (t > n || fabs(a[t][c]) < eps)
            continue;  // c列已全为0
        assert(t <= n);
        // 2.把最大行换到上面
        for (int i = c; i <= m + 1; i++) swap(a[t][i], a[r][i]);
        // 3.把当前行r的第一个数，变成1
        for (int i = m + 1; i >= c; i--) a[r][i] /= a[r][c];
        // 4.把当前列c下面的所有数，全部消成0
        for (int i = r + 1; i <= n; i++)
            if (fabs(a[i][c]) > eps)
                for (int j = m + 1; j >= c; j--)
                    a[i][j] -= a[i][c] * a[r][j];
        id[c] = r;
        rid[r] = c;
        r++;  // 从下一行开始消元下一列
    }
    //---------------------------
    // print(n, m);
    rksz = r - 1;
    freenum = m - r + 1;
    ans.resize(m + 1);
    //------------------------
    if (r <= m) {  // 说明已经提前变成梯形矩阵
        for (int i = r; i <= n; i++) {
            if (fabs(a[i][m + 1]) > eps)
                return 0;  // 左边=0，右边≠0,无解
        }
    }
    for (int i = 1; i <= m; i++) {
        if (id[i] == 0) {
            // deb(i);
            ans[i] = 1;  // 如果第i列的主元没有对应行，自由变量随机赋值
        }
    }
    // 5.唯一解，从下往上回代，得到方程的解
    for (int i = rksz; i >= 1; i--) {
        for (int j = 1; j <= m; j++) {  // 左侧自由变量残余，右侧已经算出来的以及右侧自由变量
            if (j == rid[i])
                continue;
            a[i][m + 1] -= a[i][j] * ans[j];
        }
        // deb(rid[i]);
        ans[rid[i]] = a[i][m + 1];  // 第i行的主元在rid[i]列
    }
    if (m > n) {
        return 2;
    } else {
        // m<=n;
        assert(rksz <= m);
        if (rksz == m)
            return 1;
        return 2;
    }
}
void solve() {
    int n, m;
    cin >> n >> m;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= m; j++) {
            cin >> b[i][j];
        }
    }
    for (int i = 1; i <= n; i++) cin >> b[i][m + 1];
    int t = gauss(b, n, m);  // n个方程,m个未知数
    deb(t, rksz, freenum);
    if (t == 0) {
        cout << "NO" << endl;
    } else {
        cout << "YES" << endl;
        for (int i = 1; i <= m; i++) cout << fixed << setprecision(12) << ans[i] << endl;
    }
}

n行m列的高斯消元异或

有重大越界问题，现在已经更改。就是方程数比未知数多的时候，我们需要对于提前到达最后一行的情况判断，将后面的变量全部放入自由变量

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
int rksz;
// 还没有处理n方程，m变量的情况(n>m)
int gauss(vector<vector<int>>& a, int n, int m, vector<int>& solution) {
    vector<int> freevar;           // 记录自由变量对应的列
    vector<int> pivot(n + 1, -1);  // 记录每一行的主元所在的列
    int r = 1;
    for (int c = 1; c <= m; c++) {
        int t = r;
        // 找到当前列中的主元
        for (int i = r; i <= n; i++) {
            if (a[i][c]) {
                t = i;
                break;
            }
        }
        // if (t > n)
        //     deb(t, n, c, m);
        // assert(t <= n);
        // assert(c <= m);
        if (t > n || !a[t][c]) {
            freevar.push_back(c);
            continue;  // 当前列没有主元，继续到下一列
        }
        pivot[r] = c;  // 第 r 行的主元在 c 列
        if (t != r) {  // 交换行，将主元行放在第 r 行
            for (int i = c; i <= m + 1; i++)
                swap(a[r][i], a[t][i]);
        }
        // 消去主元下方的所有行
        for (int i = r + 1; i <= n; i++) {
            if (a[i][c])
                for (int j = m + 1; j >= c; j--) a[i][j] ^= a[r][j];
        }
        r++;
    }

    // 检查是否有解
    for (int i = r; i <= n; i++) {
        if (a[i][m + 1])
            return 0;  // 无解
    }
    // int tot = 0;
    rksz = r - 1;  // 这是系数矩阵的秩
    // 自由变量根据题目要求情况去赋值
    for (auto i : freevar) solution[i] = 0;
    //------------------------------
    for (int i = rksz; i >= 1; i--) {
        int sum = a[i][m + 1];
        for (int j = 1; j <= m; j++) {
            if (j == pivot[i]) {
                continue;
            }  // 如果不是主元所在的列
            sum ^= (a[i][j] * solution[j]);  // 右边已经求出来的，左边自由变量遗留
        }
        solution[pivot[i]] = sum;  // 求解对应的主元变量
    }
    // assert(rksz <= m);
    if (rksz < m)
        return 2;  // 无穷多解
    return 1;      // 唯一解
}
// int t = gauss(b, n，m, sol);

bitset优化版本加速到O(n2m/w)O(n^2m/w)https://qoj.ac/contest/1716/problem/313

5002是最多5000个自变量，并且多留一位给常数项

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
int rksz;
// 还需要处理m方程，n变量，m>n
#define bit(x) bitset<(x)>
int gauss(vector<bit(5002)>& a, int n, int m, vector<int>& solution) {
    vector<int> freevar;           // 记录自由变量对应的列
    vector<int> pivot(n + 1, -1);  // 记录每一行的主元所在的列
    int r = 1;
    for (int c = 1; c <= m; c++) {
        int t = r;
        // 找到当前列中的主元
        for (int i = r; i <= n; i++) {
            if (a[i][c]) {
                t = i;
                break;
            }
        }
        if (t > n || !a[t][c]) {
            freevar.push_back(c);
            continue;  // 当前列没有主元，继续到下一列
        }
        pivot[r] = c;  // 第 r 行的主元在 c 列
        if (t != r) {  // 交换行，将主元行放在第 r 行
            swap(a[r], a[t]);
        }
        // 消去主元下方的所有行
        for (int i = r + 1; i <= n; i++) {
            if (a[i][c])
                a[i] ^= a[r];
        }
        r++;
    }

    // 检查是否有解
    for (int i = r; i <= n; i++) {
        if (a[i][m + 1])
            return 0;  // 无解
    }
    // int tot = 0;
    rksz = r - 1;  // 这是系数矩阵的秩
    // 自由变量根据题目要求情况去赋值
    for (auto i : freevar) solution[i] = 0;
    //------------------------------
    for (int i = rksz; i >= 1; i--) {
        int sum = a[i][m + 1];
        for (int j = 1; j <= m; j++) {
            if (j == pivot[i]) {
                continue;
            }  // 如果不是主元所在的列
            sum ^= (a[i][j] * solution[j]);  // 右边已经求出来的，左边自由变量遗留
        }
        solution[pivot[i]] = sum;  // 求解对应的主元变量
    }
    // assert(rksz <= m);
    if (rksz < m)
        return 2;  // 无穷多解
    return 1;      // 唯一解
}
// int t = gauss(b, n，m, sol);

1.Square - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

类题：只有数据范围不一样https://acm.hdu.edu.cn/showproblem.php?pid=5833

题意：给你n（<=300）个数，每个数的最大素因子均不会超过2000，每个数的范围<=1e18。每次选出一个非空子集，问你有多少种选法可以使得乘积为一个完全平方数

Sol：考虑关键性质是最大质因子不超过2000，完全平方数要求质因子的指数是偶数,也就是mod2意义下为0，问题转化为线性方程组；。我们打表发现2000以内质数只有303个，所以考虑高斯消元。b[i][j]b[i][j]表示第i个质因子在第j个数的指数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
int rksz;
int gauss(vector<vector<int>>& a, int n, int m) {
    vector<int> freevar;           // 记录自由变量对应的列
    vector<int> pivot(n + 1, -1);  // 记录每一行的主元所在的列
    int r = 1;
    for (int c = 1; c <= m; c++) {
        int t = r;
        // 找到当前列中的主元
        for (int i = r; i <= n; i++) {
            if (a[i][c]) {
                t = i;
                break;
            }
        }
        if (!a[t][c]) {
            freevar.push_back(c);
            continue;  // 当前列没有主元，继续到下一列
        }
        pivot[r] = c;  // 第 r 行的主元在 c 列
        if (t != r) {  // 交换行，将主元行放在第 r 行
            for (int i = c; i <= m + 1; i++)
                swap(a[r][i], a[t][i]);
        }
        // 消去主元下方的所有行
        for (int i = r + 1; i <= n; i++) {
            if (a[i][c])
                for (int j = m + 1; j >= c; j--) a[i][j] ^= a[r][j];
        }
        r++;
    }
    // 检查是否有解
    for (int i = r; i <= n; i++) {
        if (a[i][m + 1])
            return 0;  // 无解
    }
    rksz = r - 1;  // 这是系数矩阵的秩
    assert(rksz <= m);
    if (rksz < m)
        return 2;  // 无穷多解
    return 1;      // 唯一解
}
vector<int> primes, minp;
int cnt;
void sieve(int n) {
    minp.assign(n + 1, 0);
    primes.clear();

    for (int i = 2; i <= n; i++) {
        if (minp[i] == 0) {
            minp[i] = i;
            primes.push_back(i);
        }

        for (auto p : primes) {
            if (i * p > n) {
                break;
            }
            minp[i * p] = p;
            if (p == minp[i]) {
                break;
            }
        }
    }
    cnt = primes.size();
    deb(cnt);
}
int tt = 0;
void solve() {
    int n;
    tt++;
    cin >> n;
    vector<int> a(n + 1);
    for (int i = 1; i <= n; i++) cin >> a[i];
    vector<vector<int>> b(cnt + 1, vector<int>(n + 2));
    for (int i = 1; i <= n; i++) {
        int x = a[i];
        for (int j = 1; j <= cnt; j++) {
            int v = 0;
            while (x % primes[j - 1] == 0) {
                x /= primes[j - 1];
                v++;
            }
            b[j][i] = v % 2;
        }
    }
    int t = gauss(b, cnt, n);
    deb(rksz);
    if (t == 0)
        cout << 0 << endl;
    else {
        assert(rksz <= n);
        int ans = 1LL;
        for (int i = 1; i <= n - rksz; i++) ans = (ans * 2) % mod;
        ans = (ans + mod - 1) % mod;
        cout << "Case #" << tt << ":" << endl;
        cout << ans << endl;
    }
}

【日报配套题单】高斯消元，线性代数 - 题单 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

高斯消元模板及习题 - BigSmall_En - 博客园 (cnblogs.com)

高斯消元法(Gauss Elimination) 分析 & 题解 & 模板——czyuan原创 - Sizaif - 博客园 (cnblogs.com)