bzoj 4816: [Sdoi2017]数字表格【莫比乌斯反演+逆元】

把题意简化，就是要求
\[ \prod_{d=1}^{min(n,m)}f[d]^{\sum_{i=1}^{n}\sum_{j=1}^{m}e[gcd(i,j)==d]} \]
把幂用莫比乌斯反演转化，得到
\[ \prod_{d=1}^{min(n,m)}f[d]^{\sum_{k=1}^{min(\frac{n}{d},\frac{m}{d})}\mu(k)\left \lfloor \frac{n}{dk} \right \rfloor\left \lfloor \frac{m}{dk} \right \rfloor} \]
然后枚举q=dk
\[ \prod_{q=1}^{min(n,m)}\left ( \prod_{d|q}f[d]^{\mu(\frac{q}{d})} \right )^{\left \lfloor \frac{n}{q} \right \rfloor\left \lfloor \frac{m}{q} \right \rfloor } \]
用枚举因数的方法处理出\( f[d]^{\mu(\frac{q}{d})} \)，根据调和级数，复杂度为\( O(nlog_2n) \)，然后处理询问的时候分块，复杂度为\( O(\sqrt{n}+\sqrt{m}) \)
因为用了\( O(nlog_2n) \)的粗暴逆元求法，所以跑的比较慢…是可以线性求的。

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
const long long N=1000005,mod=1e9+7;
int T,n,m,mb
,p
,tot;
long long f
,t
,invf
,s
,invs
,ans;
bool v
;
int read()
{
int r=0,f=1;
char p=getchar();
while(p>'9'||p<'0')
{
if(p=='-')
f=-1;
p=getchar();
}
while(p>='0'&&p<='9')
{
r=r*10+p-48;
p=getchar();
}
return r*f;
}
long long inv(long long x)
{//cout<<x<<endl;
return x==1?1:(mod-mod/x)*inv(mod%x)%mod;
}
long long ksm(long long a,long long b)
{
long long r=1ll;
while(b)
{
if(b&1)
r=r*a%mod;
a=a*a%mod;
b>>=1;
}
return r;
}
int main()
{
T=read();
mb[1]=1;
for(int i=2;i<=1000000;i++)
{
if(!v[i])
{
p[++tot]=i;
mb[i]=-1;
}
for(int j=1;j<=tot&&i*p[j]<=1000000;j++)
{
int k=i*p[j];
v[k]=1;
if(i%p[j]==0)
{
mb[k]=0;
break;
}
mb[k]=-mb[i];
}
}
f[1]=1;
invf[1]=inv(1);
for(int i=2;i<=1000000;i++)
{
f[i]=(f[i-1]+f[i-2])%mod;
invf[i]=inv(f[i]);
}//cout<<"OKF"<<endl;
// for(int i=1;i<=20;i++)
// cout<<f[i]<<" "<<invf[i]<<endl;
for(int i=1;i<=1000000;i++)
t[i]=1;
for(int i=1;i<=1000000;i++)
for(int j=i;j<=1000000;j+=i)
{
if(mb[j/i]==-1)
t[j]=(long long)t[j]*(long long)invf[i]%mod;
else if(mb[j/i]==1)
t[j]=(long long)t[j]*(long long)f[i]%mod;
}//cout<<"ok"<<endl;
// for(int i=1;i<=20;i++)
// cout<<i<<" "<<t[i]<<endl;
s[0]=invs[0]=1ll;
for(int i=1;i<=1000000;i++)
{
s[i]=s[i-1]*t[i]%mod;
invs[i]=inv(s[i]);
}
while(T--)
{
n=read(),m=read();
ans=1ll;
if(n>m)
swap(n,m);
for(int i=1,la;i<=n;i=la+1)
{
int ni=n/i,mi=m/i;
la=min(n/ni,m/mi);
ans=ans*ksm(s[la]*invs[i-1]%mod,(long long)ni*mi)%mod;
}
printf("%lld\n",ans);
}
return 0;
}