BZOJ 4816: [Sdoi2017]数字表格 莫比乌斯反演

#include<bits/stdc++.h>
#define ll long long 
#define maxn 1000005 
#define N 1000003        
using namespace std;
const ll mod=1000000007;       
namespace IO
{ 
	inline void setIO(string s)
	{
		string in=s+".in"; 
		freopen(in.c_str(),"r",stdin); 
	}
	inline void shut() 
	{
		fclose(stdin); 
		fclose(stdout); 
	}
}; 
int cnt; 
int mu[maxn], vis[maxn], prime[maxn];  
ll f[maxn], inv[maxn], h[maxn];    
inline ll qpow(ll base,ll k)
{
	ll tmp=1; 
	while(k) 
	{
		if(k&1) tmp=tmp*base%mod; 
		base=base*base%mod; 
		k>>=1; 
	}
	return tmp;     
}
inline void prepare()
{  
	mu[1]=1; 
	for(int i=2;i<=N;++i) 
	{
		if(!vis[i]) prime[++cnt]=i,mu[i]=-1; 
		for(int j=1;j<=cnt&&1ll*prime[j]*i<=N;++j) 
		{
			vis[i*prime[j]]=1; 
			if(i%prime[j]==0) 
			{
				mu[i*prime[j]]=0; 
				break; 
			}
			mu[i*prime[j]]=-mu[i];   
		}
	}    
	inv[1]=1; 
	h[0]=1,h[1]=1; 
	f[0]=0,f[1]=1;   
	for(int i=2;i<=N;++i) f[i]=(f[i-1]+f[i-2])%mod, inv[i]=qpow(f[i],mod-2), h[i]=1;        
	for(int i=1;i<=N;++i) 
		if(mu[i]!=0) 
		{
			for(int j=i;j<=N;j+=i)      
			h[j]*=(mu[i]==1?f[j/i]:inv[j/i]),h[j]%=mod; 
		}
	for(int i=2;i<=N;++i) h[i]*=h[i-1], h[i]%=mod;        
}
inline ll work(int n,int m) 
{
	int i,j;  
	ll re=1,tmp; 
	for(i=1;i<=n;i=j+1)
	{
		j=min(n/(n/i), m/(m/i));    
		tmp=qpow(h[j]*qpow(h[i-1], mod-2)%mod,  1ll*(n/i)*(m/i)%(mod-1)); 
		re*=tmp, re%=mod;   
	}
	return re; 
}
int main()
{
	// IO::setIO("input");  
	prepare(); 
	int T,n,m; 
	scanf("%d",&T); 
	while(T--)
	{
		scanf("%d%d",&n,&m); 
		if(n>m) swap(n,m); 
		printf("%lld\n",work(n,m)); 
	} 
	IO::shut(); 
	return 0; 
}