BZOJ4816 || 洛谷P3704 [SDOI2017]数字表格【莫比乌斯反演】

Time Limit: 50 Sec
Memory Limit: 128 MB

Description

Doris刚刚学习了fibonacci数列。用f[i]表示数列的第i项，那么
f[0]=0,f[1]=1
f[n]=f[n-1]+f[n-2],n>=2
Doris用老师的超级计算机生成了一个n×m的表格，第i行第j列的格子中的数是f[gcd(i,j)]，其中gcd(i,j)表示i,
j的最大公约数。Doris的表格中共有n×m个数，她想知道这些数的乘积是多少。答案对10^9+7取模。

Input

有多组测试数据。
第一个一个数T，表示数据组数。
接下来T行，每行两个数n,m
T<=1000,1<=n,m<=10^6

Output

输出T行，第i行的数是第i组数据的结果

---

题目分析

令$g(d)$表示$1\le x\le n,1\le y\le m$内有多少个$(x,y)$满足$gcd(x,y)=d$
那么可以转化目标表达式
$$\prod _{i=1}^n\prod _{j=1}^{m}f(gcd(i,j))=\prod _{d=1}^{min(n,m)}f(d{)}^{g(d)}$$

这个$g(d)$几乎是莫比乌斯反演考烂了的套路了

$$g(d)=\sum _{d|k}\mu (\frac{k}{d})\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor =\sum _{t=1}^{min(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )}\mu (t)\lfloor \frac{n}{td}\rfloor \lfloor \frac{m}{td}\rfloor$$

带入原式
$$\prod _{d=1}^{min(n,m)}f(d{)}^{{\sum }_{d|k}\mu (\frac{k}{d})\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor }=\prod _{d=1}^{min(n,m)}\prod _{d|k}f(d{)}^{\mu (\frac{k}{d})\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor }$$

交换$k,d$的枚举顺序

$$\prod _{k=1}^{min(n,m)}\prod _{d|k}f(d{)}^{\mu (\frac{k}{d})\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor }$$

看着似乎有点复杂，把中间那坨抽出来表示
令$h(k)={\prod }_{d|k}f(d{)}^{\mu (\frac{k}{d})}$，那么表达式为${\prod }_{k=1}^{min(n,m)}h(k{)}^{\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor }$

预处理$h(k)$的前缀积，每个询问用整除分块优化后复杂度$O(\sqrt{n}logn)$，这里log是快速幂计算的复杂度
$h(k)$可以用枚举倍数的方法预处理，复杂度$O(nlogn)$，莫比乌斯函数预处理$O(n)$
最后总复杂度$O(nlogn+n+T\sqrt{n}logn)$

#include<iostream>
#include<vector>
#include<algorithm>
#include<queue>
#include<cstring>
#include<cstdio>
using namespace std;
typedef long long lt;
  
int read()
{
    int f=1,x=0;
    char ss=getchar();
    while(ss<'0'||ss>'9'){if(ss=='-')f=-1;ss=getchar();}
    while(ss>='0'&&ss<='9'){x=x*10+ss-'0';ss=getchar();}
    return f*x;
}
 
const lt mod=1e9+7;
const int maxn=1000010;
int T,miu[maxn];
int vis[maxn],prim[maxn],cnt;
lt f[maxn],h[maxn];
lt sum[maxn];
 
lt qpow(lt a,lt k)
{
    lt res=1;
    while(k>0){
        if(k&1) res=(res*a)%mod;
        a=(a*a)%mod; k>>=1;
    }
    return res;
}
 
void calc(int n)
{
    miu[1]=1;
    for(int i=2;i<=n;++i)
    {
        if(!vis[i]) prim[++cnt]=i,miu[i]=-1;
        for(int j=1;j<=cnt;++j)
        {
            if(i*prim[j]>n) break;
            vis[i*prim[j]]=1;
            if(i%prim[j]==0) break;
            else miu[i*prim[j]]=-miu[i];
        }
    }
     
    f[0]=0; f[1]=1;
    for(int i=2;i<=n;++i)
    f[i]=(f[i-1]+f[i-2])%mod;
     
    for(int i=0;i<=n;++i) h[i]=1;
    for(int i=1;i<=n;++i)
    {
        lt inv=1ll*qpow(f[i],mod-2);//先计算f[i]逆元，下面枚举遇到miu为-1时不用重复计算
        for(int j=i;j<=n;j+=i)
        {
            if(miu[j/i]==-1) h[j]=h[j]*inv%mod;
            else if(miu[j/i]==1) h[j]=h[j]*f[i]%mod;
        }
    }
     
    for(int i=1;i<=n;++i)
    h[i]=h[i-1]*h[i]%mod;
}
 
lt query(int n,int m)
{
    lt res=1; int lim=min(n,m);
    for(int ll=1,rr;ll<=lim;ll=rr+1)
    {
        rr=min(n/(n/ll),m/(m/ll));
        lt Mul=h[rr]*qpow(h[ll-1],mod-2)%mod , ki=1ll*(n/ll)*(m/ll);
        res=res*qpow(Mul,ki)%mod;
    }
    return res;
}
 
int main()
{
    T=read(); calc(maxn-5);
    while(T--)
    {
        int n=read(),m=read();
        printf("%lld\n",query(n,m));
    }
    return 0;
}