BZOJ 4816 [SDOI2017]数字表格 (莫比乌斯反演)

题目大意：略

洛谷传送门

第一次不看题解做出来反演题，我好菜啊

最初推的式子里是把矩阵快速幂然后欧拉降幂，然后发现弄出来的矩阵不能直接相乘..

但这依然改变不了反演玩的是套路的事实

以下默认$n \leq m$

题目让我们求$\prod\limits_{i=1}^{n} \prod\limits_{j=1}^{m}f[gcd(i,j)]$

显然直接统计$f[gcd(i,j)]$很困难，我们每次指定一个值k，求有多少个$gcd(i,j)==k$，那么$k$的贡献就是$f[k]^{ans}$

变形$\prod\limits_{k=1}^{n} f[k]^{\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{m}[gcd(i,j)==k]}$

我们只需要处理出$\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{m}[gcd(i,j)==k]$的数量即可，这是一道经典的反演题

再利用常规套路，$\varepsilon =1*\mu$变形，$\sum\limits_{i=1}^{\left \lfloor \frac{n}{k} \right \rfloor} \sum\limits_{j=1}^{\left \lfloor \frac{m}{k} \right \rfloor}\sum\limits_{d|gcd(i,j)}\mu(d)$

指数最后变成了$\sum\limits_{d=1}^{\left \lfloor \frac{n}{k} \right \rfloor}\mu(d) \left \lfloor \frac{n}{kd} \right \rfloor \left \lfloor \frac{m}{kd} \right \rfloor$

原式$\prod\limits_{k=1}^{n} f[k]^{\sum\limits_{d=1}^{\left \lfloor \frac{n}{k} \right \rfloor}\mu(d) \left \lfloor \frac{n}{kd} \right \rfloor \left \lfloor \frac{m}{kd} \right \rfloor}$

还是常规套路提取出$Q=kd$，变形$\prod\limits_{Q=1}^{n}  \prod\limits_{d|Q}f[d]^{\mu(\frac{Q}{d}) \left \lfloor \frac{n}{Q} \right \rfloor \left \lfloor \frac{m}{Q} \right \rfloor}$

令$g(Q)=\prod\limits_{d|Q} f[d]^{\mu(\frac{Q}{d})}$，显然它的前$n$项可以在$O(nlogn)$内处理出

每次询问都需要整除分块，还需要快速幂。总时间$O(T\sqrt{n}logn+nlogn)$

#include <cstdio>
#include <cstring>
#include <algorithm>
#define N1 1010000
#define ll long long
#define dd double
#define maxn 1000000
using namespace std;

const ll mod2=1000000007;
int T,n,m;
int use[N1],pr[N1],mu[N1],cnt;
int qpow(ll x,ll y,const ll &mod)
{
    int ans=1;
    while(y){
        if(y&1) ans=ans*x%mod;
        x=x*x%mod; y>>=1;
    }return ans;
}
void exgcd(ll a,ll b,ll &x,ll &y)
{
    if(!b){ x=1; y=0; return;}
    else{ exgcd(b,a%b,x,y); ll t=x; x=y; y=t-a/b*y; }
}
ll f[N1],g[N1],sg[N1];
void init()
{
    int i,j; ll x,y,inv; mu[1]=1;
    for(i=2;i<=maxn;i++)
    {
        if(!use[i]){ pr[++cnt]=i; mu[i]=-1; }
        for(j=1;j<=cnt&&i*pr[j]<=maxn;j++)
        {
            use[i*pr[j]]=1;
            if(i%pr[j]){ mu[i*pr[j]]=-mu[i]; }
            else{ mu[i*pr[j]]=0; break; }
        }
    }
    for(j=1;j<=maxn;j++) g[j]=1;
    for(f[1]=1,i=2;i<=maxn;i++)
    {
        f[i]=(f[i-1]+f[i-2])%mod2;
        exgcd(f[i],mod2,inv,y); inv=(inv%mod2+mod2)%mod2;
        for(j=1;j*i<=maxn;j++)
            if(mu[j]) g[i*j]=1ll*g[i*j]*(mu[j]>0?f[i]:inv)%mod2;
    }
    for(sg[0]=1,i=1;i<=maxn;i++) sg[i]=1ll*sg[i-1]*g[i]%mod2;
}

int main()
{
    scanf("%d",&T);
    init();
    int i,j,la,ans=0; ll inv,y;
    while(T--)
    {
        scanf("%d%d",&n,&m); if(n>m) swap(n,m); 
        for(i=1,ans=1;i<=n;i=la+1)
        {
            la=min(n/(n/i),m/(m/i));
            exgcd(sg[i-1],mod2,inv,y); inv=(inv%mod2+mod2)%mod2;
            ans=1ll*qpow(1ll*sg[la]*inv%mod2,1ll*(n/i)*(m/i),mod2)*ans%mod2;
        }
        printf("%d\n",ans);
    }
    return 0;
}