斐波那契公约数

题目描述

对于Febonacci数列,求第n项和第m项的最大公约数是多少?(n,m<=1e9>)。对于最后的结果只要输出最后的8位数字就可以了。

思路

首先注意到的是数据量巨大(1e9),数组开不下;即使使用Febonacci递推公式也会超时。不过好在Febonacci数列是有通项公式的。因为$F[n] = F[n-2] + F[n-1]$,解特征方程$X^{2}=X+1$得到解$X=(1 \pm \sqrt{5})/2$,再代入$F[n]=aX_{1}^{n} + bX_{2}^{n}$,就可以解的最终的通项公式是$F[n] = \frac{1}{\sqrt{5}} × [(\frac{1+\sqrt{5}}{2})^n-(\frac{1-\sqrt{5}}{2})^n]$。

求两个数的最大公约数可以用 辗转相除法

1
2
3
4
typedef long long LL;
LL gcd(LL a, LL b) {
return (b ? gcd(b, a % b) : a);
}

不过这里的数字会非常巨大,并且$\sqrt{5}$也不是一个整数,用快速幂的时候对double做%运算似乎会报错……

题解

正解用到了 矩阵加速一个Febonacci的结论
矩阵加速并不是必须的,但是能减少不少计算量。先看看矩阵加速的巧妙之处:

矩阵加速

如果已知:$(F[2], F[1])^T$,那么要求$(F[3], F[2])^T$,只需要乘一个矩阵:

其实这里已经能看出用矩阵的方式表达Febonacci数列的递推关系。
如果初始化矩阵为 $(F[2], F[1])^T$,那么通项公式为:

这样做就能使用 矩阵快速幂 减小复杂度到 O(logn) 算得 F[n]

再看一下这道题的 关键之处

一个Febonacci的结论

先看结论:

这个结论似乎很好记,但是不是很直接,Febonacci数列的第n项和第m项的最大公约数居然正好等于第$gcd(n,m)$项的值?

详细的推导过程是这样的:

  1. 如果$F[n]=a, F[n+1]=b$
  2. 那么$F[n+2]=a+b$, $F[n+3] = a+2b$, … ,$F[m] = F[m−n−1]a + F[m−n]b$
  3. 所以有$F[m] = F[m−n−1] ∗ F[n] + F[m−n] ∗ F[n+1]$
  4. 这样,$gcd(F[n], F[m]) = gcd(F[n]$, $F[m−n−1] ∗ F[n] + F[m−n] ∗ F[n+1])$
  5. 其中$F[m−n−1] ∗ F[n]$是$F[n]$的倍数,因此上式转为$gcd(F[n], F[m]) = gcd(F[n], F[m−n] ∗ F[n+1])$
  6. 可以通过递推证明:$gcd(F[n], F[n+1]) = 1$
  7. 因此$gcd(F[n], F[m]) = gcd(F[n], F[m−n])$,即$gcd(F[n], F[m]) = gcd(F[n], F[m%n])$

这样求解的目标就化简为求第$gcd(n,m)$项了。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include <cstring>
#include <iostream>
#define MOD 100000000

using namespace std;
typedef long long LL;

LL gcd(LL a, LL b) {
return (b ? gcd(b, a%b) : a);
}

struct Mat {
LL mat[2][2];
int h, w;
Mat(int _h=2, int _w=2) {
h = _h, w = _w;
memset(mat, 0, sizeof(mat));
}
};

Mat Multiply(const Mat &a, const Mat &b) {
Mat ret(a.h, b.w);
for (int i = 0; i < a.h; ++i)
for (int j = 0; j < b.w; ++j)
for (int k = 0; k < a.w; ++k) {
LL tmp = a.mat[i][k] * b.mat[k][j] % MOD;
ret.mat[i][j] = (ret.mat[i][j] + tmp) % MOD;
}
return ret;
}

int main() {
LL n, m;
cin >> n >> m;
LL gcd_nm = gcd(n, m);
if (gcd_nm <= 2) cout << 1 << endl;
else {
Mat raw(2, 1);
raw.mat[0][0] = raw.mat[1][0] = 1;
Mat f(2, 2);
f.mat[0][0] = f.mat[0][1] = f.mat[1][0] = 1;
gcd_nm -= 2;
while (gcd_nm) {
if (gcd_nm & 1)
raw = Multiply(f, raw);
f = Multiply(f, f);
gcd_nm >>= 1;
}
cout << raw.mat[0][0] << endl;
}
return 0;
}