DP nhân ma trận (Matrix Exponentiation)

Lũy thừa ma trận (matrix exponentiation) là kỹ thuật biểu diễn một công thức truy hồi tuyến tính dưới dạng phép nhân ma trận, rồi dùng lũy thừa nhanh (binary exponentiation) để tính số hạng thứ $N$ trong thời gian $O(k^3\log N)$ thay vì $O(N)$ như khi lặp tuần tự. Nhờ đó ta tính được những bài mà $N$ lên tới ${10}^{18}$ — quá lớn để duyệt từng bước, nhưng ${\log }_2({10}^{18})\approx 60$ nên chỉ cần khoảng 60 phép nhân ma trận.

Kỹ thuật này áp dụng cho mọi DP/truy hồi mà trạng thái tiếp theo là tổ hợp tuyến tính của một số cố định trạng thái trước đó: dãy Fibonacci, đếm số đường đi đúng $N$ bước trên đồ thị, đếm cách lát/xếp với chiều dài khổng lồ, v.v.

Ý tưởng / Trực giác

Giả sử ta gom $k$ giá trị gần nhất của dãy vào một vector trạng thái. Với Fibonacci ($F_n=F_{n-1}+F_{n-2}$), chọn vector 2 chiều $\left(\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right)$.

Câu hỏi mấu chốt: có tồn tại một ma trận cố định $M$ biến vector trạng thái tại bước $n$ thành vector tại bước $n+1$ không? Nếu truy hồi là tuyến tính thì có:

$$\left(\begin{array}{c}{F}_{n+1}\\ F_n\end{array}\right)=\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right)$$

Hàng đầu $\left(\begin{array}{cc}1& 1\end{array}\right)$ nói "$F_{n+1}=1·F_n+1·F_{n-1}$" — đúng công thức truy hồi. Hàng hai $\left(\begin{array}{cc}1& 0\end{array}\right)$ chỉ là "$F_n=F_n$" để dịch cửa sổ trượt xuống một bậc.

Vì $M$ không đổi giữa các bước, áp $M$ liên tiếp $n-1$ lần đưa vector ban đầu tới vector thứ $n$:

$$\left(\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right)=M^{n-1}\left(\begin{array}{c}{F}_1\\ F_0\end{array}\right)$$

Vì sao đúng? Phép nhân ma trận có tính kết hợp: $M·(M\cdots (Mv))=M^{n-1}v$. Mà lũy thừa $M^{n-1}$ tính được bằng kỹ thuật bình phương liên tiếp, giống lũy thừa số nguyên: $M^{13}=M^8·M^4·M^1$ (13 = 1101 nhị phân). Mỗi lần ta hoặc bình phương $M$ (tăng số mũ lên gấp đôi), hoặc nhân kết quả với $M$ hiện tại — tổng cộng chỉ $O(\log N)$ phép nhân ma trận. Đây chính là chỗ đổi $O(N)$ thành $O(\log N)$.

Trực giác tổng quát: mỗi phép nhân ma trận = "đi một bước truy hồi"; lũy thừa ma trận = "đi $N$ bước một lúc" bằng cách gộp các bước theo lũy thừa của 2.

Ví dụ chạy tay

Tính $F_5$ với $F_0=0,F_1=1$ (đáp số đúng là $5$). Dùng $M=\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)$, cần $M^4$ vì $\left(\begin{array}{c}{F}_5\\ F_4\end{array}\right)=M^4\left(\begin{array}{c}{F}_1\\ F_0\end{array}\right)$.

Số mũ $4={100}_2$. Lũy thừa nhanh duyệt bit từ thấp lên cao, giữ result (khởi tạo = ma trận đơn vị $I$) và base (khởi tạo = $M$):

p = 4 (binary 100)        result = [1 0]   base = [1 1]
                                   [0 1]          [1 0]

--- bit 0 của p = 0 -> KHÔNG nhân result ---
   base = base*base = M^2:
        [1 1][1 1]   [2 1]
        [1 0][1 0] = [1 1]
   p >>= 1  ->  p = 2 (binary 10)

   result = [1 0]   base = [2 1]
            [0 1]          [1 1]

--- bit 0 của p = 0 -> KHÔNG nhân result ---
   base = base*base = M^4:
        [2 1][2 1]   [5 3]
        [1 1][1 1] = [3 2]
   p >>= 1  ->  p = 1 (binary 1)

   result = [1 0]   base = [5 3]
            [0 1]          [3 2]

--- bit 0 của p = 1 -> NHÂN result ---  (ô đang ghi: result)
   result = result*base = I*M^4 = [5 3]
                                  [3 2]
   p >>= 1  ->  p = 0  -> DỪNG

Vậy $M^4=\left(\begin{array}{cc}5& 3\\ 3& 2\end{array}\right)$. Nhân với vector đầu:

