P3704 [SDOI2017]数字表格

洛谷链接

题意：求\(\prod_{i=1}^n\prod_{j=1}^mf(gcd(i,j))\%(1e9+7),f(i)表示斐波那契数列，其中f(0)=0,f(1)=1\)

开始化简：
\(\prod_{i=1}^n\prod_{j=1}^mf(gcd(i,j))=\)
\(\prod_{d=1}^nf(d)^{\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==d]}\)
对于幂上的式子
\(\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==d]=\)
\(\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}[gcd(i,j)==1]=\)
\(\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\sum_{x|gcd(i,j)}\mu(x)=\)
\(\sum_{x=1}^{\frac{n}{d}}\mu(x)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}} [x|gcd(i,j)]=\)
\(\sum_{x=1}^{\frac{n}{d}}\mu(x)\sum_{i=1}^{\frac{n}{dx}}\sum_{j=1}^{\frac{m}{dx}} [1|gcd(i,j)]=\)
\(\sum_{x=1}^{\frac{n}{d}}\mu(x)\frac{n}{dx}\frac{m}{dx}\)
代入原式
\(\prod_{d=1}^nf(d)^{\sum_{x=1}^{\frac{n}{d}}\mu(x)*\frac{n}{dx}*\frac{m}{dx}}\)

据以往操作，令\(T=dx\)，则式子可进一步化简为
\(\prod_{T=1}^n\prod_{d|T}f(d)^{\mu(\frac{T}{d})*\frac{n}{T}*\frac{m}{T}}\)
后面的\(\frac{n}{T}*\frac{m}{T}\)可以用整乘法在根号时间内解决，那麽我们令\(h(T)=\prod_{d|T}f(d)^{\mu(\frac{T}{d})}\),我们需要首先预处理出\(h(T)\)的乘积前缀和，然后在遇到除法时取逆元即可。

求成绩前缀和的复杂度为\(O(nlnn)\),求解复杂度为\(O(T\sqrt{n})\),题中\(n<=1e6\),可以解决。

总结：这个题是莫比乌斯反演在乘积上的应用，算是莫比乌斯反演的别样用法，同时乘法逆元的处理也很可以，展现了数学运算的奇妙。

int T;
int n,m;
ll f[maxn],g[maxn],h[maxn];

int p[maxn],cnt,mu[maxn];
bool vis[maxn];

ll qpow(ll a,ll b){
    ll sum=1;
    while(b){if(b&1)sum=1LL*sum*a%mod;a=1LL*a*a%mod;b>>=1;}
    return sum; 
}

void pre(){
    cnt=f[0]=g[0]=0;h[0]=h[1]=f[1]=g[1]=mu[1]=1;
    rep(i,2,n){
        h[i]=1LL;
        f[i]=(f[i-2]+f[i-1])%mod;
        g[i]=qpow(f[i],mod-2);
        if(!vis[i]){p[++cnt]=i;mu[i]=-1;}
        rep(j,1,cnt){
            if(1LL*p[j]*i>n) break;
            vis[p[j]*i]=1;
            if(i%p[j]) mu[i*p[j]]=-mu[i];
            else {mu[i*p[j]]=0;break;}
        }
    }
    rep(i,1,n){
        if(!mu[i])continue;
        for(int j=i;j<=n;j+=i){
            h[j]=(1LL*h[j]* ( (mu[i]==1)? f[j/i]:g[j/i]) )%mod;
        }
    }
    
    rep(i,2,n) h[i]*=h[i-1],h[i]%=mod;
}
int main(){ 
    T=read();
    n=1e6;
    pre();
    while(T--){
        n=read(),m=read();
        if(n>m)swap(n,m);
        ll ans=1,inv;
        for(ll l=1,r;l<=n;l=r+1){
            r= min(n/(n/l),m/(m/l));
            inv=h[r]*qpow(h[l-1],mod-2)%mod;
            ans = ans * qpow(inv,1LL*(n/l)*(m/l))%mod;
        }
        print((ans+mod)%mod);pts;
    }
    return 0;
}