[学习笔记] 高斯消元 - 线性代数
高斯-约旦消元
下面给两道板子
最最基础的板子，没啥哆嗦的。下面给出高斯-约旦消元解法。
#include<bits/stdc++.h>
using namespace std;
int n, dt = 1;
double eps = 1e-9, m[102][102];
int main(){
scanf("%d", &n);
for(int i=1; i<=n; ++i)
for(int j=1; j<=n+1; ++j)
scanf("%lf", &m[i][j]);
for(int i=1; i<=n; ++i){
int mx = i;
for(int j=i+1; j<=n; ++j)
if(fabs(m[mx][i]) < fabs(m[j][i])) mx = j;
if(fabs(m[mx][i]) < eps) dt = 0;
swap(m[mx], m[i]);
for(int j=n+1; j>=i; --j) m[i][j] /= m[i][i];
for(int j=1; j<=n; ++j){
if(j == i) continue;
double t = m[j][i];
for(int k=1; k<=n+1; ++k) m[j][k] -= t*m[i][k];
}
}
if(!dt){
printf("No Solution");
return 0;
}
for(int i=1; i<=n; ++i) printf("%.2f\n", m[i][n+1]);
return 0;
}
这是一道很好的题(板子)，跟上题唯一的区别就是判无解和无穷解。这题打通了说明你对高斯消元的理解也就差不多了。
搞了两天，终于把高斯-约旦消元的解调出来了。
高斯-约旦消元 和 高斯消元 的原理是一样的，只不过实现方式略有差别。高斯消元是先把主元以前的系数统统化为0，再回代求解；而高斯-约旦消元是把主元化为一，在把所有狮子的主元都化为0，让主元系数都为1并且在对角线上排列。
对于形如这样的狮子\([0\ 0\ 0\ …\ |\ k]\)， 无解的判定方法是所有系数都为0但 \(k\) 不为0，无穷解的判定方式是所有系数都为0且 \(k\) 等于0 或存在某个主元都为0。易知：判定无解在先，其次为无穷解，最后输出答。
所以，在进行高-约消元的时候，如果出现个主元的最大值为0，那么会出现两种情况，要么无解，要么无穷解，所以要查一遍这个狮子的所有系数，看是否满足上述情况。如果判不出无解，这一列主元也就没用了，为了不影响后续的消元操作，那就把这一列全部移动到右边去。选出来的这一行也移动到下边去。所以最后消完元的子矩阵应该位于完整矩阵的左上角。最后再遍历一遍整个矩阵看是否无解即可。
注意细节。
#include<bits/stdc++.h>
using namespace std;
int n, dt = 1, a, b; // a b 分别代表当前的行和列
double eps = 1e-9, m[52][52]; //eps为精度，因为涉及到小数运算，所以不太会出现完全等于0的数
int main(){
scanf("%d", &n);
a = n, b = n+1;
for(int i=1; i<=n; ++i)
for(int j=1; j<=n+1; ++j)
scanf("%lf", &m[i][j]);
for(int i=1; i<=a; ++i){
int mx = i;
for(int j=i+1; j<=b; ++j)
if(fabs(m[mx][i]) < fabs(m[j][i])) mx = j;
if(fabs(m[mx][i]) < eps){ //如果主元都为0
bool cnt = 1;
for(int j=i+1; j<=n; ++j) //查是否无解
if(fabs(m[mx][j]) >= eps) cnt = 0;
if(cnt && fabs(m[mx][n+1]) >= eps){
printf("-1");
return 0;
}
dt = 0;
swap(m[mx], m[a--]); //移动
--b;
for(int j=1; j<=n; ++j) swap(m[j][i], m[j][b]); //移动
--i;
continue; //一是防止报错 二是再消也没啥大用 毕竟主元都没了
}
swap(m[mx], m[i]);
for(int j=n+1; j>=i; --j){
if(j >= b && j <= n) continue; //避免误判
m[i][j] /= m[i][i];
}
for(int j=1; j<=n; ++j){ //这里很重要，一定要把能消的都消干净
if(j == i) continue;
double t = m[j][i];
for(int k=1; k<=n+1; ++k) m[j][k] -= t*m[i][k];
}
}
// for(int i=1; i<=n; ++i){
// for(int j=1; j<=n+1; ++j){
// printf("%.2f ", m[i][j]);
// }printf("\n");
// }
for(int i=1; i<=n; ++i){ //再扫一遍 看是否无解
bool cnt = 1;
for(int j=1; j<=n; ++j)
if(fabs(m[i][j]) >= eps) cnt = 0;
if(cnt && fabs(m[i][n+1]) >= eps){
printf("-1");
return 0;
}
}
if(!dt){ //无穷解
printf("0");
return 0;
}
for(int i=1; i<=n; ++i) printf("x%d=%.2f\n", i, m[i][n+1]);
return 0;
}
这道题没有啥太大的nan度，唯推一推狮子尔。通过题干给的公式，我们很容易能推理出每个点到球心的距离都相等这一关系，再围绕这一关系建立方程组即可。推狮子的部分就不展示了。
#include<bits/stdc++.h>
using namespace std;
int n;
double eps = 1e-9, m[15][15];
int main(){
scanf("%d", &n);
for(int i=1; i<=n+1; ++i){
for(int j=1; j<=n; ++j){
double w;
scanf("%lf", &w);
m[i-1][j] -= 2*w, m[i-1][n+1] -= w*w;
m[i][j] += 2*w, m[i][n+1] += w*w;
}
}
for(int i=1; i<=n; ++i){
int mx = i;
for(int j=i+1; j<=n; ++j)
if(fabs(m[mx][i]) < fabs(m[j][i])) mx = j;
swap(m[mx], m[i]);
for(int j=n+1; j>=i; --j) m[i][j] /= m[i][i];
for(int j=1; j<=n; ++j){
if(j == i) continue;
double t = m[j][i];
for(int k=1; k<=n+1; ++k) m[j][k] -= t*m[i][k];
}
}
for(int i=1; i<=n; ++i) printf("%.3f ", m[i][n+1]);
return 0;
}
高斯消元与异或方程组
第一个要说的就是异或方程组了。
其中：系数 \(a_{ij} = 1\ or\ 0\) ，常数项 \(k_i=1\ or\ 0\)
形如以上狮子的方程组就命名为异或方程组。他跟普通方程组没啥太大区别，甚至还要更简单。
求解异或方程组需要了解异或的性质：
简单来说就是 按位运算，同0异1。扩展性质：
于是，对于异或方程组，消元后的理想状态是矩阵的对角线上都为1，其他都为0。和普通高斯消元大同小异，事先找到一个主元系数不为零的行移动到对角线，再用这一行消掉剩下该主元系数不为零的行。如果出现形为\([0\ 0\ 0\ …\ |\ 0]\)的狮子，那就说明有无穷解。如果出现形为\([0\ 0\ 0\ …\ |\ 1]\)的狮子，那就说明无解。进一步探究，还会发现一个现象就是方程的解总为1 or 0。又由于该方程所有的常数都为0 or 1，所以很方便使用bitset或int类型进行压位操作。
下面是 OI WIKI 上给出的模板代码：
std::bitset<1010> matrix[2010]; // matrix[1~n]：增广矩阵，0 位置为常数
std::vector<bool> GaussElimination(
int n, int m) // n 为未知数个数，m 为方程个数，返回方程组的解
// （多解 / 无解返回一个空的 vector）
{
for (int i = 1; i <= n; i++) {
int cur = i;
while (cur <= m && !matrix[cur].test(i)) cur++;
if (cur > m) return std::vector<bool>(0);
if (cur != i) swap(matrix[cur], matrix[i]);
for (int j = 1; j <= m; j++)
if (i != j && matrix[j].test(i)) matrix[j] ^= matrix[i];
}
std::vector<bool> ans(n + 1, 0);
for (int i = 1; i <= n; i++) ans[i] = matrix[i].test(0);
return ans;
}
给个例题 - [poj1830] 开关问题
读完题，可以列出一个矩阵表示灯和灯之间的关系：
其中，第 \(i\) 行 第 \(j\) 列 表示操作 \(j\) 号灯能否影响 \(i\) 号灯。特别地，每个灯都能被自己影响。\(x\) 表示操作次数。于是有：
只要让等式两边的数值的个位相等即可。可以巧用异或求解。狮子形式如最上面的那个。矩阵的轶就是化成阶梯矩阵后系数全0行的行数(就是有多少不确定的解)。而解只有0或1两种类型，所以最后的答案就是2的秩次方。
#include<bits/stdc++.h>
using namespace std;
int k, n, ans;
bitset<32> m[32]; //bitset压位
int main(){
ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
cin>>k;
while(k--){
ans = 0;
for(int i=1; i<=n; ++i) m[i].reset();
cin>>n;
for(int j=1; j<=2; ++j) for(int i=1, a; i<=n; ++i) cin>>a, m[i][n+1] = m[i][n+1] ^ a, m[i][i] = 1;
int a, b, line = 1;
bool data = 0;
while(cin>>a>>b && a) m[b][a] = 1; //反向输入
for(int i=1; i<=n; ++i){
int mx = line;
while(mx <= n && !m[mx][i]) ++mx;
if(mx > n) continue;
++line;
swap(m[mx], m[i]);
for(int j=1; j<=n; ++j)
if(j != i && m[j][i]) m[j] ^= m[i];
}
for(int i=1; i<=n; ++i){
bool cnt = 1;
for(int j=1; j<=n; ++j)
if(m[i][j]) cnt = 0;
if(cnt && m[i][n+1]) data = 1;
if(cnt && !m[i][n+1]) ++ans; //统计秩
}
if(data) cout<<"Oh,it's impossible~!!\n";
else cout<<(1<<ans)<<'\n';
}
}