$$\left(\begin{array}{c}{F}_5\\ F_4\end{array}\right)=\left(\begin{array}{cc}5& 3\\ 3& 2\end{array}\right)\left(\begin{array}{c}{F}_1\\ F_0\end{array}\right)=\left(\begin{array}{cc}5& 3\\ 3& 2\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)=\left(\begin{array}{c}5\\ 3\end{array}\right)$$

Đọc ra $F_5=5$, $F_4=3$ — khớp. Chú ý phần tử $M^n[0][0]$ chính bằng $F_{n+1}$: một mẹo hay để lấy kết quả mà không cần nhân vector.

Cài đặt

Bộ khung dùng chung cho mọi bài: kiểu ma trận, phép nhân (có lấy modulo), và lũy thừa nhanh.

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

const long long MOD = 1e9 + 7;
typedef vector<vector<long long>> Matrix;

int sz; // kích thước ma trận (k x k)

// Nhân hai ma trận k x k, lấy modulo MOD
Matrix multiply(const Matrix& A, const Matrix& B) {
    Matrix C(sz, vector<long long>(sz, 0));
    for (int i = 0; i < sz; i++)
        for (int k = 0; k < sz; k++) {
            if (A[i][k] == 0) continue;        // bỏ qua phần tử 0 để tăng tốc
            for (int j = 0; j < sz; j++)
                // A[i][k], B[k][j] < MOD ~ 1e9 nên tích < 1e18, vừa khít long long
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
        }
    return C;
}

// Tính A^p bằng bình phương liên tiếp
Matrix matpow(Matrix A, long long p) {
    Matrix result(sz, vector<long long>(sz, 0));
    for (int i = 0; i < sz; i++) result[i][i] = 1; // ma trận đơn vị I
    while (p > 0) {
        if (p & 1) result = multiply(result, A);   // bit hiện tại = 1: gộp vào kết quả
        A = multiply(A, A);                          // bình phương: số mũ gấp đôi
        p >>= 1;
    }
    return result;
}

Ví dụ 1 — Fibonacci thứ $n$ ($n\le {10}^{18}$), với $F_0=0,F_1=1$:

long long fibonacci(long long n) {
    if (n == 0) return 0;
    sz = 2;
    Matrix M = {{1, 1}, {1, 0}};
    Matrix R = matpow(M, n);   // M^n[0][0] = F_{n+1}, M^n[0][1] = F_n
    return R[0][1];            // chính là F_n
}

Ví dụ 2 — Truy hồi tuyến tính bậc $k$ tổng quát $f(n)=c_{1}f(n-1)+\dots +c_{k}f(n-k)$:

// init = {f(0), f(1), ..., f(k-1)}; c = {c1, c2, ..., ck}
long long solve(long long n, vector<long long> c, vector<long long> init) {
    int k = c.size();
    sz = k;
    if (n < k) return init[n];
    // Hàng đầu là các hệ số; các hàng dưới dịch cửa sổ xuống một bậc
    Matrix M(k, vector<long long>(k, 0));
    for (int j = 0; j < k; j++) M[0][j] = c[j] % MOD;
    for (int i = 1; i < k; i++) M[i][i-1] = 1;

    // Trạng thái ban đầu V = [f(k-1), f(k-2), ..., f(0)]^T (giảm dần)
    Matrix R = matpow(M, n - (k - 1));
    long long ans = 0;
    for (int j = 0; j < k; j++)
        ans = (ans + R[0][j] * init[k - 1 - j]) % MOD;
    return ans;
}

Ví dụ 3 — Đếm số đường đi đúng $N$ bước trên đồ thị: nếu $A$ là ma trận kề ($A[i][j]$ = số cạnh từ $i$ tới $j$) thì $A^N[s][t]$ = số đường đi (theo cạnh, có thể lặp đỉnh) độ dài đúng $N$ từ $s$ tới $t$. Chỉ cần lập $A$ rồi matpow(A, N).

Độ phức tạp

Ước lượng thực tế: với $k\le 100$ và $N\le {10}^{18}$, ta có $k^3{\log }_{2}N\approx {10}^6\times 60\approx 6\times {10}^7$ — chạy trong khoảng dưới một giây. Khi $k$ lớn hơn (vd $k=200$) thì $k^3=8\times {10}^6$ và tổng lên $~5\times {10}^8$, cần tối ưu hằng số (bỏ phần tử 0, dùng mảng phẳng thay vector<vector>).

⚠️ Lỗi thường gặp

