[莫比乌斯反演] BZOJ4816: [Sdoi2017]数字表格

挺简单的。。。

∏i=1n∏j=1mfib[gcd(i,j)]=∏dfib[d]∑ni∑mj[gcd(i,j)==d]

上面这个指数是经典的反演的问题，然后就变成：

∏dfib[d]∑kμ(k)⌊ndk⌋⌊mdk⌋

套路变形一下

∏T=1(∏d|Tfib[d]μ(Td))⌊nT⌋⌊mT⌋

中间这个f(T)=∏d|Tfib[d]μ(Td)，直接暴力 nlogn 算。

∏T=1f(T)⌊nT⌋⌊mT⌋

分块求就好了。

这题我跑的很慢，于是网上看了一下发现一个小trick:(其实是我太菜之前不知道)

求一个数列的每个数的逆元和前缀积的逆元，可以先算 S−1n，然后O(n) 推。

S−1i=S−1i+1∗ai+1 a−1i=S−1i∗Si−1

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long LL;
const int maxn=1000005, N=1000000, MOD=1e9+7;
int p[maxn],mu[maxn];
bool vis[maxn];
void get_mu(){
mu[1]=1;
for(int i=2;i<=N;i++){
if(!vis[i]) p[++p[0]]=i, mu[i]=-1;
for(int j=1;j<=p[0]&&(LL)i*p[j]<=N;j++){
vis[i*p[j]]=true;
if(i%p[j]==0){ mu[i*p[j]]=0; break; }
mu[i*p[j]]=-mu[i];
}
}
}
LL Pow(LL a,LL b){
if(b==-1) return Pow(a,MOD-2);
LL res=1;
for(;b;a=a*a%MOD,b>>=1) if(b&1) res=(res*a)%MOD;
return res;
}
int _test,n,m,fib[maxn],fib_n[maxn],f[maxn],f_n[maxn];
LL ans;
int f_mul(int L,int R){ return (LL)f[R]*f_n[L-1]%MOD; }
int main(){
freopen("bzoj4816.in","r",stdin);
freopen("bzoj4816.out","w",stdout);
get_mu();
fib[0]=0; fib[1]=fib_n[1]=1; for(int i=2;i<=N;i++) fib[i]=(fib[i-1]+fib[i-2])%MOD, fib_n[i]=Pow(fib[i],MOD-2);
for(int i=1;i<=N;i++) f[i]=1;
for(int i=1;i<=N;i++)
for(int j=1;(LL)i*j<=N;j++){
int t=(mu[j]==0?1:(mu[j]==1?fib[i]:fib_n[i]));
f[i*j]=((LL)f[i*j]*t)%MOD;
}
f[0]=f_n[0]=1; for(int i=1;i<=N;i++) f[i]=((LL)f[i]*f[i-1])%MOD, f_n[i]=Pow(f[i],MOD-2);
scanf("%d",&_test);
while(_test--){
scanf("%d%d",&n,&m); int _min=min(n,m);
ans=1;
for(int i=1,nxt;i<=_min;i=nxt+1){
nxt=min(n/(n/i),m/(m/i));
ans=(ans*Pow(f_mul(i,nxt),(LL)(n/i)*(m/i)))%MOD;
}
printf("%lld\n",ans);
}
return 0;
}