Nhân ma trận (Matrix Exponentiation)

Nhân ma trận nhanh (Matrix Exponentiation) là kỹ thuật tính lũy thừa $M^n$ của một ma trận vuông $K\times K$ trong $O(K^3\log N)$ thay vì $O(N)$ phép cộng/nhân ngây thơ. Nhờ đó nó giải được công thức truy hồi tuyến tính (như Fibonacci) với $n$ lên tới ${10}^{18}$, và đếm/đo đường đi có đúng $K$ cạnh trong đồ thị — những bài mà mô phỏng từng bước là bất khả thi.

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

Mọi công thức truy hồi tuyến tính đều có thể viết lại thành một phép nhân ma trận với vector trạng thái. Lấy Fibonacci $F_n=F_{n-1}+F_{n-2}$ làm ví dụ. Gọi trạng thái là vector cột chứa hai số liên tiếp:

$$v_n=\left(\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right).$$

Ta muốn tìm một ma trận $M$ sao cho $v_n=M·v_{n-1}$. Vì $F_n=1·F_{n-1}+1·F_{n-2}$ và $F_{n-1}=1·F_{n-1}+0·F_{n-2}$:

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

Đây là chìa khoá. Áp dụng $M$ một lần là "đi một bước". Áp dụng nó $n$ lần liên tiếp:

$$v_n=M·v_{n-1}=M^2·v_{n-2}=\cdots =M^{n-1}·v_1.$$

Câu hỏi còn lại: làm sao tính $M^{n-1}$ nhanh? Nếu nhân $M$ với chính nó $n-1$ lần thì vẫn $O(n)$, chẳng lợi gì. Nhưng phép nhân ma trận có tính kết hợp, nên ta dùng đúng ý tưởng của lũy thừa nhị phân (binary exponentiation): phân tích số mũ theo cơ số 2. Ví dụ $M^{13}=M^8·M^4·M^1$ (vì $13={1101}_2$). Ta chỉ cần bình phương liên tiếp $M,M^2,M^4,M^8,\dots$ — tốn $\log n$ phép nhân ma trận — rồi nhân vào kết quả những lũy thừa ứng với bit bằng 1.

Vì sao đúng? Tính kết hợp của nhân ma trận đảm bảo $M^a·M^b=M^{a+b}$, nên gộp các "khối" lũy thừa 2 lại đúng bằng $M^n$. Một phép nhân hai ma trận $K\times K$ tốn $O(K^3)$, ta làm $O(\log n)$ lần, nên tổng cộng $O(K^3\log n)$.

Tổng quát: bất kỳ đại lượng nào ở bước $n$ được tính từ một tổ hợp tuyến tính cố định của các đại lượng ở các bước trước đều đặt được vào khuôn này. Với đồ thị, $(A^K)[i][j]$ — số đường đi $K$ cạnh từ $i$ tới $j$ — chính là một dạng "truy hồi" trên ma trận kề, vì đường đi $K$ cạnh = đi 1 cạnh rồi đi $K-1$ cạnh.

Ví dụ chạy tay

Tính $F_5$ qua $M^4·v_1$ với $v_1=\left(\begin{array}{c}{F}_1\\ F_0\end{array}\right)=\left(\begin{array}{c}1\\ 0\end{array}\right)$ và $M=\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)$.

Số mũ là $4={100}_2$. Lũy thừa nhị phân duyệt các bit của 4 từ thấp lên cao, vừa bình phương A, vừa nhân vào result khi bit = 1:

 n=4 (100b)   bit  result (ban đầu = I)        A (ban đầu = M)
 -----------  ---  --------------------------  ----------------
 bước 1       0    I (không nhân)              A = M^2 = |2 1|
                                                          |1 1|
 bước 2       0    I (không nhân)              A = M^4 = |5 3|
              ^                                          |3 2|
            bit này
 bước 3       1    result = I * M^4 = M^4       (dừng)
              ^
            bit cao nhất

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

$$M^4{v}_1=\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)=\left(\begin{array}{c}{F}_5\\ F_4\end{array}\right).$$

Kết quả $F_5=5$, đúng (dãy $0,1,1,2,3,5$). Chỉ tốn 2 phép bình phương ma trận thay vì 4 bước cộng — với $n={10}^{18}$ thì khác biệt là sống còn.

