luoguP3704 [SDOI2017]数字表格

题意

默认(nleqslant m)
所求即为：(prodlimits_{i=1}^nprodlimits_{j=1}^mf[gcd(i,j)])
枚举(gcd(i,j))变为：
(prodlimits_{k=1}^{n}f(k)^{sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=k]})
上面那个是莫比乌斯反演套路形式：
(sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=k])
(sumlimits_{i=1}^{frac{n}{k}}sumlimits_{j=1}^{frac{m}{k}}[gcd(i,j)=1])
(sumlimits_{i=1}^{frac{n}{k}}sumlimits_{j=1}^{frac{m}{k}}sumlimits_{x|gcd(i,j)}mu(x))
(sumlimits_{x=1}^{frac{n}{k}}mu(x)sumlimits_{i=1}^{frac{n}{k}}sumlimits_{j=1}^{frac{m}{k}}[x|gcd(i,j)])
(sumlimits_{x=1}^{frac{n}{k}}mu(x)sumlimits_{i=1}^{frac{n}{k*x}}sumlimits_{j=1}^{frac{m}{k*x}}1)
(sumlimits_{x=1}^{frac{n}{k}}mu(x)frac{n}{k*x}frac{m}{k*x})
代回原式：
(prodlimits_{k=1}^{n}f(k)^{sumlimits_{x=1}^{frac{n}{k}}mu(x)frac{n}{k*x}frac{m}{k*x}})
设(T=k*x)，转而枚举(T)：
(prodlimits_{T=1}^{n}prodlimits_{d|T}f(d)^{mu(frac{T}{d})frac{n}{T}frac{m}{T}})
(prodlimits_{T=1}^{n}(prodlimits_{d|T}f(d)^{mu(frac{T}{d})})^{frac{n}{T}frac{m}{T}})
显然指数部分可以除法分块，考虑如何求(prodlimits_{d|T}f(d)^{mu(frac{T}{d})})：
设(g(T)=prodlimits_{d|T}f(d)^{mu(frac{T}{d})})
在算到(f(d))时乘到(g(T))即可。
答案即为：
(prodlimits_{T=1}^{n}g(T)^{frac{n}{T}frac{m}{T}})

code:

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=1e6+10;
const ll mod=1e9+7;
int T,n,m;
ll ans=1;
ll f[maxn],g[maxn],mu[maxn],invf[maxn];
bool vis[maxn];
vector<int>prime;
inline ll power(ll x,ll k,ll mod)
{
	ll res=1;
	while(k)
	{
		if(k&1)res=res*x%mod;
		x=x*x%mod;k>>=1;
	}
	return res;
}
inline void pre_work(int n)
{
	g[0]=g[1]=1;
	vis[1]=1;mu[1]=1;
	for(int i=2;i<=n;i++)
	{
		g[i]=1;
		if(!vis[i])prime.push_back(i),mu[i]=-1;
		for(unsigned int j=0;j<prime.size()&&i*prime[j]<=n;j++)
		{
			vis[i*prime[j]]=1;
			if(i%prime[j]==0)break;
			mu[i*prime[j]]=-mu[i];	
		}	
	}
	f[0]=0,f[1]=1;invf[1]=1;
	for(int i=2;i<=n;i++)f[i]=(f[i-1]+f[i-2])%mod,invf[i]=power(f[i],mod-2,mod);
	for(int i=1;i<=n;i++)
	{
		if(!mu[i])continue;
		for(int j=i;j<=n;j+=i)
			g[j]=g[j]*(mu[i]==1?f[j/i]:invf[j/i])%mod;
	}
	for(int i=2;i<=n;i++)g[i]=g[i]*g[i-1]%mod;
}
int main()
{
	pre_work(1000000);
	scanf("%d",&T);
	while(T--)
	{
		ans=1;
		scanf("%d%d",&n,&m);
		if(n>m)swap(n,m);
		for(int l=1,r;l<=n;l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			ans=ans*power(g[r]*power(g[l-1],mod-2,mod)%mod,1ll*(n/l)*(m/l)%(mod-1),mod)%mod;
		}
		printf("%lld
",(ans+mod)%mod);
	}
	return 0;
}