【莫比乌斯反演】[SDOI2017]数字表格

原题链接

P3704 [SDOI2017]数字表格

题目大意

求：

$$\prod _{i=1}^n\prod _{j=1}^m{f}_{\gcd (i,j)}$$

其中 $f_x$ 表示斐波那契数列的第 $x$ 项，$f_0=0,f_1=1$。

题解

令 $\gcd (i,j)=d$，我们直接枚举他：

$$\prod _{d=1}^n\prod _{i=1}^n\prod _{j=1}^m{f}_d\times [\gcd (i,j)=d]$$

一个 $f_d$ 会被乘上若干次，表示如下：

$$\prod _{d=1}^n{f}_d^{\sum _{i=1}^n\sum _{j=1}^m[\gcd (i,j)=d]}$$

老套路处理一下右上角：

$$\prod _{d=1}^n{f}_d^{\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }[\gcd (i,j)=1]}$$

接下来把右上角单独拿出来看：

$$\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }[\gcd (i,j)=1]$$

使用 $[n=1]=\sum _{d\mid n}\mu (d)$ 可得：

$$\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }\sum _{k\mid \gcd (i,j)}\mu (k)$$

于是对他再改一下枚举顺序：

$$\sum _{k=1}^n\mu (k)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }[k\mid i\text{and}k\mid j]$$

这个时候直接去除内部的 $k$：

$$\sum _{k=1}^n\mu (k)\sum _{i=1}^{\lfloor \frac{n}{dk}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{dk}\rfloor }1$$

所以这一大块就是：

$$\sum _{k=1}^n\mu (k)\lfloor \frac{n}{dk}\rfloor \lfloor \frac{m}{dk}\rfloor$$

把他放回原式：

$$\prod _{d=1}^n{f}_d^{\sum _{k=1}^n\mu (k)\lfloor \frac{n}{dk}\rfloor \lfloor \frac{m}{dk}\rfloor }$$

令 $k^{\prime }=kd$，并且提出来：

$$\prod _{k=1}^n{\left(\prod _{d\mid k}{f}_d^{\mu (\frac{k}{d})}\right)}^{\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor }$$

括号内 $\mathcal{O}(n\log n)$ 预处理，外部 $\mathcal{O}(\sqrt{n})$ 解决。此外，由欧拉定理得：

$$a^{\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor }\equiv a^{\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor \mathrm{mod}\phi (m)}(\mathrm{mod}m)$$

于是我们可以把指数又缩小一点。

代码

//你 cnt=1 了吗？
#include<bits/stdc++.h>
#define int long long
#define mid ((l+r)>>1)
#define fir first
#define sec second
#define lowbit(i) (i&(-i))
using namespace std;
const int N=1e6+5;
const int inf=1e18;
struct edge{int to,nxt,l;};
inline int read(){
    char op=getchar();
    int w=0,s=1;
    while(op<'0'||op>'9'){
        if(op=='-') s=-1;
        op=getchar();
    }
    while(op>='0'&&op<='9'){
        w=(w<<1)+(w<<3)+op-'0';
        op=getchar();
    }
    return w*s;
}
//数论部分
const double pi=acos(-1);
const int mod=1e9+7;
int Mul(int a,int b){return (a*b)%mod;}
int Add(int a,int b){return (a+b)%mod;}
int Dec(int a,int b){return (a-b+mod)%mod;}
int Pow(int a,int k){
    int ans=1;
    while(k){
        if(k&1) ans=Mul(ans,a);
        a=Mul(a,a);
        k>>=1;
    }
    return ans;
}
int gcd(int x,int y){return y==0?x:gcd(y,x%y);}
int lcm(int x,int y){return x/gcd(x,y)*y;}
int inv(int x){return Pow(x,mod-2);}
void exgcd(int a,int b,int &x,int &y){
    if(b==0){
        x=1,y=0;
        return;
    }
    exgcd(b,a%b,x,y);
    int t=x;
    x=y;
    y=t-a/b*y;
}
int f[N],mu[N],invf[N],vis[N],p[N],cnt=0,ans[N];
void prework(){
    int n=N-5;
    f[0]=0,f[1]=1,mu[1]=1;
    for(register int i=2;i<=n;i++) f[i]=Add(f[i-1],f[i-2]);
    for(register int i=0;i<=n;i++) invf[i]=inv(f[i]),ans[i]=1;
    for(register int i=2;i<=n;i++){
        if(!vis[i]){
            p[++cnt]=i;
            mu[i]=-1;
        }
        for(register int j=1;j<=cnt&&p[j]*i<=n;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
    }
    for(register int i=1;i<=n;i++){
        for(register int j=i;j<=n;j+=i){
            if(mu[j/i]==1) ans[j]=Mul(ans[j],f[i]);
            if(mu[j/i]==-1) ans[j]=Mul(ans[j],invf[i]);
        }
    }
    for(register int i=1;i<=n;i++) ans[i]=Mul(ans[i-1],ans[i]);
}
int Sqrt(int n,int m){
    int tot=1;
    for(register int l=1,r;l<=min(n,m);l=r+1){
        r=min(n/(n/l),m/(m/l));
        tot=Mul(tot,Pow(Mul(ans[r],inv(ans[l-1])),(n/l)*(m/l)%(mod-1)));
    }
    return tot;
}
signed main(){
    int T=read();
    prework();
    while(T--){
        int n=read(),m=read();
        printf("%lld\n",(Sqrt(n,m)+mod)%mod);
    }
}