算法刷题系列（五十六）：矩阵快速幂——加速线性递推的多步转移

算法刷题系列（五十六）：矩阵快速幂——加速线性递推的多步转移

线性递推每一步只依赖固定数量的前项。逐项计算第 n 项需要 O(n) 步；如果递推能写成向量乘矩阵，那么连续做 n 次相同转移，就是计算转移矩阵的 n 次幂。利用二进制快速幂，可以把转移次数降到 O(log n)。

本文以斐波那契数列为例：

F(0) = 0, F(1) = 1
F(n) = F(n - 1) + F(n - 2)

一、把递推写成矩阵转移

定义状态向量：

S(n) = [F(n + 1), F(n)]ᵀ

那么：

S(n) = [[1, 1],
        [1, 0]] × S(n - 1)

记转移矩阵为 M，初始 S(0) = [1, 0]ᵀ，于是 S(n) = M^n S(0)。计算 M^n 就能得到目标项。

二、快速幂的二进制分解

指数可写为二进制位的组合。重复执行：

1. 若当前指数为奇数，把当前底数乘进答案；
2. 底数平方；
3. 指数整除 2。

矩阵乘法满足结合律，所以这些平方与累乘和整数快速幂同理。对于指数 13 = 1101₂，需要组合 M^8 × M^4 × M，而无需连续乘 13 次。

三、例子推演

M² = [[2, 1], [1, 1]]，由此 M² S(0) = [2, 1]ᵀ，对应 F(2)=1、F(3)=2。更一般地，矩阵幂会把初始状态映射到所需的后续状态。

本文实现接受非负 n，返回精确整数。Java long 的结果有范围限制；超出 long 表示范围时会发生溢出，若需大项应使用 BigInteger。Python 整数则支持任意精度，但位数本身会随结果增长。

四、Java 实现

public class MatrixFastPower {
    private static long[][] multiply(long[][] a, long[][] b) {
        int n = a.length;
        long[][] result = new long[n][n];
        for (int i = 0; i < n; i++) {
            for (int k = 0; k < n; k++) {
                for (int j = 0; j < n; j++) {
                    result[i][j] += a[i][k] * b[k][j];
                }
            }
        }
        return result;
    }

    private static long[][] power(long[][] base, long exponent) {
        int n = base.length;
        long[][] result = new long[n][n];
        for (int i = 0; i < n; i++) {
            if (base[i].length != n) {
                throw new IllegalArgumentException("matrix must be square");
            }
            result[i][i] = 1;
        }

        long[][] factor = base;
        long remaining = exponent;
        while (remaining > 0) {
            if ((remaining & 1L) != 0) {
                result = multiply(result, factor);
            }
            remaining >>= 1;
            if (remaining > 0) {
                factor = multiply(factor, factor);
            }
        }
        return result;
    }

    public static long fibonacci(long n) {
        if (n < 0) {
            throw new IllegalArgumentException("n must be non-negative");
        }
        long[][] transition = {{1, 1}, {1, 0}};
        long[][] result = power(transition, n);
        return result[1][0];
    }

    public static void main(String[] args) {
        if (fibonacci(0) != 0 || fibonacci(1) != 1
                || fibonacci(10) != 55 || fibonacci(20) != 6765) {
            throw new AssertionError("matrix power test failed");
        }
    }
}

对于本例的 2×2 矩阵，每次矩阵乘法是常数级；更一般的 d×d 矩阵使用朴素乘法时，单次代价为 O(d³)。

五、Python 实现

def multiply(a, b):
    size = len(a)
    result = [[0] * size for _ in range(size)]
    for row in range(size):
        for middle in range(size):
            for column in range(size):
                result[row][column] += a[row][middle] * b[middle][column]
    return result

def matrix_power(base, exponent):
    if exponent < 0:
        raise ValueError("exponent must be non-negative")
    size = len(base)
    if any(len(row) != size for row in base):
        raise ValueError("matrix must be square")

    result = [[int(row == column) for column in range(size)]
              for row in range(size)]
    factor = [row[:] for row in base]
    remaining = exponent
    while remaining > 0:
        if remaining & 1:
            result = multiply(result, factor)
        remaining >>= 1
        if remaining:
            factor = multiply(factor, factor)
    return result

def fibonacci(n):
    if n < 0:
        raise ValueError("n must be non-negative")
    transition = [[1, 1], [1, 0]]
    return matrix_power(transition, n)[1][0]

assert fibonacci(0) == 0
assert fibonacci(1) == 1
assert fibonacci(10) == 55
assert fibonacci(100) == 354224848179261915075

六、正确性与复杂度

快速幂循环维护两个量：result 是已经处理的指数位所对应的矩阵乘积，factor 是当前二进制位的幂。指数为奇数时将该幂乘入答案；随后把底数平方、指数右移。循环结束时所有指数位都已处理，因此答案为原矩阵的指定次幂。由于矩阵乘法满足结合律，分组顺序不改变结果。

对 d×d 矩阵，朴素相乘耗时 O(d³)；快速幂做 O(log n) 次乘法，总时间 O(d³ log n)，矩阵存储空间 O(d²)。若采用取模版本，每次乘加后要按题意取模，并注意 Java 中间乘积也可能在取模前溢出。

七、边界与常见错误

1. 指数 0 的结果：矩阵的 0 次幂是单位矩阵，初始答案必须是单位矩阵。
2. 矩阵方向混淆：行向量与列向量写法会导致转移矩阵转置；本文固定列向量。
3. 斐波那契下标错一位：状态定义是 [F(n+1), F(n)]，所以返回矩阵幂结果的 [1][0]。
4. 只把底数平方却漏掉答案累乘：指数的每个二进制位都要处理。
5. 不检查方阵：矩阵幂仅对方阵定义。
6. 整数溢出：Java long 不是任意精度；改用取模或 BigInteger 需同步替换矩阵运算。
7. 递推不是线性固定阶：若转移依赖输入变化或非线性条件，不能直接套用固定矩阵幂。

小结

矩阵快速幂将固定线性转移的重复应用转为矩阵乘方，再通过指数二分把步数降到对数级。先选能表达递推历史的状态向量，再推导转移矩阵，最后核对初始向量与目标下标。