[SDOI2017]数字表格【莫比乌斯反演】

[SDOI2017]数字表格
求$$\prod_{i=1}^n\prod_{j=1}mf(gcd(i, j))$$,其中\(f\)为\(fibonacci\)数列
\(T\le 1000, 1\le n, m\le {10}^6\)

\[\begin{aligned}ans&=\prod_{i=1}^n\prod_{j=1}^mf(gcd(i, j)) \\ &=\prod_{d=1}^nf(d)^{\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j) == d]} \\ &=\prod_{d=1}^nf(d)^{\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)\lfloor\frac{n}{xd}\rfloor\lfloor\frac{m}{xd}\rfloor}\\ &=\prod_{T=1}^n(\prod_{d|T}f(d)^{\mu(\frac{T}{d})})^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor} \end{aligned} \]

先暴力筛\(G(T)=\prod_{d|T}f(d)^{\mu(\frac{T}{d})}\)
就可以\(O(\sqrt{N}logN)\)的时间算$$\prod_{T=1}^nG(T)\rfloor\lfloor\frac{m}{T}\rfloor}$$

ll pow(ll x, ll k){ll res=1; while(k){if(k&1) res=res*x%p; x=x*x%p; k >>= 1;} return res;}
void init(){
	miu[1]=1; f[1]=invf[1]=1; g[1]=g[0]=1;
	for(int i=2; i < Maxn; i++){
		if(!pr[i]) pr[++ptot]=i, miu[i]=-1;
		f[i]=(f[i-1]+f[i-2])%p; invf[i]=pow(f[i], p-2); g[i]=1;
		for(int j=1, x; j <= ptot && (x=pr[j]*i) < Maxn; j++){
			pr[x]=1; if(i%pr[j] == 0) break; miu[x]=-miu[i];
		}
	}
	for(int i=1; i < Maxn; i++) if(miu[i]) for(int j=i; j < Maxn; j+=i) 
			g[j]=g[j]*(miu[i] > 0 ? f[j/i] : invf[j/i])%p;
	for(int i=2; i < Maxn; i++) g[i]=g[i-1]*g[i]%p;
}
void solve(){
	init(); int T=read(); 
	while(T--){
		n=read(), m=read(); if(n > m) swap(n, m); ll ans=1;
		for(ll l=1, r=0; r < n; l=r+1){
			r=min(n/(n/l), m/(m/l)); ll delta=g[r]*pow(g[l-1], p-2)%p;
			ans=ans*pow(delta, (n/l)*(m/l))%p;
		}
		printf("%lld\n", ans);
	}
}