Matrix Multiplication

Một bước truy hồi là một phép nhân ma trận

Gom hai số liền nhau thành một cột hai số, gọi là vectơ trạng thái, $w_i = \begin{pmatrix} F_{i+1} \\ F_i \end{pmatrix}$. Một bước truy hồi cho $F_{i+2} = 1 \cdot F_{i+1} + 1 \cdot F_i$ và $F_{i+1} = 1 \cdot F_{i+1} + 0 \cdot F_i$: mỗi tọa độ mới là một tổ hợp tuyến tính của các tọa độ cũ, tức tổng của các tọa độ cũ, mỗi tọa độ nhân với một hệ số cố định. Viết hệ số của từng tổ hợp thành một hàng:

Sau $n$ bước $w_n = Q^n w_0$. Chỗ lợi nằm ở một quan sát: hai bước liên tiếp cũng là một phép biến đổi cố định, gói được thành $Q \cdot Q = Q^2$; gói được $Q^2$ thì gói tiếp $Q^4$, $Q^8$, và để tới $Q^n$ chỉ cần chừng $\log_2 n$ lần gói.

Với hai ma trận vuông cỡ $k \times k$, tích $C = AB$ có $C_{ij} = \sum_{t} A_{it} B_{tj}$: ô $(i, j)$ là hàng $i$ của $A$ nhân từng ô với cột $j$ của $B$ rồi cộng lại. Định nghĩa này được chọn để $(AB)\,v = A\,(B\,v)$. Ba tính chất cần dùng: có ma trận đơn vị $I$ (đường chéo bằng 1, các ô khác bằng 0), chính là $A^0$; kết hợp, $(AB)C = A(BC)$, nên lũy thừa bằng bình phương chạy đúng như với số; và không giao hoán, $AB$ thường khác $BA$. Các lũy thừa của cùng một ma trận thì giao hoán, nhưng khi nhân nhiều ma trận khác nhau, thứ tự là một phần của đáp án.

Viết $n$ trong hệ nhị phân: $A^n$ là tích của các $A^{2^b}$ ứng với những bit $b$ bằng 1, và mỗi $A^{2^{b+1}}$ có từ $A^{2^b}$ bằng một lần bình phương. Với $n = 10 = 1010_2$: $Q^2 = \begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix}$, $Q^4 = \begin{pmatrix} 5 & 3 \\ 3 & 2 \end{pmatrix}$, $Q^8 = \begin{pmatrix} 34 & 21 \\ 21 & 13 \end{pmatrix}$, rồi $Q^{10} = Q^8 Q^2 = \begin{pmatrix} 89 & 55 \\ 55 & 34 \end{pmatrix}$.

Số 55 nằm ngoài đường chéo không phải ngẫu nhiên: quy nạp cho $Q^n = \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix}$ với mọi $n \ge 1$, vì nhân $\begin{pmatrix} a & b \\ c & d \end{pmatrix}$ với $Q$ bên phải cho $\begin{pmatrix} a + b & a \\ c + d & c \end{pmatrix}$. Với $n = 0$, $Q^0 = I$ có ô ngoài đường chéo bằng $0 = F_0$.

Mỗi phép nhân hai ma trận $k \times k$ tốn $k^3$ phép nhân-cộng; lũy thừa số mũ $n$ làm khoảng $\log_2 n$ lần bình phương và tối đa chừng ấy lần nhân vào kết quả, tổng $O(k^3 \log n)$. Về số học: mọi ô nằm trong $[0, 10^9 + 6]$, nên tích hai ô nhỏ hơn $2^{60}$ và lấy % MOD ngay sau mỗi phép nhân-cộng là an toàn. Modulo lớn hơn cỡ $3 \cdot 10^9$ thì cần __int128.

Ví dụChỉ Số Fibonacci

Đề. Cho $0 \le n \le 10^{18}$, in $F_n$ theo modulo $10^9 + 7$.

Lời giải. Như đã thấy, $Q^n = \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix}$ với $Q = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}$, nên chỉ cần lũy thừa $Q$ lên $n$ rồi đọc ô ngoài đường chéo.

Bước 1: nhân hai ma trận theo modulo. Ba vòng lặp, mỗi tích lấy dư ngay; hai số dưới $10^9 + 7$ nhân với nhau còn vừa ll.