• Tràn số khi nhân. Tích $A[i][k]\times B[k][j]$ với hai thừa số cỡ ${10}^9$ cho ra $~{10}^{18}$, sát trần long long ($\approx 9.2\times {10}^{18}$). Nếu cộng dồn nhiều tích trước khi lấy modulo, hoặc dùng int, sẽ tràn. Cách tránh: ép kiểu long long, lấy % MOD ngay sau mỗi lần cộng một tích.
• Khởi tạo sai phần tử trung tính. matpow phải bắt đầu result = ma trận đơn vị $I$ (đường chéo 1, còn lại 0), không phải ma trận 0. Lấy ma trận 0 thì mọi tích về 0 và kết quả luôn sai. (Trong biến thể min-plus, phần tử trung tính lại là ma trận có 0 trên đường chéo, $+\infty$ ngoài đường chéo — đừng nhầm.)
• Lệch số mũ (off-by-one). Phải xác định rõ vector ban đầu ứng với chỉ số nào. Với $\left(\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right)=M^{n-1}\left(\begin{array}{c}{F}_1\\ F_0\end{array}\right)$ thì số mũ là $n-1$, không phải $n$. Sai một bậc cho ra $F_{n\pm 1}$. Luôn kiểm bằng một giá trị nhỏ tự tính tay (như $F_5=5$ ở trên).
• Quên xử lý $N$ nhỏ. Khi $N<k$ (vd hỏi $F_0$ hay $f(1)$ với truy hồi bậc 3), số mũ $n-(k-1)$ có thể âm, làm matpow chạy sai. Phải trả thẳng init[n] cho các trường hợp cơ sở trước khi gọi lũy thừa.
• Modulo âm. Khi truy hồi có hệ số trừ (vd $f(n)=f(n-1)-f(n-2)$), kết quả % MOD có thể âm. Sửa bằng ((x % MOD) + MOD) % MOD.
• Tưởng nhanh nhưng vẫn $O(N)$ ẩn. Nếu lỡ tính $M^N$ bằng vòng lặp nhân $N$ lần thay vì bình phương liên tiếp, độ phức tạp tụt về $O(k^{3}N)$ — vô dụng với $N={10}^{18}$. Đảm bảo p >>= 1 chứ không phải p--.

Biến thể / Mở rộng

• Thêm hằng số / đa thức vào truy hồi. Với $f(n)=f(n-1)+f(n-2)+n+c$, đưa thêm các biến phụ ($n$ và hằng $1$) vào vector trạng thái để giữ tính tuyến tính. Ví dụ vector $[f(n),f(n-1),n,1{]}^T$ và ma trận tương ứng tự cập nhật $n\to n+1$.
• Min-plus (tropical) matrix exponentiation. Thay phép $(\times ,+)$ trong nhân ma trận bằng $(+,min)$: $C[i][j]={min}_k(A[i][k]+B[k][j])$. Khi đó $A^k[s][t]$ cho đường đi ngắn nhất dùng đúng $k$ cạnh (phần tử trung tính là $+\infty$ ngoài đường chéo, $0$ trên đường chéo). Đây là kỹ thuật cho bài graphpath2.
• Lũy thừa ma trận trên số mũ. Khi dãy nhân (vd $f_x=c^{2x-6}{f}_{x-1}{f}_{x-2}{f}_{x-3}$), lấy logarit/quan tâm số mũ của từng thừa số rồi dùng matrix exp trên các số mũ (theo modulo $p-1$ nhờ định lý Fermat nhỏ). Đây là bài prodrec.
• Liên hệ. Kỹ thuật bình phương liên tiếp giống hệt lũy thừa nhanh số nguyên; ý tưởng "gộp nhiều bước truy hồi" cũng xuất hiện trong nhiều bài quy hoạch động.

Bài tập luyện

• Số Fibonacci (fibonum) — (Intermediate) Tính $F_n$ với $n\le {10}^{18}$: bài mẫu kinh điển để dựng ma trận $2\times 2$ và lũy thừa nhanh.
• Đường đi trong đồ thị I (graphpath1) — (Advanced) Đếm số đường đi đúng $k$ cạnh từ đỉnh $1$ tới $n$: chính là $A^k[1][n]$ với $A$ là ma trận kề.
• Đường đi trong đồ thị II (graphpath2) — (Advanced) Đường đi ngắn nhất đúng $k$ cạnh: luyện biến thể min-plus (tropical) của lũy thừa ma trận.
• Truy hồi tích luỹ (prodrec) — (Expert) Dãy nhân $f_x=c^{2x-6}{f}_{x-1}{f}_{x-2}{f}_{x-3}$: luyện lũy thừa ma trận trên số mũ (modulo $p-1$).
• Ngôn ngữ COWBASIC (cowbasic) — (Expert) Mô phỏng vòng lặp lồng nhau tới ${10}^5$ lần: biểu diễn mỗi vòng lặp thành ma trận biến đổi tuyến tính rồi lũy thừa.