Description
Doris刚刚学习了(fibonacci)数列。用(f[i])表示数列的第(i)项，那么
[egin{align}f[0]&=0\f[1]&=1\f[n]&=f[n-1]+f[n-2],n>=2end{align}
]
Doris用老师的超级计算机生成了一个(n imes m)的表格，第(i)行第(j)列的格子中的数是(f[mathrm{gcd}(i,j)])，其中(mathrm{gcd}(i,j))表示(i),(j)的最大公约数。Doris的表格中共有(n imes m)个数，她想知道这些数的乘积是多少。答案对(10^9+7)取模。
Input
有多组测试数据。第一个一个数(T)，表示数据组数。
接下来(T)行，每行两个数(n,m)
(T leq 1000),(1 leq n),(m leq 10^6)
Output
输出(T)行，第(i)行的数是第(i)组数据的结果
Sample Input
3
2 3
4 5
6 7
Sample Output
1
6
960
题解
要求的东西是
[ans=prod_{i=1}^nprod_{j=1}^mf[mathrm{gcd}(i,j)]
]
假设(n<m)，根据常见套路，枚举约数
[ans=prod_{d=1}^nf(d)^{sum_{i=1}^nsum_{j=1}^m[mathrm{gcd}(i,j)==d]}
]
把指数拿出来化简
[egin{align}
&sum_{i=1}^nsum_{j=1}^m[mathrm{gcd}(i,j)==d]\
=&sum_{i=1}^{left lfloor frac{n}{d}
ight
floor}sum_{j=1}^{left lfloor frac{m}{d}
ight
floor}[mathrm{gcd}(i,j)==1]\
=&sum_{i=1}^{left lfloor frac{n}{d}
ight
floor}sum_{j=1}^{left lfloor frac{m}{d}
ight
floor}sum_{k|mathrm{gcd}(i,j)}mu(k)\
=&sum_{k=1}^{left lfloor frac{n}{d}
ight
floor}mu(k)left lfloor frac{n}{kd}
ight
floorleft lfloor frac{m}{kd}
ight
floor
end{align}
]
令(T=kd)，代入上式
[sum_{k=1}^{left lfloor frac{n}{d}
ight
floor}mu(k)left lfloor frac{n}{kd}
ight
floorleft lfloor frac{m}{kd}
ight
floor=sum_{T=1}^nmu(frac{T}{d})left lfloor frac{n}{T}
ight
floorleft lfloor frac{m}{T}
ight
floor\
]
把(T)放到外面枚举
[ans=prod_{T=1}^nprod_{d|T}f(d)^{mu(frac{T}{d})left lfloor frac{n}{T}
ight
floorleft lfloor frac{m}{T}
ight
floor}
]
令函数(g(x)=prod_{d|x}f(d)^{mu(frac{x}{d})})
[ans=prod_{T=1}^ng(T)^{left lfloor frac{n}{T}
ight
floorleft lfloor frac{m}{T}
ight
floor}
]
只要我们求出了(g(x))的前缀积以及对应的逆元就可以整除分块了。
只要预处理出莫比乌斯函数和(f)以及(f)的逆元，(g(x))就可以通过枚举每个数的筛它的倍数求出，复杂度是调和级数。
问题就是如何求所有(f)以及对应的逆元(invf)。
这里我们需要一个辅助的东西(mulf)，表示(f)的前缀积，即(mulf(x)=prod_{i=1}^x f(i)).
假设我们知道了(mulf(i+1))的逆元(invm(i+1))，那么有
[egin{align}
invm(i+1) &equiv invm(i) imes invf(i+1) (mathrm{mod} P)\
invm(i)&equiv invm(i+1) imes f(i+1)(mathrm{mod} P)
end{align}
]
而(invf(i)equiv invm(i)*mulf(i-1)(mathrm{mod} P)),这样我们就可以求出(f)以及(invf)了。
用同样的方法，我们也可以求出(g(x))的前缀积以及对应的逆元。
代码
#include<bits/stdc++.h>
#define MAXN 1000010
#define P 1000000007
#define LL long long
namespace IO{
char buf[1<<15],*fs,*ft;
inline char gc(){return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;}
inline int qr(){
int x=0,rev=0,ch=gc();
while(ch<'0'||ch>'9'){if(ch=='-')rev=1;ch=gc();}
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=gc();}
return rev?-x:x;}
}using namespace IO;
using namespace std;
LL f[MAXN],invm[MAXN],invf[MAXN],g[MAXN],invg[MAXN],mulf[MAXN],mulg[MAXN],ans;
int p[MAXN],cnt,u[MAXN],T,N,M;
bool np[MAXN];
inline LL qp(LL x,LL y){
LL ret=1;
while(y){
if(y&1)ret=ret*x%P;
x=x*x%P;y>>=1;
}
return ret;
}
int main(){
#ifndef ONLINE_JUDGE
freopen("product.in","r",stdin);
freopen("product.out","w",stdout);
#endif
f[0]=0;
f[1]=u[1]=np[1]=g[1]=mulf[1]=mulg[1]=invm[0]=1;
for(int i=2;i<=1000000;i++){
f[i]=(f[i-1]+f[i-2])%P;
mulf[i]=mulf[i-1]*f[i]%P;
g[i]=1;
if(!np[i])p[++cnt]=i,u[i]=-1;
for(int j=1;j<=cnt;j++){
int t=p[j]*i;
if(t>1000000)break;
np[t]=1;
if(i%p[j]==0){u[t]=0;break;}
u[t]=-u[i];
}
}
invm[1000000]=qp(mulf[1000000],P-2);
invf[1000000]=qp(f[1000000],P-2);
for(int i=999999;i>=1;i--){
invm[i]=invm[i+1]*f[i+1]%P;
invf[i]=mulf[i-1]*invm[i]%P;
}
for(int i=2;i<=1000000;i++){
for(int j=1,k;(k=j*i)<=1000000;j++){
if(u[j]==1)g[k]=g[k]*f[i]%P;
else if(u[j]==-1)g[k]=g[k]*invf[i]%P;
}
mulg[i]=mulg[i-1]*g[i]%P;
}
invg[1000000]=qp(mulg[1000000],P-2);
for(int i=999999;i>=0;i--)invg[i]=invg[i+1]*g[i+1]%P;
T=qr();
while(T--){
N=qr();M=qr();ans=1;
if(N>M)swap(N,M);
for(int i=1,j;i<=N;i=j+1){
j=min(N/(N/i),M/(M/i));
ans=ans*qp(mulg[j]*invg[i-1]%P,1LL*(N/i)*(M/i)%(P-1))%P;
}
printf("%lld
",ans);
}
return 0;
}