mat mul(mat a, mat b) {
    ll k = a.size();
    mat c(k, vector<ll>(k, 0));
    for (ll i = 0; i < k; i++)
        for (ll t = 0; t < k; t++)
            for (ll j = 0; j < k; j++) c[i][j] = (c[i][j] + a[i][t] * b[t][j]) % MOD;
    return c;
}

Bước 2: lũy thừa nhanh. Kết quả bắt đầu bằng ma trận đơn vị. Duyệt các bit của $e$ từ thấp lên cao: a lần lượt là $A^1, A^2, A^4, \dots$, bit nào bằng 1 thì nhân a vào kết quả. Đây là lũy thừa nhanh của bài 08, áp cho ma trận.

mat power(mat a, ll e) {
    ll k = a.size();
    mat r(k, vector<ll>(k, 0));
    for (ll i = 0; i < k; i++) r[i][i] = 1;
    for (; e > 0; e >>= 1) {
        if (e & 1) r = mul(r, a);
        a = mul(a, a);
    }
    return r;
}

Bước 3: đọc đáp án. Ô [1][0] của $Q^n$ là $F_n$. Với $n = 0$ kết quả là ma trận đơn vị và ô đó cho đúng 0.

mat q = {{1, 1}, {1, 0}};
cout << power(q, n)[1][0] << '\n';

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

const ll MOD = 1e9 + 7;

mat mul(mat a, mat b) {
    ll k = a.size();
    mat c(k, vector<ll>(k, 0));
    for (ll i = 0; i < k; i++)
        for (ll t = 0; t < k; t++)
            for (ll j = 0; j < k; j++) c[i][j] = (c[i][j] + a[i][t] * b[t][j]) % MOD;
    return c;
}

mat power(mat a, ll e) {
    ll k = a.size();
    mat r(k, vector<ll>(k, 0));
    for (ll i = 0; i < k; i++) r[i][i] = 1;
    for (; e > 0; e >>= 1) {
        if (e & 1) r = mul(r, a);
        a = mul(a, a);
    }
    return r;
}

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

    ll n;
    cin >> n;
    mat q = {{1, 1}, {1, 0}};
    cout << power(q, n)[1][0] << '\n';
}

Input10

Output55

Khoảng 60 lần nhân ma trận $2 \times 2$, cỡ $10^3$ phép tính.

Truy hồi tuyến tính bậc k

Một truy hồi tuyến tính bậc $k$ cho $k$ số hạng đầu $a_0, \dots, a_{k-1}$ và $a_i = c_1 a_{i-1} + c_2 a_{i-2} + \dots + c_k a_{i-k}$ với $i \ge k$. Lấy trạng thái là cửa sổ gồm $k$ số hạng liên tiếp, số mới nhất ở trên cùng, $w_i = (a_{i+k-1}, \dots, a_{i})^T$. Tọa độ đầu của $w_{i+1}$ là $a_{i+k} = c_1 a_{i+k-1} + \dots + c_k a_i$, một tổ hợp với hệ số $c_1, \dots, c_k$; mỗi tọa độ còn lại chỉ là tọa độ liền trên của $w_i$ dời xuống một ô. Vậy

và $a_n$ là tọa độ cuối của $M^n w_0$ với $w_0 = (a_{k-1}, \dots, a_0)^T$. Với $a_i = a_{i-1} + 2 a_{i-3}$ và $a_0 = a_1 = a_2 = 1$: $(c_1, c_2, c_3) = (1, 0, 2)$, hệ số bằng 0 vẫn chiếm đúng một ô trong hàng đầu. Một bước cho $w_1 = (3, 1, 1)^T$, bước nữa cho $w_2 = (5, 3, 1)^T$; đó là $a_3 = 3$ và $a_4 = 5$.

Chương trình tổng quát đọc $k$, các hệ số $c_1, \dots, c_k$, các số hạng đầu $a_0, \dots, a_{k-1}$ rồi $n$; tọa độ cuối của $M^n w_0$ là tích của hàng cuối của $M^n$ với $w_0$.

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

const ll MOD = 1e9 + 7;

mat mul(mat a, mat b) {
    ll k = a.size();
    mat c(k, vector<ll>(k, 0));
    for (ll i = 0; i < k; i++)
        for (ll t = 0; t < k; t++)
            for (ll j = 0; j < k; j++) c[i][j] = (c[i][j] + a[i][t] * b[t][j]) % MOD;
    return c;
}

