Luogu3704 数字表格 | Arahc's Home

题面

简化题意

多组数据，每组数据给定 $n,m$ ，求：

\prod_{i=1}^n \prod_{j=1}^m fib_{\gcd(i,j)}

其中 $fib_i$ 表示斐波那契第 $i$ 项， $fib_1=fib_2=1$ 。

$n,m\leq 10^6,T\leq 10^3$ 。

题解

不失一般性，设 $n<m$ 。

不妨考虑数列第 $i$ 项能产生多少贡献，答案可以写为：

\prod_{k=1}^n fib_i^{\sum_{i=0}^n \sum_{j=0}^m [\gcd(i,j)=k]}

先把指数提出来处理一下：

\sum_{i=0}^n \sum_{j=0}^m [\gcd(i,j)=k]

莫反老套路了吧，中间过程省略：

\prod_{d=1}^{\lfloor \frac{n}{k} \rfloor} \mu(d) \left\lfloor \frac{n}{kd} \right\rfloor \left\lfloor \frac{m}{kd} \right\rfloor

带回去：

\prod_{k=1}^n fib_k^{\sum_{d=1}^{\lfloor \frac{n}{k} \rfloor} \mu(d) \lfloor \frac{n}{kd} \rfloor \lfloor \frac{m}{kd} \rfloor}

然后考虑把指数整理一下：

\prod_{k=1}^n \left( fib_k^{\sum_{d=1}^{\lfloor \frac{n}{k} \rfloor} \mu(d)} \right)^{\lfloor \frac{n}{kd} \rfloor \lfloor \frac{m}{kd} \rfloor}

里面的内容，实际上就是 $\sum_{d=1}^{\lfloor \frac{n}{k} \rfloor} \mu(d)$ 这么多个 $fib_k$ 相乘，莫反老套路，设 $p=kd$ ：

\prod_{k=1}^n \left( \prod_{k\mid p} fib_k^{\mu(\frac{p}{k})} \right)^{\lfloor \frac{n}{p} \rfloor \lfloor \frac{m}{p} \rfloor}

分别考虑括号内和括号外，对于括号外，数论分块即可。

括号内的东西，设：

f(n)=\prod_{k\mid n} fib_k^{\mu(\frac{n}{k})}

直接暴力预处理，用埃式筛即可。

代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int max_n=1000006,N=1000000,mod=1000000007;
inline int read(){
    int x=0;bool w=0;char c=getchar();
    while(c<'0' || c>'9') w|=c=='-',c=getchar();
    while(c>='0' && c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return w?-x:x;
}
inline void write(int x){
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10^48);
}
inline int mi(int a,int p=mod-2){
    int res=1;
    while(p){
        if(p&1) res*=a,res%=mod;
        a*=a,a%=mod,p>>=1;
    }
    return res;
}

int n,m,f[max_n],fib[max_n],iv[max_n];

bool isp[max_n];
int pr[max_n],cnt,mu[max_n];
inline void init(int n){
    memset(isp,1,sizeof(isp)),isp[1]=0,mu[1]=1;
    fib[1]=fib[2]=iv[1]=iv[2]=1;
    for(register int i=3;i<=n;++i)
        fib[i]=(fib[i-1]+fib[i-2])%mod,
        iv[i]=mi(fib[i]);
    for(register int i=2;i<=n;++i){
        if(isp[i]){
            pr[++cnt]=i,
            mu[i]=-1;
        }
        for(register int j=1;j<=cnt && i*pr[j]<=n;++j){
            isp[i*pr[j]]=0;
            if(i%pr[j]==0){
                mu[i*pr[j]]=0;
                break;
            }
            mu[i*pr[j]]=-mu[i];
        }
    }
    for(register int i=0;i<=n;++i)
        f[i]=1;
    for(register int i=1;i<=n;++i) if(mu[i])
        for(register int j=i;j<=n;j+=i)
            f[j]=f[j]*(mu[i]>0?fib[j/i]:iv[j/i])%mod;
    for(register int i=2;i<=n;++i)
        f[i]=f[i]*f[i-1]%mod;
}

signed main(){
    init(N);
    for(register int T=read();T>0;--T){
        n=read(),m=read();
        if(n>m) n^=m^=n^=m;
        int ans=1;
        for(register int l=1,r;l<=n;l=r+1){
            r=min(n/(n/l),m/(m/l));
            int tmp=f[r]*mi(f[l-1])%mod;
            ans=ans*mi(tmp,(n/l)*(m/l)%(mod-1))%mod;
        }
        write((ans+mod)%mod),putchar('\n');
    }
    return 0;
}