高斯消元
title: 高斯消元
categories:
- ICPC
tags:
- null
abbrlink: 55f2dbb7
date: 2023-10-29 00:00:00
高斯消元
1.高斯消元解线性方程组。这个模板只处理n方程n个未知数的，分别处理无解，无穷多解，以及唯一解的情况。结果存在最后一列中（一开始的常数项），
- 无解返回0
- 唯一解返回1
- 无穷多解返回2
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
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=B$问题的模板：处理了无穷多解
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列的高斯消元异或
有重大越界问题，现在已经更改。就是方程数比未知数多的时候，我们需要对于提前到达最后一行的情况判断，将后面的变量全部放入自由变量
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(n^2m/w)$https://qoj.ac/contest/1716/problem/313
5002是最多5000个自变量，并且多留一位给常数项
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]$表示第i个质因子在第j个数的指数。
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)
本博客所有文章除特别声明外，均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来源 爱飞鱼的blog！