mat power(mat a, ll e) {
    ll k = a.size();
    mat r(k, vector<ll>(k, 0));
    for (ll i = 0; i < k; i++) r[i][i] = 1;
    for (; e > 0; e >>= 1) {
        if (e & 1) r = mul(r, a);
        a = mul(a, a);
    }
    return r;
}

ll k, n, c[105], a[105];

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

    cin >> k;
    for (ll i = 0; i < k; i++) cin >> c[i];
    for (ll i = 0; i < k; i++) cin >> a[i];
    cin >> n;
    mat m(k, vector<ll>(k, 0));
    for (ll j = 0; j < k; j++) m[0][j] = c[j] % MOD;
    for (ll i = 1; i < k; i++) m[i][i - 1] = 1;
    mat p = power(m, n);
    ll ans = 0;
    for (ll j = 0; j < k; j++) ans = (ans + p[k - 1][j] * (a[k - 1 - j] % MOD)) % MOD;
    cout << ans << '\n';
}

Input3
1 0 2
1 1 1
6

Output13

Dãy là $1, 1, 1, 3, 5, 7, 13, \dots$ nên $a_6 = 13$. Chi phí $O(k^3 \log n)$: với $n = 10^{18}$ và $k = 10$ chỉ khoảng $1{,}2 \cdot 10^5$ phép, với $k = 100$ khoảng $1{,}2 \cdot 10^8$, còn $k = 300$ đã cỡ $3 \cdot 10^9$. Khi $n$ chỉ cỡ $10^7$ mà $k$ nhỏ, cứ tính từng số hạng bằng $O(nk)$ cho gọn. Chỗ hay sai nhất là thứ tự: hàng đầu của $M$ phải khớp với thứ tự các tọa độ của $w_i$, và các hệ số bằng 0 không được “nén” khỏi hàng.

Chạy nhanh hơn. Chi phí tăng theo lập phương của $k$, nên việc đầu tiên là giữ $k$ nhỏ. Ngoài ra có mấy cách tiết kiệm. Thường ta chỉ cần $M^n w$ cho một vectơ $w$, không cần cả $M^n$: vì các $M^{2^b}$ giao hoán với nhau, nhân dần từng $M^{2^b}$ vào vectơ, mỗi lần $O(k^2)$ thay vì $O(k^3)$, và số phép nhân ma trận với ma trận chỉ còn các lần bình phương. Khi có nhiều truy vấn trên cùng một ma trận, dựng sẵn các $M^{2^b}$ một lần rồi mỗi truy vấn tốn $O(k^2 \log n)$. Ma trận chuyển thường rất thưa; một dòng if (!a[i][t]) continue; trong vòng lặp ngoài của phép nhân bỏ qua cả vòng trong. Riêng Fibonacci còn có fast doubling: $F_{2m+1} = F_{m+1}^2 + F_m^2$ và $F_{2m} = F_m(2F_{m+1} - F_m)$, chỉ theo dõi một cặp số. Khi $k$ lên vài trăm, người ta tính $x^n$ modulo đa thức đặc trưng (kỹ thuật Kitamasa) với $O(k^2 \log n)$; đó là chủ đề nâng cao hơn.

Ví dụNhững Viên Ngọc Ma Thuật

Đề. Xếp một hàng gồm số viên ngọc tùy ý, mỗi viên dạng thường dài 1 hoặc dạng mở rộng dài $m$ ($2 \le m \le 100$); thứ tự các dạng được phân biệt. Đếm số cấu hình có tổng độ dài đúng $n \le 10^{18}$, theo modulo $10^9 + 7$.

Lời giải. Gọi $f(i)$ là số cấu hình dài $i$, $f(0) = 1$. Viên đầu là dạng thường thì phần còn lại dài $i - 1$, là dạng mở rộng thì dài $i - m$; hai trường hợp rời nhau và phủ hết, nên $f(i) = f(i - 1) + f(i - m)$. Đây là truy hồi bậc $m$, và $f(i) = 1$ với $0 \le i < m$ (chỉ xếp được toàn viên thường). Với $n = 4$, $m = 3$: các cấu hình là $1+1+1+1$, $1+3$, $3+1$.

Bước 1: ma trận dời cửa sổ. Trạng thái là cửa sổ $m$ giá trị liên tiếp $(f(i), f(i - 1), \dots, f(i - m + 1))$. Hàng đầu của ma trận tính giá trị mới $f(i + 1) = f(i) + f(i - m + 1)$, nên có số 1 ở cột 0 và cột $m - 1$; các hàng sau dời cửa sổ xuống một ô.

