CC++ & Algorithm

矩阵快速幂求斐波那契数列

极难2
语言版本：通用
概述：将斐波那契递推转化为矩阵乘法，用矩阵快速幂在 O(log n) 时间内求出第 n 项。

用矩阵快速幂秒杀斐波那契数列——从递推到指数爆炸

开头：为什么你需要这个技巧

假设你想知道斐波那契数列的第100项，用循环加法很快就算出来了。但假如有人问你第 10^18 项呢？循环一亿次都嫌慢，更别提 10^18 次了，太阳爆炸都算不完。这时你需要一个神技：矩阵快速幂。它能让你在 O(log n) 的时间里求出巨大的斐波那契数，哪怕 n 是 10^18，也只需要大约 60 次计算。

斐波那契数列你一定不陌生：0, 1, 1, 2, 3, 5, 8, 13, 21, ... 定义是
F(0)=0, F(1)=1, F(n)=F(n-1)+F(n-2)。
如果让你求第 100 项，用递归或循环 O(n) 可以，但如果 n = 10^18 呢？循环上亿次肯定不行。这时候矩阵快速幂就登场了，它可以把时间复杂度降到 O(log n)。怎么做到的呢？让我们把递推关系写成矩阵形式。

从递推到矩阵：把加法变成乘法

递推的矩阵表示

观察递推式：F(n) = F(n-1) + F(n-2)。
如果我们能构造一个矩阵，使得下一个状态 [F(n), F(n-1)] 由上一步状态 [F(n-1), F(n-2)] 乘以某个矩阵得到，那该多好。试试看：

想要
[F(n), F(n-1)] = [F(n-1)+F(n-2), F(n-1)]
左边是 1×2 向量，右边是 1×2 向量。但矩阵乘法需要左边的向量是行向量？其实，我们可以用列向量更方便。
令向量 V(n) = [F(n); F(n-1)]（2×1 列向量）。那么：

V(n) = [F(n) ] = [1F(n-1) + 1F(n-2)]
[F(n-1) ] [1F(n-1) + 0F(n-2)]

写成矩阵乘法：V(n) = M * V(n-1)，其中
M = [[1, 1], [1, 0]]。

举一个生活中的例子：假设你每天攒钱，第一天有 1 元，第二天有 0 元，之后每天的钱 = 前一天的钱 + 前两天的钱（很奇怪但类似斐波那契）。那么 V(n) 就是第 n 天和第 n-1 天的钱数。矩阵 M 就像一个“攒钱规则”，每次把状态 V(n-1) 变成 V(n)。只要一直乘以这个矩阵，就能算出任何一天的钱数。

验证：
M * [F(n-1); F(n-2)] = [1F(n-1)+1F(n-2), 1F(n-1)+0F(n-2)] = [F(n), F(n-1)]，正好是 V(n)。

所以 V(n) = M * V(n-1)。
继续：V(n-1) = M * V(n-2)，所以 V(n) = M^2 * V(n-2)。
一直递推下去，得到：
V(n) = M^n * V(0)，其中 V(0) = [F(0); F(-1)]？
注意边界：通常 F(0)=0, F(1)=1，那么 V(1)=[F(1); F(0)]=[1;0]。
所以 V(n) = M^(n-1) * V(1) for n>=1。
也可以写成：M^n 乘以 [1;0] 得到 [F(n+1); F(n)]。
但为了简单，我们让 n 对应 F(n+1)。常见做法：首先定义初始向量 V0 = [F(1), F(0)] = [1,0]（写成列向量）。
则 [F(n+1); F(n)] = M^n * V0。
因此 F(n) 就是 M^(n-1) * V0 的第二个元素（对于 n>=1）。或者直接定义 M 的 n 次方的左上角就是 F(n)。

事实上，用归纳法可以证明：
M^n = [[F(n+1), F(n)], [F(n), F(n-1)]]，其中 F(0)=0, F(1)=1。
所以 M^n 的左上角就是 F(n+1)，右上角是 F(n)。因此，要得到 F(n)，只需计算 M^(n-1) 的右上角元素（或 M^n 的右上角？需要小心）。最好采用标准方法：
令 M = [[1,1],[1,0]]，则 M^n 的左上角 = F(n+1)，右上角 = F(n)，左下角 = F(n)，右下角 = F(n-1)。
所以 F(n) = M^n 的右上角元素。那直接计算 M^n 取右上角即可，其中 F(0) 可以通过 n=0 得到（M^0 = I，右上角是0）。

这样，求 F(n) 就变成了计算矩阵 M 的 n 次幂，然后取右上角元素。

矩阵快速幂：如何用 log 时间求次幂？

普通求一个数的 n 次幂可以用快速幂（二进制分解）。同样，矩阵也可以快速幂：

  • 把指数 n 写成二进制，例如 n=13 = 1101_2。
  • 从低位到高位：如果当前位是1，就把当前的累乘结果乘上相应幂次的矩阵。
  • 每一步把底数平方，前进到下一个二进制位。

因为矩阵乘法满足结合律，所以这种方法正确。每次乘法是常数时间（2×2 矩阵），复杂度 O(log n)。

生活类比：你有一张优惠券，每 2 天可以翻倍一次。如果你有 100 天，需要翻 50 次才能拿到全部优惠？不，如果你每次把剩余天数减半，同时将翻倍次数进行平方，就能很快算出结果。这就是快速幂的思路。

算法步骤

  1. 定义 2×2 矩阵 M = [[1,1],[1,0]]。
  2. 若 n=0，直接返回 0。
  3. 计算 M^(n-1) 或 M^n（根据约定）。为了统一，我们用 M = [[1,1],[1,0]]，求 M^n，取右上角元素作为 F(n)。对于 n=0，M^0 = I，右上角是0，正确。对于 n=1，M^1 右上角是1，正确。所以直接计算 M^n 即可。
  4. 使用矩阵快速幂计算 M^n，然后返回结果矩阵的 (0,1) 元素（第0行第1列）。

注意：结果可能很大，可以取模。

新手容易犯的常见错误

  1. n 的索引混乱：不同资料对 F(0)=0 或 F(1)=1 的约定可能不同。检查你的输出：当 n=0 应该返回 0，n=1 返回 1，n=2 返回 1 等等。
  2. 忘记处理 n=0：直接计算 M^0 右上角是 0，没问题，但如果你用 M^(n-1) 就会出问题。统一用 M^n 更安全。
  3. 取模只对加法和乘法取：矩阵乘法中要每一步都取模，否则中间结果可能爆炸。
  4. 用 long long 却忘了取模溢出：虽然矩阵元素不大，但两个 long long 相乘会溢出，建议在乘法后立即取模。
  5. 把行和列搞反：我们用的是 2×2 矩阵，右上角索引是 [0][1]。

完整可运行的代码示例

C++ 代码

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

// 2x2 矩阵乘法，可带模数
vector<vector<long long>> matMul2x2(const vector<vector<long long>>& A,
                                    const vector<vector<long long>>& B,
                                    long long mod = 0) {
    vector<vector<long long>> C(2, vector<long long>(2, 0));  // 结果矩阵初始化为0
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            long long sum = 0;  // 临时累加器
            for (int k = 0; k < 2; k++) {
                sum += A[i][k] * B[k][j];  // 矩阵乘法
                if (mod) sum %= mod;       // 取模
            }
            C[i][j] = sum;
        }
    }
    return C;
}

// 2x2 矩阵快速幂
vector<vector<long long>> matPow2x2(vector<vector<long long>> M, long long n, long long mod = 0) {
    vector<vector<long long>> res = {{1, 0}, {0, 1}}; // 单位矩阵，相当于1
    while (n > 0) {
        if (n & 1LL) {                    // 如果当前二进制位是1
            res = matMul2x2(res, M, mod); // 乘上当前幂次的矩阵
        }
        M = matMul2x2(M, M, mod);        // 底数平方
        n >>= 1;                         // 右移一位
    }
    return res;
}

// 求斐波那契数列第 n 项，返回 F(n)
long long fibonacci(long long n, long long mod = 0) {
    if (n == 0) return 0;
    // 基矩阵 M = [[1,1],[1,0]]
    vector<vector<long long>> M = {{1, 1}, {1, 0}};
    auto Mn = matPow2x2(M, n, mod);  // 计算 M^n
    // 右上角元素 (0,1) 对应 F(n)
    long long ans = Mn[0][1];
    if (mod) ans %= mod;  // 再取一次模确保结果正确
    return ans;
}

int main() {
    long long n = 50;  // 第50项斐波那契数
    cout << "F(" << n << ") = " << fibonacci(n) << endl;
    // 输出 12586269025
    // 如果取模，例如 mod = 1000000007
    long long mod = 1000000007;
    cout << "F(" << n << ") mod " << mod << " = " << fibonacci(n, mod) << endl;
    return 0;
}

Python 代码

def mat_mul_2x2(A, B, mod=None):
    """2x2 矩阵乘法，可选模"""
    C = [[0, 0], [0, 0]]  # 结果矩阵
    for i in range(2):
        for j in range(2):
            s = A[i][0] * B[0][j] + A[i][1] * B[1][j]  # 2x2矩阵乘法展开
            if mod is not None:
                s %= mod
            C[i][j] = s
    return C

def mat_pow_2x2(M, n, mod=None):
    """2x2 矩阵快速幂"""
    # 单位矩阵
    res = [[1, 0], [0, 1]]
    base = M
    while n > 0:
        if n & 1:                     # 二进制位为1
            res = mat_mul_2x2(res, base, mod)
        base = mat_mul_2x2(base, base, mod)  # 平方
        n >>= 1                       # 右移
    return res

