刚开始看到“斐波那契公约数”这个标题时，我以为是道纯数论题：给两个斐波那契数求个 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 和 0，gcd(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_{k+1}) = 1。

证明很直接。利用 F_{k+1} - F_k = F_{k-1}，可以得到 gcd(F_k, F_{k+1}) = gcd(F_k, F_{k-1})。不断把这个式子递归应用下去，最终会落在 gcd(F_1, F_0) 上。如果约定 F_1 = F_2 = 1，F_0 = 0，那么 gcd(1, 0) = 1，引理成立。

第二个引理：加法恒等式 F_{m+n} = F_{m+1}F_n + F_mF_{n-1}。

这个等式可以通过矩阵乘法快速得到：把斐波那契递推写成矩阵 A = [[1,1],[1,0]]，那么 A^{m+n} = A^m × A^n，取结果的某个分量就能推出上述恒等式。当然也可以用数学归纳法硬证，但矩阵视角更直观，也顺便为后面的矩阵快速幂做了铺垫。

第三个引理：若 m 整除 n，则 F_m 整除 F_n。

设 n = qm，对 q 做归纳。q = 1 时显然成立。若 F_m | F_{qm} 已成立，用加法恒等式展开 F_{(q+1)m} = F_{qm+m} = F_{qm+1}F_m + F_{qm}F_{m-1}，右边两项都含有 F_m 因子，所以 F_m | F_{(q+1)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_n：

F_n = F_{qm+r} = F_{qm+1}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_{qm+1}F_r, F_m)

又因为 F_{qm+1} 与 F_{qm} 互质（这是第一个引理），而 F_m 整除 F_{qm}，所以 F_{qm+1} 与 F_m 必然互质。把 F_{qm+1} 从 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 等于 1，F_1 = 1，对应答案就是 1。这些边界在代码里最好显式处理，虽然矩阵快速幂在 k=0 时返回单位矩阵，取第 0 行第 1 列恰好也会得到 0，但写清楚总归更稳妥。

3. 用矩阵快速幂求第 k 个斐波那契数

3.1 矩阵构造的核心思路

现在问题变成了：给定 k = gcd(n, m)，求 F_k mod 100000000。

矩阵快速幂可以做到 O(log k)。原理是把递推关系写成矩阵形式：

[F_{n+1}] [1 1] [F_n] [F_n ] = [1 0] [F_{n-1}]

所以如果定义 A = [[1,1],[1,0]]，那么 A^n 乘向量 [F_1, F_0]^T 后，得到的向量是 [F_{n+1}, F_n]^T。也就是说，A^0 的第 0 行第 1 列，或者第 1 行第 0 列，就是 F_n。

为什么常用 A = [[1,1],[1,0]] 而不是别的形式？因为这个矩阵恰好把“一步递推”翻译成了“一次线性变换”。第一个分量是新的 F_{n+1}，第二个分量是旧的 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_{k+1})：

F_{2k} = F_k × (2 × F_{k+1} - F_k) F_{2k+1} = F_k² + F_{k+1}²

利用这两个公式，可以按 k 的二进制拆分递归降低指数。代码也很短：

pair<ll, 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_{k+1}。递归到 k=0 时返回 (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 = 55，F_15 = 610，真实 gcd 是 5。如果先对 12 取模，55 mod 12 = 7，610 mod 12 = 10，gcd(7, 10) = 1，完全对不上。

这说明一个道理：不要试图绕过数学性质去“模拟”题面。竞赛题给出大范围的数据，往往就是暗示你要找到结构上的捷径，而不是硬算。

4.3 一组可以手动验证的边界数据

我建议写完代码后先用小数据做对照，下面这组是我当时拿来自测用的，全部都是手工可算的范围。

其中 12 和 18 这一组值得多说一句：F_12 = 144，F_18 = 2584，gcd(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 性能对比

我们用一张表直观对比几种常见做法的复杂度：

实际运行中，k 即使是 10^9 量级，二进制位数也不过 30，2×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) 去用，建议花十分钟把这个带余除法的证明过程亲手推一遍。推完之后你会发现，矩阵快速幂只是收尾工具，而那个漂亮的恒等式才是这道题真正的灵魂。