mat t(m, vector<ll>(m, 0));
t[0][0] = t[0][m - 1] = 1;
for (ll i = 1; i < m; i++) t[i][i - 1] = 1;
vector<ll> w(m, 1);

Cửa sổ đầu là $(f(m - 1), \dots, f(0))$, toàn số 1.

Bước 2: nhân vectơ thay vì nhân hai ma trận. Ta chỉ cần $T^n w$, nên giữ vectơ w và duyệt bit của $n$: bit 1 thì w = T * w (chỉ $m^2$ phép), luôn bình phương T. Hàm mul bỏ qua các ô 0 của a, giúp ma trận thưa chạy nhanh hơn.

for (; n > 0; n >>= 1) {
    if (n & 1) w = mulvec(t, w);
    t = mul(t, t);
}

Bước 3: đọc đáp án. Sau $n$ bước, cửa sổ là $(f(n + m - 1), \dots, f(n))$, nên $f(n)$ nằm ở tọa độ cuối w[m - 1].

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

const ll MOD = 1e9 + 7;

mat mul(mat a, mat b) {
    ll k = a.size();
    mat c(k, vector<ll>(k, 0));
    for (ll i = 0; i < k; i++)
        for (ll t = 0; t < k; t++) {
            if (!a[i][t]) continue;
            for (ll j = 0; j < k; j++) c[i][j] = (c[i][j] + a[i][t] * b[t][j]) % MOD;
        }
    return c;
}

vector<ll> mulvec(mat a, vector<ll> v) {
    ll k = a.size();
    vector<ll> w(k, 0);
    for (ll i = 0; i < k; i++)
        for (ll j = 0; j < k; j++) w[i] = (w[i] + a[i][j] * v[j]) % MOD;
    return w;
}

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

    ll n, m;
    cin >> n >> m;
    mat t(m, vector<ll>(m, 0));
    t[0][0] = t[0][m - 1] = 1;
    for (ll i = 1; i < m; i++) t[i][i - 1] = 1;
    vector<ll> w(m, 1);
    for (; n > 0; n >>= 1) {
        if (n & 1) w = mulvec(t, w);
        t = mul(t, t);
    }
    cout << w[m - 1] << '\n';
}

Input4 3

Output3

Ma trận $100 \times 100$: 60 lần bình phương, mỗi lần $10^6$ phép, khoảng $6 \cdot 10^7$, vừa một giây; phần nhân vectơ chỉ $6 \cdot 10^5$.

Thêm tọa độ phụ: hằng số, tổng và chỉ số

Trạng thái không nhất thiết chỉ gồm các số hạng của dãy. Tọa độ nào cũng được đưa vào, miễn giá trị mới của nó tính được từ các tọa độ hiện có bằng một tổ hợp tuyến tính. Ba trường hợp dùng nhiều nhất:

Hằng số. Xét $a_i = a_{i-1} + a_{i-2} + 1$. Số $+1$ không phải là tổ hợp của $a$, nên ta thêm một tọa độ luôn bằng 1 với hàng $(0, 0, 1)$ để nó tự giữ nguyên: trạng thái $(a_{i+1}, a_i, 1)^T$ và ma trận có các hàng $(1, 1, 1)$, $(1, 0, 0)$, $(0, 0, 1)$. Cũng cách đó biến hàm bậc nhất $t \mapsto At + B$ (biến $t$ thành $At + B$) thành ma trận $\begin{pmatrix} A & B \\ 0 & 1 \end{pmatrix}$ tác động lên $(t, 1)^T$.

Tổng tiền tố. Muốn có $s_i = F_0 + \dots + F_{i-1}$, để ý $s_{i+1} = s_i + F_i$. Thêm $s_i$ vào trạng thái: $(F_{i+1}, F_i, s_i)^T$ với ma trận có các hàng $(1, 1, 0)$, $(1, 0, 0)$, $(0, 1, 1)$. Tổng trên đoạn $F_l + \dots + F_r$ là $s_{r+1} - s_l$, đúng ý của tổng tiền tố ở bài 02 nhưng bây giờ $r$ có thể tới $10^{18}$.

