第一次见到这种形式，推了两步就被卡住了。没想到还能这么迁移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; }