[BZOJ4816][莫比乌斯反演]SDOI2017：数字表格

BZOJ4816
都是套路
$$\prod _{i=1}^n\prod _{j=1}^{m}fib[gcd(i,j)]$$
$$\prod _{k=1}^{min(n,m)}{fib[k]}^{\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }gcd(i,j)==k}$$
$$\prod _{k=1}^{min(n,m)}{fib[k]}^{\sum _{d=1}^{\lfloor n/d\rfloor }\mu (d)\lfloor \frac{n}{kd}\rfloor \lfloor \frac{m}{kd}\rfloor }$$
令$T=kd$
$$\prod _{T=1}^{min(n,m)}\prod _{d|T}{fib[d]}^{\mu (\frac{T}{d})\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor }$$
预处理${\prod }_{d|T}{fib[d]}^{\mu (\frac{T}{d})}$是$O(nlogn)$的，然后就完了

Code：

#include<bits/stdc++.h>
#define ll long long
#define mod 1000000007
using namespace std;
inline int read(){
	int res=0,f=1;char ch=getchar();
	while(!isdigit(ch)) {if(ch=='-') f=-f;ch=getchar();}
	while(isdigit(ch)) {res=(res<<1)+(res<<3)+(ch^48);ch=getchar();}
	return res*f;
}
inline int add(int x,int y){x+=y;if(x>=mod) x-=mod;return x;}
inline void inc(int &x,int y){x+=y;if(x>=mod) x-=mod;}
inline int mul(int x,int y){return (ll)x*y%mod;}
inline void Mul(int &x,int y){x=(ll)x*y%mod;}
inline int ksm(int a,int b){int res=1;for(;b;b>>=1,a=mul(a,a)) if(b&1) res=mul(res,a);return res;}
const int N=1e6+5;
int pri[N],cnt,mu[N],g[N],fib[N],inv[N];
int pt[N];
inline void init(int n){
	mu[1]=1;
	for(int i=2;i<=n;i++){
		if(!pt[i]) pri[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt && i*pri[j]<=n;j++){
			pt[i*pri[j]]=1;
			if(i%pri[j]==0) break;
			mu[i*pri[j]]=-mu[i];
		}
	}
	fib[0]=0,fib[1]=1;
	for(int i=2;i<=n;i++) fib[i]=add(fib[i-1],fib[i-2]);
	for(int i=0;i<=n;i++) g[i]=1;
	for(int i=1;i<=n;i++) inv[i]=ksm(fib[i],mod-2);
	for(int i=1;i<=n;i++){
		if(!mu[i]) continue;
		for(int j=i;j<=n;j+=i)
			Mul(g[j],mu[i]==1?fib[j/i]:inv[j/i]);
	}
	for(int i=1;i<=n;i++) Mul(g[i],g[i-1]);
}
int main(){
	init(1000000);
	int t=read();
	while(t--){
		int n=read(),m=read(),ans=1;
		if(n>m) swap(n,m);
		for(int i=1,j;i<=n;i=j+1){
			j=min(n/(n/i),m/(m/i));
			Mul(ans,ksm(mul(g[j],ksm(g[i-1],mod-2)),1ll*(n/i)*(m/i)%(mod-1)));
		}
		cout<<ans<<"\n";
	}
	return 0;
}