Chỉ số $i$ và các lũy thừa của nó. Nếu truy hồi có số hạng dạng $i$ hay $i^2$, thêm các tọa độ $i$, $i^2$ và hằng $1$: khi $i$ tăng một, $i + 1$ và $(i+1)^2 = i^2 + 2i + 1$ đều là tổ hợp tuyến tính của $(i^2, i, 1)$.

Điều ta không làm được là các số hạng phi tuyến như $a_{i-1} \cdot a_{i-2}$. Khi truy hồi chỉ gồm các phép nhân, ta chuyển sang theo dõi số mũ: nếu $b_i = b_{i-1} b_{i-2}$ thì mỗi $b_i$ là tích các lũy thừa của $b_0, b_1$ với số mũ theo một truy hồi cộng kiểu Fibonacci.

Ví dụLặp Hàm Tuyến Tính

Đề. Cho $f(t) = At + B$ với $1 \le A, B, x \le 10^9$. Bắt đầu từ $x$, áp $f$ đúng $n \le 10^{18}$ lần; in kết quả theo modulo $10^9 + 7$.

Lời giải. Thêm tọa độ hằng 1: trạng thái $(t, 1)^T$ và một lần áp $f$ là nhân với $\begin{pmatrix} A & B \\ 0 & 1 \end{pmatrix}$, vì $At + B \cdot 1$ là tọa độ trên và $1$ tự giữ nguyên. Áp $n$ lần là nhân với ma trận ấy mũ $n$. Có công thức đóng $A^n x + B \frac{A^n - 1}{A - 1}$ nhưng nó phải chia, và vỡ khi $A \equiv 1$; ma trận thì không phải bận tâm. Với ví dụ $3 \to 11 \to 27 \to 59 \to 123$.

Bước 1: nhân và lũy thừa ma trận. Hai hàm mul và power giống hệt ví dụ Fibonacci.

Bước 2: dựng ma trận rồi đọc tọa độ trên. Lấy $A$, $B$, $x$ theo modulo trước khi nhân để tích không tràn. Gọi $R$ là ma trận mũ $n$, thì $R \cdot (x, 1)^T$ có tọa độ trên $R_{00} x + R_{01}$.

mat r = power({{A % MOD, B % MOD}, {0, 1}}, n);
cout << (r[0][0] * (x % MOD) + r[0][1]) % MOD << '\n';

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

const ll MOD = 1e9 + 7;

mat mul(mat a, mat b) {
    ll k = a.size();
    mat c(k, vector<ll>(k, 0));
    for (ll i = 0; i < k; i++)
        for (ll t = 0; t < k; t++)
            for (ll j = 0; j < k; j++) c[i][j] = (c[i][j] + a[i][t] * b[t][j]) % MOD;
    return c;
}

mat power(mat a, ll e) {
    ll k = a.size();
    mat r(k, vector<ll>(k, 0));
    for (ll i = 0; i < k; i++) r[i][i] = 1;
    for (; e > 0; e >>= 1) {
        if (e & 1) r = mul(r, a);
        a = mul(a, a);
    }
    return r;
}

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

    ll A, B, n, x;
    cin >> A >> B >> n >> x;
    mat r = power({{A % MOD, B % MOD}, {0, 1}}, n);
    cout << (r[0][0] * (x % MOD) + r[0][1]) % MOD << '\n';
}

Input2 5 4 3

Output123

Khoảng 60 lần nhân ma trận $2 \times 2$.

Đếm đường đi và quy hoạch động bằng ma trận kề

Đồ thị có hướng gồm các đỉnh (vẽ thành chấm tròn, đánh số $1..n$) và các cạnh (mũi tên đi từ một đỉnh tới một đỉnh); chương 3 học kỹ về đồ thị, ở đây chỉ cần chừng đó. Xét đồ thị có hướng với $n$ đỉnh. Ma trận kề $A$ có $A_{uv}$ bằng số cạnh từ $u$ tới $v$. Một đường đi (walk) có $k$ cạnh là dãy đỉnh $u_0, u_1, \dots, u_k$ mà mỗi $u_{j-1} \to u_j$ là một cạnh; đỉnh và cạnh được lặp lại. Điều cần nhớ: $(A^k)_{uv}$ là số đường đi đúng $k$ cạnh từ $u$ tới $v$. Quy nạp theo $k$: phân loại đường đi $k + 1$ cạnh theo đỉnh $x$ ngay trước đỉnh cuối, phần đầu có $(A^k)_{ux}$ cách, rồi cạnh $x \to v$ có $A_{xv}$ cách; cộng theo $x$ ra đúng công thức nhân ma trận.

