求
[Ans=prod_{i=1}^Nprod_{j=1}^MFib[(i,j)]
]
连乘的反演，其实并没有什么不一样
我们把套路柿子拿出来
[F(n)=sum_{i=1}^Nsum_{j=1}^M[n|(i,j)]=left lfloor frac{N}{n}
ight
floor imes left lfloor frac{M}{n}
ight
floor=sum_{n|d}f(d)
]
[f(n)=sum_{i=1}^Nsum_{j=1}^M[n=(i,j)]=sum_{n|d}mu(frac{d}{n})left lfloor frac{N}{d}
ight
floor left lfloor frac{M}{d}
ight
floor
]
我们要求的就是
[Ans=prod_{i=1}^NFib(i)^{f(i)}
]
把它化开
[Ans=prod_{i=1}^NFib(i)^{sum_{i|d}mu(frac{d}{i})left lfloor frac{N}{d}
ight
floor left lfloor frac{M}{d}
ight
floor}
]
非常显然的就是
[Ans=prod_{d=1}^N(prod_{i|d}Fib(i)^{mu(frac{d}{i})})^{left lfloor frac{N}{d}
ight
floor left lfloor frac{M}{d}
ight
floor}
]
利用调和级数在(O(nlogn))的时间内处理出(prod_{i|d}Fib(i)^{mu(frac{d}{i})})的值，做一个前缀积就好了，之后整除分块和快速幂一起上就好了
代码
#include<iostream>
#include<cstring>
#include<cstdio>
#define re register
#define LL long long
#define maxn 1000005
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
const LL mod=1e9+7;
inline int read()
{
char c=getchar();int x=0;
while(c<'0'||c>'9') c=getchar();
while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
}
LL fib[maxn],pre[maxn];
int T,N[1005],M[1005],U;
int f[maxn],p[maxn],mu[maxn];
LL exgcd(LL a,LL b,LL &x,LL &y) {if(!b) return x=1,y=0,a;LL r=exgcd(b,a%b,y,x);y-=a/b*x;return r;}
inline LL quick(LL a,LL b) {LL S=1;while(b) {if(b&1ll) S=S*a%mod;b>>=1ll;a=a*a%mod;} return S;}
inline LL inv(LL a){LL x,y;LL r=exgcd(a,mod,x,y);return (x%mod+mod)%mod;}
inline LL solve(LL a,int b) {if(!b) return 1;if(b==1) return a;if(b==-1) return inv(a);}
int main()
{
T=read();
for(re int i=1;i<=T;i++) N[i]=read(),M[i]=read();
for(re int i=1;i<=T;i++) if(N[i]>M[i]) std::swap(N[i],M[i]);
for(re int i=1;i<=T;i++) U=max(U,N[i]);
mu[1]=1,f[1]=1,pre[1]=1,pre[0]=1;
for(re int i=2;i<=U;i++)
{
pre[i]=1;
if(!f[i]) p[++p[0]]=i,mu[i]=-1;
for(re int j=1;j<=p[0]&&p[j]*i<=U;j++) {f[p[j]*i]=1;if(i%p[j]==0) break;mu[p[j]*i]=-1*mu[i];}
}
fib[1]=fib[2]=1;
for(re int i=3;i<=U;i++) fib[i]=fib[i-1]+fib[i-2],fib[i]%=mod;
for(re int i=1;i<=U;i++)
for(re int j=1;j*i<=U;j++) pre[i*j]*=solve(fib[i],mu[j]),pre[i*j]%=mod;
for(re int i=1;i<=U;i++) pre[i]*=pre[i-1],pre[i]%=mod;
for(re int t=1;t<=T;t++)
{
int n=N[t],m=M[t];
LL ans=1;
for(re LL l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans*=quick(pre[r]*inv(pre[l-1])%mod,(n/l)*(m/l)%(mod-1));ans%=mod;
}
printf("%lld
",ans);
}
return 0;
}