100006. 【SDOI2017】数字表格

Description

Doris 刚刚学习了 fibnacci 数列。用 \(f[i]\) 表示数列的第 \(i\) 项, 那么,

\[\begin{aligned} f[0] &=0 \\ f[1] &=1 \\ f[n] &=f[n-1]+f[n-2], n \geq 2 \end{aligned} \]

Doris 用老师的超级计算机生成了一个 \(n \times m\) 的表格, 第 \(i\) 行第 \(j\) 列的格子中的数是 \(f[\operatorname{gcd}(i, j)]\), 其中 \(\operatorname{gcd}(i, j)\) 表示 \(i\) 与 \(j\) 的最大公约数。
Doris 的表格中共有 \(n \times m\) 个数, 她想知道这些数的乘积是多少。 这些数的乘积实在是太大了, 所以 Doris 只想知道乘积对 1000000007 取模后的结果。

Solution

喜闻乐见推式子。

令 \(n\le m\)。

\(\prod\limits_{i=1}^{n}\prod\limits_{j=1}^{m}f[\gcd(i,j)]=\prod\limits_{k=1}^nf_k^{\sum^{n}_{i=1}\sum_{j=1}^m[\gcd(i,j)=k]}\)

指数是典型的莫反式子。

\(\begin{aligned}\sum\limits^{n}_{i=1}\sum\limits_{j=1}^m[\gcd(i,j)=k]&=\sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{k}\rfloor}[gcd(i,j)=1]\\&=\sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum\limits_{d|\gcd(i,j)}\mu(d)\\&=\sum\limits_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor\end{aligned}\)

​

带回原式。枚举 \(T=kd\)。

\(\begin{aligned}\prod\limits_{k=1}^nf_k^{\sum^{n}_{i=1}\sum_{j=1}^m[\gcd(i,j)=k]}&=\prod\limits_{k=1}^{n}f_k^{\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor}\\&=\prod\limits_{T=1}^n(\prod_{k|T}f_k^{\mu(\frac{T}{k})})^{\lfloor\dfrac{n}{T}\rfloor\lfloor\dfrac{m}{T}\rfloor}\end{aligned}\)

令 \(F_n=\prod_{d|n}f_d^{\mu(\frac{n}{d})}\)。

则

\(\prod\limits_{T=1}^n(\prod_{k|T}f_k^{\mu(\frac{T}{k})})^{\lfloor\dfrac{n}{T}\rfloor\lfloor\dfrac{m}{T}\rfloor}=\prod_{T=1}^n F_T^{\lfloor\dfrac{n}{T}\rfloor\lfloor\dfrac{m}{T}\rfloor}\)

\(F\) 可以暴力预处理，枚举 \(d\) 和倍数，得到 \(n\)。

处理完 \(F\) 后，答案用整除分块求。

Code

#include<cstdio>
#include<algorithm>
#define N 1000000
#define mod 1000000007
#define ll long long
using namespace std;
int cnt,T,n,m,mu[N+5],p[N+5];
ll fa,fb,ans,f[N+5],fv[N+5];
bool bj[N+5];
ll ksm(ll x,ll y)
{
    ll res=1;
    while (y)
    {
        if (y&1) res=res*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return res;
}
int main()
{
    freopen("product.in","r",stdin);
    freopen("product.out","w",stdout);
    mu[1]=1;
    for (int i=2;i<=N;++i)
    {
        if (!bj[i]) p[++cnt]=i,mu[i]=-1;
        for (int j=1;j<=cnt&&i*p[j]<=N;++j)
        {
            bj[i*p[j]]=true;
            if (i%p[j]==0)
            {
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
    }
    for (int i=1;i<=N;++i)
        f[i]=fv[i]=1;
    fa=1;fb=0;
    for (int i=1;i<=N;++i)
    {
        fb=(fa+fb)%mod;
        fa=(fb-fa+mod)%mod;//计算当前斐波那契值
        ll x[3];
        x[0]=ksm(fb,mod-2);x[1]=1;x[2]=fb;
        for (int j=i,k=1;j<=N;j+=i,++k)
        {
            f[j]=f[j]*x[1+mu[k]]%mod;
            fv[j]=fv[j]*x[1-mu[k]]%mod;//记录逆元
        }
    }
    f[0]=fv[0]=1;
    for (int i=1;i<=N;++i)
    {
        f[i]=f[i-1]*f[i]%mod;
        fv[i]=fv[i-1]*fv[i]%mod;
    }
    scanf("%d",&T);
    while (T--)
    {
        scanf("%d%d",&n,&m);
        if (n>m) swap(n,m);
        ans=1;
        for (int l=1,r;l<=n;l=r+1)
        {
            r=min(n/(n/l),m/(m/l));
            ans=ans*ksm(f[r]*fv[l-1]%mod,(ll)(n/l)*(ll)(m/l))%mod;
        }
        printf("%lld\n",ans);
    }
    return 0;
}