题目链接
题解
满满的反演的套路的味道
[Ans= prod_{i=1}^{n} prod_{j=1}^{m} f[gcd(i,j)]
]
常规操作枚举约数
[prod_{d=1}^nprod_{i=1}^nprod_{j=1}^m gcd(i,j)==d ? f[gcd(i,j)]
]
上面的式子可以化为
[prod_{d=1}^{n}f[d]^{sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)==d]}
]
对于右上角的式子
[sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)==d]
]
很眼熟...推下式子
[=sum_{i=1}^{⌊ dfrac{n}{d}⌋}sum_{j=1}^{⌊dfrac{m}{d}⌋}[gcd(i,j)==1]
]
[=sum_{i=1}^{min(dfrac{n}{d},dfrac{m}{d})} mu(i)sum_{i|p,pleq n/d }sum_{i|j,jleq m/d}1
]
[=sum_{i=1}^{n/d}mu(i)left lfloor frac{n}{id}
ight
floor left lfloorfrac{m}{id}
ight
floor
]
令(p=id)那么原式变成了
[prod_{p=1}^{min(n,m)} prod_{d|p} f[d]^{left lfloor frac{n}{p}
ight
floor left lfloor frac{m}{p}
ight
floor mu(frac{p}{d})} = prod_{p=1}^{min(n,m)} (prod_{d|p} f[d]^{mu(frac{p}{d})})^{left lfloor frac{n}{p}
ight
floor left lfloor frac{m}{p}
ight
floor}
]
对于每个(p),令(F[p]=prod_{d|p} f[d]^{mu(frac{p}{d})}),(F[p])的值是固定的，用筛法求出(F[p])，做前缀积，对与(left lfloor frac{n}{p}
ight
floor left lfloor frac{m}{p}
ight
floor)数论分块一下
复杂度为(O((n+Tsqrt{n})log(1e9+7)))
代码
#include<cstdio>
#include<algorithm>
const int maxn = 1000007;
const int mod = 1e9+7;
int f[maxn+7],invf[maxn+7],F[maxn+7];
bool isprime[maxn+7];int mu[maxn+7],prime[maxn+7],num;
int pow(int a,int p) {
int ret=1;
for(;p;p>>=1,a=(1LL*a*a)%mod){
if(p&1) ret=(1LL*ret*a)%mod;
}
return ret;
}
void init() {
f[1]=1;
for(int i=2;i<maxn;++i) f[i]=(f[i-1]+f[i-2])%mod;
for(int i=1;i<maxn;++i) {
invf[i]=pow(f[i],mod-2)%mod;
//printf("%d %d
",invf[i],f[i]);
}
isprime[1]=1;mu[1]=1;F[0]=F[1]=1;
//printf("%d
",mu[1]);
for(int i=2;i<=maxn-7;++i) {
F[i]=1;
if(!isprime[i]) prime[++num]=i,mu[i]=-1;
for(int j=1;j<=num&&prime[j]*i<=maxn-7;++j) {
isprime[i*prime[j]]=1;
if(i%prime[j]==0) {mu[i*prime[j]]=0;break;}
else mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=maxn-7;++i) {
if(mu[i]==0)continue;
for(int j=i;j<=maxn-7;j+=i) {
if(mu[i]==1) F[j]=(1ll*F[j]*f[j/i])%mod;
if(mu[i]==-1) F[j]=(1ll*F[j]*invf[j/i])%mod;
}
}
for(int i=1;i<=maxn-7;++i) F[i]=(1ll*F[i]*F[i-1])%mod;//,printf("%d
",F[i]);
}
int query(int a,int b) {
int n=std::min(a,b);long long ans=1;
for(int nxt,i=1;i<=n;i=nxt+1) {
nxt=std::min(a/(a/i),b/(b/i));
long long tmp=1ll*F[nxt]*pow(F[i-1],mod-2)%mod;
ans=(1ll*ans*pow(tmp,1ll*(a/i)*(b/i)%(mod-1)))%mod;
}
return (ans+mod)%mod;
}
int main() {
//freopen("001.out","w",stdout);
init();
int T;scanf("%d",&T);
for(int a,b;T--;) {
scanf("%d%d",&a,&b);
printf("%d
",query(a,b));
}
return 0;
}