书海扬帆的博客

题目描述

对于$Fibonacci$数列:$1,1,2,3,5,8,13……$大家应该很熟悉吧~~~
但是现在有一个很“简单”问题:第$n$项和第$m$项的最大公约数是多少?

输入输出格式

输入格式:

两个正整数$n$和$m$。$(n,m<=10^9)$
注意:数据很大

输出格式:

${F_n}$和${F_m}$的最大公约数。
由于看了大数字就头晕,所以只要输出最后的8位数字就可以了。

输入输出样例

输入样例#1:
4 7
输出样例#1:
1

说明

用递归&递推会超时
用通项公式也会超时


解题思路

起先看到这个题,递推不行,通项也不行,那只能使用矩阵加速来解决了。(矩阵加速是啥?)
那么怎样运用矩阵加速求解呢?显然,一个比较朴素的方法就是直接算出${F_n}$和${F_m}$。但是考虑到数据范围是$(n,m<=10^9)$,因此该算法会超时!
正在冥思苦想之时,看到一个有趣的结论:

在$Fibonacci$数列中,对于任意两项${F_n}$和${F_m}$,均有$gcd(F_n,F_m)={F_{gcd(n,m)}}$成立

接下来,我们只需证明该结论的正确性。

证明

若$n=m$,显然成立,下证$n≠m$时的情况.
不失一般性,
不妨设n>m.

引理1

【证明】
由辗转相除法可得,

$gcd(f[n+1],f[n])=gcd(f[n+1]−f[n],f[n])$

$=gcd(f[n−1],f[n])$

$=gcd(f[n−1],f[n]−f[n−1])$

$=gcd(f[n−1],f[n−2])$

$=……$

$=gcd(f[1],f[2])=gcd(1,1)=1.$

引理2

【证明】

$f[m+n]=f[m+n−1]+f[m+n−2]$

$=2∗f[m+n−2]+f[m+n−3]$

$=……$

设$f[m+n]=a[x]∗f[m+n−x]+b[x]∗f[m+n−x−1]$

$=a[x]∗(f[m+n−x−1]+f[m+n−x−2])+b[x]∗f[m+n−x−1]$

$=(a[x]+b[x])∗f[m+n−x−1]+a[x]∗f[m+n−x−2]$

当$x=1$时,$f[m+n]=a[1]∗f[m+n−1]+b[1]∗f[m+n−1−1]$,

有$a[1]=f[1]=f[2]=1,b[1]=f[1]=1$;

当$x=2$时,$f[m+n]=a[2]∗f[m+n−2]+b[2]∗f[m+n−2−1]$,

有$a[2]=f[1]+f[2]=f[3]=2,b[2]=a[1]=f[2]=1$;

假设当$x=k$时,$a[k]=f[k+1],b[k]=a[k−1]=f[k]$,

那么当$x=k+1$时,$a[k+1]=a[k]+b[k]=f[k+1]+f[k]=f[k+2]b[k+1]=a[k]=f[k+1]$;

所以当$x=n$时有

$f[m+n]=a[n]f[m]+b[n]f[n−1]$

$=f[m−1]f[n]+f[m]f[n+1]$.

证毕。

引理3

【证明】

由引理2得

$gcd(f[m+n],f[n])=gcd(f[m−1]f[n]+f[m]f[n+1],f[n])$

$=gcd(f[m]f[n+1],f[n])$

$=gcd(f[m],f[n])∗gcd(f[n+1,f[n])$

由引理1,

gcd(f[n+1],f[n])=1,

因此原式

$=gcd(f[m],f[n]).$

证毕。

由引理3证明结论

由引理3可得,

$gcd(f[m],f[n])=gcd(f[m%n],f[n]),$

记$m1=m%n$,则

$gcd(f[m],f[n])=gcd(f[m1],f[n%m1])$

……

$=gcd(0,f[gcd(m,n)])$

我们已知用计算机求解最大公约数的算法如下:

1
2
3
4
int gcd(int a,int b)
{
return b==0?a:gcd(b,a%b);
}

不难看出,在递归过程中,

证毕。

运用结论结合矩阵加速求解

斐波那契数列的初始矩阵为
单位矩阵为

对于第$gcd(n,m)$项,只需快速幂$gcd(n,m)−2$次即可,结果$\mod 1e8$.

代码

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
53
54
55
56
57
58
59
60
61
62
63
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std;
typedef long long ll;
ll n,m;
ll gcd(ll a,ll b)
{
if(b==0) return a;
return gcd(b,a%b);
}
const ll mod=1e8;
struct matrix
{
ll n,m;
ll z[105][105];
matrix()
{
n=m=0;
memset(z,0,sizeof(z));
}
}m1,e;
matrix operator*(const matrix &m1,const matrix &m2) {
matrix m3;
int n=m1.n,m=m1.m,k=m2.m;
m3.n=n;m3.m=k;
for (int a=1;a<=n;a++)
for (int b=1;b<=m;b++)
for (int c=1;c<=k;c++)
m3.z[a][c] = m3.z[a][c]%mod+m1.z[a][b]*m2.z[b][c]%mod;
return m3;
}
inline matrix Pow(matrix &m1, ll k)
{
matrix ans=e;
while(k)
{
if(k&1) ans=ans*m1;
m1=m1*m1;
k>>=1;
}
return ans;
}
int main()
{
scanf("%lld%lld",&n,&m);
ll cur=gcd(n,m);
if(cur==1 || cur==2)
{
printf("1");
return 0;
}
e.n=e.m=2;
e.z[1][1]=e.z[1][2]=e.z[2][1]=1;
m1.n=1,m1.m=2;
m1.z[1][1]=m1.z[1][2]=1;
//给定A1 A2 乘n次得到第N+2项
matrix ans=Pow(e,cur-3);
ans=m1*ans;
printf("%lld",ans.z[1][1]%mod);
return 0;
}


本站使用 Material-X 作为主题 , 总访问量为 次 。
载入天数...载入时分秒... 字数统计:725.9k