Gaussian Elimination

Gauss-Jordan

Ma trận mở rộng $[A \mid b]$ cỡ $n \times (n + 1)$. Duyệt từng cột $c$ với con trỏ hàng $r$:

1. Trong các hàng $r..n-1$, chọn hàng có $|a_{ic}|$ lớn nhất làm hàng trụ. Nếu nó nhỏ hơn $\varepsilon$, cột $c$ không có trụ (biến tự do), sang cột sau.
2. Đổi hàng trụ lên vị trí $r$, rồi khử cột $c$ ở mọi hàng khác (Gauss-Jordan), tăng $r$.

Sau cùng: nếu một hàng $i \ge r$ có vế phải khác 0 thì hệ vô nghiệm (hàng $0 = b_i$); ngược lại nếu $r < n$ thì vô số nghiệm; còn lại nghiệm duy nhất $x_i = b_i / a_{ii}$. Kiểm vô nghiệm trước, vì hệ có thể vừa có biến tự do vừa mâu thuẫn.

Chọn trụ lớn nhất giữ sai số nhỏ, vì chia cho số gần 0 khuếch đại sai số làm tròn. Ngưỡng $\varepsilon$ khoảng $10^{-8}$ với double, nhỏ hơn với long double. Khi in, đổi giá trị có trị tuyệt đối dưới nửa đơn vị cuối thành 0 để khỏi in -0.00.

Trong hình, ô cam là phần tử trụ, hàng xanh là hàng trụ. Sau mỗi lần khử, cột trụ chỉ còn một ô khác 0.

Ví dụBa Trạng Thái Của Hệ

Đề. $n \le 50$ phương trình $n$ ẩn, hệ số nguyên nhỏ. In $-1$ nếu vô nghiệm, $0$ nếu vô số nghiệm, ngược lại in từng xi= giá trị với hai chữ số thập phân.

Lời giải. Khử Gauss-Jordan trên ma trận mở rộng $n \times (n + 1)$. Cột không có trụ (mọi phần tử gần 0) là một ẩn tự do. Sau khi khử, một hàng toàn 0 bên trái mà vế phải khác 0 nghĩa là vô nghiệm; không vô nghiệm mà có ẩn tự do thì vô số nghiệm; còn lại nghiệm duy nhất. Tính bằng long double với $\varepsilon = 10^{-11}$, để giá trị đúng bằng nửa đơn vị làm tròn (như $0{,}125$) không bị sai số đẩy sang phía kia.

Bước 1: chọn trụ lớn nhất. r là số hàng trụ đã dùng. Với mỗi cột c, tìm hàng có |a[p][c]| lớn nhất từ hàng r trở đi. Trụ lớn nhất làm hệ số f không vượt 1 về trị tuyệt đối, sai số không bị khuếch đại. Trụ gần 0 thì bỏ qua cột, r không tăng.

int p = r;
for (int i = r; i < n; i++)
    if (fabsl(a[i][c]) > fabsl(a[p][c])) p = i;
if (fabsl(a[p][c]) < EPS) continue;
swap(a[p], a[r]);

Bước 2: khử mọi hàng khác. Khử cột c ở mọi hàng khác, cả phía trên: đây là Gauss-Jordan, xong là đọc được nghiệm luôn.

for (int i = 0; i < n; i++)
    if (i != r) {
        ld f = a[i][c] / a[r][c];
        for (int j = c; j <= n; j++) a[i][j] -= f * a[r][j];
    }
r++;

Bước 3: phân loại. Các hàng từ r trở đi có vế trái toàn 0; vế phải khác 0 ở hàng nào là vô nghiệm. r < n mà không vô nghiệm là vô số nghiệm. r == n thì $x_c = a_{cn} / a_{cc}$; số rất gần 0 đưa về 0 để khỏi in -0.00.

for (int i = r; i < n; i++)
    if (fabsl(a[i][n]) > EPS) {
        cout << -1 << '\n';
        return 0;
    }
if (r < n) {
    cout << 0 << '\n';
    return 0;
}

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define ld long double

