资讯动态

斐波那契公约数:从gcd性质到矩阵快速幂的算法简析

发布时间:2026/10/17 18:29:21 来源:尧图企业网站定制
刚开始看到“斐波那契公约数”这个标题时我以为是道纯数论题给两个斐波那契数求个 gcd 不就行了但真正动手之后才发现n 和 m 最大到 10^9别说存下完整的斐波那契数列就是拿高精度都顶不住。这道题表面上是“公约数”实际上是一道“先做数学化简再做快速幂”的组合题而那个化简性质就是斐波那契数列里非常经典的 gcd 性质。这篇文章我把自己从读题、推到矩阵快速幂、再到提交排坑的完整过程整理出来适合已经有矩阵快速幂基础、但又好奇“为什么能把下标直接 gcd”的算法竞赛选手。1. 题目表面是“公约数”实际考的是“把下标做 gcd”1.1 直接求 F[n] 和 F[m] 会撞上几堵墙先明确题目要什么给定 n、m求 gcd(F[n], F[m]) 对 100000000 取模。第一时间想到的做法是分别算出 F[n] 和 F[m]再求最大公约数。但问题来了。第一堵墙是数值规模。斐波那契数列是指数增长的n 到 10^9 时F[n] 的位数大概是 n 乘以 log10(黄金比例)约等于 0.208987 乘以 n也就是超过两亿位十进制数字。这种规模的数字别说普通程序跑不动连“输出”都是一件灾难性的事。第二堵墙是“边算边取模再 gcd”行不通。很多人会想那我分别算 F[n] mod 100000000 和 F[m] mod 100000000然后对两个余数求 gcd答案不就是题目要的吗这个思路是错的因为取模运算和最大公约数运算不可交换。举个简单例子gcd(6, 10) 2但两个数都对 5 取模后得到 1 和 0gcd(1, 0) 1和原来完全不同。第三堵墙是即使你想用快速幂去算第 n 个斐波那契数的真实值你也根本做不到因为快速幂得到的是模某个数之后的余数而不是完整整数。所以题目给的模数 100000000 不是在帮你“简化输出”而是在提示你答案本身可以只保留模意义下的结果但你必须找到一个不依赖完整数值的解法。1.2 数据范围给出的路线图n, m 最大到 10^9意味着任何 O(n) 或 O(m) 的递推都不可行。哪怕 n 和 m 是互质的最坏也要线性扫到接近 10^9。按普通评测机一秒能跑大约 10^8 次简单运算来估算纯递推算 F[10^9] 肯定超时。所以题目里真正的计算瓶颈不是 gcd而是“求第 k 个斐波那契数”这个子问题。gcd(n, m) 本身的复杂度是 O(log min(n, m))连塞牙缝都不够。换句话说只要能把题目化简成“求某个下标对应的斐波那契数”剩下的事就是快速幂干的事。顺着这个方向唯一合理的路线就是先证明一个数学性质把 gcd(F[n], F[m]) 转化成 F[gcd(n, m)]然后利用矩阵快速幂或快速双倍公式在 O(log k) 时间内求出第 k 项。整道题的代码量并不大难的是你是否能想到这个性质以及能不能把证明看透。1.3 模数 100000000 的特殊性这道题的模数是 100000000也就是 10^8分解质因数是 2^8 × 5^8。它既不是质数也不是常见的 1e9 7。这个细节很值得停下来想一下。如果某个解法里需要用到费马小定理求逆元比如用组合数公式或通项公式去算斐波那契数那就必须要求模数是质数。100000000 显然不满足。但矩阵快速幂和快速双倍公式从头到尾只用了加法和乘法不涉及除法更不涉及求逆元所以即使模数是合数也能正常工作。另外模数决定了乘法过程中数值范围的上限。两个小于 10^8 的数相乘结果小于 10^16远在 C 的 long long 可表示范围约 9.22 × 10^18之内。因此在这个模数下用 long long 做中间变量是安全的。2. gcd(F_n, F_m) F_gcd(n,m)这道题最值得背的性质2.1 三个引理要完整理解 gcd(F_n, F_m) F_gcd(n,m) 这个结论先铺三个引理。第一个引理相邻斐波那契数互质即 gcd(F_k, F_{k1}) 1。证明很直接。利用 F_{k1} - F_k F_{k-1}可以得到 gcd(F_k, F_{k1}) gcd(F_k, F_{k-1})。不断把这个式子递归应用下去最终会落在 gcd(F_1, F_0) 上。如果约定 F_1 F_2 1F_0 0那么 gcd(1, 0) 1引理成立。第二个引理加法恒等式 F_{mn} F_{m1}F_n F_mF_{n-1}。这个等式可以通过矩阵乘法快速得到把斐波那契递推写成矩阵 A [[1,1],[1,0]]那么 A^{mn} A^m × A^n取结果的某个分量就能推出上述恒等式。当然也可以用数学归纳法硬证但矩阵视角更直观也顺便为后面的矩阵快速幂做了铺垫。第三个引理若 m 整除 n则 F_m 整除 F_n。设 n qm对 q 做归纳。q 1 时显然成立。若 F_m | F_{qm} 已成立用加法恒等式展开 F_{(q1)m} F_{qmm} F_{qm1}F_m F_{qm}F_{m-1}右边两项都含有 F_m 因子所以 F_m | F_{(q1)m}。归纳完成。2.2 核心证明一次带余除法化归有了三个引理主性质就不难证了。设 n qm r其中 0 ≤ r m。分两种情况。如果 r 0说明 m | n由第三个引理直接得到 F_m | F_n因此 gcd(F_n, F_m) F_m F_gcd(n,m)结论成立。如果 r 0用加法恒等式展开 F_nF_n F_{qmr} F_{qm1}F_r F_{qm}F_{r-1}对 F_m 取最大公约数时注意 F_{qm} 能被 F_m 整除所以 F_{qm}F_{r-1} 这一整项在模 F_m 的意义下等于 0。于是gcd(F_n, F_m) gcd(F_{qm1}F_r, F_m)又因为 F_{qm1} 与 F_{qm} 互质这是第一个引理而 F_m 整除 F_{qm}所以 F_{qm1} 与 F_m 必然互质。把 F_{qm1} 从 gcd 里提出来不影响结果于是gcd(F_n, F_m) gcd(F_r, F_m)这一步非常关键它告诉我们对斐波那契数列的下标做“带余除法”之后gcd 的结果也随之缩小成 gcd(F_r, F_m)。这正是整数欧几里得算法的同款过程。不断对下标对做辗转相除最终下标对会变成 (d, 0) 或 (0, d)其中 d gcd(n, m)。最后一步gcd(F_d, F_0) F_d因为 F_0 0而任何数与 0 的最大公约数等于它本身。整个证明到此闭环。2.3 下标约定和边界情况证明里依赖 F_0 0、F_1 1 这个下标约定。洛谷这一系列题通常习惯 F_1 F_2 1于是 F_0 0 是自然补出来的。当 gcd(n, m) 的结果特别小的时候要注意如果 gcd 等于 0第 0 项斐波那契数按定义就是 0。如果 gcd 等于 1F_1 1对应答案就是 1。这些边界在代码里最好显式处理虽然矩阵快速幂在 k0 时返回单位矩阵取第 0 行第 1 列恰好也会得到 0但写清楚总归更稳妥。3. 用矩阵快速幂求第 k 个斐波那契数3.1 矩阵构造的核心思路现在问题变成了给定 k gcd(n, m)求 F_k mod 100000000。矩阵快速幂可以做到 O(log k)。原理是把递推关系写成矩阵形式[F_{n1}] [1 1] [F_n] [F_n ] [1 0] [F_{n-1}]所以如果定义 A [[1,1],[1,0]]那么 A^n 乘向量 [F_1, F_0]^T 后得到的向量是 [F_{n1}, F_n]^T。也就是说A^0 的第 0 行第 1 列或者第 1 行第 0 列就是 F_n。为什么常用 A [[1,1],[1,0]] 而不是别的形式因为这个矩阵恰好把“一步递推”翻译成了“一次线性变换”。第一个分量是新的 F_{n1}第二个分量是旧的 F_n矩阵乘法天然保留了递推所需的两个相邻项信息。理解了这一点后面就容不容易把矩阵方向写反。3.2 C 实现模板下面给出我实际提交用的版本主要部分是 2×2 矩阵的乘法与快速幂。因为矩阵只有 2×2没必要用 vector 套 vector直接写成结构体数组即可。#include bits/stdc.h using namespace std; using ll long long; const ll MOD 100000000LL; struct Mat { ll a[2][2]; }; Mat mul(const Mat A, const Mat B, ll mod) { Mat C; C.a[0][0] C.a[0][1] 0; C.a[1][0] C.a[1][1] 0; for (int i 0; i 2; i) { for (int k 0; k 2; k) { if (A.a[i][k] 0) continue; for (int j 0; j 2; j) { C.a[i][j] (C.a[i][j] A.a[i][k] * B.a[k][j]) % mod; } } } return C; } Mat mpow(Mat base, ll p, ll mod) { Mat res{}; res.a[0][0] res.a[1][1] 1; res.a[0][1] res.a[1][0] 0; while (p 0) { if (p 1) res mul(res, base, mod); base mul(base, base, mod); p 1; } return res; } ll fib(ll k) { if (k 0) return 0; Mat base{}; base.a[0][0] 1; base.a[0][1] 1; base.a[1][0] 1; base.a[1][1] 0; Mat P mpow(base, k, MOD); return P.a[0][1] % MOD; } int main() { ll n, m; cin n m; ll g gcd(n, m); cout fib(g) \n; return 0; }这里有几个容易写错的地方一是单位矩阵必须初始化为左上角和右下角为 1如果忘了初始化 res 或全初始化为 0任何矩阵乘完都变成零矩阵结果恒为 0二是乘法里要记得取模虽然 2×2 矩阵很小但中间乘积可能超过 int 范围必须用 long long三是在乘法循环里可以加一个 if 判断跳过为 0 的元素这是个小优化不影响正确性。3.3 不建矩阵的快速双倍公式如果不想背矩阵还有一套快速双倍公式只需要维护一对相邻项 (F_k, F_{k1})F_{2k} F_k × (2 × F_{k1} - F_k) F_{2k1} F_k² F_{k1}²利用这两个公式可以按 k 的二进制拆分递归降低指数。代码也很短pairll, ll fib2(ll k) { if (k 0) return {0, 1}; auto [f, fnext] fib2(k / 2); ll c f * ((2 * fnext % MOD - f MOD) % MOD) % MOD; ll d (f * f % MOD fnext * fnext % MOD) % MOD; if (k 1) return {d, (c d) % MOD}; else return {c, d}; }这个版本返回的 pair 中第一个元素是 F_k第二个元素是 F_{k1}。递归到 k0 时返回 (F_0, F_1) (0, 1)。当 k 是偶数时c 就是 F_k当 k 是奇数时需要额外加上 F_{k-1} 才能凑出 F_k。两种方法时间复杂度都是 O(log k)。矩阵法写起来规律性强快速双倍法常数小一点我个人的建议是至少把其中一种练熟另一种留着互相验证。4. 提交过程踩过的坑与边界数据复盘4.1 模数抄错的教训我第一次写这道题时顺手在全局定义 const ll MOD 1000000007。当时还觉得自己写得很顺畅直到提交后拿到一个 WA才意识到题目要的模数是 100000000少了一个 7差得很远。这个错误非常隐蔽因为程序本身逻辑、边界全对只是常量写错了。后来我会在敲代码之前先把题面上的模数单独一行注释出来而不是靠记忆。4.2 取模后的余数直接做 gcd 为什么错前文已经说过取模和 gcd 不可交换。我在测试阶段专门构造过反例F_10 55F_15 610真实 gcd 是 5。如果先对 12 取模55 mod 12 7610 mod 12 10gcd(7, 10) 1完全对不上。这说明一个道理不要试图绕过数学性质去“模拟”题面。竞赛题给出大范围的数据往往就是暗示你要找到结构上的捷径而不是硬算。4.3 一组可以手动验证的边界数据我建议写完代码后先用小数据做对照下面这组是我当时拿来自测用的全部都是手工可算的范围。nmgcd(n,m)正确答案371F_1 1693F_3 210155F_5 512186F_6 88124F_4 3110000000001F_1 1其中 12 和 18 这一组值得多说一句F_12 144F_18 2584gcd(144, 2584) 8正好等于 F_6因为 gcd(12, 18) 6。如果你直接手算验证会发现这个组合非常有说服力。4.4 溢出和负数处理的细节模数 100000000 不算大但快速双倍公式里的中间表达式 2 * fnext - f 可能会产生负数。所以要写成 2 * fnext % MOD - f MOD 然后再整体取模确保结果非负。在矩阵乘法里A.a[i][k] * B.a[k][j] 最大约 10^16远没到 long long 上限所以这里用 long long 就足够。但如果哪天真遇到更大的模数比如 10^18 级别就得准备 __int128 或者做乘法拆分这是后话。4.5 性能对比我们用一张表直观对比几种常见做法的复杂度做法时间复杂度能否通过高精度暴力递推 F[n]O(n) 且数值位数爆炸不可能边算边取模递推 F[n] 和 F[m]O(n)10^9 超时先算 F[n] mod MOD 再 gcdO(n)10^9 超时且思路错误数学化简 矩阵快速幂O(log gcd(n,m))轻松通过数学化简 快速双倍公式O(log gcd(n,m))轻松通过实际运行中k 即使是 10^9 量级二进制位数也不过 302×2 矩阵快速幂只需要几十次矩阵乘法。评测机上单组数据跑完基本是 0 毫秒级别瓶颈反而在读入输出和 gcd 的常数开销。4.6 一道题能带出来的思维习惯做完 P1306 之后我最大的收获不是背熟了一个矩阵快速幂模板而是记住了“gcd 下标化”这个套路。以后只要看到形如 gcd(f(n), f(m)) 的结构第一步永远是先考虑 f 这个函数是否存在类似“f 的 gcd 等于下标 gcd 的 f”的性质。斐波那契有这样的性质很多线性递推数列也有只是系数条件不同。做竞赛题最常见的卡题原因往往不是不会写代码而是脑袋里少了那条通往性质的路。我自己在实际操作中的体会是写题解和写代码不一样。代码只需要最终正确但真正把一个性质吃透还是要回到证明本身。如果你也曾经从别人题解里直接抄下 gcd(F_n, F_m) F_gcd(n,m) 去用建议花十分钟把这个带余除法的证明过程亲手推一遍。推完之后你会发现矩阵快速幂只是收尾工具而那个漂亮的恒等式才是这道题真正的灵魂。

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