[SDOI2017] 数字表格

求 \(\prod\limits_{i=1}^{n} \prod\limits_{j=1}^{m} f_{\gcd(i,j)}\) 的值，对 \(10^9 + 7\) 取模，\(f\) 是斐波那契数列。

\(T \leq 10^3\) 组数据，\(1 \leq n,m \leq 10^6\)。

---

\[\begin{aligned} & \prod\limits_{i=1}^{n} \prod\limits_{j=1}^{m} f_{\gcd(i,j)} \\ = & \prod\limits_{d=1}^{\min(n,m)} f_d^{\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum\limits_{j=1}^{\lfloor \frac{m}{d} \rfloor} [\gcd(i,j)=1]} \\ = & \prod\limits_{d=1}^{\min(n,m)} f_d^{\sum\limits_{t=1}^{\lfloor \frac{\min(n,m)}{d} \rfloor} \mu(t) \lfloor \frac{n}{dt} \rfloor \lfloor \frac{m}{dt} \rfloor} \\ = & \prod\limits_{k=1}^{\min(n,m)} \prod\limits_{d \mid k} f_d^{\mu(\frac{k}{d}) \lfloor \frac{n}{k} \rfloor \lfloor \frac{m}{k} \rfloor} \\ = & \prod\limits_{k=1}^{\min(n,m)} (\prod\limits_{d \mid k} f_d^{\mu(\frac{k}{d})})^{\lfloor \frac{n}{k} \rfloor \lfloor \frac{m}{k} \rfloor} \\ = & \prod\limits_{k=1}^{\min(n,m)} g(k)^{\lfloor \frac{n}{k} \rfloor \lfloor \frac{m}{k} \rfloor} \end{aligned} \]

\(g(k)\) 容易在 \(O(n \log n)\) 预处理，剩余的部分整除分块即可。

// Ad astra per aspera
#include<iostream>
#include<cstdio>
#include<vector>
using namespace std;
const long long mod=1000000007;
const int N=1000000;
int mu[N+10];
long long f[N+10],inv_f[N+10],g[N+10],pre_g[N+10],inv_g[N+10];
long long pow_mod(long long num1,long long num2){
	num2%=mod-1;
	long long num3=1;
	while(num2){
		if(num2&1){
			num3=num3*num1%mod;
		}
		num1=num1*num1%mod;
		num2>>=1;
	}
	return num3;
} 
int main(){
	f[1]=1;
	f[2]=1;
	for(int i=3;i<=N;i++){
		f[i]=(f[i-2]+f[i-1])%mod;
	}
	for(int i=1;i<=N;i++){
		inv_f[i]=pow_mod(f[i],mod-2);
		g[i]=1;
	}
	mu[1]=1;
	for(int i=1;i<=N;i++){
		for(int j=2;j<=N/i;j++){
			mu[i*j]-=mu[i];
		}
	}
	for(int i=1;i<=N;i++){
		for(int j=1;j<=N/i;j++){
			if(mu[j]==-1){
				g[i*j]=g[i*j]*inv_f[i]%mod;
			}
			if(mu[j]==1){
				g[i*j]=g[i*j]*f[i]%mod;
			}
		}
	}
	pre_g[0]=inv_g[0]=1;
	for(int i=1;i<=N;i++){
		pre_g[i]=pre_g[i-1]*g[i]%mod;
		inv_g[i]=inv_g[i-1]*pow_mod(g[i],mod-2)%mod;
	}
	int T;
	scanf("%d",&T);
	while(T--){
		int n,m;
		scanf("%d %d",&n,&m);
		int l=1;
		long long ans=1;
		while(l<=min(n,m)){
			int r=min(n/(n/l),m/(m/l));
			ans=ans*pow_mod(pre_g[r]*inv_g[l-1]%mod,(long long)(n/l)*(m/l))%mod;
			l=r+1;
		}
		printf("%lld\n",ans);
	}
    return 0;
}