题意
(T) 组询问，每组询问给定 (n,m)，求 (prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}) 对 (10^9+7) 取模的值，(F) 是 Fibonacci 数列。
( exttt{Data Range:}1leq Tleq 10^3,1leq n,mleq 10^6)
题解
( exttt{Y})( exttt{LWang}) 是神！
考虑套路地枚举一发 (gcd)，有
[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{d=1}^{n}prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}[gcd(i,j)=d]
]
再考虑套路地把 (gcd) 提出来
[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{d=1}^{n}prodlimits_{i=1}^{leftlfloorfrac{n}{d}
ight
floor}prodlimits_{j=1}^{leftlfloorfrac{m}{d}
ight
floor}F_{d}[gcd(i,j)=1]
]
然后我们可以将这个式子写成 (F_d) 的幂的形式
[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{d=1}^{n}F_d^{sumlimits_{i=1}^{leftlfloorfrac{n}{d}
ight
floor}sumlimits_{j=1}^{leftlfloorfrac{m}{d}
ight
floor}[gcd(i,j)=1]}
]
注意到上面其实就是一个套路反演
[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{d=1}^{n}F_d^{sumlimits_{k=1}^{leftlfloorfrac{n}{d}
ight
floor}mu(k)leftlfloorfrac{n}{dk}
ight
floorleftlfloorfrac{m}{dk}
ight
floor}
]
然后给你写回来
[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{d=1}^{n}prodlimits_{k=1}^{leftlfloorfrac{n}{d}
ight
floor}F_d^{mu(k)leftlfloorfrac{n}{dk}
ight
floorleftlfloorfrac{m}{dk}
ight
floor}
]
现在做代换 (T=dk)
[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{T=1}^{n}prodlimits_{dmid T}F_d^{muleft(frac{T}{d}
ight)leftlfloorfrac{n}{T}
ight
floorleftlfloorfrac{m}{T}
ight
floor}
]
见证奇迹的时刻到了
[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{T=1}^{n}left(prodlimits_{dmid T}F_d^{muleft(frac{T}{d}
ight)}
ight)^{leftlfloorfrac{n}{T}
ight
floorleftlfloorfrac{m}{T}
ight
floor}
]
括号里面的东西可以预处理出来，括号外面整除分块即可，由于可能要求很多次逆元，所以建议写个线性求逆。
代码
#include<bits/stdc++.h>
using namespace std;
typedef int ll;
typedef long long int li;
const ll MAXN=1e6+51,MOD=1e9+7;
ll test,ptot,n,m,x;
ll np[MAXN],prime[MAXN],mu[MAXN],fib[MAXN],f[MAXN],prod[MAXN];
ll pr[MAXN],invf[MAXN],invp[MAXN];
inline ll read()
{
register ll num=0,neg=1;
register char ch=getchar();
while(!isdigit(ch)&&ch!='-')
{
ch=getchar();
}
if(ch=='-')
{
neg=-1;
ch=getchar();
}
while(isdigit(ch))
{
num=(num<<3)+(num<<1)+(ch-'0');
ch=getchar();
}
return num*neg;
}
inline ll qpow(ll base,ll exponent)
{
ll res=1;
while(exponent)
{
if(exponent&1)
{
res=(li)res*base%MOD;
}
base=(li)base*base%MOD,exponent>>=1;
}
return res;
}
inline void sieve(ll limit)
{
np[1]=mu[1]=f[1]=fib[1]=1;
for(register int i=2;i<=limit;i++)
{
fib[i]=(fib[i-1]+fib[i-2])%MOD,f[i]=1;
if(!np[i])
{
prime[++ptot]=i,mu[i]=-1;
}
for(register int j=1;i*prime[j]<=limit;j++)
{
np[i*prime[j]]=1;
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
}
inline ll calc(ll n,ll m)
{
ll res=1,cur;
for(register int l=1,r;l<=min(n,m);l=r+1)
{
r=min(n/(n/l),m/(m/l)),cur=1,cur=(li)prod[r]%MOD*invp[l-1]%MOD;
res=(li)res*qpow(cur,(li)(n/l)*(m/l)%(MOD-1))%MOD;
}
return res;
}
int main()
{
test=read(),sieve(1e6+10),prod[0]=pr[0]=1;
for(register int i=1;i<=1e6;i++)
{
pr[i]=(li)pr[i-1]*fib[i]%MOD;
}
x=qpow(pr[1000000],MOD-2),invf[1]=1;
for(register int i=1e6-1;i;i--)
{
invf[i+1]=(li)pr[i]*x%MOD,x=(li)x*fib[i+1]%MOD;
}
for(register int i=1;i<=1e6;i++)
{
for(register int j=1;i*j<=1e6;j++)
{
f[i*j]=(li)f[i*j]*(mu[j]==1?fib[i]:mu[j]==-1?invf[i]:1)%MOD;
}
}
for(register int i=1;i<=1e6;i++)
{
prod[i]=(li)prod[i-1]*f[i]%MOD,pr[i]=(li)pr[i-1]*prod[i]%MOD;
}
x=qpow(pr[1000000],MOD-2),invp[0]=invp[1]=1;
for(register int i=1e6-1;i;i--)
{
invp[i+1]=(li)pr[i]*x%MOD,x=(li)x*prod[i+1]%MOD;
}
for(register int i=0;i<test;i++)
{
n=read(),m=read(),printf("%d
",calc(n,m));
}
}