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[N],p[N],tot;
long long f[N],t[N],invf[N],s[N],invs[N],ans;
bool v[N];
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;
}