矩阵快速幂与线性递推

前言：递推太慢了，怎么办？

斐波那契数列你一定见过：$1, 1, 2, 3, 5, 8, 13, \ldots$，每一项是前两项之和。如果要算第 100 项，循环 99 次就行——很快。

但如果要算第 $10^{18}$ 项呢？即使每秒算 $10^9$ 次，也需要 30 年。

等等——快速幂不是能把 $n$ 次操作压缩到 $\log n$ 次吗？能不能用类似的思路加速递推？

可以，但需要一个中间步骤。快速幂能加速的是 $a^n$ 这种同一个东西重复相乘的操作。而斐波那契递推是 $f_n = f_{n-1} + f_{n-2}$，每一步是加法，不是乘法。

解决方法是把递推关系写成矩阵乘法。矩阵乘法本质上就是一组线性运算的打包——加权和。一旦递推变成了”矩阵的幂”，快速幂就能直接用了。

这篇笔记不需要你事先学过线性代数。我会从矩阵乘法的定义开始讲（只需要中学数学就能理解），然后展示怎么把递推”翻译”成矩阵形式，最后用快速幂加速。

一、线性递推的问题

斐波那契数列：$f_0 = 0$，$f_1 = 1$，$f_n = f_{n-1} + f_{n-2}$。

求第 $n$ 项，朴素递推 $O(n)$。$n = 10^{18}$ 呢？不行。

关键洞察：递推关系是线性的，每一步只做加权和。线性变换可以用矩阵乘法表示，而矩阵的幂可以用快速幂加速。

二、矩阵乘法基础

在构造转移矩阵之前，我们需要先了解矩阵乘法。如果你没学过线性代数，不用怕——算法中用到的矩阵乘法只需要一个规则：“行乘列”。

2.1 定义

$A$ 是 $m \times p$ 矩阵，$B$ 是 $p \times n$ 矩阵，乘积 $C = AB$ 是 $m \times n$ 矩阵：

$$c_{ij} = \sum_{k=1}^{p} a_{ik} \cdot b_{kj}$$

2.2 性质

矩阵乘法最关键的性质是结合律——$(AB)C = A(BC)$。这意味着连续乘很多次可以任意加括号，这正是快速幂能用在矩阵上的原因。

其他性质：

• 结合律：$(AB)C = A(BC)$ — 这是一切的基础，使得快速幂适用
• 不满足交换律：$AB \ne BA$（一般情况）
• 单位矩阵 $I$：$AI = IA = A$

typedef vector<vector<long long>> Mat;

Mat mat_mul(const Mat& A, const Mat& B, long long mod) {
    int m = A.size(), p = A[0].size(), n = B[0].size();
    Mat C(m, vector<long long>(n, 0));
    for (int i = 0; i < m; i++)
        for (int k = 0; k < p; k++)
            for (int j = 0; j < n; j++)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % mod;
    return C;
}

Mat mat_pow(Mat A, long long n, long long mod) {
    int sz = A.size();
    Mat res(sz, vector<long long>(sz, 0));
    for (int i = 0; i < sz; i++) res[i][i] = 1;  // 单位矩阵
    while (n) {
        if (n & 1) res = mat_mul(res, A, mod);
        A = mat_mul(A, A, mod);
        n >>= 1;
    }
    return res;
}

注意三重循环的顺序 i-k-j：把 k 放中间，内层 j 循环连续访问内存，比 i-j-k 快得多。

三、构造转移矩阵

3.1 斐波那契数列

递推：$f_n = f_{n-1} + f_{n-2}$。

把状态写成向量 $\begin{pmatrix} f_n \\ f_{n-1} \end{pmatrix}$，目标是从 $\begin{pmatrix} f_{n-1} \\ f_{n-2} \end{pmatrix}$ 转移过来：

$$\begin{pmatrix} f_n \\ f_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} f_{n-1} \\ f_{n-2} \end{pmatrix}$$

验证：第一行 $1 \cdot f_{n-1} + 1 \cdot f_{n-2} = f_n$ ✓，第二行 $1 \cdot f_{n-1} + 0 \cdot f_{n-2} = f_{n-1}$ ✓。

