第一次见到这种形式，推了两步就被卡住了。没想到还能这么迁移QAQ，真强！
题目设 $f(i)$为第$i$项斐波那契数列
要求计算
$$ans=\prod_{i=1}^n \prod_{j=1}^m f(gcd(i,j))$$
枚举$d=gcd$得到
$$ans=\prod_{d=1}^n f(d)^{\sum_{i=1}^n \sum_{j=1}^m[gcd(i,j)==d]}$$
然后反演得到
$$ans=\prod_{d=1}^n f(d)^{\sum_{T=1}^{\lfloor \frac{n}{d} \rfloor}\mu(T) \lfloor \frac{n}{dT} \rfloor \lfloor \frac{m}{dT} \rfloor}$$
然后呢，令$p=dT$，得到
$$ans=\left( \prod_{d=1}^nf(d)^{\mu(\frac{p}{d})} \right) ^{\sum_{p=1}^n \lfloor \frac{n}{p} \rfloor \lfloor \frac{m}{p} \rfloor}$$
这样我们就发现我们可以对$p$进行分块去搞了
令
$$g(p)=\prod_{d|p} f(d)^{\mu(\frac{p}{d})}$$
然后我们就可以通过预处理$g(i)$和它的前缀积，然后枚举$p$分块去搞，用$last$前缀积除以$i-1$的前缀积就可以得到这一段的乘积，然后将它们对 $\lfloor \frac{n}{i} \rfloor \lfloor \frac{m}{i} \rfloor$做快速幂即可
我们怎么预处理$g(i)$呢？
我们可以枚举$d$，然后统计它对它的倍数的贡献
观察式子$f(d)^{\mu(\frac{p}{d})}$，$\mu$的取值只有三种。为$0$时不用管，为$1$时需要乘$f(i)$，为$-1$时需要乘$f(i)^{-1}$
关键问题就是预处理出$f(i)^{-1}$
我们线性推就好了。设$$s[x]=\prod_{i=1}^x f(x)$$ $$rs[x]=s[x]^{-1}$$ $$ni[x]=f(x)^{-1}$$
于是$$rs[x-1]=rs[x]*f[x]$$ $$ni[x]=rs[x]*s[x-1]$$我们只需要求出$rs[n]$然后倒着推一遍就行了（撒花！）
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#define N 1000100
#define pos(i,a,b) for(LL i=(a);i<=(b);i++)
#define pos2(i,a,b) for(LL i=(a);i>=(b);i--)
using namespace std;
#define LL long long
const int p=(LL)1e9+7;
int t;
LL n,m;
int notprime[N],prime[N],mu[N];
LL x,y;
LL s[N],rs[N],f[N],ni[N];
LL g[N];
void kuogcd(LL a,LL b,LL &x,LL &y){
if(!b){
x=1;y=0;return;
}
kuogcd(b,a%b,x,y);
LL temp=x;x=y;y=temp-a/b*y;
}
LL qpow(LL a,LL k){
LL res=1;
while(k>0){
if(k&1){
res=(res*a)%p;
}
a=(a*a)%p;k>>=1;
}
return res;
}
void get_pre(){
notprime[1]=mu[1]=1;
pos(i,2,N-10){
if(!notprime[i]){
prime[++prime[0]]=i;
mu[i]=-1;
}
for(int j=1;j<=prime[0]&&prime[j]*i<=N-10;j++){
notprime[i*prime[j]]=1;
if(i%prime[j]==0){
mu[i*prime[j]]=0;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
f[1]=1;
pos(i,2,N-10) f[i]=(f[i-1]+f[i-2])%p;
s[0]=1;
pos(i,1,N-10){
s[i]=s[i-1]*f[i]%p;
}
kuogcd(s[N-10],p,x,y);
while(x<0) x+=p;
rs[N-10]=x%p;
pos2(i,N-10,1){
rs[i-1]=rs[i]*f[i]%p;
ni[i]=rs[i]*s[i-1]%p;
}
pos(i,1,N-10) g[i]=1;
pos(i,1,N-10){
pos(j,1,N-10){
if(i*j>N-10) break;
int pp=i*j;
if(mu[j]==1) g[pp]=(g[pp]*f[i])%p;
else if(mu[j]==-1) g[pp]=(g[pp]*ni[i])%p;
}
}
g[0]=1;
pos(i,1,N-10) g[i]=(g[i]*g[i-1])%p;
}
LL get_ans(){
LL res(1);
LL last(0);
for(LL i=1;i<=n;i=last+1){
last=min(n/(n/i),m/(m/i));
kuogcd(g[i-1],p,x,y);
while(x<0) x+=p;
res=res*qpow(g[last]*x%p,(n/i)*(m/i))%p;
}
return res%p;
}
int main(){
get_pre();
scanf("%d",&t);
while(t--){
scanf("%lld%lld",&n,&m);if(n>m) n^=m,m^=n,n^=m;
printf("%lld\n",get_ans());
}
return 0;
}
