bzoj4816[SDOI2017]数字表格

题目链接
以下除法均指除后下取整

\[\prod_{i=1}^n \prod_{j=1}^mf(gcd(i,j))\\ =\prod_{x}f(x)^{\sum_i \sum_j [gcd(i,j)=x] } \\ =\prod_{x}f(x)^{ \sum_{x|d} \mu(\frac{d}{x})\frac{n}{d} \frac{m}{d} }\\ =\prod_{d}(\prod_{x|d} f(x)^{\mu(\frac{d}{x})} )^{\frac{n}{d}\frac{m}{d}}\]

然后把中间那个括号里的东西看成\(g(d)\)预处理,即可做到\(O(nlogn)\)预处理\(O(\sqrt n *logn)\)询问

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<vector>
#include<algorithm>
#include<cmath>
#define P puts("lala")
#define cp cerr<<"lala"<<endl
#define ln putchar('\n')
#define pb push_back
#define fi first
#define se second
#define mkp make_pair
using namespace std;
inline int read()
{
    char ch=getchar();int g=1,re=0;
    while(ch<'0'||ch>'9') {if(ch=='-')g=-1;ch=getchar();}
    while(ch<='9'&&ch>='0') re=(re<<1)+(re<<3)+(ch^48),ch=getchar();
    return re*g;
}
typedef long long ll;
typedef pair<int,int> pii;

const int N=1000050;
const int inf=0x3f3f3f3f;
const int mod=1e9+7;
inline ll qpow(ll a,ll n)
{
	ll ans=1;
	for(;n;n>>=1,a=a*a%mod) if(n&1) ans=ans*a%mod;
	return ans;
}

int f[N],finv[N],g[N],ginv[N],prime[N],cnt=0,mu[N];
bool isnotprime[N];

void init()
{
	int n=N-50;
	f[1]=1; g[1]=1;
	for(int i=2;i<=n;++i) f[i]=(f[i-1]+f[i-2])%mod,g[i]=1;
	for(int i=1;i<=n;++i) finv[i]=qpow(f[i],mod-2);

	isnotprime[1]=1; mu[1]=1;
	for(int i=2;i<=n;++i)
	{
		if(!isnotprime[i]) prime[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*prime[j]<=n;++j)
		{
			isnotprime[i*prime[j]]=1;
			if(!(i%prime[j])) break;
			mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=1;i<=n;++i) for(int j=i;j<=n;j+=i)
	{
		if(mu[j/i]>0) g[j]=1ll*g[j]*f[i]%mod;
		else if(!mu[j/i]) ;
		else g[j]=1ll*g[j]*finv[i]%mod;
	}
	g[0]=1; ginv[0]=1;
	for(int i=1;i<=n;++i) g[i]=1ll*g[i]*g[i-1]%mod;
	for(int i=1;i<=n;++i) ginv[i]=qpow(g[i],mod-2);
}

int main()
{
#ifndef ONLINE_JUDGE
	freopen("1.in","r",stdin);freopen("1.out","w",stdout);
#endif
	init();
	int T=read();
	for(int cas=1;cas<=T;++cas)
	{
		int n=read(),m=read();
		int M=min(n,m),ans=1;
		for(int i=1,j;i<=M;i=j+1)
		{
			j=min(n/(n/i),m/(m/i));
			ans=1ll*ans*qpow(1ll*g[j]*ginv[i-1]%mod,1ll*(n/i)*(m/i))%mod;
		}
		printf("%d\n",ans);
	}
	return 0;
}