Trong hình, ô hàng 1, cột 4 của $A^2$ bằng $A_{12} A_{24} + A_{13} A_{34} = 2$, ứng với hai đường đi $1 \to 2 \to 4$ và $1 \to 3 \to 4$.

Đừng nhầm đường đi này với đường đi đơn: nó cho đi lại đỉnh cũ, nên một chu trình chỉ có ba đỉnh vẫn có đường đi dài 100 cạnh. Từ đỉnh $s$ tới $t$: lấy ô $(A^k)_{st}$. Tổng số đường đi $k$ cạnh: cộng mọi ô của $A^k$. Đường đi không quá $k$ cạnh: thêm một đỉnh mới $z$ có khuyên $z \to z$ và cạnh $v \to z$ từ mọi đỉnh $v$; mỗi đường đi $j \le k$ cạnh khớp với đúng một đường đi $k + 1$ cạnh kết thúc ở $z$.

Quy hoạch động theo lớp. Nhiều bài đếm xâu quy về đếm đường đi trong một đồ thị trạng thái nhỏ, trong đó đỉnh là thông tin cần nhớ về phần xâu đã viết và cạnh là ký tự được thêm vào. Đếm xâu nhị phân độ dài $n$ không có ba chữ số 1 liên tiếp chỉ cần nhớ số chữ số 1 ở cuối xâu, tức ba trạng thái $0, 1, 2$: thêm 0 đưa mọi trạng thái về 0, thêm 1 đưa $s$ lên $s + 1$ nếu $s < 2$. Với vectơ cột số xâu ở từng trạng thái,

với cạnh $s \to s'$ nằm ở hàng $s'$, cột $s$. Các vectơ $w_1, \dots, w_4$ có tổng $2, 4, 7, 13$; tám xâu độ dài 3 trừ xâu $111$ đúng bằng 7. Mọi quy hoạch động theo lớp mà bảng chuyển giống nhau ở mọi lớp và chỉ cộng các giá trị của lớp trước đều lũy thừa hóa được như vậy.

Ví dụĐếm Hành Trình Đúng K Cạnh

Đề. Đồ thị có hướng đơn $n \le 100$ đỉnh. Đếm tổng số đường đi (walk) đúng $k \le 10^9$ cạnh với đầu, cuối tùy ý, theo modulo $10^9 + 7$.

Lời giải. Ô $(u, v)$ của $A^k$ là số đường đi $k$ cạnh từ $u$ tới $v$, nên đáp án là tổng mọi ô của $A^k$. Ta không cần cả ma trận: với $\mathbf{1}$ là vectơ toàn số 1, $(A^k \mathbf{1})_u = \sum_v (A^k)_{uv}$ là số đường đi $k$ cạnh bắt đầu tại $u$, và đáp án là tổng các tọa độ. Quy hoạch động theo từng độ dài thì tốn $k \cdot n^2 = 10^{13}$.

Bước 1: ma trận kề. Cạnh $u \to v$ đặt ở hàng $u$, cột $v$ (đổi về đánh số từ 0). Khi đó (A * w)[u] cộng w[v] qua mọi cạnh $u \to v$: số đường đi bắt đầu tại $u$ bằng tổng qua đỉnh kế tiếp.

mat a(n, vector<ll>(n, 0));
while (m--) {
    ll u, v;
    cin >> u >> v;
    a[u - 1][v - 1] = 1;
}

Bước 2: nhân vectơ theo các bit của $k$. w bắt đầu là $\mathbf{1}$; bit 1 thì w = a * w, luôn bình phương a. Thứ tự nhân không quan trọng vì mọi lũy thừa của cùng một ma trận giao hoán với nhau.

vector<ll> w(n, 1);
for (; k > 0; k >>= 1) {
    if (k & 1) w = mulvec(a, w);
    a = mul(a, a);
}

Bước 3: cộng các tọa độ.

ll ans = 0;
for (ll x : w) ans = (ans + x) % MOD;

Với ví dụ, từ mỗi đỉnh chỉ có một đường đi 3 cạnh, tổng 2.

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

const ll MOD = 1e9 + 7;

