bzoj4816: [Sdoi2017]数字表格
题目链接
题解
满满的反演的套路的味道
\[Ans= \prod_{i=1}^{n} \prod_{j=1}^{m} f[gcd(i,j)]
\]
常规操作枚举约数
\[\prod_{d=1}^n\prod_{i=1}^n\prod_{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,p\leq n/d }\sum_{i|j,j\leq m/d}1
\]
\[=\sum_{i=1}^{n/d}\mu(i)\left \lfloor \frac{n}{id}\right \rfloor \left \lfloor\frac{m}{id} \right \rfloor
\]
令\(p=id\)那么原式变成了
\[\prod_{p=1}^{min(n,m)} \prod_{d|p} f[d]^{\left \lfloor \frac{n}{p} \right \rfloor \left \lfloor \frac{m}{p} \right \rfloor \mu(\frac{p}{d})} = \prod_{p=1}^{min(n,m)} (\prod_{d|p} f[d]^{\mu(\frac{p}{d})})^{\left \lfloor \frac{n}{p} \right \rfloor \left \lfloor \frac{m}{p} \right \rfloor}
\]
对于每个\(p\),令\(F[p]=\prod_{d|p} f[d]^{\mu(\frac{p}{d})}\),\(F[p]\)的值是固定的，用筛法求出\(F[p]\)，做前缀积，对与\(\left \lfloor \frac{n}{p} \right \rfloor \left \lfloor \frac{m}{p} \right \rfloor\)数论分块一下
复杂度为\(O((n+T\sqrt{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\n",invf[i],f[i]);
}
isprime[1]=1;mu[1]=1;F[0]=F[1]=1;
//printf("%d\n",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\n",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\n",query(a,b));
}
return 0;
}