const ld EPS = 1e-11;
int n;
ld a[55][55];

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    cin >> n;
    for (int i = 0; i < n; i++)
        for (int j = 0; j <= n; j++) cin >> a[i][j];
    int r = 0;
    for (int c = 0; c < n && r < n; c++) {
        int p = r;
        for (int i = r; i < n; i++)
            if (fabsl(a[i][c]) > fabsl(a[p][c])) p = i;
        if (fabsl(a[p][c]) < EPS) continue;
        swap(a[p], a[r]);
        for (int i = 0; i < n; i++)
            if (i != r) {
                ld f = a[i][c] / a[r][c];
                for (int j = c; j <= n; j++) a[i][j] -= f * a[r][j];
            }
        r++;
    }
    for (int i = r; i < n; i++)
        if (fabsl(a[i][n]) > EPS) {
            cout << -1 << '\n';
            return 0;
        }
    if (r < n) {
        cout << 0 << '\n';
        return 0;
    }
    cout << fixed << setprecision(2);
    for (int i = 0; i < n; i++) {
        ld x = a[i][n] / a[i][i];
        if (fabsl(x) < 0.0005) x = 0;
        cout << "x" << i + 1 << "=" << (double)x << '\n';
    }
}

ChạyInput2
1 1 5
2 -1 1

Outputx1=2.00
x2=3.00

Độ phức tạp: $O(n^3)$.

Ma trận nghịch đảo

Theo mô-đun nguyên tố $p$, mọi số khác 0 có nghịch đảo, nên khử Gauss chạy chính xác: hàng trụ là hàng bất kỳ có phần tử khác 0, chia bằng nhân nghịch đảo (lũy thừa nhanh), không cần $\varepsilon$.

Ma trận nghịch đảo: khử Gauss-Jordan trên $[A \mid I]$. Khi bên trái thành $I$, bên phải là $A^{-1}$, vì dãy phép biến đổi hàng chính là nhân trái với $A^{-1}$. Gặp cột không tìm được trụ thì $A$ không khả nghịch. Ở cột $c$, hàng trụ đã bằng 0 ở các cột trước $c$ nên chỉ cần khử từ cột $c$ trở đi.

Ví dụMa Trận Đảo

Đề. Ma trận $n \times n$, $n \le 400$, phần tử theo mô-đun $10^9 + 7$. In $A^{-1}$ hoặc No Solution.

Lời giải. Ghép ma trận đơn vị vào bên phải thành $[A \mid I]$ cỡ $n \times 2n$, khử Gauss-Jordan theo mô-đun. Mỗi phép biến đổi hàng là nhân trái với một ma trận, nên khi bên trái thành $I$ thì tích các phép đó là $A^{-1}$, và bên phải chính là $A^{-1}$. Theo mô-đun nguyên tố thì mọi phần tử khác 0 đều có nghịch đảo, nên trụ là bất kỳ phần tử khác 0 nào; không tìm được trụ ở một cột là $A$ suy biến.

Bước 1: ma trận mở rộng. a[i][n + i] = 1 đặt ma trận đơn vị bên phải.

Bước 2: tìm trụ. Hàng p từ c trở xuống có phần tử khác 0, đổi lên hàng c.

int p = c;
while (p < n && !a[p][c]) p++;
if (p == n) {
    cout << "No Solution\n";
    return 0;
}
swap(a[p], a[c]);

Bước 3: chuẩn hóa và khử. Nhân cả hàng trụ với nghịch đảo của trụ để trụ thành 1, rồi trừ f lần hàng trụ khỏi mọi hàng khác để cột c chỉ còn số 1 ở hàng c. Vòng j chạy từ c vì các cột trước của hàng trụ đã bằng 0. Phép trừ cộng thêm MOD để kết quả trong $[0, p)$.

ll inv = power(a[c][c], MOD - 2);
for (int j = c; j < 2 * n; j++) a[c][j] = a[c][j] * inv % MOD;
for (int i = 0; i < n; i++)
    if (i != c && a[i][c]) {
        ll f = a[i][c];
        for (int j = c; j < 2 * n; j++) a[i][j] = (a[i][j] - f * a[c][j] % MOD + MOD) % MOD;
    }

#include <bits/stdc++.h>
using namespace std;
#define ll long long

const ll MOD = 1e9 + 7;
int n;
ll a[405][810];