Kiểm tra nhanh một phép bình phương để thấy cơ chế:

$$M^2=\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)=\left(\begin{array}{cc}1·1+1·1& 1·1+1·0\\ 1·1+0·1& 1·1+0·0\end{array}\right)=\left(\begin{array}{cc}2& 1\\ 1& 1\end{array}\right).$$

Cài đặt

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

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

// Nhân hai ma trận A (n x p) và B (p x m), kết quả lấy theo MOD.
Matrix multiply(const Matrix& A, const Matrix& B) {
    int n = A.size(), p = B.size(), m = B[0].size();
    Matrix C(n, vector<long long>(m, 0));
    for (int i = 0; i < n; i++)
        for (int k = 0; k < p; k++)
            if (A[i][k])                       // bỏ qua phần tử 0 để tăng tốc
                for (int j = 0; j < m; j++)
                    // A[i][k] đã < MOD nên A[i][k]*B[k][j] < 1e18, an toàn với long long
                    C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
    return C;
}

// Tính A^n bằng lũy thừa nhị phân: O(K^3 log n).
Matrix mat_pow(Matrix A, long long n) {
    int sz = A.size();
    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 (n > 0) {
        if (n & 1) result = multiply(result, A);    // bit thấp nhất = 1 -> nhân vào
        A = multiply(A, A);                          // bình phương để "lên một bậc lũy thừa 2"
        n >>= 1;
    }
    return result;
}

// Fibonacci thứ n (F_0 = 0, F_1 = 1) theo MOD.
long long fibonacci(long long n) {
    if (n == 0) return 0;
    Matrix M = {{1, 1}, {1, 0}};
    Matrix R = mat_pow(M, n);     // R = M^n; khi đó R[0][1] = F_n
    return R[0][1];
}

int main() {
    long long n; 
    cin >> n;
    cout << fibonacci(n) << "\n";
}

Với truy hồi tuyến tính bậc $k$ tổng quát $a_n=c_1{a}_{n-1}+\cdots +c_k{a}_{n-k}$, ma trận chuyển là ma trận đồng hành (companion matrix):

// init = {a_0, a_1, ..., a_{k-1}};  c = {c_1, ..., c_k}
long long linear_recurrence(vector<long long> c, vector<long long> init, long long n) {
    int k = c.size();
    if (n < k) return init[n] % MOD;
    Matrix M(k, vector<long long>(k, 0));
    for (int j = 0; j < k; j++) M[0][j] = c[j];          // hàng đầu là các hệ số
    for (int i = 1; i < k; i++) M[i][i-1] = 1;           // dịch trạng thái xuống dưới
    Matrix R = mat_pow(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)) % MOD;
    return ans % MOD;
}

Với đếm đường đi đúng $K$ cạnh trong đồ thị $N$ đỉnh: lấy ma trận kề $A$ (với $A[i][j]=1$ nếu có cạnh $i\to j$), thì $(A^K)[i][j]$ là số đường đi $K$ cạnh từ $i$ tới $j$.

Matrix A_k = mat_pow(adj, K);   // adj là ma trận kề N x N
long long paths = A_k[s][t];    // số đường đi đúng K cạnh từ s tới t

Độ phức tạp

Vì sao $O(K^3\log n)$? Một phép nhân hai ma trận $K\times K$ có ba vòng lặp lồng nhau nên tốn $O(K^3)$. Lũy thừa nhị phân thực hiện $O(\log n)$ phép nhân (mỗi bit của $n$ tốn nhiều nhất 2 phép: một bình phương, một nhân vào kết quả). Nhân hai con số đó ra $O(K^3\log n)$.

Bộ nhớ $O(K^2)$: ta chỉ lưu vài ma trận $K\times K$ (kết quả, ma trận đang bình phương, ma trận tạm trong multiply). Không cần lưu lịch sử các bước, vì lũy thừa nhị phân chỉ giữ trạng thái hiện tại.

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

