luoguP3704 [SDOI2017]数字表格

题意

默认\(n\leqslant m\)
所求即为：\(\prod\limits_{i=1}^n\prod\limits_{j=1}^mf[\gcd(i,j)]\)
枚举\(\gcd(i,j)\)变为：
\(\prod\limits_{k=1}^{n}f(k)^{\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=k]}\)
上面那个是莫比乌斯反演套路形式：
\(\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=k]\)
\(\sum\limits_{i=1}^{\frac{n}{k}}\sum\limits_{j=1}^{\frac{m}{k}}[\gcd(i,j)=1]\)
\(\sum\limits_{i=1}^{\frac{n}{k}}\sum\limits_{j=1}^{\frac{m}{k}}\sum\limits_{x|\gcd(i,j)}\mu(x)\)
\(\sum\limits_{x=1}^{\frac{n}{k}}\mu(x)\sum\limits_{i=1}^{\frac{n}{k}}\sum\limits_{j=1}^{\frac{m}{k}}[x|\gcd(i,j)]\)
\(\sum\limits_{x=1}^{\frac{n}{k}}\mu(x)\sum\limits_{i=1}^{\frac{n}{k*x}}\sum\limits_{j=1}^{\frac{m}{k*x}}1\)
\(\sum\limits_{x=1}^{\frac{n}{k}}\mu(x)\frac{n}{k*x}\frac{m}{k*x}\)
代回原式：
\(\prod\limits_{k=1}^{n}f(k)^{\sum\limits_{x=1}^{\frac{n}{k}}\mu(x)\frac{n}{k*x}\frac{m}{k*x}}\)
设\(T=k*x\)，转而枚举\(T\)：
\(\prod\limits_{T=1}^{n}\prod\limits_{d|T}f(d)^{\mu(\frac{T}{d})\frac{n}{T}\frac{m}{T}}\)
\(\prod\limits_{T=1}^{n}(\prod\limits_{d|T}f(d)^{\mu(\frac{T}{d})})^{\frac{n}{T}\frac{m}{T}}\)
显然指数部分可以除法分块，考虑如何求\(\prod\limits_{d|T}f(d)^{\mu(\frac{T}{d})}\)：
设\(g(T)=\prod\limits_{d|T}f(d)^{\mu(\frac{T}{d})}\)
在算到\(f(d)\)时乘到\(g(T)\)即可。
答案即为：
\(\prod\limits_{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\n",(ans+mod)%mod);
	}
	return 0;
}