ll power(ll x, ll k) {
    ll r = 1;
    x %= MOD;
    while (k) {
        if (k & 1) r = r * x % MOD;
        x = x * x % MOD;
        k >>= 1;
    }
    return r;
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    cin >> n;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) cin >> a[i][j];
        a[i][n + i] = 1;
    }
    for (int c = 0; c < n; c++) {
        int p = c;
        while (p < n && !a[p][c]) p++;
        if (p == n) {
            cout << "No Solution\n";
            return 0;
        }
        swap(a[p], a[c]);
        ll inv = power(a[c][c], MOD - 2);
        for (int j = c; j < 2 * n; j++) a[c][j] = a[c][j] * inv % MOD;
        for (int i = 0; i < n; i++)
            if (i != c && a[i][c]) {
                ll f = a[i][c];
                for (int j = c; j < 2 * n; j++) a[i][j] = (a[i][j] - f * a[c][j] % MOD + MOD) % MOD;
            }
    }
    for (int i = 0; i < n; i++) {
        for (int j = n; j < 2 * n; j++) cout << a[i][j] << ' ';
        cout << '\n';
    }
}

ChạyInput2
1 2
3 5

Output1000000002 2
3 1000000006

Ước lượng phép tính: $n \cdot n \cdot 2n \approx 1{,}3 \cdot 10^8$ phép nhân theo mô-đun.

GF(2) và bitset

Phương trình XOR (bật tắt đèn, kiểm chẵn lẻ) là hệ trên GF(2): cộng là XOR, nhân là AND, mọi hệ số khác 0 đều bằng 1. Một hàng gồm $n$ hệ số và vế phải là một bitset, nên khử một hàng tốn $O(n / 64)$, cả hệ $O(n^3 / 64)$.

bitset nhanh hơn mảng bool khoảng 13 lần: mỗi phép XOR xử lý 64 hệ số cùng lúc. Quy tắc: hệ mod 2 thì hàng luôn là bitset.

Khi phương trình đến lần lượt và cần biết lúc nào đủ để xác định nghiệm, chèn từng phương trình như chèn vào cơ sở XOR: duyệt biến từ nhỏ tới lớn, gặp biến đã có hàng trụ thì XOR với hàng đó, gặp biến chưa có thì nhận phương trình làm hàng trụ của biến ấy. Hạng đạt $n$ là lúc nghiệm duy nhất; sau đó khử ngược từ biến lớn về nhỏ để đọc nghiệm.

Ví dụMốc Đủ Dữ Kiện

Đề. $N \le 1000$ bit chưa biết, $M \le 2000$ phép đo theo thứ tự, mỗi phép đo cho XOR của một tập bit (dữ liệu nhất quán). Tìm $K$ nhỏ nhất để $K$ phép đo đầu xác định duy nhất mọi bit và in các bit (?y7M# cho 1, Earth cho 0), hoặc Cannot Determine.

Lời giải. Mỗi phép đo là một phương trình trên GF(2), lưu thành bitset gồm $N$ hệ số và vế phải ở bit $N$. Chèn lần lượt như chèn vào cơ sở XOR; phép đo đầu tiên làm hạng đạt $N$ là $K$. Sau đó khử ngược để mỗi hàng chỉ còn biến trụ của nó, và vế phải chính là giá trị của biến.

Bước 1: chèn một phép đo. row[j] là hàng trụ của biến $j$: bit thấp nhất bật của nó là $j$. Đi $j$ từ 0 lên: bit $j$ tắt thì bỏ qua; bật mà chưa có trụ thì nhận làm trụ và dừng; bật mà đã có trụ thì XOR để tắt bit $j$, rồi đi tiếp. Phép đo đi hết mà không thành trụ là phụ thuộc các phép đo trước, không thêm thông tin. Dữ liệu nhất quán nên không cần kiểm mâu thuẫn.

for (int j = 0; j < n; j++) {
    if (!v[j]) continue;
    if (!has[j]) {
        has[j] = 1, row[j] = v, rk++;
        break;
    }
    v ^= row[j];
}
if (rk == n) K = k;

Vòng đọc có điều kiện !K nên dừng ngay khi đủ hạng; các phép đo còn lại không cần đọc.

Bước 2: khử ngược. $j$ từ $N - 1$ về 0. Lúc xử lý $j$, các hàng $t > j$ đã chỉ còn bit $t$, nên XOR với chúng xóa sạch các bit lớn hơn $j$ trong row[j]. Xong thì row[j][n] là giá trị của biến $j$.

for (int j = n - 1; j >= 0; j--)
    for (int t = j + 1; t < n; t++)
        if (row[j][t]) row[j] ^= row[t];

#include <bits/stdc++.h>
using namespace std;
#define ll long long

const int N = 1005;
int n, m, rk;
bitset<N> row[N];
bool has[N];

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    cin >> n >> m;
    int K = 0;
    for (int k = 1; k <= m && !K; k++) {
        string s;
        int b;
        cin >> s >> b;
        bitset<N> v;
        for (int j = 0; j < n; j++) v[j] = s[j] == '1';
        v[n] = b;
        for (int j = 0; j < n; j++) {
            if (!v[j]) continue;
            if (!has[j]) {
                has[j] = 1, row[j] = v, rk++;
                break;
            }
            v ^= row[j];
        }
        if (rk == n) K = k;
    }
    if (!K) {
        cout << "Cannot Determine\n";
        return 0;
    }
    for (int j = n - 1; j >= 0; j--)
        for (int t = j + 1; t < n; t++)
            if (row[j][t]) row[j] ^= row[t];
    cout << K << '\n';
    for (int j = 0; j < n; j++) cout << (row[j][n] ? "?y7M#" : "Earth") << '\n';
}