mat mul(mat a, mat b) {
    ll k = a.size();
    mat c(k, vector<ll>(k, 0));
    for (ll i = 0; i < k; i++)
        for (ll t = 0; t < k; t++) {
            if (!a[i][t]) continue;
            for (ll j = 0; j < k; j++) c[i][j] = (c[i][j] + a[i][t] * b[t][j]) % MOD;
        }
    return c;
}

vector<ll> mulvec(mat a, vector<ll> v) {
    ll k = a.size();
    vector<ll> w(k, 0);
    for (ll i = 0; i < k; i++)
        for (ll j = 0; j < k; j++) w[i] = (w[i] + a[i][j] * v[j]) % MOD;
    return w;
}

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

    ll n, m, k;
    cin >> n >> m >> k;
    mat a(n, vector<ll>(n, 0));
    while (m--) {
        ll u, v;
        cin >> u >> v;
        a[u - 1][v - 1] = 1;
    }
    vector<ll> w(n, 1);
    for (; k > 0; k >>= 1) {
        if (k & 1) w = mulvec(a, w);
        a = mul(a, a);
    }
    ll ans = 0;
    for (ll x : w) ans = (ans + x) % MOD;
    cout << ans << '\n';
}

Input2 2 3
1 2
2 1

Output2

30 lần bình phương ma trận $100 \times 100$, mỗi lần $10^6$ phép.

Đổi phép cộng và phép nhân

Trong chứng minh $(AB)C = A(BC)$ ta chỉ dùng: phép cộng giao hoán, kết hợp và có phần tử không; phép nhân kết hợp, có phần tử đơn vị, và phân phối với phép cộng. Cấu trúc nào có đủ các tính chất ấy (nửa vành, semiring) thì phép nhân ma trận và lũy thừa bằng bình phương dùng nguyên được. Ba biến thể hay gặp:

• $(\min, +)$: thay “cộng” bằng $\min$ và “nhân” bằng $+$, $C_{ij} = \min_t (A_{it} + B_{tj})$. Khi $A_{uv}$ là trọng số cạnh (vô cực nếu không có cạnh), $(A^k)_{uv}$ là độ dài đường đi ngắn nhất đúng $k$ cạnh. Cùng khuôn với $\max$ cho đường đi dài nhất.
• $(\vee, \wedge)$: “có hay không” thay cho “bao nhiêu”, để biết $v$ có tới được từ $u$ bằng đúng $k$ cạnh không.
• $(\oplus, \wedge)$: XOR thay cộng, AND thay nhân, cho các truy hồi trên bit.

Hai phần tử đặc biệt đổi theo cấu trúc. Phần tử không của $\min$ là $+\infty$, phần tử đơn vị của $+$ là 0, nên ma trận đơn vị của $(\min, +)$ có $0$ trên đường chéo và $\infty$ ở nơi còn lại, không phải ma trận 1/0 quen thuộc. Chương trình sau in độ dài đường đi ngắn nhất đúng $k$ cạnh từ đỉnh 1 tới đỉnh $n$ (hoặc $-1$). Bỏ qua các ô $\infty$ trước khi cộng để $\infty + 5$ không bị nhầm thành giá trị hữu hạn.

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

const ll INF = 1e18;

mat mul(mat a, mat b) {
    ll k = a.size();
    mat c(k, vector<ll>(k, INF));
    for (ll i = 0; i < k; i++)
        for (ll t = 0; t < k; t++) {
            if (a[i][t] == INF) continue;
            for (ll j = 0; j < k; j++)
                if (b[t][j] < INF) c[i][j] = min(c[i][j], a[i][t] + b[t][j]);
        }
    return c;
}

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

    ll n, m, k;
    cin >> n >> m >> k;
    mat a(n, vector<ll>(n, INF));
    while (m--) {
        ll u, v, w;
        cin >> u >> v >> w;
        a[u - 1][v - 1] = min(a[u - 1][v - 1], w);
    }
    mat r(n, vector<ll>(n, INF));
    for (ll i = 0; i < n; i++) r[i][i] = 0;
    for (; k > 0; k >>= 1) {
        if (k & 1) r = mul(r, a);
        a = mul(a, a);
    }
    cout << (r[0][n - 1] == INF ? -1 : r[0][n - 1]) << '\n';
}

Input4 6 4
1 2 4
1 3 1
3 2 1
2 4 5
3 4 6
2 3 1

Output9

