把题意简化，就是要求
[prod_{d=1}^{min(n,m)}f[d]^{sum_{i=1}^{n}sum_{j=1}^{m}e[gcd(i,j)==d]}
]
把幂用莫比乌斯反演转化，得到
[prod_{d=1}^{min(n,m)}f[d]^{sum_{k=1}^{min(frac{n}{d},frac{m}{d})}mu(k)left lfloor frac{n}{dk}
ight
floorleft lfloor frac{m}{dk}
ight
floor}
]
然后枚举q=dk
[prod_{q=1}^{min(n,m)}left ( prod_{d|q}f[d]^{mu(frac{q}{d})}
ight )^{left lfloor frac{n}{q}
ight
floorleft lfloor frac{m}{q}
ight
floor }
]
用枚举因数的方法处理出( f[d]^{mu(frac{q}{d})} )，根据调和级数，复杂度为( O(nlog_2n) )，然后处理询问的时候分块，复杂度为( O(sqrt{n}+sqrt{m}) )
因为用了( O(nlog_2n) )的粗暴逆元求法，所以跑的比较慢…是可以线性求的。
#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
const long long N=1000005,mod=1e9+7;
int T,n,m,mb[N],p[N],tot;
long long f[N],t[N],invf[N],s[N],invs[N],ans;
bool v[N];
int read()
{
int r=0,f=1;
char p=getchar();
while(p>'9'||p<'0')
{
if(p=='-')
f=-1;
p=getchar();
}
while(p>='0'&&p<='9')
{
r=r*10+p-48;
p=getchar();
}
return r*f;
}
long long inv(long long x)
{//cout<<x<<endl;
return x==1?1:(mod-mod/x)*inv(mod%x)%mod;
}
long long ksm(long long a,long long b)
{
long long r=1ll;
while(b)
{
if(b&1)
r=r*a%mod;
a=a*a%mod;
b>>=1;
}
return r;
}
int main()
{
T=read();
mb[1]=1;
for(int i=2;i<=1000000;i++)
{
if(!v[i])
{
p[++tot]=i;
mb[i]=-1;
}
for(int j=1;j<=tot&&i*p[j]<=1000000;j++)
{
int k=i*p[j];
v[k]=1;
if(i%p[j]==0)
{
mb[k]=0;
break;
}
mb[k]=-mb[i];
}
}
f[1]=1;
invf[1]=inv(1);
for(int i=2;i<=1000000;i++)
{
f[i]=(f[i-1]+f[i-2])%mod;
invf[i]=inv(f[i]);
}//cout<<"OKF"<<endl;
// for(int i=1;i<=20;i++)
// cout<<f[i]<<" "<<invf[i]<<endl;
for(int i=1;i<=1000000;i++)
t[i]=1;
for(int i=1;i<=1000000;i++)
for(int j=i;j<=1000000;j+=i)
{
if(mb[j/i]==-1)
t[j]=(long long)t[j]*(long long)invf[i]%mod;
else if(mb[j/i]==1)
t[j]=(long long)t[j]*(long long)f[i]%mod;
}//cout<<"ok"<<endl;
// for(int i=1;i<=20;i++)
// cout<<i<<" "<<t[i]<<endl;
s[0]=invs[0]=1ll;
for(int i=1;i<=1000000;i++)
{
s[i]=s[i-1]*t[i]%mod;
invs[i]=inv(s[i]);
}
while(T--)
{
n=read(),m=read();
ans=1ll;
if(n>m)
swap(n,m);
for(int i=1,la;i<=n;i=la+1)
{
int ni=n/i,mi=m/i;
la=min(n/ni,m/mi);
ans=ans*ksm(s[la]*invs[i-1]%mod,(long long)ni*mi)%mod;
}
printf("%lld
",ans);
}
return 0;
}