51nod 1355 斐波那契的最小公倍数 （数论+莫比乌斯反演）

题目描述

传送门

题目大意：斐波那契数列定义如下：

F(0) = 0 F(1) = 1

F(n) = F(n-1) + F(n-2)

给出n个正整数a1, a2,…… an，求对应的斐波那契数的最小公倍数，由于数字很大，输出Mod 1000000007的结果即可。

题解

首先需要知道斐波那契数列的一个性质gcd(f[a],f[b])=f[gcd(a,b)]

这个怎么证明？先证明gcd(f[a],f[a+1])=1

f[0]=1,f[1]=1,gcd(f[0],f[1])=1

利用更相减损，gcd(f[a],f[a+1])=gcd(f[a],f[a]+f[a−1])=gcd(f[a],f[a−1])，最终得到gcd(f[a],f[a+1])=gcd(f[0],f[1])

还需要用到一个式子f[n+m]=f[n+1]f[m]+f[n]f[m−1]

这道题有一种有趣的组合的证明方法（摘自小火车的论文）





gcd(f[a],f[b])=gcd(f[a],f[a]f[b−a+1]+f[a−1]f[b−a])=gcd(f[a],f[a−1]f[b−a])

因为gcd(f[a],f[a−1])=1，那么gcd(f[a],f[b])=gcd(f[a],f[b−a]),可以递归，就是辗转相除嘛。那么上面的问题得证。

lcm(f[a],f[b])=f[a]f[b]f[gcd(a,b)]

设读入的下标构成的集合为S，那么lcm=∏T∈S,T≠∅f[gcd(T)](−1)|T|−1,其中gcd(T)表示的是集合T中所有数的最大公约数

然后考虑枚举gcd(T)=i

lcm=∏i=1maxf[i]∑T∈S,T≠∅,gcd(T)=i(−1)|T|−1

式子的瓶颈在于gcd(T)=i,我们可以利用莫比乌斯反演进行转化。

设g(i)=∑T∈S,T≠∅,gcd(T)=i(−1)|T|−1,h(i)=∑T∈S,T≠∅,i|gcd(T)(−1)|T|−1

T≠∅比较讨厌，我们发现空集比较讨厌，考虑直接提出来。然后再把指数中的−1提出来，式子可以化简成h(i)=1−∑T∈S,i|gcd(T)(−1)|T|

∑T∈S,i|gcd(T)(−1)|T|=∑|S|i=0(−1)i∗C(|S|,i),C表示组合数

其实就是杨辉三角中某一行的偶数列的和-奇数列的和，在|S|>0是式子的值恒等于0

所以h(i)=[下标中有i的倍数]

根据莫比乌斯反演g(n)=∑d|nh(d)∗μ(dn)

那么lcm=∏i=1maxf[i]g[i]可以用O(nlogn)的时间预处理g[i]数组，然后再用快速幂求解。

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 1000003
#define p 1000000007
#define LL long long
using namespace std;
int mu
,a
,cnt
,h
,g
,n,pd
,prime
,f
;
LL quickpow(LL num,int x)
{
LL ans=1,base=num%p;
x=(x%(p-1)+p-1)%(p-1);
while (x) {
if (x&1) ans=ans*base%p;
x>>=1;
base=base*base%p;
}
return ans;
}
void init(int n)
{
mu[1]=1;
for (int i=2;i<=n;i++) {
if (!pd[i]) {
prime[++prime[0]]=i;
mu[i]=-1;
}
for (int j=1;j<=prime[0];j++) {
if (prime[j]*i>n) break;
pd[i*prime[j]]=1;
if (i%prime[j]==0) break;
mu[i*prime[j]]=-mu[i];
}
}
f[0]=0; f[1]=1;
for (int i=2;i<=n;i++) f[i]=(f[i-1]+f[i-2])%p;
for (int i=1;i<=n;i++) {
if (!h[i]) continue;
for (int j=i;j<=n;j+=i) g[i]+=h[j]*mu[j/i];
}
}
int main()
{
freopen("a.in","r",stdin);
scanf("%d",&n);
int mx=0;
for (int i=1;i<=n;i++) scanf("%d",&a[i]),mx=max(a[i],mx);
for (int i=1;i<=n;i++) {
for (int x=1;x*x<=a[i];x++)
if (a[i]%x==0) h[x]=1,h[a[i]/x]=1;
}
init(mx);
LL ans=1;
for (int i=1;i<=mx;i++)
ans=ans*quickpow(f[i],g[i])%p;
printf("%lld\n",ans);
}