所以：

$$\begin{pmatrix} f_n \\ f_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{n-1} \begin{pmatrix} f_1 \\ f_0 \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{n-1} \begin{pmatrix} 1 \\ 0 \end{pmatrix}$$

矩阵的 $n-1$ 次幂用快速幂算，$O(\log n)$ 次矩阵乘法。矩阵大小固定为 $2 \times 2$，每次乘法 $O(1)$，总计 $O(\log n)$。

3.2 一般的 $k$ 阶线性递推

斐波那契是 2 阶递推（每一项依赖前 2 项）。一般地，如果每一项依赖前 $k$ 项：

$$f_n = c_1 f_{n-1} + c_2 f_{n-2} + \cdots + c_k f_{n-k}$$

状态向量 $\begin{pmatrix} f_n \\ f_{n-1} \\ \vdots \\ f_{n-k+1} \end{pmatrix}$，转移矩阵：

$$T = \begin{pmatrix} c_1 & c_2 & \cdots & c_{k-1} & c_k \\ 1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & \cdots & 0 & 0 \\ \vdots & & \ddots & & \vdots \\ 0 & 0 & \cdots & 1 & 0 \end{pmatrix}$$

第一行是递推系数，下面是单位阵的”平移”——把 $f_{n-i}$ 移到下一个位置。矩阵大小 $k \times k$，每次乘法 $O(k^3)$，总时间 $O(k^3 \log n)$。

3.3 带常数项或额外项

递推 $f_n = f_{n-1} + f_{n-2} + 3$（带常数）？在状态向量末尾加一个恒为 1 的分量：

$$\begin{pmatrix} f_n \\ f_{n-1} \\ 1 \end{pmatrix} = \begin{pmatrix} 1 & 1 & 3 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} f_{n-1} \\ f_{n-2} \\ 1 \end{pmatrix}$$

恒为 1 的分量通过矩阵最后一行的 $(0, 0, 1)$ 保持不变，常数 3 被加到递推中。

四、加速原理的本质

为什么矩阵快速幂能加速？因为结合律：

$$T^n = \begin{cases} (T^{n/2})^2, & n \text{ 偶} \\ T \cdot T^{n-1}, & n \text{ 奇} \end{cases}$$

和快速幂一模一样——只不过运算对象从整数变成了矩阵。结合律保证 $T \cdot T \cdot T \cdots$（$n$ 个）可以任意加括号，对半折叠。

五、实战：斐波那契第 $n$ 项模 $p$

求 $f_n \bmod p$，$0 \le n \le 10^{18}$，$p = 10^9 + 7$。

#include <cstdio>
#include <vector>
using namespace std;

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

Mat mat_mul(const Mat& A, const Mat& B) {
    int m = A.size(), n = B[0].size(), p = A[0].size();
    Mat C(m, vector<long long>(n, 0));
    for (int i = 0; i < m; i++)
        for (int k = 0; k < p; k++)
            for (int j = 0; j < n; j++)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
    return C;
}

Mat mat_pow(Mat A, long long n) {
    Mat res(2, vector<long long>(2, 0));
    res[0][0] = res[1][1] = 1;
    while (n) {
        if (n & 1) res = mat_mul(res, A);
        A = mat_mul(A, A);
        n >>= 1;
    }
    return res;
}

int main() {
    long long n;
    scanf("%lld", &n);

    if (n == 0) { printf("0\n"); return 0; }

    Mat T = {{1, 1}, {1, 0}};
    Mat R = mat_pow(T, n - 1);

    // R * (f_1, f_0)^T，f_1 = 1, f_0 = 0，所以结果是 R[0][0]
    printf("%lld\n", R[0][0]);
    return 0;
}

解析：$2 \times 2$ 矩阵的快速幂，$\log_2(10^{18}) \approx 60$ 次迭代，每次 8 次乘法，总共约 480 次乘法——对比朴素递推的 $10^{18}$ 次，加速了 $2 \times 10^{15}$ 倍。这就是线性代数的力量：把递推的”时序依赖”压缩为矩阵的”空间组合”，然后利用结合律折叠。