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

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

$$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}$$

$$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}$$

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

#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 阅读(...) 评论(...) 编辑 收藏