题目大意：略
第一次不看题解做出来反演题，我好菜啊
最初推的式子里是把矩阵快速幂然后欧拉降幂，然后发现弄出来的矩阵不能直接相乘..
但这依然改变不了反演玩的是套路的事实
以下默认$n leq m$
题目让我们求$prodlimits_{i=1}^{n} prodlimits_{j=1}^{m}f[gcd(i,j)]$
显然直接统计$f[gcd(i,j)]$很困难，我们每次指定一个值k，求有多少个$gcd(i,j)==k$，那么$k$的贡献就是$f[k]^{ans}$
变形$prodlimits_{k=1}^{n} f[k]^{sumlimits_{i=1}^{n} sumlimits_{j=1}^{m}[gcd(i,j)==k]}$
我们只需要处理出$sumlimits_{i=1}^{n} sumlimits_{j=1}^{m}[gcd(i,j)==k]$的数量即可，这是一道经典的反演题
再利用常规套路，$varepsilon =1*mu$变形，$sumlimits_{i=1}^{left lfloor frac{n}{k} ight floor} sumlimits_{j=1}^{left lfloor frac{m}{k} ight floor}sumlimits_{d|gcd(i,j)}mu(d)$
指数最后变成了$sumlimits_{d=1}^{left lfloor frac{n}{k} ight floor}mu(d) left lfloor frac{n}{kd} ight floor left lfloor frac{m}{kd} ight floor$
原式$prodlimits_{k=1}^{n} f[k]^{sumlimits_{d=1}^{left lfloor frac{n}{k} ight floor}mu(d) left lfloor frac{n}{kd} ight floor left lfloor frac{m}{kd} ight floor}$
还是常规套路提取出$Q=kd$，变形$prodlimits_{Q=1}^{n} prodlimits_{d|Q}f[d]^{mu(frac{Q}{d}) left lfloor frac{n}{Q} ight floor left lfloor frac{m}{Q} ight floor}$
令$g(Q)=prodlimits_{d|Q} f[d]^{mu(frac{Q}{d})}$，显然它的前$n$项可以在$O(nlogn)$内处理出
每次询问都需要整除分块，还需要快速幂。总时间$O(Tsqrt{n}logn+nlogn)$
1 #include <cstdio> 2 #include <cstring> 3 #include <algorithm> 4 #define N1 1010000 5 #define ll long long 6 #define dd double 7 #define maxn 1000000 8 using namespace std; 9 10 const ll mod2=1000000007; 11 int T,n,m; 12 int use[N1],pr[N1],mu[N1],cnt; 13 int qpow(ll x,ll y,const ll &mod) 14 { 15 int ans=1; 16 while(y){ 17 if(y&1) ans=ans*x%mod; 18 x=x*x%mod; y>>=1; 19 }return ans; 20 } 21 void exgcd(ll a,ll b,ll &x,ll &y) 22 { 23 if(!b){ x=1; y=0; return;} 24 else{ exgcd(b,a%b,x,y); ll t=x; x=y; y=t-a/b*y; } 25 } 26 ll f[N1],g[N1],sg[N1]; 27 void init() 28 { 29 int i,j; ll x,y,inv; mu[1]=1; 30 for(i=2;i<=maxn;i++) 31 { 32 if(!use[i]){ pr[++cnt]=i; mu[i]=-1; } 33 for(j=1;j<=cnt&&i*pr[j]<=maxn;j++) 34 { 35 use[i*pr[j]]=1; 36 if(i%pr[j]){ mu[i*pr[j]]=-mu[i]; } 37 else{ mu[i*pr[j]]=0; break; } 38 } 39 } 40 for(j=1;j<=maxn;j++) g[j]=1; 41 for(f[1]=1,i=2;i<=maxn;i++) 42 { 43 f[i]=(f[i-1]+f[i-2])%mod2; 44 exgcd(f[i],mod2,inv,y); inv=(inv%mod2+mod2)%mod2; 45 for(j=1;j*i<=maxn;j++) 46 if(mu[j]) g[i*j]=1ll*g[i*j]*(mu[j]>0?f[i]:inv)%mod2; 47 } 48 for(sg[0]=1,i=1;i<=maxn;i++) sg[i]=1ll*sg[i-1]*g[i]%mod2; 49 } 50 51 int main() 52 { 53 scanf("%d",&T); 54 init(); 55 int i,j,la,ans=0; ll inv,y; 56 while(T--) 57 { 58 scanf("%d%d",&n,&m); if(n>m) swap(n,m); 59 for(i=1,ans=1;i<=n;i=la+1) 60 { 61 la=min(n/(n/i),m/(m/i)); 62 exgcd(sg[i-1],mod2,inv,y); inv=(inv%mod2+mod2)%mod2; 63 ans=1ll*qpow(1ll*sg[la]*inv%mod2,1ll*(n/i)*(m/i),mod2)*ans%mod2; 64 } 65 printf("%d ",ans); 66 } 67 return 0; 68 }