Đường đi $1 \to 3 \to 2 \to 3 \to 4$ có đúng 4 cạnh và độ dài 9; đường $1 \to 3 \to 4$ chỉ dài 7 nhưng có hai cạnh nên không được tính. Cũng vì phép nhân ma trận không giao hoán, khi nhân nhiều ma trận khác nhau (mỗi phần tử của dãy một ma trận, hỏi tích trên một đoạn) ta phải giữ đúng thứ tự; cây phân đoạn của bài 17 làm được việc này.

Ví dụTruy Hồi Bit

Đề. Cho $K \le 100$ số hạng đầu $A_1..A_K$ và $K$ hệ số $C_1..C_K$, đều là số không dấu 32 bit. Với $i > K$, $A_i = (C_1 \wedge A_{i-1}) \oplus \dots \oplus (C_K \wedge A_{i-K})$. In $A_M$ với $M \le 10^9$.

Lời giải. Đúng truy hồi bậc $K$, chỉ đổi cặp phép toán sang $(\oplus, \wedge)$: “cộng” là XOR, “nhân” là AND. Lũy thừa ma trận chỉ cần phép cộng giao hoán, kết hợp, có phần tử không, và phép nhân kết hợp, phân phối với phép cộng. XOR và AND thỏa hết, theo từng bit: phần tử không của XOR là 0, phần tử đơn vị của AND là số toàn bit 1 (~0u). Kiểu unsigned 32 bit giữ đúng mọi phép bit.

Bước 1: nhân ma trận với cặp phép mới. Như hàm mul quen thuộc, thay + bằng ^ và * bằng &, không cần modulo.

mat mul(mat a, mat b) {
    ll k = a.size();
    mat c(k, vector<unsigned>(k, 0));
    for (ll i = 0; i < k; i++)
        for (ll t = 0; t < k; t++) {
            if (!a[i][t]) continue;
            for (ll j = 0; j < k; j++) c[i][j] ^= a[i][t] & b[t][j];
        }
    return c;
}

Bước 2: ma trận dời cửa sổ. Trạng thái là cửa sổ $(A_K, \dots, A_1)$. Hàng đầu chứa các hệ số $C_1..C_K$; các ô “1” dời cửa sổ xuống phải là phần tử đơn vị của AND, tức ~0u, chứ không phải số 1.

mat t(k, vector<unsigned>(k, 0));
for (ll j = 0; j < k; j++) t[0][j] = c[j];
for (ll i = 1; i < k; i++) t[i][i - 1] = ~0u;
vector<unsigned> w(k);
for (ll i = 0; i < k; i++) w[i] = a[k - 1 - i];

Bước 3: đi $M - K$ bước. Nếu $M \le K$ thì in thẳng $A_M$. Ngược lại nhân vectơ theo các bit của $M - K$, rồi đọc tọa độ trên w[0].

for (ll e = m - k; e > 0; e >>= 1) {
    if (e & 1) w = mulvec(t, w);
    t = mul(t, t);
}

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using mat = vector<vector<unsigned>>;

mat mul(mat a, mat b) {
    ll k = a.size();
    mat c(k, vector<unsigned>(k, 0));
    for (ll i = 0; i < k; i++)
        for (ll t = 0; t < k; t++) {
            if (!a[i][t]) continue;
            for (ll j = 0; j < k; j++) c[i][j] ^= a[i][t] & b[t][j];
        }
    return c;
}

vector<unsigned> mulvec(mat a, vector<unsigned> v) {
    ll k = a.size();
    vector<unsigned> w(k, 0);
    for (ll i = 0; i < k; i++)
        for (ll j = 0; j < k; j++) w[i] ^= a[i][j] & v[j];
    return w;
}

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

    ll k, m;
    cin >> k >> m;
    vector<unsigned> a(k), c(k);
    for (ll i = 0; i < k; i++) cin >> a[i];
    for (ll i = 0; i < k; i++) cin >> c[i];
    if (m <= k) {
        cout << a[m - 1] << '\n';
        return 0;
    }
    mat t(k, vector<unsigned>(k, 0));
    for (ll j = 0; j < k; j++) t[0][j] = c[j];
    for (ll i = 1; i < k; i++) t[i][i - 1] = ~0u;
    vector<unsigned> w(k);
    for (ll i = 0; i < k; i++) w[i] = a[k - 1 - i];
    for (ll e = m - k; e > 0; e >>= 1) {
        if (e & 1) w = mulvec(t, w);
        t = mul(t, t);
    }
    cout << w[0] << '\n';
}