高斯 - 约当消元法
定义
高斯消元方法
给定 \(n\) 元一次方程组
\[\begin{cases}
a_{1,1}x_1+a_{1,2}x_2+\cdots+a_{1,n}x_n=b_1\\
a_{2,1}x_1+a_{2,2}x_2+\cdots+a_{2,n}x_n=b_2\\
\cdots\\
a_{n,1}x_1+a_{n,2}x_2+\cdots+a_{n,n}x_n=b_n\\
\end{cases}
\]
请求出方程组的解的情况：
-
无解；
-
无穷多解；
-
唯一解。
对于这样的问题，我们可以使用 高斯消元法 进行求解，当然高斯消元法有一个回代的过程，代码略长，而且精度较低。
所以我们隆重推出 高斯-约旦消元法 ！！！
回顾一下我们是怎么手算的，一般用的都是 加减消元法，普通高斯和高斯-约旦用的都是加减消元。
在此之前，我们需要了解一下矩阵初等变换。
在线性代数中，矩阵初等行变换 是指以下三种变换类型 ：
-
交换矩阵的两行；
-
用一个非零数 \(k\) 乘矩阵的某一行所有元素；
-
把矩阵的某一行所有元素乘以一个数 \(k\) 后加到另一行对应的同一列的元素上；
类似地，把以上的 行 改为 列 便得到 矩阵初等列变换 的定义。
矩阵初等行变换与初等列变换合称为 矩阵初等变换。
若矩阵 \(A\) 经过有限次的初等行变换变为矩阵 \(B\)，则矩阵 \(A\) 与矩阵 \(B\) 行等价；若矩阵 \(A\) 经过有限次的初等列变换变为矩阵 \(B\)，则矩阵 \(A\) 与矩阵 \(B\) 列等价；若矩阵 \(A\) 经过有限次的初等变换变为矩阵 \(B\)，则矩阵 \(A\) 与矩阵 \(B\) 等价。
当然列的用不着
首先有一个由系数构成的 \(n\times n\) 的矩阵
\[\begin{bmatrix}
a_{1,1}&a_{1,2}&\cdots&a_{1,n}\\
a_{2,1}&a_{2,2}&\cdots&a_{2,n}\\
\vdots&\vdots&\ddots&\vdots\\
a_{n,1}&a_{n,2}&\cdots&a_{n,n}\\
\end{bmatrix}
\]
然后是一个由常数构成的 \(n\times 1\) 的列向量
\[\begin{bmatrix}
b_1\\
b_2\\
\vdots\\
b_n
\end{bmatrix}
\]
把它们放在一起构成一个 \(n\times(n+1)\) 的增广矩阵：
\[\begin{bmatrix}
a_{1,1}&a_{1,2}&\cdots&a_{1,n}&\mid&b_1\\
a_{2,1}&a_{2,2}&\cdots&a_{2,n}&\mid&b_2\\
\vdots&\vdots&\ddots&\vdots&\mid&\vdots\\
a_{n,1}&a_{n,2}&\cdots&a_{n,n}&\mid&b_n\\
\end{bmatrix}
\]
我们遍历每一列，对于第 \(i\) 列选出最大的 未处理过的行 上的数作为主元，意味着加减消元后除了主元这一行外其他行的第 \(i\) 列都为 \(0\)。
选最大的作为主元的原因是避免精度误差。
然后就是加减消元了。
举个例子
\[\begin{cases}
2x-y+z=1\\
4x+y-z=5\\
x+y+z=0
\end{cases}
\]
增广矩阵就是
\[\begin{bmatrix}
2&-1&1&\mid&1\\
4&1&-1&\mid&5\\
1&1&1&\mid&0
\end{bmatrix}
\]
先是第 \(1\) 列，选 \(4\) 为主元。
\(4\) 在第 \(2\) 行，将第 \(2\) 行与第 \(1\) 行交换。
\[\begin{bmatrix}
4&1&-1&\mid&5\\
2&-1&1&\mid&1\\
1&1&1&\mid&0
\end{bmatrix}
\]
对于第 \(2\) 行，第一列上的 \(2\) 是 \(4\) 的 \(\dfrac{1}{2}\)。
\[\begin{bmatrix}
4&1&-1&\mid&5\\
2-4\times\dfrac{1}{2}&-1-1\times\dfrac{1}{2}&1-(-1)\times\dfrac{1}{2}&\mid&1-5\times\dfrac{1}{2}\\
1&1&1&\mid&0
\end{bmatrix}
\]
化简得
\[\begin{bmatrix}
4&1&-1&\mid&5\\
0&-\dfrac{3}{2}&\dfrac{3}{2}&\mid&-\dfrac{3}{2}\\
1&1&1&\mid&0
\end{bmatrix}
\]
对第 \(3\) 行的处理同理
\[\begin{bmatrix}
4&1&-1&\mid&5\\
0&-\dfrac{3}{2}&\dfrac{3}{2}&\mid&-\dfrac{3}{2}\\
0&\dfrac{3}{4}&\dfrac{5}{4}&\mid&-\dfrac{5}{4}
\end{bmatrix}
\]
现在到了第 \(2\) 列，未处理过的是 \(2,3\) 行，选最大的 \(\dfrac{3}{4}\) 为主元。
将第 \(3\) 行与第 \(2\) 行交换
\[\begin{bmatrix}
4&1&-1&\mid&5\\
0&\dfrac{3}{4}&\dfrac{5}{4}&\mid&-\dfrac{5}{4}\\
0&-\dfrac{3}{2}&\dfrac{3}{2}&\mid&-\dfrac{3}{2}
\end{bmatrix}
\]
消元得
\[\begin{bmatrix}
4&0&-\dfrac{8}{3}&\mid&\dfrac{20}{3}\\
0&\dfrac{3}{4}&\dfrac{5}{4}&\mid&-\dfrac{5}{4}\\
0&0&4&\mid&-4
\end{bmatrix}
\]
第 \(3\) 列，未处理过的是第 \(3\) 行，选 \(4\) 作主元。
\[\begin{bmatrix}
4&0&0&\mid&4\\
0&\dfrac{3}{4}&0&\mid&0\\
0&0&4&\mid&-4
\end{bmatrix}
\]
这样就把原来的矩阵通过初等变换，使得系数矩阵只有主对角线位置的元素非零，即一个对角矩阵。
上面那个矩阵的意思是
\[\begin{cases}
4x=4\\
\dfrac{3}{4}y=0\\
4z=-4
\end{cases}
\]
所以再把系数除过去就得到
\[\begin{cases}
x=1\\
y=0\\
z=-1
\end{cases}
\]
请自行检验。
时间复杂度为 \(\mathcal{O}(n^3)\)。
当然方程还可能是无解或无穷多解。
考虑一个一元一次方程 \(ax=b\) 的解的情况：
- 无解：\(a=0,b\ne0\)；
- 无穷多解：\(a=b=0\)；
- 唯一解：\(a\ne0\)。
- 所以当发现某一列的主元为 \(0\) 时，因为主元是最大的，所以意味着这一列全都是 \(0\)，那么要么无解，要么有无穷多解。
#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 10, mod = 1e9 + 7;
double eps = 1e-7;
signed main()
{
int n; cin >> n;
vector<vector<double>> a(n + 5, vector<double> (n + 6));
for(int i = 1; i <= n; i++){
for(int j = 1; j <= n + 1; j++){
cin >> a[i][j];
}
}
for(int j = 1; j <= n; j++){ // 枚举列
int max = j;
for(int k = j + 1; k <= n; k++){
if(fabs(a[k][j]) > fabs(a[max][j])){
max = k; // 得到 j 列的最大系数所在的行
}
}
// 最大元系数调换到第 j 行
swap(a[j], a[max]);
if(fabs(a[j][j]) < eps){
return cout << "No Solution" << '\n', 0;
}
// 主元系数变为 1
for(int k = n + 1; k >= 1; k--){
a[j][k] = a[j][k] / a[j][j];
}
// 其余行系数变成 0
for(int i = 1; i <= n; i++){
if(i == j) continue;
double tmp = a[i][j];
for(int k = 1; k <= n + 1; k++){
a[i][k] = a[i][k] - tmp * a[j][k];
}
}
}
for(int i = 1; i <= n; i++) printf("%.2lf\n", a[i][n + 1]);
return 0;
}
关于判断无解和无穷多解
关于判断无解和无穷多解
这下要具体到到底是无解还是无穷多解了。
这就是高斯-约旦的一个缺点：相比于普通高斯，它更难判断无解和无穷多解。
但是也是可以处理的。
具体地，我们用 \(r\) 来记录当前行，如果主元为 \(0\)，那么 \(\operatorname{continue}\) 到下一列，但 \(r\) 不变；否则消元后令 \(r\gets r+1\)。
遍历完所有列后：
- 若 \(r=n\)，说明有唯一解；
- 若 \(r<n\)，说明第 \(r+1\sim n\) 行系数矩阵全都是 \(0\)，若列向量全是 \(0\)，说明有无穷多解，否则无解。
#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 10, mod = 1e9 + 7;
double eps = 1e-7;
signed main()
{
int n; cin >> n;
vector<vector<double>> a(n + 5, vector<double> (n + 6));
for(int i = 1; i <= n; i++){
for(int j = 1; j <= n + 1; j++){
cin >> a[i][j];
}
}
int r = 1;
for(int j = 1; j <= n; j++){ // 枚举列
int max = r;
for(int k = r + 1; k <= n; k++){
if(fabs(a[k][j]) > fabs(a[max][j])){
max = k; // 得到 j 列的最大系数所在的行
}
}
// 最大元系数调换到第 r 行
swap(a[r], a[max]);
if(fabs(a[r][j]) < eps) continue;
// 主元系数变为 1
for(int k = n + 1; k >= 1; k--){
a[r][k] = a[r][k] / a[r][j];
}
// 其余行系数变成 0
for(int i = 1; i <= n; i++){
if(i == r) continue;
double tmp = a[i][j];
for(int k = 1; k <= n + 1; k++){
a[i][k] = a[i][k] - tmp * a[r][k];
}
}
r ++;
}
r --;
bool vis1 = false, vis2 = false;
for(int i = 1; i <= n; i++){
bool ok = false;
for(int j = 1; j <= n; j++){
if(fabs(a[i][j]) > eps) ok = true;
}
if(!ok && fabs(a[i][n + 1]) < eps) vis1 = true;
if(!ok && fabs(a[i][n + 1]) > eps) vis2 = true;
}
if(vis2) return cout << "-1" << '\n', 0;
if(vis1) return cout << "0" << '\n', 0;
for(int i = 1; i <= n; i++){
cout << "x" << i;
printf("=%.2lf\n", a[i][n + 1]);
}
return 0;
}
关于多项异或线性方程式求解
XOR Neighbors
观察每个点的性质, 发现对于任何一个点 \(u\) 来说, 设其邻点为 \(v\), 则有 \(\bigoplus_{v_i \in \text{Adj}(u)} x_{v_i} = 0\), 那么这样其实就符合一个异或方程组, 对于 \(n\) 个点就有 \(n\) 个方程组, 那么直接使用高斯消元求, 但是异或方程式与普通的方程式并不一样, 这里介绍一下异或方程主要特征:
- 在消除当前行的其余元时, 应逆序在这一行上异或上主元乘当前元, 这样可以使当前行的元变为 \(0\)
- 异或方程式无法像一般线性方程式那样可以化简成行最简形式, 经过高斯消元之后会变成上三角的样式
- 在求解时, 应从最后一行一步一步向上求解, 由于不是行最简, 所以要把所有的非零元代入方程求解
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e6 + 10, mod = 1e9 + 7;
signed main()
{
std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
int n, m; cin >> n >> m;
vector<vector<int>> a(n + 10, vector<int> (n + 10));
for(int i = 1; i <= m; i++){
int u, v; cin >> u >> v;
a[u][v] = a[v][u] = 1;
}
int r = 1;
for(int j = 1; j <= n; j++){
int mx = r;
for(int i = r + 1; i <= n; i++){
if(a[i][j] > a[mx][j]){
mx = i;
}
}
swap(a[r], a[mx]);
for(int i = r + 1; i <= n; i++){
if(i == r || a[i][j] == 0) continue;
for(int k = n + 1; k >= j; k--){
a[i][k] = a[i][k] ^ (a[i][j] * a[r][k]);
}
}
r++;
}
int tot = 0;
for(int i = n; i >= 1; i--){
for(int j = i + 1; j <= n; j++){
a[i][n + 1] = a[i][n + 1] ^ (a[i][j] * a[j][n + 1]);
}
if((a[i][i] == 0 && a[i][n + 1] != 0) || (a[i][i] != 0 && a[i][n + 1] == 0)){
cout << "No" << '\n';
return 0;
}
if(a[i][i] == 0 && a[i][n + 1] == 0) a[i][n + 1] = (1 << tot++);
}
cout << "Yes" << '\n';
for(int i = 1; i <= n; i++) cout << a[i][n + 1] << ' ';
return 0;
}
例题
题目简述: 有一个由 n*n 个小方砖组成的正方形墙。一些砖块是白色的，而一些是黄色的。
鲍勃是一名画家，他想把所有的砖块都涂成黄色。但鲍勃的刷子有点问题，一旦他用这个刷子涂上砖块 (i, j)，
那么 (i, j)、(i-1, j)、(i+1, j)、(i, j-1) 和 (i, j+1) 的砖块颜色都会发生变化。
你的任务是找到鲍勃应该涂的最少砖块数量，以使所有的砖块都变成黄色。
本题乍一看很难想到使用高斯消元，下面通过建模转化成线性方程组问题。
共有 \(k = n * n\) 个方砖，那么设 \(x_i\) 为当前方砖刷的次数，明显的，每个方砖刷的次数最多为 \(1\)。
定义 \(a_{i,j}\) 表示两块砖的关系，若刷了第 \(j\) 块砖头对 \(i\) 有影响那么就是 \(1\)，特别的 \(a_{ii} = 1\)。
对于方砖 \(i\) 的颜色，设黄色为 \(0\)，白色为 \(1\)，那么最后颜色取决于初始颜色 \(s_i\) 和所有其\(a_{ij} = 1\) 的异或操作情况
故有如下方程组。
\[\begin{cases}
s_1\oplus{a_{1,1}}\&\&{x_1}\oplus{a_{1,2}}\&\&{x_2}\oplus\cdots\oplus{a_{1,k}}\&\&{x_k}=0\\
s_2\oplus{a_{2,1}}\&\&{x_1}\oplus{a_{2,2}}\&\&{x_2}\oplus\cdots\oplus{a_{2,k}}\&\&{x_k}=0\\
\cdots\\
s_k\oplus{a_{k,1}}\&\&{x_1}\oplus{a_{k,2}}\&\&{x_2}\oplus\cdots\oplus{a_{k,k}}\&\&{x_k}=0\\
\end{cases}
\]
那么 \(x_i = 1\) 就是被刷了一次，那么最后的答案就是 \(\sum\limits_{i=1}^{k}{x_i}\)
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string.h>
#include<algorithm>
#include<map>
#include<queue>
#include<vector>
#include<cmath>
#include<stdlib.h>
#include<time.h>
using namespace std;
const int N = 300, mod = 1e9 + 7;
int a[N][N], dx[] = {0, 0, 1, -1}, dy[] = {1, -1, 0, 0};
void solve(){
int n; cin >> n;
memset(a, 0, sizeof a);
for(int i = 1; i <= n * n; i ++){
char x; cin >> x;
if(x == 'y') a[i][n * n + 1] = 0;
else if(x == 'w') a[i][n * n + 1] = 1;
}
for(int i = 1; i <= n; i ++){
for(int j = 1; j <= n; j++){
a[n * (i - 1) + j][n * (i - 1) + j] = 1;
for(int k = 0; k < 4; k ++){
int xx = i + dx[k], yy = j + dy[k];
if(xx < 1 || xx > n || yy < 1 || yy > n) continue;
a[n * (i - 1) + j][n * (xx - 1) + yy] = 1;
}
}
}
int r = 1;
for(int j = 1; j <= n * n; j ++){
int max = r;
for(int k = r + 1; k <= n * n; k ++){
if(abs(a[k][j]) > abs(a[max][j])) max = k;
}
if(a[max][j] == 0) continue;
for(int k = 1; k <= n * n + 1; k ++) swap(a[r][k], a[max][k]);
for(int i = 1; i <= n * n; i ++){
if(i == r || a[i][j] == 0) continue;
for(int k = 1; k <= n * n + 1; k ++) a[i][k] ^= a[r][k];
}
r ++;
}
for(int i = r; i <= n * n; i++){
if(a[i][n * n + 1]) return cout << "inf" << '\n', void();
}
int res = 0;
vector<int> x(N);
for(int i = 1; i <= n * n; i ++){
x[i] = a[i][n * n + 1];
for(int j = i + 1; j <= n * n; j++) x[i] ^= (a[i][j] && x[j]);
res += x[i];
}
cout << res << '\n';
}
signed main()
{
std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
int t; cin >> t; while(t--) solve();
return 0;
}
开关问题
与上述的异或方程式同理, 不再赘述, 每个灯都能影响其他的灯, 所以可以列出对所有对该灯有影响的异或方程式: \(\bigoplus_ {A_{i,j}*x_i} = end_i\), 对于一个自由变量, 那么我们就有开或不开两种, 乘\(2\)即可
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e6 + 10, mod = 1e9 + 7;
void solve(){
int n; cin >> n;
vector<vector<int>> a(n + 10, vector<int> (n + 10));
for(int i = 1; i <= n; i++) cin >> a[i][n + 1];
for(int i = 1; i <= n; i++){
int x; cin >> x;
a[i][n + 1] ^= x;
a[i][i] = 1;
}
int x, y;
while(cin >> x >> y && (x || y)){
a[y][x] = 1;
}
int r = 1;
for(int j = 1; j <= n; j++){
int mx = r;
for(int i = r + 1; i <= n; i++){
if(a[i][j] > a[mx][j]){
mx = i;
}
}
swap(a[r], a[mx]);
for(int i = r + 1; i <= n; i++){
if(i == r || a[i][j] == 0) continue;
for(int k = n + 1; k >= j; k--){
a[i][k] = a[i][k] ^ (a[i][j] * a[r][k]);
}
}
r++;
}
int res = 1;
for(int i = n; i >= 1; i--){
for(int j = i + 1; j <= n; j++){
a[i][n + 1] = a[i][n + 1] ^ (a[i][j] * a[j][n + 1]);
}
if(a[i][i] == 0 && a[i][n + 1] != 0){
return cout << "Oh,it's impossible~!!" << '\n', void();
}
if(a[i][i] == 0 && a[i][n + 1] == 0) res *= 2;
}
cout << res << '\n';
}
signed main(){
std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int t;cin>>t;while(t--)solve();
}