• Tràn số khi nhân (overflow). $A[i][k]·B[k][j]$ với hai số cỡ ${10}^9$ cho tích cỡ ${10}^{18}$ — vượt int (giới hạn ~$2,1\times {10}^9$) và sát trần long long ($\approx 9,2\times {10}^{18}$). Dùng long long cho tích, và luôn lấy modulo trước để mọi phần tử ma trận đều $<MOD$, đảm bảo tích từng cặp $<{10}^{18}$. Cộng dồn nhiều tích vẫn an toàn vì ta % MOD ngay sau mỗi lần cộng.

• Quên khởi tạo ma trận đơn vị $I$. Nếu result khởi tạo bằng ma trận 0 (hoặc bằng $A$), kết quả sẽ sai hoàn toàn. Phần tử trung hoà của phép nhân ma trận là $I$ (đường chéo 1, còn lại 0), giống vai trò của số 1 trong lũy thừa số học.

• Sai số mũ (off-by-one). Với Fibonacci $v_n=M^{n-1}{v}_1$ chứ không phải $M^n{v}_1$; tuỳ cách định nghĩa trạng thái mà số mũ lệch 1. Hãy kiểm tra bằng tay với $n$ nhỏ ($F_5=5$) trước khi tin vào công thức.

• Quên xử lý biên $n$ nhỏ. Với truy hồi bậc $k$, khi $n<k$ thì đáp án nằm thẳng trong giá trị khởi tạo init[n], không được đưa vào mat_pow (số mũ âm sẽ sai hoặc vòng while bỏ qua, trả về $I$ cho kết quả lệch). Kiểm tra if (n < k) return init[n]; ngay đầu hàm.

• Modulo âm khi bài toán có phép trừ. Nếu công thức truy hồi có hệ số âm (vd $a_n=a_{n-1}-a_{n-2}$), kết quả % MOD có thể âm trong C++. Sửa bằng ((x % MOD) + MOD) % MOD. Tương tự khi đưa hệ số âm vào ma trận, hãy chuẩn hoá về $[0,MOD)$ trước.

• Modulo trên bài "đường đi ngắn nhất K cạnh". Bài tối ưu (min/max) dùng đại số min-plus: phép "nhân" thay bằng ${min}_k(A[i][k]+B[k][j])$ và "cộng" thay bằng $min$. Đừng lấy modulo và đừng dùng phép cộng/nhân thường — phần tử trung hoà lúc này là ma trận toàn $+\infty$ với đường chéo 0, không phải $I$ số học.

Biến thể / Mở rộng

• Đại số min-plus / max-plus (tropical semiring): thay $(+,\times )$ bằng $(min,+)$ để tìm đường đi ngắn nhất có đúng $K$ cạnh, hoặc bằng $(max,+)$ cho bài tối đa. Khung lũy thừa nhị phân giữ nguyên, chỉ đổi định nghĩa multiply và ma trận đơn vị.
• Truy hồi có hằng số (affine): với $a_n=c_1{a}_{n-1}+\cdots +d$ (có số hạng cộng cố định $d$), thêm một thành phần hằng $1$ vào vector trạng thái và một hàng/cột tương ứng trong ma trận để "mang" hằng số đó. Đây là mấu chốt để biến vòng lặp gán tuyến tính thành lũy thừa ma trận.
• Liên hệ: kỹ thuật này dựa trực tiếp trên lũy thừa nhị phân áp dụng lên ma trận thay vì số.

Bài tập luyện

• Số Fibonacci (fibonum) — (Intermediate) tính $F_n$ với $n\le {10}^{18}$, bài nhập môn kinh điển: dựng ma trận $2\times 2$ rồi 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) tìm đường đi ngắn nhất đúng $k$ cạnh, luyện biến thể min-plus: thay phép nhân ma trận bằng $min(A[i][k]+B[k][j])$.
• Truy hồi tích luỹ (prodrec) — (Expert) truy hồi nhân $f_x=c^{2x-6}{f}_{x-1}{f}_{x-2}{f}_{x-3}$: lấy log/đếm số mũ để đưa về truy hồi tuyến tính trên mũ rồi lũy thừa ma trận (kết hợp với định lý Fermat nhỏ).
• Ngôn Ngữ COWBASIC (cowbasic) — (Expert) mỗi vòng lặp MOO lặp tới ${10}^5$ lần là một phép gán affine: gộp thân vòng thành một ma trận chuyển trạng thái rồi lũy thừa để mô phỏng cả vòng trong $O(\log n)$.