「题解」「SDOI2017」数字表格 product

SDOI2017 Round 1 Day 1 T1，一道 莫比乌斯反演 的题目。

题目链接：Luogu P3704/BZOJ 4816/LibreOJ 2000/SDOI2017 R1D1T1。

题目

题目描述

$T$ 组数据，每组数据有两个数 $n,m$，求
$$\prod^n _ {i=1}\prod^m _ {j=1}f _ {\gcd(i,j)}\pmod p$$

其中 $f _ i$ 表示斐波那契数列的第 $i$ 项，$p=10^9+7$。

数据范围

时空限制

题解

思路

显然看出
$$
\begin{aligned}
\texttt{ans}&\equiv\prod^n _ {i=1}\prod^m _ {j=1}f _ {\gcd(i,j)}\pmod p\\
&\equiv\prod _ {x=1}\prod^n _ {i=1}\prod _ {\gcd(i,j)=x,j\leq m}f _ {\gcd(i,j)}\pmod p\\
&\equiv\prod _ {x=1}f _ x^{\sum^n _ {i=1}\sum^m _ {j=1}[\gcd(i,j)=x]}\pmod p\\
&\equiv\prod _ {x=1}f _ x^{\sum^{\lfloor\frac{n}{x}\rfloor} _ {i=1}\sum^{\lfloor\frac{m}{x}\rfloor} _ {j=1}[\gcd(i,j)=1]}\pmod p\\
\end{aligned}
$$

然后又
$$
\sum^{\lfloor\frac{n}{x}\rfloor} _ {i=1}\sum^{\lfloor\frac{m}{x}\rfloor} _ {j=1}[\gcd(i,j)=1]=\sum^{\lfloor\frac{\max\{n,m\}}{d}\rfloor} _ {i=1}\mu(i)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor
$$

继而
$$
\texttt{ans}\equiv\prod^{\max\{n,m\}} _ {T=1}(\prod _ {d\mid T}f^{\mu(\frac{T}{d})} _ d)^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}\pmod p
$$

然后剩下部分可以直接算了。

代码

#include<bits/stdc++.h>
using namespace std;
#define reg register
#define MOD 1000000007
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++)
static char buf[100000],*p1=buf,*p2=buf;
inline int read(void){
    reg char ch=getchar();
    reg int res=0;
    while(ch<'0'||'9'<ch)ch=getchar();
    while('0'<=ch&&ch<='9')res=10*res+ch-'0',ch=getchar();
    return res;
}

const int MAXN=1000000+5;

inline int Add(reg int a,reg int b){
    reg int sum=a+b;
    return sum>=MOD?sum-MOD:sum;
}

inline int Mul(reg int a,reg int b){
    return 1ll*a*b%MOD;
}

inline int pow(reg int x,reg int exp,reg int mod){
    reg int res=1;
    while(exp){
        if(exp&1)
            res=1ll*res*x%mod;
        x=1ll*x*x%mod;
        exp>>=1;
    }
    return res;
}

inline int inv(reg int x,reg int mod){
    return pow(x,mod-2,mod);
}

bool vis[MAXN];
int tot,prime[MAXN];
int mu[MAXN];

int f[MAXN],invf[MAXN];
int F[MAXN];

inline void Init(reg int n){
    f[1]=invf[1]=F[0]=F[1]=1;
    mu[1]=1;
    for(reg int i=2;i<=n;++i){
        f[i]=(f[i-1]+f[i-2])%MOD;
        invf[i]=inv(f[i],MOD);
        F[i]=1;
        if(!vis[i]){
            mu[i]=-1;
            prime[++tot]=i;
        }
        for(reg int j=1;j<=tot&&i*prime[j]<=n;++j){
            vis[i*prime[j]]=true;
            if(i%prime[j]==0)
                break;
            mu[i*prime[j]]=-mu[i];
        }
    }
    for(reg int i=1;i<=n;++i){
        if(!mu[i])
            continue;
        for(reg int j=i;j<=n;j+=i)
            F[j]=Mul(F[j],mu[i]==1?f[j/i]:invf[j/i]);
    }
    for(reg int i=2;i<=n;++i)
        F[i]=Mul(F[i],F[i-1]);
    return;
}

int n,m;

int main(void){
    Init(1e6);
    reg int T=read();
    while(T--){
        n=read(),m=read();
        if(n>m)
            swap(n,m);
        reg int ans=1;
        for(reg int l=1,r;l<=n;l=r+1){
            r=min(n/(n/l),m/(m/l));
            reg int Inv=Mul(F[r],inv(F[l-1],MOD));
            ans=Mul(ans,pow(Inv,1ll*(n/l)*(m/l)%(MOD-1),MOD));
        }
        printf("%d\n",ans);
    }
    return 0;
}