Skip to content

Instantly share code, notes, and snippets.

@GoBigorGoHome
Last active April 23, 2019 13:04
Show Gist options
  • Save GoBigorGoHome/11def4b30b10bd347c56485b75a0d500 to your computer and use it in GitHub Desktop.
Save GoBigorGoHome/11def4b30b10bd347c56485b75a0d500 to your computer and use it in GitHub Desktop.
F_p 上的高斯消元。
/* a tutorial on Gauss-Jordan elimination: https://cp-algorithms.com/linear_algebra/linear-system-gauss.html
From the article:
Strictly speaking, the method described below should be called "Gauss-Jordan", or Gauss-Jordan elimination,
because it is a variation of the Gauss method, described by Jordan in 1887.
*/
using num = modnum<1000003>;
using mat = vv<num>;
using vec = vector<num>;
// precondition: 系数矩阵不含0,否则需要选主元。
// a 是增广矩阵,n 行 n+1 列。
vec gauss_jordan(mat a, int n) {
for (int i = 0; i < n; i++) {
// num t = a[i][i].inv();
for (int j = i + 1; j <= n; j++) //ERROR-PRONE
a[i][j] /= a[i][i];
for (int j = 0; j < n; j++)
if (i != j) {
for (int k = i + 1; k <= n; k++) {
a[j][k] -= a[j][i] * a[i][k];
}
}
}
vec x(n);
rng (i, 0, n) x[i] = a[i][n];
return x;
}
// 带选主元的版本 (select a pivoting row)
using num = modnum<1000003>;
using mat = vv<num>;
using vec = vector<num>;
// a 是增广矩阵,n 行 n+1 列。
vec gauss_jordan(mat a, int n) {
for (int i = 0; i < n; i++) {
// select a pivoting row
int pivot = i;
while (pivot < n && a[pivot][i] == 0) ++pivot;
if (pivot == n) // 无解或解不唯一
return vec();
swap(a[i], a[pivot]);
for (int j = i + 1; j <= n; j++) //ERROR-PRONE
a[i][j] /= a[i][i];
for (int j = 0; j < n; j++) {
if (i != j) {
for (int k = i + 1; k <= n; k++) {
a[j][k] -= a[j][i] * a[i][k];
}
}
}
}
vec x(n);
rng (i, 0, n) x[i] = a[i][n];
return x;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment