很多刷《信息学奥赛一本通》提高篇的同学，看到 1644 这题都会有点发怵。题目名字叫“佳佳的 Fibonacci”，看似只是求斐波那契相关的和，但 n 的范围给到 10^18，普通的 for 循环连边都摸不到。第一次做的时候我也被这个 n 吓了一跳，后来才意识到这题本质上是矩阵快速幂的经典入门例题，考的是能不能把“带权斐波那契前缀和”转化为矩阵状态转移。今天把这题的完整思路、推导过程和代码一次性说透，尤其适合卡在数学推导上的朋友。
1. 题目到底在问什么：读懂“佳佳的 Fibonacci”的核心要求
1.1 题干拆解与数据范围解析
这道题的实际要求是这样的：已知斐波那契数列 F(1)=1，F(2)=1，对于 n≥3 有 F(n)=F(n-1)+F(n-2)。现在定义
T(n) = 1×F(1) + 2×F(2) + 3×F(3) + ... + n×F(n)
输入一行两个整数 n 和 m，要求输出 T(n) mod m 的值。
注意这个“mod m”非常关键。m 没有给任何额外限制，也就是说 m 不一定是质数，甚至可能很小。这个细节直接决定了我们不能往数论逆元、费马小定理那类方向想，因为那些方法大多要求模数是质数。题目老老实实只给加法和乘法，那就是在提示：请用只依赖加法和乘法的算法来解决，矩阵快速幂正好符合这个条件。
数据范围方面，n 可以大到 10^18，这基本堵死了递推求和的路。哪怕你从 F(1) 一路推到 F(n)，光求斐波那契数列本身用 O(n) 也完全扛不住。所以需要一种能把 O(n) 压缩到 O(log n) 的工具，矩阵快速幂就是这个工具。
1.2 为什么不能用暴力：从 O(n) 到 O(log n) 的动机
有些同学会觉得，斐波那契数列第 n 项不是有通项公式吗？确实有，通项公式长这样：
F(n) = (1/√5) × [((1+√5)/2)^n - ((1-√5)/2)^n]
但这个公式里带着无理数和除法，放到整数取模的场景里极其难受，尤其是 m 还是任意整数的时候。你想用快速幂去算无理数的高次幂，再对任意 m 取模，中间必然出现浮点误差，结果根本不可靠。
再看 T(n) 的定义，它比单纯的 F(n) 更麻烦，因为每一项都乘了一个位置系数 i。这个系数 i 会随着递推不断变化，没办法像常量一样提到外面。不过换个角度想：n 虽然大，但如果我们把“n 乘以某个斐波那契项”也当成一种状态，就能构造出一个线性递推关系，然后用矩阵乘法一步到位。这就是这道题整个解法的灵魂。
2. 解这类题的通用武器：矩阵快速幂是什么、为什么能用
2.1 线性递推与矩阵乘法的对应关系
先回忆一下最简单的斐波那契递推。F(n+1) = F(n) + F(n-1)，F(n) = F(n)，这个二元组可以写成矩阵形式：
[ F(n+1) ] = [ 1 1 ] [ F(n) ] [ F(n) ] [ 1 0 ] [ F(n-1) ]
每次乘一次这个 2×2 矩阵，二元组就从“第 n 项的状态”推进到“第 n+1 项的状态”。这个思想可以推广：只要下一个时刻的每个新状态都能表示成当前时刻各个状态的线性组合，那么整个递推过程就是一次矩阵乘法。
为什么这个形式有用？因为矩阵乘法有结合律。计算 A 的 n 次方，可以像普通整数快速幂一样，把指数 n 拆成二进制，每次把矩阵平方，复杂度从 O(n) 降到 O(log n)。对一个二进制位数为 60 的数来说，log 级别的计算量几乎可以忽略不计。
2.2 快速幂思想：把 n 次递推压缩到 log n 次矩阵乘法
快速幂的本质是“二分加速”。想要求矩阵 A 的 n 次方，先看 n 的二进制表示。比如 n=13，二进制是 1101，也就是 13=8+4+1。那我们只需要预先算好 A^1、A^2、A^4、A^8，再把这四个乘起来就行了。每一次乘法都是 O(矩阵大小^3)，矩阵大小在这里最多也就是 5 阶，所以整体复杂度非常低。
实际写代码时用的是迭代式快速幂，维护一个结果矩阵 res，初始为单位矩阵，然后不断右移指数 b。如果当前最低位是 1，就把 res 乘上当前的 A；每一轮都把 A 自乘一次。单位矩阵在这里相当于整数快速幂里的初始值 1，乘任何矩阵都保持不变，所以结果一定对。
2.3 本题最关键的一步：加权系数 n 为什么会成为“状态变量”
把 T(n) 与 T(n-1) 做差，立刻得到：
T(n) - T(n-1) = n×F(n)
这个式子告诉我们，只要能快速求出每一个 n×F(n)，再加起来就能得到 T(n)。难点在于 n×F(n) 的递推不是那么直接。n 每次加 1，前面乘的系数也加 1，这就不是固定系数递推了。
但仔细展开一下：
(n+1)×F(n+1) = (n+1)×(F(n) + F(n-1)) = n×F(n) + n×F(n-1) + F(n) + F(n-1)
看到没有？右边的四样东西分别是 n×F(n)、n×F(n-1)、F(n)、F(n-1)，这些全部都可以成为状态，而且下一时刻的状态是当前状态的线性组合。这就是“把动态系数变成新增状态”的典型操作。原本让人头疼的 n，被并进了状态向量里，整个递推重新变成了线性齐次递推。
3. 从零推导 5 维矩阵：逐步构造状态转移
3.1 状态向量的选择：F(n), F(n-1), nF(n), nF(n-1), T(n)
基于上面的分析，我构造的 5 维状态向量是：
v(n) = [ F(n), F(n-1), n×F(n), n×F(n-1), T(n) ] 的转置
五个分量分别对应五个“量”：
- x1 = F(n)：当前斐波那契项
- x2 = F(n-1)：前一项，斐波那契二项递推需要它
- x3 = n×F(n)：位置系数乘以当前项，直接关系到 T 的增量
- x4 = n×F(n-1)：位置系数乘以前一项，用于构造 x3 的下一时刻
- x5 = T(n)：我们最终要求的答案
之所以要 x4，是因为要算 (n+1)×F(n+1)，展开后会出现 n×F(n-1)，它不在别的状态里，必须有专门的位置存。
3.2 逐行推导转移矩阵，为什么每一行这样写
现在从 v(n) 推到 v(n+1)，设转移矩阵为 M，即 v(n+1) = M × v(n)。分别看五个新状态怎么由旧状态表示。
第一行，新的第一个状态是 F(n+1)： F(n+1) = F(n) + F(n-1) = x1 + x2 所以 M 的第一行是 [1, 1, 0, 0, 0]。
第二行，新的第二个状态是 F(n)： 它就是旧的 x1，所以第二行是 [1, 0, 0, 0, 0]。
第三行，新的第三个状态是 (n+1)×F(n+1)。刚才推导过： (n+1)×F(n+1) = n×F(n) + n×F(n-1) + F(n) + F(n-1) = x3 + x4 + x1 + x2 所以第三行是 [1, 1, 1, 1, 0]。
第四行，新的第四个状态是 (n+1)×F(n)： (n+1)×F(n) = n×F(n) + F(n) = x3 + x1 所以第四行是 [1, 0, 1, 0, 0]。
第五行，新的第五个状态是 T(n+1)： T(n+1) = T(n) + (n+1)×F(n+1) = x5 + (x1 + x2 + x3 + x4) 所以第五行是 [1, 1, 1, 1, 1]。
整理成矩阵：
M = [ 1 1 0 0 0 ] [ 1 0 0 0 0 ] [ 1 1 1 1 0 ] [ 1 0 1 0 0 ] [ 1 1 1 1 1 ]
做到这步，最核心的数学部分已经完成。剩下的就是把初始状态放进去，然后乘上 M 的 n-1 次方。
3.3 初始向量与答案提取：v(1) 与 T(n) 的关系
接下来确定初始向量 v(1)。注意递推里用到了 F(0)，为了满足 F(2)=F(1)+F(0)=1，取 F(0)=0 是合理的。
v(1) = [ F(1), F(0), 1×F(1), 1×F(0), T(1) ] = [ 1, 0, 1, 0, 1 ]
这里 T(1) = 1×F(1) = 1，所以第五个分量是 1。
因此最终答案就是矩阵乘法的第五个分量：
answer = (M^(n-1) × v(1)) 的第 5 个元素
n=1 时，M^0 是单位矩阵，乘完 v(1) 后第五个分量依然是 1，符合 T(1)=1 的要求，所以不需要另加特判。
4. 完整代码实现与关键细节
4.1 C++ 矩阵快速幂模板（针对本题）
下面是完整的 C++ 实现。矩阵大小固定为 5，代码里直接用 5×5 数组，清晰直接：
#include <iostream> #include <cstring> using namespace std; typedef long long ll; const int SZ = 5; ll n, m; struct Mat { ll a[SZ][SZ]; Mat(bool flag = false) { memset(a, 0, sizeof(a)); if (flag) { for (int i = 0; i < SZ; i++) a[i][i] = 1; } } }; Mat mul(const Mat& A, const Mat& B) { Mat C; for (int i = 0; i < SZ; i++) { for (int k = 0; k < SZ; k++) { if (A.a[i][k] == 0) continue; for (int j = 0; j < SZ; j++) { C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j]) % m; } } } return C; } Mat pow_mat(Mat A, ll p) { Mat res(true); while (p) { if (p & 1) res = mul(res, A); A = mul(A, A); p >>= 1; } return res; } int main() { cin >> n >> m; Mat M; M.a[0][0] = 1; M.a[0][1] = 1; M.a[1][0] = 1; M.a[2][0] = 1; M.a[2][1] = 1; M.a[2][2] = 1; M.a[2][3] = 1; M.a[3][0] = 1; M.a[3][2] = 1; M.a[4][0] = 1; M.a[4][1] = 1; M.a[4][2] = 1; M.a[4][3] = 1; M.a[4][4] = 1; Mat P = pow_mat(M, n - 1); ll init[SZ] = {1, 0, 1, 0, 1}; ll ans = 0; for (int i = 0; i < SZ; i++) { ans = (ans + P.a[4][i] * init[i]) % m; } cout << ans << endl; return 0; }代码里最核心的乘法我做了个小优化：第二层循环只用 k 来过滤，如果 A.a[i][k] 是 0 就跳过。这个优化对矩阵快速幂很有效，因为矩阵很多元素是 1，但过滤掉 0 能让常数变小，实际跑起来速度快不少。
4.2 边界情况的处理与取模陷阱
第一个边界情况是 n=1。上面代码里 pow_mat 的 p=n-1=0，res 初始为单元矩阵，乘完还是单元矩阵。ans 算出来就是 init[4]=1，所以输出 1 mod m。m=1 时会输出 0，也符合模运算的定义，不需要特判。
第二个容易踩的坑是 m 很小甚至等于 1。取模运算在矩阵乘法里每一步都做，m=1 时所有结果全变成 0，输出就是 0，这是对的。但如果你只在最后取模，中间乘法溢出，结果就完全错了。所以务必在矩阵乘法的内层循环里就取模，不要等到最后统一取。
第三个坑是乘法的溢出。矩阵元素的值都在 [0, m-1] 之间，两个元素相乘最大约是 (m-1)^2。如果 m 给到 10^9，平方接近 10^18，long long 刚好还能扛住（上限约 9.22×10^18）。但如果 m 给到 10^18 量级，直接乘就会溢出。保守做法是写一个快速乘函数，把乘法改成二进制加法；或者用 GCC 的 __int128 类型。竞赛环境一般支持 __int128，我建议评测范围不明时直接上 __int128，省心。
// 如果想更保险，可以把乘法函数改成： Mat mul(const Mat& A, const Mat& B) { Mat C; for (int i = 0; i < SZ; i++) for (int k = 0; k < SZ; k++) { if (A.a[i][k] == 0) continue; for (int j = 0; j < SZ; j++) { C.a[i][j] = (C.a[i][j] + (__int128)A.a[i][k] * B.a[k][j]) % m; } } return C; }4.3 一个更快的写法：状态向量直接继续乘矩阵
上面的做法是算 M^(n-1) 之后再去乘初始向量。更常见的竞赛写法是先把初始向量看成 5×1 的矩阵，直接把两者一起放进快速幂循环里，省一次矩阵乘法的常数。不过这种写法对初学者来说不太直观，我在上面保留了“先算矩阵幂再乘向量”的方案，逻辑更清楚。真正卡时间时再考虑合并优化即可。
5. 常见问题与实战排查
5.1 矩阵乘法顺序搞反导致答案错误
快速幂里 res = mul(res, A) 和 res = mul(A, res) 是两种写法，结果多数情况下不同，因为矩阵乘法不满足交换律。很多人习惯从整数快速幂模板迁移，在整数里乘的顺序无所谓，到了矩阵就翻车。
解决方案是固定使用“res 乘在左边、A 乘在右边”或始终“把本轮矩阵 A 乘到 res 右侧”，并配合初始单位矩阵验证。验证方法：把指数 p 设成 1，看结果是否等于 A；设成 2，看是否等于 A×A。如果顺序反了，p=2 时很可能得到不同答案。
5.2 单位矩阵初始化错误
单位矩阵写错是另一个常见低级错误。单位矩阵是对角线全 1，其他地方全 0，它在矩阵乘法里扮演“1”的角色。有些人会把单位矩阵初始成全 1 矩阵，结果快速幂结果被整体放大，T(n) 自然不对。
给一个简单的自测：写一个函数打印矩阵，把 pow_mat 的指数先设成 0 和 1，肉眼检查输出是否分别为单位矩阵和原矩阵。
5.3 斐波那契下标不一致导致初始向量错误
题目给定的 F(1)=1，F(2)=1。为了方便矩阵递推，我会引入 F(0)=0。这个扩展是合法的，因为 F(2)=F(1)+F(0)=1+0=1。但有些同学习惯用 F(0)=1、F(1)=1 的版本（从 F(0) 开始定义），如果不统一，初始向量就会错。
建议把初始向量写成表格，逐个分量核对：
| 分量 | 含义 | 初始值 |
|---|---|---|
| 第1个 | F(1) | 1 |
| 第2个 | F(0) | 0 |
| 第3个 | 1×F(1) | 1 |
| 第4个 | 1×F(0) | 0 |
| 第5个 | T(1) | 1 |
5.4 大 n 下用 long long 的问题
n 本身是 10^18 量级，读入要用 long long，不能用 int。还有矩阵快速幂里的 p 参数也要是 long long。别小看这个，很多人把 pow_mat 的参数写成 int，n=10^18 时直接溢出成负数或截断值，结果完全不可预测。我见过不少同学在这上面卡了很久，最后发现只是函数签名写错了。
5.5 验证技巧：小数据暴力对拍
代码写完后不要急着提交，先用小数据验证。自己另外写一个暴力程序，对 n=1 到 n=20 的所有结果逐一对比。暴力程序直接按递推式算 F 和 T，逻辑非常简单。对拍结果全对之后再提交，能帮你过滤掉绝大部分推导和实现错误。
真实对拍示例，比如：
- n=2, m=1000：T(2)=1×1+2×1=3，输出 3
- n=3, m=1000：T(3)=1+2+3×2=9，输出 9
- n=4, m=1000：T(4)=1+2+6+4×3=21，输出 21
这几个小值都和矩阵代码跑出来的结果一致，基本可以放心提交。
6. 从这道题延伸开：矩阵快速幂的更多玩法
6.1 同一思路的变形题目
理解了“把带位置系数的量放进状态向量”这个技巧之后，很多题都可以秒解。比如说求 S(n)=F(1)+F(2)+...+F(n)，只需要在二元组后面加一个前缀和状态；求带平方系数 Σ i^2×F(i)，就再把 n^2×F(n)、n^2×F(n-1)、n×F(n)、n×F(n-1) 一起纳入状态向量。每多一类“动态系数”，就多几个状态，但方法论完全一样。
这也是为什么很多老师会说矩阵快速幂是“套路题”：只要你能设计出状态向量，剩下的就是机械地构造矩阵。
6.2 斐波那契数列相关的经典矩阵构造汇总
| 需求 | 状态向量 | 说明 |
|---|---|---|
| F(n) | [F(n), F(n-1)] | 最基本的二元递推 |
| 前缀和 S(n) | [F(n), F(n-1), S(n)] | 加一个累加状态 |
| 加权和 T(n) | [F(n), F(n-1), nF(n), nF(n-1), T(n)] | 本题方法 |
| 两个斐波那契项相乘 | [F(n)F(n), F(n)F(n-1), F(n-1)F(n-1)] | 利用乘法分配律展开 |
这些构造共同的核心思想是：每一个“看起来不齐次”的项，都想办法用已有状态的不同倍数组合表示出来，只要能做到，就一定能写成矩阵。
6.3 为什么这个知识点在竞赛里如此重要
矩阵快速幂涉及的远不只是斐波那契。任何线性递推，包括常系数齐次线性递推、部分含有简单多项式系数的非齐次递推，只要项数不多，都可以转成矩阵幂。再加上快速幂本身是 O(log n)，对题目里动不动给到 10^18 的数据范围有着天然优势。所以它在信息学奥赛里几乎是必考考点，从普及组到提高组再到省选，都会以各种形式出现。
我个人在做这类题时最大的体会是：矩阵推导阶段宁可慢一点，把每一行转移都写在纸上，也不要直接在脑子里空想。把 v(n) 的五个分量列出来，再写 v(n+1) 的五个分量一一对应，十分钟就能保证矩阵构造不出错；省掉这一步直接拍代码，调试反而可能花掉一小时。题解看懂只是第一步，自己独立推一遍、再拿小数据对拍一次，这个矩阵你才能算真正会用了。后面再遇到“佳佳的 Fibonacci”这类带权前缀和的题目，看到 n=10^18 就不会慌，直接条件反射地往状态向量里塞系数项，这就是一道送分题了。