def fibonacci(n, mod=None):
    """返回第 n 项斐波那契数"""
    if n == 0:
        return 0
    M = [[1, 1], [1, 0]]
    Mn = mat_pow_2x2(M, n, mod)   # M^n
    return Mn[0][1]  # 右上角

# 测试
if __name__ == "__main__":
    n = 50
    print(f"F({n}) = {fibonacci(n)}")  # 12586269025
    mod = 1000000007
    print(f"F({n}) mod {mod} = {fibonacci(n, mod)}")

为什么这么快？

因为矩阵快速幂只用 O(log n) 次矩阵乘法，每次乘法运算量固定（2×2 矩阵乘法只有8次乘法和4次加法）。所以对于 n=10^18，只需约60次乘，瞬间就能算出。而普通循环需要 10^18 次，地球都毁灭了还没有结果。

更多练习与拓展

练习题

  1. 基础验证：用矩阵快速幂求 F(100) 并验证与循环结果一致（F(100) = 354224848179261915075）。
  2. 批量输出：修改代码，求 F(1) 到 F(10) 并打印。
  3. 取模应用：如果要求 F(n) 模一个很大的质数（如 1000000007），怎么改？已经实现了。你可以试试 n=10^18 模 1e9+7 的结果。
  4. 挑战：尝试用矩阵快速幂求“卢卡斯数列”（Lucas numbers），它们的递推公式也是 L(n) = L(n-1)+L(n-2)，但初始值不同：L(0)=2, L(1)=1。请修改矩阵或初始向量来实现。

与真实世界的联系

  • 你在银行存钱，复利计算可以用快速幂；
  • 某些游戏中的“每日签到奖励”如果按斐波那契增长，可以用这个快速算出第100天的奖励；
  • 编程竞赛中，很多线性递推（比如泰波那契数列、佩尔数列）都能用矩阵快速幂解决。只要递推式是常系数线性齐次的，就可以构造矩阵。

总结：从递推到矩阵的魔法

将线性递推转化为矩阵乘法是一种强大的技巧。对于斐波那契数列，我们只需要一个简单的 2×2 矩阵，就能在 log 时间内求出任意项。这种思想可以推广到任意常系数线性递推（如高阶斐波那契等）。矩阵快速幂让原本 O(n) 的问题变成 O(log n)，是算法竞赛中必会的神技。

相关指引

  • 如果你想深入理解矩阵快速幂在其他问题中的应用，可以学习：
    • “常系数齐次线性递推”的矩阵解法
    • 用矩阵快速幂求“阿格拉沃尔数列”
    • 矩阵快速幂在图论中计算路径数量（如“恰好走 k 步”的方案数）
  • 更高级的内容：使用 特征多项式 优化到 O(k^2 log n) 甚至 O(k log k log n)，但那是另一座高山了。

现在，你可以用矩阵快速幂去挑战那些天文数字般的斐波那契项了！

例题精讲

1单选题

使用矩阵快速幂计算斐波那契数列的第n项（n≥0），算法的时间复杂度是？

AO(n)
BO(log n)
CO(n²)
DO(2ⁿ)
2单选题

已知斐波那契数列递推关系：F(n)=F(n-1)+F(n-2)，采用矩阵快速幂，对应的递推矩阵是？

A[[1,1],[1,0]]
B[[1,0],[1,1]]
C[[1,1],[0,1]]
D[[1,0],[0,1]]
3判断题

定义F(0)=0, F(1)=1，递推矩阵M=[[1,1],[1,0]]，则Mⁿ的第一行第一列元素等于F(n+1)。

4填空题
以下是一段用矩阵快速幂求斐波那契数列第n项（模MOD）的代码，补全矩阵乘法函数中的空缺。

int MOD = 1000000007;
struct Matrix { long long a[2][2]; };
Matrix mul(Matrix A, Matrix B) {
    Matrix C = {0};
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 2; j++)
            for (int k = 0; k < 2; k++)
                C.a[i][j] = (C.a[i][j] + ___) % MOD;
    return C;
}
Matrix pow(Matrix M, int n) {
    Matrix res = {1,0,0,1}; // 单位矩阵
    while (n) {
        if (n & 1) res = mul(res, M);
        M = mul(M, M);
        n >>= 1;
    }
    return res;
}
5单选题

在矩阵快速幂求解斐波那契数列时，若n极大（如10¹⁸）且未对中间结果取模，以下说法正确的是？

A结果依然精确，因为矩阵快速幂数值稳定
B结果可能溢出，但取模后计算仍可得到正确余数
C必须使用高精度整数存储，否则无法计算
D矩阵快速幂无法处理n超过10⁹的情况