ChạyInput3 4
110 1
011 1
101 0
001 1

Output4
?y7M#
Earth
?y7M#

Độ phức tạp: mỗi lần XOR hai bitset dài 1001 bit tốn khoảng 16 phép trên số 64 bit: tổng $O(MN^2 / 64)$.

Định thức

Định thức là tích các phần tử trên đường chéo sau khi đưa về tam giác trên, mỗi lần đổi hai hàng thì đổi dấu. Theo mô-đun nguyên tố thì khử như thường.

Khi mô-đun không nguyên tố, phần tử trụ có thể không có nghịch đảo, nên khử kiểu Euclid: với hai hàng $i$ và $j$, trừ $\lfloor a_{ii} / a_{ji} \rfloor$ lần hàng $j$ khỏi hàng $i$ rồi đổi chỗ, lặp tới khi $a_{ji} = 0$. Các phép này chỉ trừ bội nguyên, không chia. Trong một cột, giá trị trụ chỉ giảm qua các lần gọi Euclid, nên tổng số bước của cột là $O(n + \log p)$ và cả thuật toán $O(n^2 (n + \log p))$.

Ví dụĐịnh Thức Trên Vành

Đề. Ma trận $n \times n$, $n \le 600$, phần tử dưới $p \le 10^9 + 7$, $p$ có thể là hợp số. In $\det A \bmod p$.

Lời giải. $p$ có thể là hợp số nên không chia được. Đưa ma trận về tam giác trên bằng cách làm như thuật toán Euclid: với cột $i$ và mỗi hàng $j > i$, trừ $\lfloor a_{ii} / a_{ji} \rfloor$ lần hàng $j$ khỏi hàng $i$ rồi đổi chỗ hai hàng, tới khi $a_{ji} = 0$. Trừ bội của một hàng không đổi định thức, đổi hai hàng thì đổi dấu. Cuối cùng định thức là tích đường chéo, nhân thêm dấu.

Bước 1: hàng là con trỏ. row[i] trỏ tới hàng $i$; swap(row[i], row[j]) đổi hai hàng trong $O(1)$ thay vì chép 600 số.

Bước 2: một bước Euclid. Trừ t lần hàng j khỏi hàng i (chỉ từ cột i, các cột trước đã bằng 0), rồi đổi chỗ và đổi dấu. Phép trừ có thể ra số âm nên cộng p. Giống Euclid, số bước là $O(\log p)$.

while (row[j][i]) {
    ll t = row[i][i] / row[j][i];
    if (t)
        for (ll k = i; k < n; k++) {
            row[i][k] = (row[i][k] - t * row[j][k]) % p;
            if (row[i][k] < 0) row[i][k] += p;
        }
    swap(row[i], row[j]);
    det = (p - det) % p;
}

Bước 3: tích đường chéo. det = 1 % p ở đầu để đúng cả khi $p = 1$; sau mỗi cột nhân thêm row[i][i].

#include <bits/stdc++.h>
using namespace std;
#define ll long long

ll n, p, a[605][605], *row[605];

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    cin >> n >> p;
    for (ll i = 0; i < n; i++) {
        row[i] = a[i];
        for (ll j = 0; j < n; j++) cin >> a[i][j];
    }
    ll det = 1 % p;
    for (ll i = 0; i < n; i++) {
        for (ll j = i + 1; j < n; j++)
            while (row[j][i]) {
                ll t = row[i][i] / row[j][i];
                if (t)
                    for (ll k = i; k < n; k++) {
                        row[i][k] = (row[i][k] - t * row[j][k]) % p;
                        if (row[i][k] < 0) row[i][k] += p;
                    }
                swap(row[i], row[j]);
                det = (p - det) % p;
            }
        det = det * row[i][i] % p;
    }
    cout << det << '\n';
}

