高斯消元
别随便做初等列变换！要不然解出来解的顺序会变！
系数化为1的时候，要把主元系数单独存起来，否则就覆盖掉了
高斯-约旦消元
用 \(i \in [1,n]\) 标定当前求解第几个未知数（第几列），\(line\) 存储当前处理到第几行，每次进行：
- 初等行变换，找到第\(i\)个数（主元系数）绝对值最大的一行，若不为0，使用swap提上来；若为0，则\(line\)不变，\(i\)加一，跳过该元。
- 对角化，枚举每一行消去当前元。
- 把所有系数不为0的消干净后，用剩下几行（系数此时应该全为0）来判断无解/无数解。
- 若剩下几行常数项均为0，则说明其为自由元，可取任意值（最好取0或模数）；否则为无解。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52 | #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\)行，也可以用：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76 | // 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。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21 | 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;
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 | 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\)，求解其逆矩阵的方法如下：
- 构造 \(n \times 2n\) 的矩阵 \((A, I_n)\)，拼在右边；
- 用高斯消元法将其化简为最简形 \((I_n, A^{-1})\)，即可得到 \(A\) 的逆矩阵 \(A^{-1}\)。如果最终最简形的左半部分不是单位矩阵 \(I_n\)，则矩阵 \(A\) 不可逆。
求线性基
线性基
求异或方程组
待完成{: #待完成 .hash}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 | 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;
}
|