[学习笔记] 高斯消元 - 线性代数

合集 - 学习笔记(14)

1.[学习笔记] MST(最小生成树) - 图论2024-01-132.[学习笔记] LCA - 图论2024-04-113.[学习笔记] 树上差分 - 图论2024-04-13

4.[学习笔记] 高斯消元 - 线性代数2024-04-17

5.[学习笔记] 丢番图方程 & 同余 & 逆元 - 数论2024-04-216.[学习笔记] 质数与唯一分解定理 - 数论2024-05-027.[学习笔记] 乘性函数 - 数论2024-05-068.[学习笔记] hash & kmp & Trie树 - 字符串2024-05-109.[学习笔记] 单调队列优化DP - DP2024-05-3110.[学习笔记] 斜率优化DP - DP2024-06-1611.[学习笔记] 树链剖分(重链剖分) - 图论 & 数据结构2024-06-1912.[学习笔记] 动态开点权值线段树合并 - 数据结构2024-07-0313.[学习笔记] 长链剖分 - 图论2024-07-1214.[学习笔记] Splay & Treap 平衡树 - 数据结构2024-08-29

收起

高斯-约旦消元

下面给两道板子

【模板】高斯消元法

最最基础的板子，没啥哆嗦的。下面给出高斯-约旦消元解法。

#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;
}

[SDOI2006] 线性方程组

这是一道很好的题(板子)，跟上题唯一的区别就是判无解和无穷解。这题打通了说明你对高斯消元的理解也就差不多了。

搞了两天，终于把高斯-约旦消元的解调出来了。

高斯-约旦消元 和 高斯消元 的原理是一样的，只不过实现方式略有差别。高斯消元是先把主元以前的系数统统化为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;
}

[JSOI2008] 球形空间产生器

这道题没有啥太大的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;
}

高斯消元与异或方程组

第一个要说的就是异或方程组了。

$$\left\{\begin{matrix} a_{11}*x_1\oplus a_{12}*x_2\oplus a_{13}*x_3…=k_1\\ a_{21}*x_1\oplus a_{22}*x_2\oplus a_{23}*x_3…=k_2\\ a_{31}*x_1\oplus a_{32}*x_2\oplus a_{33}*x_3…=k_3\\ \end{matrix}\right.$$

其中：系数 $a_{ij} = 1\ or\ 0$ ，常数项 $k_i=1\ or\ 0$

形如以上狮子的方程组就命名为异或方程组。他跟普通方程组没啥太大区别，甚至还要更简单。

求解异或方程组需要了解异或的性质：

$$\begin{split} 1\oplus 1=0\\ 0\oplus 0=0\\ 1\oplus 0=1\\ 0\oplus 1=1 \end{split}$$

简单来说就是 按位运算，同0异1。扩展性质：

$$\begin{split} A\oplus A=0\\ \end{split}$$

于是，对于异或方程组，消元后的理想状态是矩阵的对角线上都为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] 开关问题

读完题，可以列出一个矩阵表示灯和灯之间的关系：

$$\begin{bmatrix} 1 & 0 &1 & 1\\ 0 &1 &1 &1 \\ 1 &0 &1 &0 \\ 0 &0 & 0 &1 \end{bmatrix}$$

其中，第 $i$ 行 第 $j$ 列 表示操作 $j$ 号灯能否影响 $i$ 号灯。特别地，每个灯都能被自己影响。$x$ 表示操作次数。于是有：

$$\begin{bmatrix} 1 & 0 &1 & 1\\ 0 &1 &1 &1 \\ 1 &0 &1 &0 \\ 0 &0 & 0 &1 \end{bmatrix}* \begin{bmatrix} x_1\\ x_2\\ x_3\\ x_4 \end{bmatrix}= \begin{bmatrix} x_1+0+x_1+x_1&=K_1\\ 0+x_2+x_2+x_2&=K_2\\ x_3+0+x_3+0&=K_3\\ 0+0+0+x_4&=K_4 \end{bmatrix}$$

只要让等式两边的数值的个位相等即可。可以巧用异或求解。狮子形式如最上面的那个。矩阵的轶就是化成阶梯矩阵后系数全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';
	}
}

高斯消元与同余方程组