高斯消元
高斯消元用于在 \(O(n^3)\) 时间内实现各种矩阵相关的计算，如矩阵行列式，矩阵求逆，解线性方程组等。
高斯消元的核心是：将原矩阵通过初等行变换，转换为行最简形。
行最简形
若行阶梯形矩阵满足：
- 非零行的首非零元为 \(1\)；
- 首非零元所在的列的其他元均为 \(0\)；
则称此矩阵为行最简形矩阵。例如：
\[
A=
\left(
\begin{matrix}
1& 0& -1& 0& 4\\
0& 1& -1& 0& 3\\
0& 0& 0& 1& -3\\
0& 0& 0& 0& 0\\
\end{matrix}
\right)
\]
对于给矩阵 \(A\)，我们从左到右考虑 \(A\) 的每一列，并将前 \(i\) 列消元的结果存储在前 \(cur\) 行。显然有 \(cur\le i\)，当且仅当前 \(i\) 列满秩，等号成立。
- 如果第 \(i\) 列全都为 \(0\)，则 \(cur\) 不变，直接跳过该列；
- 否则：
- 选取 \([cur,n]\) 行中该列绝对值最大的一行，将它交换到第 \(cur\) 行（主元法，为了减小误差）；
- 将第 \(cur\) 行除以 \(a_{cur,i}\)，使 \(a_{cur,i}\) 变为 \(1\)；
- 用第 \(cur\) 行消去其它行的第 \(i\) 列；
- \(cur\leftarrow cur+1\)；
经过这些消元之后，矩阵将会变为行最简形。当矩阵满秩时，其左侧是一个单位矩阵。若矩阵不满秩，则所有空行将堆积在矩阵的下部。
关于主元法
如果主元 \(a_{cur,i}\) 的绝对值太小，后续消元步骤中需要将其他元素除以该主元，导致乘数极大，这种放大效应会使原有数据中的微小舍入误差被显著放大。主元法通过选择较大的主元，确保除数不会过小，从而控制乘数的量级。
由 deepseek 回答。
int cur = 1;
for(int i = 1; i <= n; i++) {
int mx = cur;
for(int j = cur + 1; j <= n; j++) {
if(abs(a[j][i]) > abs(a[mx][i])) mx = j;
}
if(abs(a[mx][i]) < eps) continue;
for(int j = 1; j <= n + 1; j++) {
swap(a[cur][j], a[mx][j]);
}
for(int j = n + 1; j >= i; j--) {
a[cur][j] /= a[cur][i];
}
for(int j = 1; j <= n; j++) {
if(j == cur) continue;
for(int k = n + 1; k >= i; k--) {
a[j][k] -= a[j][i] * a[cur][k];
}
}
++cur;
}
求解线性方程组
题意
已知 \(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}\]
请根据输入的数据，编程输出方程组的解的情况。
如果有唯一解，则输出解。你的结果被认为正确，当且仅当对于每一个 \(x_i\) 而言结果值与标准答案值的绝对误差不超过 \(0.01\)。
如果方程组无解输出 \(-1\)；
如果有无穷多实数解，输出 \(0\)；
对于 \(100\%\) 的数据，\(1 \le n \le 50\)。对于 \(\forall 1\le i, j \le n\)，有 \(|a_{i, j}| \le 100\)，\(|b_i| \le 300\)。
将所有未知数的系数排成一个矩阵，然后将常数项放在右边，得到的矩阵叫做增广矩阵：
\[
\left(
\begin{matrix}
a_{1,1}& a_{1,2}& a_{1,3}& \cdots& a_{1,m}\\
a_{2,1}& a_{2,2}& a_{2,3}& \cdots& a_{2,m}\\
a_{3,1}& a_{3,2}& a_{3,3}& \cdots& a_{3,m}\\
& \vdots & & \ddots& \vdots\\
a_{n,1}& a_{n,2}& a_{n,3}& \cdots& a_{n,m}\\
\end{matrix}
\ \middle|\ \begin{matrix}
b_1\\
b_2\\
b_3\\
\vdots\\
b_n
\end{matrix}
\right)
\]
我们希望增广矩阵的左侧变为一个单位矩阵，这样右侧的常数项就是方程的解：
\[
\left(
\begin{matrix}
1& 0& 0& \cdots& 0& 0\\
0& 1& 0& \cdots& 0& 0\\
0& 0& 1& \cdots& 0& 0\\
& \vdots & & \ddots& \vdots& \vdots\\
0& 0& 0& \cdots& 1& 0\\
0& 0& 0& \cdots& 0& 1\\
\end{matrix}
\ \middle|\ \begin{matrix}
x_1\\
x_2\\
x_3\\
\vdots\\
x_{n-1}\\
x_n
\end{matrix}
\right)
\]
我们只需要对原始的增广矩阵进行高斯消元，使其左半部分变为行最简形即可。
对于无穷解的情况，矩阵下半部分会有空行；对于无解的情况，将会存在一行，系数项全都为 \(0\)，而常数项不为 \(0\)（\(0\ne 0\)）。特判这两种情况即可。
代码
| #include<iostream>
#include<iomanip>
#include<cstdio>
#include<cmath>
#define ld long double
using namespace std;
const int N = 55;
const ld eps = 1e-6l;
int n;
ld a[N][N];
int main() {
cin >> n;
for(int i = 1; i <= n; i++) {
for(int j = 1; j <= n + 1; j++) {
cin >> a[i][j];
}
}
int cur = 1;
for(int i = 1; i <= n; i++) {
int mx = cur;
for(int j = cur + 1; j <= n; j++) {
if(abs(a[j][i]) > abs(a[mx][i])) mx = j;
}
if(abs(a[mx][i]) < eps) continue;
for(int j = 1; j <= n + 1; j++) {
swap(a[cur][j], a[mx][j]);
}
for(int j = n + 1; j >= i; j--) {
a[cur][j] /= a[cur][i];
}
for(int j = 1; j <= n; j++) {
if(j == cur) continue;
for(int k = n + 1; k >= i; k--) {
a[j][k] -= a[j][i] * a[cur][k];
}
}
++cur;
}
for(int i = 1; i <= n; i++) {
if(abs(a[i][n + 1]) < eps) continue;
bool flag = true;
for(int j = i; j <= n; j++) {
if(abs(a[i][j]) > eps) {
flag = false;
break;
}
}
if(flag) {
cout << -1 << endl;
return 0;
}
}
for(int i = 1; i <= n; i++) {
if(abs(a[i][i]) < eps) {
cout << 0 << endl;
return 0;
}
}
for(int i = 1; i <= n; i++) {
printf("x%d=%.2Lf\n", i, a[i][n + 1]);
}
return 0;
}
/*
4
0 0 1 1 1
0 0 1 1 2
0 0 0 0 0
0 0 0 0 0
10
1 1 2 5 4 3 4 6 2 1 1
5 4 8 1 3 6 4 1 4 3 2
6 4 2 5 1 3 1 4 2 6 3
1 9 5 6 7 3 2 8 4 2 4
2 9 4 8 6 7 5 3 1 5 5
5 6 8 4 7 1 6 3 4 2 6
8 1 4 7 6 5 2 6 9 3 7
9 4 7 5 6 3 2 1 4 7 8
6 4 8 5 7 3 2 1 9 4 9
10 3 8 17 14 11 10 18 13 5 9
10
1 1 2 5 4 3 4 0 2 1 1
5 4 8 1 3 6 4 0 4 3 2
6 4 2 5 1 3 1 0 2 6 3
1 9 5 6 7 3 2 0 4 2 4
2 9 4 8 6 7 5 0 1 5 5
5 6 8 4 7 1 6 0 4 2 6
8 1 4 7 6 5 2 0 9 3 7
9 4 7 5 6 3 2 0 4 7 8
6 4 8 5 7 3 2 0 9 4 9
9 3 6 3 9 0 1 0 5 1 10
2 2 4 10 8 6 8 12 4 2 2
8 1 4 7 6 5 2 6 9 3 7
10 3 8 17 14 11 10 18 13 5 9
*/
|
在一些题目中，各个变量的关系较为复杂，没有简单的递推式。我们就可以使用高斯消元，解出答案。
见概率与期望
求解行列式 & 数论高消
题意
给定一个 \(n\) 阶行列式 \(A\) 和模数 \(p\)，求 \(|A|\bmod p\)。不保证 \(p\) 是质数。
我们只需要把 \(A\) 消成上三角矩阵即可，考虑如何在没有逆元的情况下求出精确解。我们可以使用欧几里德算法，使用模运算代替有理数除法。我们依次枚举每一列 \(i\)，然后尝试把该列的 \([i+1,n]\) 行都消为 \(0\)。对于每一个 \(k\in [i+1,n]\)，我们只需要将 \(a_{i,i}\) 和 \(a_{k,i}\) 进行辗转相除（所在的整行也同时进行运算），一定能把其中一项消成 \(0\)，然后把第 \(i\) 列非零的那一行交换到第 \(i\) 行即可。最终得到一个上三角矩阵。
容易发现时间消耗不超过 \(O(n^3\log V)\)。然而我们可以进一步证明时间复杂度为 \(O(n^3)\)。
时间复杂度分析
能够证明，使用欧几里德算法求出 \(n\) 个数的 \(\gcd\)，均摊进行 \(O(n+\log V)\) 次取模操作（证明详见欧几里德算法）。
每对 \(a_{i,i}\) 和 \(a_{k,i}\) 进行一次取模操作，都会对 \(i,k\) 两行进行变换，产生 \(O(n)\) 的时间复杂度。因此总时间为 \(O(n^2\log V+n^3)\)。
数论高消的核心是：用取模运算（本质是加减）代替有理数除法，可以在相同的时间复杂度内将原矩阵消为上三角矩阵。
代码
| #include<iostream>
#define int long long
using namespace std;
const int N = 610;
int n, ans = 1, MOD;
int a[N][N];
int op[2][2];
void md(int &x) {
x = (x % MOD + MOD) % MOD;
}
signed main() {
cin >> n >> MOD;
for(int i = 1; i <= n; i++) {
for(int j = 1; j <= n; j++) {
cin >> a[i][j];
}
}
for(int i = 1; i <= n; i++) {
for(int j = i + 1; j <= n; j++) {
if(a[i][i] < a[j][i]) {
for(int k = n; k >= i; k--) {
swap(a[i][k], a[j][k]);
}
ans = MOD - ans;
}
if(!a[i][i]) continue;
while(a[j][i]) {
for(int k = n; k >= i + 1; k--) {
int x = a[i][k], y = a[j][k];
a[i][k] = y;
a[j][k] = x - (a[i][i] / a[j][i]) * y % MOD;
md(a[j][k]);
}
a[i][i] %= a[j][i];
swap(a[i][i], a[j][i]);
ans = MOD - ans;
}
}
}
for(int i = 1; i <= n; i++) ans = ans * a[i][i] % MOD;
cout << ans << endl;
return 0;
}
|