PUSH NOTE : 高斯消元.md

docs/数学/线性代数/高斯消元.md

@@ -0,0 +1,260 @@
+---
+tags:
+  - 高斯消元
+  - 数学
+date: 2023-08-26
+publish: true
+---
+
+
+**别随便做初等列变换！要不然解出来解的顺序会变！**
+
+**系数化为1的时候，要把主元系数单独存起来，否则就覆盖掉了**
+
+## 高斯-约旦消元
+
+用 $i \in [1,n]$ 标定当前求解第几个未知数（第几列），$line$ 存储当前处理到第几行，每次进行：
+
+1. 初等行变换，找到第$i$个数（主元系数）绝对值最大的一行，若不为0，使用swap提上来；若为0，则$line$不变，$i$加一，跳过该元。
+2. 对角化，枚举每一行消去当前元。
+3. 把所有系数不为0的消干净后，用剩下几行（系数此时应该全为0）来判断无解/无数解。
+4. 若剩下几行常数项均为0，则说明其为自由元，可取任意值（最好取0或模数）；否则为无解。
+
+```cpp
+#include <bits/stdc++.h>
+using std::cin;
+const int maxn = 55;
+int n;
+double a[maxn][maxn];
+const double eps = 1e-200;
+
+int Gauss() {
+    int line = 1;                  // 枚举行
+    for (int i = 1; i <= n; i++) { // 枚举列
+        int maxl = line;           // 找绝对值最大主元系数所在行
+        // 不是绝对值的话 0卡在这个中间！！！
+        for (int k = line + 1; k <= n; k++) {
+            if (fabs(a[k][i]) > fabs(a[maxl][i])) maxl = k;
+        }
+        if (fabs(a[maxl][i]) < eps) continue; // 该列主元系数都为0 先跳过该元（列）
+        if (maxl != line) std::swap(a[maxl], a[line]);
+        // swap可以挪一整行 把系数不为0的提上来
+
+        // 对角化 枚举每个式子消去当前元
+        for (int k = 1; k <= n; k++) {
+            if (k == line) continue; // 不能把自己消了
+            double t = a[k][i] / a[line][i]; // 倍数
+            for (int j = 1; j <= n + 1; j++) {
+                a[k][j] -= t * a[line][j];
+            }
+        }
+        line++;
+    }
+    // 如果一切正常，line应该为n+1
+    line--;
+    if (line < n) {
+        for (int i = line + 1; i <= n; i++)
+            if (fabs(a[i][n + 1]) > eps) return -1;
+        return 0;
+    }
+    for (int i = 1; i <= n; i++) {
+        a[i][n + 1] /= a[i][i];
+        if (fabs(a[i][n + 1]) < eps) a[i][n + 1] = 0;
+    }
+    return 1;
+}
+
+int main() {
+    cin >> n;
+    for (int i = 1; i <= n; i++)
+        for (int j = 1; j <= n + 1; j++) cin >> a[i][j];
+    int flag = Gauss();
+    if (flag == 1) {
+        for (int i = 1; i <= n; i++) printf("x%d=%.2lf\n", i, a[i][n + 1]);
+    } else printf("%d", flag);
+}
+```
+
+**注意，如果无数解要得到一组可行解（且按顺序输出解）， 别忘了再带回去消元**
+
+巧妙的做法：如果不取模，令其等于0；如果是取模的方程，直接令其等于模数或0
+
+如果增广矩阵不是$n$行，也可以用：
+
+```cpp
+// cnt::有多少个方程
+// m  ::未知数个数
+int Gauss(){
+	// 只负责消元
+	int line = 1;
+	for (int i = 1; i <= m; i++) { // 当前在求第几个未知数
+		int maxl = line;
+		for (int j = line + 1; j <= cnt; j++)
+			if (a[j][i]) maxl = j;
+		if (a[maxl][i]==0) continue; // 该列主元系数都为0 先跳过该元（列）
+		std::swap(a[line], a[maxl]);
+		// 对角化
+		int tt = inv(a[line][i]); // 另存起来
+		for (int j = i; j <= m + 1; j++) a[line][j] = a[line][j]*tt % mod; // 别忘了取模
+		// 系数化为1，j从i或者1都行
+		for (int j = 1; j <= cnt; j++) {
+			if (j!=line && a[j][i]!=0) { // 别把自己消了 跳过0 小小**优化**一下
+				int t = a[j][i]; // 倍数
+				for (int k = i; k <= m + 1; k++)
+					a[j][k] = ((a[j][k] - t * a[line][k]) % mod + mod) % mod;
+				// 消当前i元，k从i或者1都行
+			}
+		}
+		line++;
+	}
+	return line-1; //返回有解的个数
+}
+
+//////////////////////////////////////////////////////////////////////////////
+
+// 如果不取模，令其等于0
+// 如果是取模的方程，直接令其等于模数或0
+int ans[maxm];
+void GetAns(int line){
+	for (int i = line + 1; i <= cnt; i++) {
+		if (a[i][m+1]) {
+			cout << -1 << "\n"; // 无解
+			return;
+		}
+	}
+	// 输出
+	// 有line个数的解是确定的
+	for (int i = 1; i <= line; i++) {
+		int p = i;
+		while (!a[i][p]) p++; // 找到系数不为0的一列
+		ans[p] = a[i][m+1];
+	}
+	for (int i = 1; i <= m; i++) cout << (ans[i] ? ans[i] : mod) << " ";
+	cout << "\n";
+}
+
+//////////////////////////////////////////////////////////////////////////////
+
+// 这玩意儿一点用都没有。。。除非你想指定解
+void GetAns2(int line) {
+	for (int i = line + 1; i <= cnt; i++) {
+		if (a[i][m+1]) {
+			cout << -1 << "\n"; // 无解
+			return;
+		}
+	}
+	// 带入特殊解
+	int p=1;
+	for (int i = 1; i <= line; i++,p++) {
+		while (a[i][p]==0) { // 该列主元系数都为0 追加
+			cnt++;
+			a[cnt][p]=1,a[cnt][m+1]=mod; // here...随意取
+			p++;
+		}
+	}
+	Gauss();
+	for (int i = 1; i <= m; i++) {
+		cout <<(a[i][m+1]?a[i][m+1]:mod)<< ' ';
+	}
+	cout << '\n';
+}
+```
+
+## 求行列式
+
+![[行列式#定义|行列式 > 定义]]
+
+方法是$O(n!)$，显然不可，我们可以化为$O(n^3)$。
+
+高斯约旦消元时，**在系数化为1之前记录下来系数**。行列式为对角线元素之积，符号由交换行的数量来确定（如果为奇数，则行列式的符号应颠倒）。
+
+**如果在某个时候，我们在当前列中找不到非零单元，则算法应停止并返回 0。**
+
+```cpp
+const double EPS = 1E-9;
+int n;
+double det = 1,a[maxn][maxn];
+for (int i = 0; i < n; ++i) {
+  int k = i;
+  for (int j = i + 1; j < n; ++j)
+    if (abs(a[j][i]) > abs(a[k][i])) k = j;
+  if (abs(a[k][i]) < EPS) {
+    det = 0;
+    break;
+  }
+  swap(a[i], a[k]);
+  if (i != k) det = -det;
+  det *= a[i][i];
+  for (int j = i + 1; j < n; ++j) a[i][j] /= a[i][i];
+  for (int j = 0; j < n; ++j)
+    if (j != i && abs(a[j][i]) > EPS)
+      for (int k = i + 1; k < n; ++k) a[j][k] -= a[i][k] * a[j][i];
+}
+
+cout << det;
+```
+
+```cpp
+LL a[K][K];
+void gauss(int n) {
+	LL det = 1;
+	for (int i = 0; i < n; ++i) {
+		int k = i;
+		for (int j = i + 1; j < n; ++j)
+			if (abs(a[j][i]) > abs(a[k][i])) k = j;
+		// 提上来最大的一行
+		if (a[k][i] == 0) {
+			det = 0;
+			break;
+		}
+		std::swap(a[i], a[k]);
+		// 变号
+		if (i != k) det = -det;
+		det = det*a[i][i]%p;
+		// 系数化为1
+		for (int j = i + 1; j < n; ++j) a[i][j] = a[i][j]*inv(a[i][i])%p;
+		for (int j = 0; j < n; ++j)
+			if (j != i && a[j][i]!=0)
+				for (int k = i + 1; k < n; ++k) a[j][k] = ((a[j][k] - a[i][k] * a[j][i]%p)%p+p)%p;
+	}
+	cout << (det%p+p)%p << '\n';
+}
+```
+
+
+## 矩阵求逆
+
+两个相乘为单位矩阵的矩阵，互为逆矩阵。给出 $n$ 阶方阵 $A$，求解其逆矩阵的方法如下：
+
+1. 构造 $n \times 2n$ 的矩阵 $(A, I_n)$，拼在右边；
+2. 用高斯消元法将其化简为最简形 $(I_n, A^{-1})$，即可得到 $A$ 的逆矩阵 $A^{-1}$。如果最终最简形的左半部分不是单位矩阵 $I_n$，则矩阵 $A$ 不可逆。
+
+## 求线性基
+
+[[线性基|线性基]]
+
+
+## 求异或方程组
+
+ #待完成
+
+```cpp
+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;
+}
+```