ChạyInput2 8
1 2
3 4

Output6

Độ phức tạp: $O(n^3 + n^2 \log p)$.

Lập hệ kỳ vọng

Bước đi ngẫu nhiên trên đồ thị có chu trình không DP theo thứ tự được, nhưng kỳ vọng số lần ghé mỗi đỉnh thỏa một hệ tuyến tính. Gọi $x_v$ là kỳ vọng số lần quá trình đứng ở $v$:

Hệ có nghiệm duy nhất khi quá trình chắc chắn dừng (có xác suất dừng dương). Xác suất dừng tại $v$ bằng $x_v$ nhân xác suất dừng mỗi lần ghé. Với $n$ vài trăm, Gauss $O(n^3)$ là đủ.

Ví dụĐiểm Dừng Của Hành Trình

Đề. Đồ thị vô hướng $N \le 300$ đỉnh. Bắt đầu tại 1; ở mỗi đỉnh, với xác suất $r = P/Q$ thì dừng, ngược lại đi theo một cạnh kề ngẫu nhiên. In xác suất dừng tại từng đỉnh.

Lời giải. Gọi $x_v$ là kỳ vọng số lần đứng ở $v$. Lần đứng ở $v$ hoặc là lúc bắt đầu ($v = 1$), hoặc đến từ một hàng xóm $u$: đứng ở $u$, không dừng (xác suất $1 - r$), và chọn đúng cạnh tới $v$ (xác suất $1 / \deg u$). Vậy

Mỗi lần đứng ở $v$ dừng với xác suất $r$, nên xác suất dừng tại $v$ là $r \cdot x_v$.

Bước 1: dựng hệ. Hàng $v$ có 1 ở cột $v$, trừ $(1 - r) / \deg u$ ở cột của mỗi hàng xóm $u$, vế phải là 1 khi $v$ là đỉnh xuất phát.

for (int v = 0; v < n; v++) {
    a[v][v] += 1;
    for (int u : adj[v]) a[v][u] -= (1 - r) / deg[u];
    a[v][n] = v == 0;
}

Bước 2: khử. Gauss-Jordan với trụ lớn nhất như ví dụ đầu bài; f == 0 thì bỏ qua cả hàng, vì ma trận thưa ở những bước đầu. Đỉnh không tới được từ 1 có phương trình không có vế phải và hệ cho $x_v = 0$, không cần xử lý riêng.

Bước 3: in. $r \cdot x_v$; kết quả gần 0 đưa về 0 để khỏi in -0.000000000.

double x = r * a[v][n] / a[v][v];
if (fabs(x) < 5e-10) x = 0;

#include <bits/stdc++.h>
using namespace std;
#define ll long long

const double EPS = 1e-12;
int n, m, deg[305];
double P, Q, a[305][305];
vector<int> adj[305];

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    cin >> n >> m >> P >> Q;
    for (int i = 0; i < m; i++) {
        int u, v;
        cin >> u >> v;
        u--, v--;
        adj[u].push_back(v);
        adj[v].push_back(u);
        deg[u]++, deg[v]++;
    }
    double r = P / Q;
    for (int v = 0; v < n; v++) {
        a[v][v] += 1;
        for (int u : adj[v]) a[v][u] -= (1 - r) / deg[u];
        a[v][n] = v == 0;
    }
    for (int c = 0; c < n; c++) {
        int p = c;
        for (int i = c; i < n; i++)
            if (fabs(a[i][c]) > fabs(a[p][c])) p = i;
        swap(a[p], a[c]);
        if (fabs(a[c][c]) < EPS) continue;
        for (int i = 0; i < n; i++)
            if (i != c) {
                double f = a[i][c] / a[c][c];
                if (f == 0) continue;
                for (int j = c; j <= n; j++) a[i][j] -= f * a[c][j];
            }
    }
    cout << fixed << setprecision(9);
    for (int v = 0; v < n; v++) {
        double x = r * a[v][n] / a[v][v];
        if (fabs(x) < 5e-10) x = 0;
        cout << x << '\n';
    }
}

ChạyInput2 1 1 1
1 2

Output1.000000000
0.000000000

Độ phức tạp: $O(N^3)$.