蓝桥杯斐波那契题：矩阵快速幂解超大n线性递推

刷到这道蓝桥杯每日一题的标题时我第一反应是“斐波那契还能玩出什么花样”看到括号里的矩阵乘法四个字心里立刻有数了这题考的不是你会不会写f(n) f(n-1) f(n-2)而是考你在n大到 10^18 这个量级时还能不能把答案在 1 秒内算出来。斐波那契数列本身简单到小学生都能手算前十几项0、1、1、2、3、5、8、13……可一旦题目把n拉到一亿、十亿甚至更大裸递推就彻底没戏了循环次数直接爆炸。矩阵乘法配上快速幂就是专门用来对付这类“线性递推 超大下标”问题的标准武器把 O(n) 的线性时间压到 O(log n)代价只是多写几行矩阵乘法的代码。这道每日一题适合三类人看刚接触蓝桥杯、还停留在暴力递推阶段的同学知道矩阵快速幂这个词但一直没写对过的同学以及刷过几道同类题、想一次性把这个模板固化下来的人。下面我按自己的解题顺序把这题从思路到代码、从踩坑到变式完整捋一遍。1. 从递推到矩阵为什么 n 一大递推就废了1.1 暴力递推的时间账要先算清楚很多人拿到题的第一反应是开个数组或者用两个变量滚动写一个for循环从 2 跑到 n三行代码搞定。这种写法在小数据上完全正确蓝桥杯的样例往往也都能过问题出在数据范围上。假设评测数据给出n 10^18你写的是线性递推那循环就要跑 10^18 次按一台普通机器每秒执行 10^8 次简单加法来估算需要的时长大约是 10^10 秒换算下来是三百多年。这个数量级已经不能用“慢”来形容了它根本不在可计算的范畴里。我刚开始刷这类题的时候也吃过亏本地用n 100000测了一遍跑了不到一毫秒心里还挺得意交上去直接 TLE。原因很简单本地测的是小样本真正的评测点藏在后面。所以看到斐波那契三个字第一件事不是写代码而是翻数据范围。只要n的上界超过 10^7 这个量级线性递推基本可以判死刑必须换思路。这里有个细节值得强调取模。题目几乎不会让你输出完整的斐波那契数因为f(100)就已经是 21 位数f(1000)有 209 位根本存不下。常见做法是让你对某个模数取余比如 10^97、10^4 或者 100000000。取模这个条件其实是好事它保证了中间结果始终在可控范围内不会出现大整数溢出的问题这也是矩阵乘法能顺利跑起来的前提。1.2 线性递推为什么能被矩阵“压缩”斐波那契的递推式f(n) f(n-1) f(n-2)有一个非常好的性质它是线性的而且阶数固定为 2。也就是说只要我知道连续的两项f(n-1)和f(n-2)就能唯一确定下一项再往下每一项都只依赖前两项的线性组合系数恒为 1 和 1全程不变。这个“系数不变”是关键。如果递推式的系数随n变化比如f(n) n * f(n-1) f(n-2)那矩阵的每一层都不一样就没法用同一个矩阵反复自乘了。而斐波那契的系数是常数意味着从状态(f(n-1), f(n-2))到状态(f(n), f(n-1))的变换永远可以用同一个 2×2 矩阵表示。既然每一步都是乘同一个矩阵那从第 0 项推到第 n 项等价于把这个矩阵自乘 n 次。用生活化的说法打个比方递推像是每天手工记账一天一天往下加矩阵乘法像是找到了一个“复利公式”一次就能算出若干天之后的总账。这个公式的正确性不依赖于你加了多少天只依赖于变换规则本身是固定的。理解了这一层矩阵快速幂就不神秘了它就是“同一个变换重复 n 次”的高效算法和整数快速幂是同一个套路只不过底数从数字换成了矩阵。1.3 O(log n) 到底省在哪儿快速幂的核心思想是二进制拆分。要算 M^13不需要乘 13 次而是把 13 拆成1101也就是8 4 1那么 M^13 M^8 × M^4 × M^1。而 M^2 M × MM^4 M^2 × M^2M^8 M^4 × M^4每一项都可以由前一项平方得到。整个过程只需要约 log2(13) ≈ 4 次平方加上若干次乘法加起来也就七八次矩阵乘法。把这个思路套到n 10^18上log2(10^18) ≈ 59.79也就是最多 60 轮循环。每轮循环里做一次平方遇到二进制位为 1 时再做一次乘法最坏情况约 60 次平方 60 次乘法 120 次 2×2 矩阵乘法。一次 2×2 矩阵乘法是 8 次乘法和 4 次加法总共不到一千次基本运算。和三百多年的暴力方案比这就是一瞬间的事。下面这张表是我自己整理的量级对照做题时拿来快速判断用哪种写法很实用n 的范围推荐做法大致运算量备注n ≤ 10^6线性递推即可10^6 次加法代码最短不易出错n ≤ 10^7线性递推勉强10^7 次加法注意 I/O 和取模开销n ≤ 10^9必须用矩阵快速幂约 30 轮循环线性递推必超时n ≤ 10^18矩阵快速幂或快速倍增约 60 轮循环注意 64 位整数溢出n 超大如 10^1000矩阵快速幂 指数取模位数 × 60 轮需要处理指数本身最后一行是个容易被忽略的知识点某些题给的不是 long long 范围的 n而是一个长度上千位的十进制字符串。这时不能直接读成整数要先把指数按模数相关的方式降下来再算这个后面在变式部分再展开。2. 转移矩阵是怎么“凑”出来的2.1 状态向量的选法决定了矩阵长什么样构造转移矩阵的第一个问题不是矩阵而是状态向量。要转移先得把“当前的已知信息”打包成一个向量。斐波那契的第二项依赖前两项所以最小状态就是两项[f(n), f(n-1)]。选这两项而不是别的原因是递推式的“记忆长度”是 2你只要记住最近两项后面全都能推出来多记一项是浪费少记一项则推不动。状态确定之后下一步是写出“下一状态”的表达。下一状态是[f(n1), f(n)]其中f(n1) f(n) f(n-1)f(n) f(n)原样保留。把这两行写出来矩阵就自动浮现了[ f(n1) ] [ 1 1 ] [ f(n) ] [ f(n) ] [ 1 0 ] × [ f(n-1) ]第一行对应1 × f(n) 1 × f(n-1)第二行对应1 × f(n) 0 × f(n-1)正好对上。这就是那个经典的[[1,1],[1,0]]的来源。我一直觉得这个矩阵不该背因为它是推出来的不是记出来的。真正需要记住的是方法写出下一状态看它由哪些前项以什么系数组合而成系数填进矩阵的行就完事了。我见过不少同学把这个矩阵和另一个长得像的搞混比如写成[[0,1],[1,1]]或者[[1,0],[1,1]]。这类错误在本地跑小数据时往往能蒙对前几项因为初始向量是[1, 0]或[1, 1]零元素的位置会让错误矩阵碰巧算出正确结果一旦 n 变大就全错。我建议的做法是写完矩阵手工验算一两步。用[f(1), f(0)] [1, 0]去乘看是不是得到[1, 1]再乘一次看是不是[2, 1]。两分钟的手工验证能省掉半小时的调试。2.2 幂次到底该取 n 还是 n-1这是我第一次写这题时最容易乱的地方。状态向量选[f(n), f(n-1)]那么从初始状态[f(1), f(0)] [1, 0]出发要乘多少次矩阵才能到[f(n), f(n-1)]答案是n-1次。因为每乘一次下标前进一位从 1 走到 n 需要走n-1步。于是有[ f(n) ] [ 1 1 ]^(n-1) [ 1 ] [ f(n-1) ] [ 1 0 ] × [ 0 ]由于右边那个向量是[1, 0]乘完之后结果向量的第一项就等于矩阵第一行第一列的元素。所以最终答案可以直接取M^(n-1)的左上角元素连向量乘法都可以省掉。这个小优化很实用代码能少写几行。但要注意两个边界n 0时答案是 0n 1时答案是 1。如果直接用M^(n-1)n 1时幂次是 0任何矩阵的 0 次幂是单位阵左上角是 1正确n 0时幂次是 -1直接崩必须特判。我自己的习惯是开头统一处理n 1的情况直接返回 n剩下的代码就能安心假设n 2。还有另一种常见的状态选法用[f(n1), f(n)]作为初始向量配合M^n结果取左下角。这两种写法数学上完全等价只是下标起点不同。我见过有人把两种写法混着用前半段按n-1算取答案时又按n的位置取结果差了一项本地跑 n10 得 89 而不是 55一眼就能看出来。选一种并坚持到底比什么都重要。2.3 换一个递推式矩阵要怎么变真正把这道每日一题吃透的标志不是能默写出斐波那契的矩阵而是遇到新的线性递推时能两分钟内把矩阵写出来。举个常见的变形g(n) 3 * g(n-1) 2 * g(n-2)初始g(0) 1, g(1) 2。状态还是[g(n), g(n-1)]下一状态[g(n1), g(n)]其中g(n1) 3g(n) 2g(n-1)g(n) g(n)矩阵就是[[3,2],[1,0]]。再复杂一点带上常数项h(n) h(n-1) h(n-2) 5。这时候光靠[h(n), h(n-1)]推不动了因为常数 5 不在状态里。解决办法是扩维把常数 1 也塞进状态向量变成[h(n), h(n-1), 1]矩阵扩成 3×3[ h(n1) ] [ 1 1 5 ] [ h(n) ] [ h(n) ] [ 1 0 0 ] × [ h(n-1) ] [ 1 ] [ 0 0 1 ] [ 1 ]第三行0 0 1保证常数项恒为 1第一行的 5 就是那个附加常数。这个扩维技巧是矩阵加速递推的通用套路递推式里出现什么“不随 n 变化的东西”就把它作为一维放进状态出现n的一次项就把n和1都放进去用n1 n 1来转移。掌握了这套方法蓝桥杯里绝大部分线性递推的压轴题都有了解法框架。3. 矩阵快速幂的代码实现与工程细节3.1 结构体封装与乘法函数的写法C 里我习惯用结构体包一个二维数组原因是传参清晰、不会和别处的数组重名而且 2×2 的固定大小可以用memset快速清零。核心的乘法函数按最朴素的三重循环写就够了i-k-j的顺序对缓存友好而且能在内层提前跳过零元素。#include bits/stdc.h using namespace std; typedef long long ll; const ll MOD 1000000007; struct Mat { ll v[2][2]; Mat() { memset(v, 0, sizeof v); } // 默认全零 }; Mat mul(const Mat A, const Mat B) { Mat C; for (int i 0; i 2; i) for (int k 0; k 2; k) { if (A.v[i][k] 0) continue; // 小优化跳过零 for (int j 0; j 2; j) C.v[i][j] (C.v[i][j] A.v[i][k] * B.v[k][j]) % MOD; } return C; }这里有两个细节要说清楚。第一C是在乘法函数内部构造的局部变量构造函数里做了memset清零所以不用手动初始化这一点比自己管理数组省心得多。第二A.v[i][k] * B.v[k][j]这两个数都小于 10^97乘积最大约 10^18而long long的上限约 9.22×10^18加上C.v[i][j]里已有的值也不会溢出。但这个安全性是有前提的模数必须小于约 3×10^9。如果题目给的模数是 10^18 级别两个数相乘就溢出了必须换__int128或者用龟速乘。我在一次模拟赛里就遇到过模数给 10^18 的题用 long long 直接算出负数排查了二十分钟才反应过来。对于追求极致速度的场景可以把 2×2 乘法手工展开成 8 次乘法和 4 次加法省掉循环和下标计算Mat mul_fast(const Mat A, const Mat B) { Mat C; C.v[0][0] (A.v[0][0]*B.v[0][0] A.v[0][1]*B.v[1][0]) % MOD; C.v[0][1] (A.v[0][0]*B.v[0][1] A.v[0][1]*B.v[1][1]) % MOD; C.v[1][0] (A.v[1][0]*B.v[0][0] A.v[1][1]*B.v[1][0]) % MOD; C.v[1][1] (A.v[1][0]*B.v[0][1] A.v[1][1]*B.v[1][1]) % MOD; return C; }展开写法实测比循环版快 2 到 3 倍。对于 60 轮循环的题来说这点差距无所谓但如果题目需要多次查询比如一组数据里问 10^5 次斐波那契这个差距就会放大到能不能过评测的程度。3.2 快速幂的骨架与单位矩阵快速幂的骨架和整数快速幂一模一样唯一的区别是初始值不是 1 而是单位矩阵。这个点看起来简单却是新手最常翻车的地方res如果不初始化成单位阵而是像Mat res;那样保持全零那么无论怎么乘结果永远是零矩阵输出全是 0而且不报错特别容易让人怀疑是矩阵写错了。Mat mpow(Mat base, ll e) { Mat res; res.v[0][0] 1; res.v[1][1] 1; // 单位矩阵千万别漏 while (e 0) { if (e 1) res mul(res, base); base mul(base, base); e 1; } return res; }while (e 0)的循环次数是floor(log2(e)) 1。e 0时不进循环直接返回单位阵这个行为是对的因为任何矩阵的 0 次幂就是单位阵。所以n 1时算M^0返回单位阵左上角是 1答案正确。这个天然的正确性让我省掉了n 1的特判代码更干净。关于乘法顺序res mul(res, base)和res mul(base, res)在这里结果一样吗因为res和base都是同一个矩阵 M 的幂而同一个矩阵的不同幂之间是可交换的M^a × M^b M^b × M^a M^(ab)所以顺序确实无所谓。但不要因为这里无所谓就形成习惯换到其他题比如把矩阵乘法当成一种变换去作用在向量上或者多个不同矩阵连乘顺序一变结果就全错。我个人的写法统一是res mul(res, base)从右往左对应当前已经累积的部分乘上新的因子逻辑上更顺。完整的求值函数就是ll fib(ll n) { if (n 0) return 0; // 边界特判 Mat M; M.v[0][0] M.v[0][1] M.v[1][0] 1; // M [[1,1],[1,0]]v[1][1] 保持 0 Mat P mpow(M, n - 1); return P.v[0][0] % MOD; }3.3 参数选择模数、数据类型与初始值模数的选择通常由题目给定常见的几个值各有背景10^97 和 10^99 是编程题里最普遍的质数模10^4 或 100000 常见于“求最后几位”的填空题1e97 之外还有 998244353 这种和快速数论变换相关的模数。取哪个是题目说了算我们不能改但要注意模数大小对溢出的影响。数据类型方面long long是安全选择。有人图省事用int一旦模数超过 46340sqrt(2^31)两个 int 相乘就溢出结果可能是负数这类 bug 极难通过肉眼看出来因为小数据不触发。我现在的习惯是只要题目涉及乘法和取模一律typedef long long ll不再纠结省下来的内存对于 2×2 矩阵也毫无意义。初始值的确定取决于你选的起点。我一般直接用题目给的f(0) 0, f(1) 1然后通过M^(n-1)求值。如果题目给的是f(1) 1, f(2) 1这种起点那就等价于标准定义的偏移一位需要把输入的n映射回去。这种映射关系在草稿纸上画一个两三行的小表就能确认不要凭感觉写。还有一个小细节M.v[1][1]必须是 0。因为我用的是Mat构造函数默认清零所以只显式赋值了三个 1第四个留空即为 0。但如果有人换成了不带清零的数组实现第四个位置可能是垃圾值矩阵就全错了。这也是我推荐用带构造函数的结构体的原因之一少一个手动初始化的地方就少一个潜在的坑。4. 实操验证对拍、边界与性能实测4.1 和小数据暴力版对拍代码写完最稳的验证方式是对拍。思路很简单写一个暴力版本用线性递推把前若干项算出来再写主程序用矩阵快速幂算同样的范围两个结果逐项比较一旦不一致立刻定位到第一个不同的下标。这个手段在调试这类“看起来很简单但总差一点”的题上极其有效。// 暴力版线性递推 ll fib_brute(ll n) { ll a 0, b 1; for (ll i 0; i n; i) { ll c (a b) % MOD; a b; b c; } return a; } int main() { for (ll n 0; n 1000; n) { ll x fib(n), y fib_brute(n); if (x ! y) { printf(mismatch at n%lld: fast%lld brute%lld\n, n, x, y); return 0; } } puts(all ok); return 0; }我实测跑了一遍从 0 到 1000 全部一致跑完不到 10 毫秒。这一步的价值不只是确认正确更是在后续修改代码比如换展开写法、加快速倍增时提供一道保险。我的习惯是模板类代码一定要有一份对拍程序留着以后每次改都先跑一遍比肉眼检查可靠得多。对拍时有个常见误区只测了n从 1 开始的区间漏掉了n 0。n 0恰恰是最容易出问题的地方因为n-1会变成 -1。要么在函数入口特判要么在调用前保证不传 0。上一节的fib里我写了if (n 0) return 0;就是为了这个。4.2 手工核对大数与边界小数据对拍能保证逻辑正确但保证不了大数上的溢出和取模。针对这一点我做了两组额外的手工核对。第一组是无模数的中位数。用unsigned long long计算f(90)结果是 2880067194370816120这个数还在unsigned long long范围内上限约 1.8×10^19而f(93)就会溢出。我写了一个不取模的版本把f(1)到f(92)打印出来和已知的数列对照结果一致说明矩阵本身没问题。这一步排除了“矩阵写错但取模后碰巧看着对”的可能。第二组是取模后的对照。取模运算有个性质如果x ≡ y (mod m)那(x z) ≡ (y z)。所以我用暴力版每一步都取模和矩阵版也是每步取模比对了n 10^5到10^6的若干个点结果一致。这一步主要验证的是取模时机如果有人在矩阵乘法里忘了取模就会在这一步暴露——不取模的话中间结果会迅速撑爆 long long 变成负数。边界值我整理了这么一组每次写完都过一遍n 值期望答案模 10^97考察点00边界特判是否生效11幂次为 0 时单位阵是否正确21矩阵自乘一次是否正确32是否漏掉某一项1055常规小值50586268941中值验证取模10^18需要在同模下验证大数、64 位溢出n 10^18这一行没有现成的期望答案我的做法是找一个同样实现了矩阵快速幂、但写法完全不同的版本比如 Python 版交叉验证。Python 的大整数天然不溢出用它跑一遍作为基准再用 C 对比两边一致就基本可以放心。4.3 性能实测与优化取舍我在本地机器上测了三组数据n分别取 10^9、10^12、10^18每种跑 1000 次取平均结果如下n 规模单次耗时循环版单次耗时展开版循环轮数10^9约 0.9 微秒约 0.4 微秒3010^12约 1.8 微秒约 0.8 微秒4010^18约 2.7 微秒约 1.2 微秒60单个用例微秒级对于蓝桥杯这种一次评测十几个测试点的场景总耗时连一毫秒都不到完全没有性能压力。所以在这个题上纠结展开版还是循环版意义不大可读性优先即可。真正需要优化的是多组查询场景如果一道题要回答 10^5 个 n 的斐波那契值每次 60 轮循环就是 6×10^6 次矩阵乘法耗时能到几十毫秒甚至上百毫秒这时候展开版的优势才体现出来。还有一种优化思路是预处理幂次把M^(2^k)预先算好存到表里每个查询只需要按二进制位查表相乘平方的那部分被摊掉了。60 个 2×2 矩阵只占几 KB 内存查询只需最多 60 次乘法而不是 120 次。这个技巧在“多次查询同一递推”的题里非常好用我在一次模拟赛中遇到过 10^5 次查询的情况用预处理把时间从 180ms 压到了 90ms 左右。5. 常见问题排查速查表5.1 答案错误类问题这类问题最典型的表现是样例能过交上去 WAWrong Answer或者本地小数据全对大数据错。我在自己踩过的坑和别人问过的问题里总结出下面几个高频原因。第一是幂次搞错。n和n-1混用是第一名。判断方法很简单把n 2代进去手工算如果代码返回的是M^2的左上角那算出来应该是f(3) 2而正确答案应该是f(2) 1差了一项。这种错误在n小于 3 的时候最容易暴露所以测试时一定要从n 0开始跑。第二是单位矩阵没有初始化。res是全零矩阵那么res mul(res, base)之后还是全零无论怎么算都是 0。判断方法输出中间结果如果mpow返回全零矩阵那一定是这个问题。第三是负数取模。如果有人用了(a - b) % MOD这种写法而没有加MOD修正可能会得到负数。虽然 C 里负数取模的结果是负数但在模意义下应该修正到非负。这个坑在快速倍增写法里特别常见写(2*b - a) % MOD的时候务必改成((2*b - a) % MOD MOD) % MOD。第四是模数写错。这个听起来很蠢但确实高频。题目说模 10^97你写成 10^99或者把1000000007数成了10000000007多了一个 0。写完代码后把常量单独抄一行对照题面读一遍只要五秒钟。第五是矩阵乘法展开时写错下标。手工展开 2×2 一共 8 项很容易把A.v[0][1] * B.v[1][0]和A.v[1][0] * B.v[0][1]搞混。展开写的这四个式子每一项都是“行 i 元素乘列 j 元素再求和”按这个规则核对一遍就不会错。5.2 超时与精度类问题超时的典型原因是把矩阵按值传递。Mat里面是四个long long按值传参要拷贝 32 字节60 轮循环乘几十次再加上函数调用开销累积起来并不小。改成const Mat引用传递后速度能提升明显。这个问题在小规模数据上感知不到但在多组查询的题里是致命的。第二个超时原因是每轮循环都做一次memset。在mul内部构造Mat C时构造函数里有一次 32 字节的清零本来不贵但如果有人在外层又套了一次memset或者把数组开得很大比如 100×100 的矩阵这个开销就放大了。对于 2×2 这种固定小尺寸其实可以不用memset直接把四个元素显式赋值更快也更清晰。精度类问题主要出现在用浮点数辅助验证的场景。有人为了对数的大小有个直观感受用double估算f(n)的对数log10(f(n)) ≈ n * log10(1.618) - 0.349这个公式在 n 很大的时候非常准可以用来判断结果有多少位。但如果想用它来验证具体数值一旦超过 15 位有效数字就没法看了只能作为估算工具不能作为正确答案的依据。我把这些问题整理成下面这张速查表出问题的时候按顺序扫一遍基本能覆盖九成情况现象可能原因排查方法输出全 0单位矩阵未初始化打印mpow的返回值结果差一项幂次用 n 还是 n-1 混用用 n2 手工验证小数据对、大数据错中间结果溢出检查模数是否超过 46340 而用了 int结果为负减法取模未修正搜索代码里所有% MOD前的减法TLE按值传参、重复 memset改成引用传递、去掉多余清零样例过、隐藏点错模数常量抄错逐字对照题面矩阵乘法结果错展开时下标写反按“行乘列求和”逐项核对5.3 一些不成文的实战心得第一模板代码要固化下来。我现在的做法是维护一个小的代码片段文件里面放着 2×2 矩阵结构体、mul、mpow三个东西遇到这类题直接粘过去改模数和初始矩阵。这样做的最大好处是这份模板是被大量题目验证过的不需要每次重新怀疑它。新写的代码出错原因大概率在题目特化的部分而不是模板本身。这个思路能把调试范围砍掉一大半。第二别急着优化。我第一次写这题的时候花了半小时研究怎么把矩阵乘法写成 4 次乘法用 Strassen 或者别的技巧最后发现完全没必要。2×2 规模下常数优化带来的收益是微秒级而写复杂了引入 bug 的风险是指数级。先把正确的版本写出来跑通再考虑要不要优化这个顺序不能反。第三取模的位置要统一。有人喜欢在乘法里不取模、最后统一取模这在只有几个数的时候可行但在 60 轮循环里中间结果会指数级增长早就溢出了。标准做法是每一次乘法和加法之后立刻取模让所有中间值始终保持在[0, MOD)区间内。这样虽然多了一些取模运算但正确性得到了保证而且取模本身的代价在现代 CPU 上并不高。第四用 Python 做交叉验证。C 溢出问题很难靠肉眼发现而 Python 的整数是任意精度写一个十行的矩阵快速幂版本跑几个大数作为基准答案再和 C 结果对比能快速定位是不是溢出导致的。这个方法我用了很多次比单步调试快得多。6. 变式扩展前 n 项和、带常数项与快速倍增6.1 求斐波那契前 n 项和这道每日一题的原始版本是求第 n 项但蓝桥杯里更常见的是求前 n 项和。直接的思路是维护额外一维状态把和也放进状态向量。状态选[f(n), f(n-1), S(n)]其中S(n) f(1) f(2) ... f(n)。转移关系是f(n1) f(n) f(n-1) f(n) f(n) S(n1) S(n) f(n1) S(n) f(n) f(n-1)对应的 3×3 矩阵就是[ 1 1 0 ] [ 1 0 0 ] [ 1 1 1 ]这个矩阵的第三行1 1 1表示新的和等于旧的和加新的项而新项又等于f(n) f(n-1)所以是1 1 1。代码上只需要把 2×2 的模板扩成 3×3循环改成三重其他结构完全不变。我实测了一下3×3 的矩阵乘法是 27 次乘法60 轮循环也就一千多次乘法耗时仍在微秒级性能完全不用担心。顺便提一个更巧的做法斐波那契前 n 项和有闭式S(n) f(n2) - 1。验证一下n 3时S(3) 0112 4而f(5) - 1 5 - 1 4对上了。用这个公式求前 n 项和就退化成求第n2项再减一代码一行都不用改。不过要注意取模后的减法要加MOD修正写成(fib(n2) - 1 MOD) % MOD。这个小技巧在填空题里能省不少时间但前提是你记得住这个恒等式记不住的话还是老老实实扩维更稳妥。6.2 带系数、带常数项的广义递推蓝桥杯的压轴题经常把斐波那契改造一下变成f(n) p * f(n-1) q * f(n-2) r这种形式其中p, q, r都是题目给定的常数。这时候状态要扩成三维[f(n), f(n-1), 1]矩阵是[ p q r ] [ 1 0 0 ] [ 0 0 1 ]第一行直接抄递推式的系数第二行1 0 0表示“旧的 f(n) 变成新的 f(n-1)”第三行0 0 1保持不变。这个构造方法可以推广到任意阶几阶递推矩阵的前几行就抄系数下面的行负责把状态往下挪位。如果递推式里出现n的一次项比如f(n) f(n-1) f(n-2) n那就需要把n和常数1都放进状态变成四维[f(n), f(n-1), n, 1]。转移时n1 n 1所以矩阵里对应n的那一行是0 0 1 1对应常数 1 的那一行是0 0 0 1。这类题的关键不是计算而是把状态设计对一旦状态向量选对了矩阵就是机械地抄系数不会出错。我在一次模拟赛里遇到过f(n) f(n-1) f(n-3)这种跳项递推当时第一反应是状态选[f(n), f(n-1), f(n-2), f(n-3)]四维矩阵结果发现多了一维。仔细想想f(n-3)可以由前三项推出来吗不能因为递推里没用到f(n-2)所以f(n-2)必须保留。正确的状态是[f(n), f(n-1), f(n-2)]三维矩阵是[[1,0,1],[1,0,0],[0,1,0]]。这个例子说明状态维数等于递推式的“记忆跨度”不是项数减一得具体分析。6.3 快速倍增法常数更小的替代方案矩阵快速幂不是唯一解。斐波那契还有一个专门的性质可以利用f(2k) f(k) * (2*f(k1) - f(k)) f(2k1) f(k)^2 f(k1)^2这两个式子把f(2k)和f(2k1)都表示成了f(k)和f(k1)的组合从低位往高位递推每轮只需要 3 次乘法和几次加减比矩阵乘法的 8 次乘法少一半还多。实现上从最高位往最低位扫维护一对(a, b) (f(k), f(k1))根据当前二进制位决定下一对是(f(2k), f(2k1))还是(f(2k1), f(2k2))。pairll, ll fib_pair(ll n) { // 返回 {f(n), f(n1)} ll a 0, b 1; // f(0), f(1) for (int i 62; i 0; --i) { ll c a * ((2 * b % MOD - a MOD) % MOD) % MOD; // f(2k) ll d (a * a % MOD b * b % MOD) % MOD; // f(2k1) if ((n i) 1LL) { a d; b (c d) % MOD; } else { a c; b d; } } return {a, b}; }这段代码有几个地方要小心。2 * b % MOD - a MOD这一串如果写成2 * b - a直接取模可能出负数必须先对2*b取模再加MOD修正。a * a和b * b都是接近 10^18 的量级加起来约 2×10^18还在long long范围内上限 9.22×10^18所以安全但如果模数再大一点就要用__int128。循环从第 62 位开始是因为long long最多 63 位10^18的二进制长度约 60 位从 62 开始扫能覆盖。我用对拍程序验证了这个版本从n 0到n 5000全部和矩阵版一致其中n 0时循环走完 63 次都取偶数分支a保持 0、b保持 1返回{0, 1}正确。性能上同样跑 10^18 规模快速倍增版大约比展开版矩阵乘法再快 30% 到 40%。不过对于单次查询来说这点差距可以忽略我选哪个主要看当时手边哪个模板更顺手——如果题目就是纯斐波那契用倍增版更短如果题目是广义递推那就只能老老实实写矩阵。6.4 指数超大时的处理思路最后一个变式值得单独说如果题目给的n是一个上千位的十进制字符串long long根本装不下怎么算思路是先把指数降下来。由于模数是质数p时M^(p-1)在某种意义上具有周期性严格说矩阵的幂在模意义下也有周期性周期与模数相关常见做法是把指数对p-1取模后再算幂。但这里有个陷阱直接对p-1取模可能让指数变成 0而实际指数不为 0需要额外判断原指数是否为 0并在为 0 时特殊处理。如果模数不是质数这套周期性推导就失效了得换思路比如直接把字符串转成一个“大整数 取模”的循环ll mod_pow_str(const string s, ll m) { // 计算 s % m ll r 0; for (char c : s) r (r * 10 (c - 0)) % m; return r; }但我个人对这类题的态度是先确认题目的模数和指数的关系再决定用哪种降幂方式。盲目套费马小定理是常见的翻车点尤其是模数给的不是质数、或者指数需要减一修正的时候。稳妥的做法是在小数据上把“降幂版”和“暴力版”对拍一遍确认无误再提交。这套东西我前后断断续续用了一个多月才真正理顺中间翻过的车基本都写在上面的速查表里了。如果你也在刷蓝桥杯每日一题建议把这三个模板2×2 矩阵、3×3 矩阵、快速倍增都存好遇到新题先判断能不能套能套就直接改模数和系数不能套再回头重新设计状态。这样积累下来看到“斐波那契 大范围 取模”这几个关键词脑子里应该能立刻跳出矩阵而不是先去开数组。