列主元GaussJordan消元法
struct Matrix {
static const int MAXN = 505;
static const int MAXM = 505;
int n, m;
long double A[MAXN][MAXM];
long double x[MAXM];
int pos[MAXM];
Matrix() {
ms(A);
}
void show() {
printf("Matrix n=%d m=%d
", n, m);
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= m; ++j)
printf("%6.2f ", A[i][j]);
printf("| %6.2f
", A[i][m + 1]);
}
}
// 列主元GaussJordan消元法
int GaussJordan() {
// 把矩阵化简为行最简形矩阵
int row = 0; // 当前已有的非零行的数量
for(int col = 1; col <= m && row <= n; ++col) {
// 在剩余行中，寻找列主元所在行
int pivot = row + 1;
for(int i = pivot + 1; i <= n; ++i) {
if(abs(A[i][col]) > abs(A[pivot][col]))
pivot = i;
}
if(abs(A[pivot][col]) <= EPS) {
pos[col] = -1; // 在剩余行中，列主元为0，该变量是自由变量
continue;
}
pos[col] = ++row;
if(pivot != row) {
// 行交换，行列式的值变成相反数
for (int j = m + 1; j >= col; --j)
swap(A[pivot][j], A[row][j]);
}
// 列主元置1
for (int j = m + 1; j >= col; --j)
A[row][j] /= A[row][col];
// 在所有其他行中，消去该变量
for(int i = 1; i <= n; ++i) {
if(i == row || abs(A[i][col]) <= EPS)
continue; // 跳过当前行和不存在该变量的行
for (int j = m + 1; j >= col; --j)
A[i][j] -= A[row][j] * A[i][col];
}
}
// 检查无解的情形
for(int i = 1; i <= n; ++i) {
if(abs(A[i][m + 1]) <= EPS)
continue; // 空行/非空行 = 0值，均不影响
bool ok = 0;
for(int j = 1; j <= m; ++j) {
if(abs(A[i][j]) <= EPS)
continue;
ok = 1;
break;
}
if(!ok)
return -1; // 失败：空行 = 非0值，返回-1
}
// 回代
for(int j = 1; j <= m; ++j) {
if(pos[j] == -1) {
x[j] = 0; // 自由变量，可以自由取值，默认置0
continue;
}
x[j] = A[pos[j]][m + 1]; // 列主元置1，自由变量均置0的情形
}
return row; //成功：返回矩阵的秩
}
} mat;
Gauss转换为行最简形，传入原本增广矩阵的左侧A（记得m是A的列数，不是增广矩阵的列数），传出一个位置向量，表示原本矩阵的新矩阵的第i行应原本矩阵的第pos[i]行，传出一个系数向量，表示在交换位置之后，位于这个新位置的行所乘的k。
Solve传入一个位置向量，和原本增广矩阵的右侧b，对这个增广矩阵进行求解，解出所有约束变量的取值（自由变量是可以任意取值的）。
【小贴士】
- 未受到方程组约束的自由变量，可以任意取值。对于自由变量的任意一组取值，都可以对约束变量进行适当赋值得到一组解。常见于统计模2意义下的高斯消元的解的个数，注意首先要保证方程组有解。
- 高斯消元法分为两个步骤，首先先把 (Ax=b) 的 (A) 转换成行最简形，复杂度为 (O(n^3)) ，然后代入正确位置的 (b) 进行求解，复杂度为 (O(n)) ，在 CF113D 中， (A) 是不变的，可以利用同一个 (A) 对不同的 (b) 分别解出 (x) 。（也可以使用LU分解法）
两个步骤高斯消元法：
struct Matrix {
static const int MAXN = 505;
static const int MAXM = 505;
int n, m;
long double A[MAXN][MAXM];
long double x[MAXM];
int pos[MAXM];
Matrix() {
ms(A);
}
void show() {
printf("Matrix n=%d m=%d
", n, m);
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= m; ++j)
printf("%6.2f ", (double)A[i][j]);
printf("| %6.2f
", (double)A[i][m + 1]);
}
}
struct Operation {
int type;
int x, y;
long double k;
} op;
vector<Operation> vec;
// 列主元GaussJordan消元法
int GaussJordan1() {
vec.clear();
// 把矩阵化简为行最简形矩阵
int row = 0; // 当前已有的非零行的数量
for(int col = 1; col <= m && row <= n; ++col) {
// 在剩余行中，寻找列主元所在行
int pivot = row + 1;
for(int i = pivot + 1; i <= n; ++i) {
if(abs(A[i][col]) > abs(A[pivot][col]))
pivot = i;
}
if(abs(A[pivot][col]) <= EPS) {
pos[col] = -1; // 在剩余行中，列主元为0，该变量是自由变量
continue;
}
pos[col] = ++row;
if(pivot != row) {
vec.pb({1, pivot, row, 0});
// 行交换，行列式的值变成相反数
for (int j = m + 1; j >= col; --j)
swap(A[pivot][j], A[row][j]);
}
vec.pb({2, row, 0, A[row][col]});
// 列主元置1
for (int j = m + 1; j >= col; --j)
A[row][j] /= A[row][col];
// 在所有其他行中，消去该变量
for(int i = 1; i <= n; ++i) {
if(i == row || abs(A[i][col]) <= EPS)
continue; // 跳过当前行和不存在该变量的行
vec.pb({3, i, row, A[i][col]});
for (int j = m + 1; j >= col; --j)
A[i][j] -= A[row][j] * A[i][col];
}
}
return row; //成功：返回矩阵的秩
}
int GaussJordan2() {
for(const Operation &op : vec) {
int type = op.type, x = op.x, y = op.y;
long double k = op.k;
if(type == 1) {
swap(A[x][m + 1], A[y][m + 1]);
continue;
}
if(type == 2) {
A[x][m + 1] /= k;
continue;
}
if(type == 3) {
A[x][m + 1] -= A[y][m + 1] * k;
continue;
}
}
// 检查无解的情形
for(int i = 1; i <= n; ++i) {
if(abs(A[i][m + 1]) <= EPS)
continue; // 空行/非空行 = 0值，均不影响
bool ok = 0;
for(int j = 1; j <= m; ++j) {
if(abs(A[i][j]) <= EPS)
continue;
ok = 1;
break;
}
if(!ok)
return -1; // 失败：空行 = 非0值，返回-1
}
// 回代
for(int j = 1; j <= m; ++j) {
if(pos[j] == -1) {
x[j] = 0; // 自由变量，可以自由取值，默认置0
continue;
}
x[j] = A[pos[j]][m + 1]; // 列主元置1，自由变量均置0的情形
}
return 1; //成功
}
} mat;
求解行列式：
const double EPS = 1E-9;
int n;
vector<vector<double> > a(n, vector<double>(n));
double 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 (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;
- 无向图的生成树计数：矩阵树定理 在操作中一般直接计算上面这个行列式的值。