【BZOJ4816】数字表格（SDOI2017）-莫比乌斯反演+数论分块

测试地址：数字表格
做法：本题需要用到莫比乌斯反演+数论分块。
本题要求的是：
ans=∏ni=1∏mj=1f(gcd(i,j))$ans=\prod _{i=1}^n\prod _{j=1}^{m}f(gcd(i,j))$
不妨设n<m$n<m$，照常转成枚举因数d$d$的形式：
ans=∏nd=1f(d)∑ni=1∑mj=1[gcd(i,j)=d]$ans=\prod _{d=1}^{n}f(d{)}^{\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=d]}$
上面那个幂中的式子是不是看起来非常熟悉，我们在很多题目中都推过这个式子，由莫比乌斯反演定理的第二种形式得：
∑ni=1∑mj=1[gcd(i,j)=d]=∑d|iμ(id)⌊ni⌋⌊mi⌋$\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=d]=\sum _{d|i}\mu (\frac{i}{d})\lfloor \frac{n}{i}\rfloor \lfloor \frac{m}{i}\rfloor$
把这个式子带回ans$ans$的式子中，有：
ans=∏nd=1f(d)∑d|iμ(id)⌊ni⌋⌊mi⌋$ans=\prod _{d=1}^{n}f(d{)}^{\sum _{d|i}\mu (\frac{i}{d})\lfloor \frac{n}{i}\rfloor \lfloor \frac{m}{i}\rfloor }$
=∏nd=1∏d|if(d)μ(id)⌊ni⌋⌊mi⌋$=\prod _{d=1}^n\prod _{d|i}f(d{)}^{\mu (\frac{i}{d})\lfloor \frac{n}{i}\rfloor \lfloor \frac{m}{i}\rfloor }$
交换d,i$d,i$的位置，有：
ans=∏ni=1∏d|if(d)μ(id)⌊ni⌋⌊mi⌋$ans=\prod _{i=1}^n\prod _{d|i}f(d{)}^{\mu (\frac{i}{d})\lfloor \frac{n}{i}\rfloor \lfloor \frac{m}{i}\rfloor }$
显然可以加个括号：
ans=∏ni=1(∏d|if(d)μ(id))⌊ni⌋⌊mi⌋$ans=\prod _{i=1}^n(\prod _{d|i}f(d{)}^{\mu (\frac{i}{d})}{)}^{\lfloor \frac{n}{i}\rfloor \lfloor \frac{m}{i}\rfloor }$
那么我们只要预处理出括号内的东西，就可以数论分块了，注意到这个式子显然可以O(nlogn)$O(n\log n)$求出（暴力枚举每个数的倍数更新贡献），而数论分块的时间复杂度为O(Tn−−√logn)$O(T\sqrt{n}\log n)$，所以可以通过此题。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1000000007;
int T;
ll n[1010],m[1010],maxn,f[1000010],inv[1000010];
ll mu[1000010],prime[1000010],F[1000010],Finv[1000010],prod[1000010],prodinv[1000010];
bool vis[1000010]={0};

void calc_mu()
{
    mu[1]=1;
    prime[0]=0;
    for(ll i=2;i<=maxn;i++)
    {
        if (!vis[i])
        {
            prime[++prime[0]]=i;
            mu[i]=-1;
        }
        for(ll j=1;j<=prime[0]&&i*prime[j]<=maxn;j++)
        {
            vis[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
}

ll power(ll a,ll b)
{
    ll s=1,ss=a;
    while(b)
    {
        if (b&1) s=s*ss%mod;
        ss=ss*ss%mod;b>>=1;
    }
    return s;
}

void init()
{
    scanf("%d",&T);
    maxn=0;
    for(int i=1;i<=T;i++)
    {
        scanf("%lld%lld",&n[i],&m[i]);
        if (n[i]>m[i]) swap(n[i],m[i]);
        maxn=max(maxn,n[i]);
    }

    calc_mu();
    f[0]=0,f[1]=1;
    inv[1]=1;
    for(ll i=2;i<=maxn;i++)
    {
        f[i]=(f[i-2]+f[i-1])%mod;
        inv[i]=power(f[i],mod-2);
    }

    for(ll i=1;i<=maxn;i++)
        F[i]=Finv[i]=1;
    prod[0]=prodinv[0]=1;
    for(ll i=1;i<=maxn;i++)
    {
        for(ll j=1;i*j<=maxn;j++)
        {
            if (mu[j]==1) F[i*j]=F[i*j]*f[i]%mod,Finv[i*j]=Finv[i*j]*inv[i]%mod;
            if (mu[j]==-1) F[i*j]=F[i*j]*inv[i]%mod,Finv[i*j]=Finv[i*j]*f[i]%mod;
        }
        prod[i]=prod[i-1]*F[i]%mod;
        prodinv[i]=prodinv[i-1]*Finv[i]%mod;
    }
}

void work()
{
    for(int i=1;i<=T;i++)
    {
        ll ans=1;
        for(ll j=n[i];j>=1;j=max(n[i]/(n[i]/j+1),m[i]/(m[i]/j+1)))
        {
            ll l=max(n[i]/(n[i]/j+1),m[i]/(m[i]/j+1))+1,r=j;
            ans=ans*power(prod[r]*prodinv[l-1]%mod,(n[i]/j)*(m[i]/j))%mod;
        }
        printf("%lld\n",ans);
    }
}

int main()
{
    init();
    work();

    return 0;
}