[BZOJ 4816]数字表格 莫比乌斯反演

第一次见到这种形式,推了两步就被卡住了。没想到还能这么迁移QAQ,真强!

题目设 $f(i)$为第$i$项斐波那契数列

要求计算

$$ans=\prod_{i=1}^n \prod_{j=1}^m f(gcd(i,j))$$

枚举$d=gcd$得到

$$ans=\prod_{d=1}^n f(d)^{\sum_{i=1}^n \sum_{j=1}^m[gcd(i,j)==d]}$$

然后反演得到

$$ans=\prod_{d=1}^n f(d)^{\sum_{T=1}^{\lfloor \frac{n}{d} \rfloor}\mu(T) \lfloor \frac{n}{dT} \rfloor \lfloor \frac{m}{dT} \rfloor}$$

然后呢,令$p=dT$,得到

$$ans=\left( \prod_{d=1}^nf(d)^{\mu(\frac{p}{d})} \right) ^{\sum_{p=1}^n  \lfloor \frac{n}{p} \rfloor \lfloor \frac{m}{p} \rfloor}$$

这样我们就发现我们可以对$p$进行分块去搞了

$$g(p)=\prod_{d|p} f(d)^{\mu(\frac{p}{d})}$$

然后我们就可以通过预处理$g(i)$和它的前缀积,然后枚举$p$分块去搞,用$last$前缀积除以$i-1$的前缀积就可以得到这一段的乘积,然后将它们对 $\lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor$做快速幂即可

我们怎么预处理$g(i)$呢?

我们可以枚举$d$,然后统计它对它的倍数的贡献

观察式子$f(d)^{\mu(\frac{p}{d})}$,$\mu$的取值只有三种。为$0$时不用管,为$1$时需要乘$f(i)$,为$-1$时需要乘$f(i)^{-1}$

关键问题就是预处理出$f(i)^{-1}$

我们线性推就好了。设$$s[x]=\prod_{i=1}^x f(x)$$ $$rs[x]=s[x]^{-1}$$ $$ni[x]=f(x)^{-1}$$

于是$$rs[x-1]=rs[x]*f[x]$$ $$ni[x]=rs[x]*s[x-1]$$我们只需要求出$rs[n]$然后倒着推一遍就行了(撒花!)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#define N 1000100
#define pos(i,a,b) for(LL i=(a);i<=(b);i++)
#define pos2(i,a,b) for(LL i=(a);i>=(b);i--)
using namespace std;
#define LL long long
const int p=(LL)1e9+7;
int t;
LL n,m;
int notprime[N],prime[N],mu[N];
LL x,y;
LL s[N],rs[N],f[N],ni[N];
LL g[N];
void kuogcd(LL a,LL b,LL &x,LL &y){
	if(!b){
		x=1;y=0;return;
	}
	kuogcd(b,a%b,x,y);
	LL temp=x;x=y;y=temp-a/b*y;
}
LL qpow(LL a,LL k){
	LL res=1;
	while(k>0){
		if(k&1){
			res=(res*a)%p;
		}
		a=(a*a)%p;k>>=1;
	}
	return res;
}
void get_pre(){
	notprime[1]=mu[1]=1;
	pos(i,2,N-10){
		if(!notprime[i]){
			prime[++prime[0]]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=prime[0]&&prime[j]*i<=N-10;j++){
			notprime[i*prime[j]]=1;
			if(i%prime[j]==0){
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=-mu[i];
		}
	}
	f[1]=1;
	pos(i,2,N-10) f[i]=(f[i-1]+f[i-2])%p;
	s[0]=1;
	pos(i,1,N-10){
		s[i]=s[i-1]*f[i]%p;
	}
	kuogcd(s[N-10],p,x,y);
	while(x<0) x+=p;
	rs[N-10]=x%p;
	pos2(i,N-10,1){
		rs[i-1]=rs[i]*f[i]%p;
		ni[i]=rs[i]*s[i-1]%p;
	}
	pos(i,1,N-10) g[i]=1;
	pos(i,1,N-10){
		pos(j,1,N-10){
			if(i*j>N-10) break;
			int pp=i*j;
			if(mu[j]==1) g[pp]=(g[pp]*f[i])%p;	
			else if(mu[j]==-1) g[pp]=(g[pp]*ni[i])%p;
		}
	}
	g[0]=1;
	pos(i,1,N-10) g[i]=(g[i]*g[i-1])%p;
}
LL get_ans(){
	LL res(1);
	LL last(0);
	for(LL i=1;i<=n;i=last+1){
		last=min(n/(n/i),m/(m/i));
		kuogcd(g[i-1],p,x,y);
		while(x<0) x+=p;
		res=res*qpow(g[last]*x%p,(n/i)*(m/i))%p;
	}
	return res%p;
}
int main(){
	get_pre();
	scanf("%d",&t);
	while(t--){
		scanf("%lld%lld",&n,&m);if(n>m) n^=m,m^=n,n^=m;
		printf("%lld\n",get_ans());
	}
	return 0;
}

  

posted @ 2017-12-10 21:16 Hallmeow 阅读(...) 评论(...) 编辑 收藏