BZOJ4816: [Sdoi2017]数字表格

【传送门：BZOJ4816】

---

简要题意：

　　设f[i]为斐波那契数列的第i项（f[0]=0,f[1]=1）

　　有一个n*m的表格，第i行第j列的格子上的数为f[gcd(i,j)]，求出n*m的所有格子的乘积

---

题解：

　　莫比乌斯反演

　　实际上就是求$\prod_{i=1}^{n}\prod_{j=1}^{m}f[gcd(i,j)]$

　　按照套路先枚举d=gcd(i,j)，然后根据反演公式得到$$\prod_{d=1}^{n}f[d]^{\sum_{i=1}^{\frac{n}{d}}\mu(i)*\frac{n}{d*i}*\frac{m}{d*i}}$$

　　还不能直接做，时间复杂度有点高，我们设x=di，设$g(d)=\prod_{d|x}f[d]^{\mu(\frac{x}{d})}$，得到原式等于$\prod_{x=1}^{n}g(x)^{\frac{n}{x}*\frac{m}{x}}$

　　求g(d)就用更新倍数的方法求，然后直接整除分块就行了

　　PS：BZOJ毒瘤卡常，求逆元时需要O(n)求逆元

---

参考代码：

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long LL;
LL Mod=1e9+7;
LL p_mod(LL a,LL b)
{
    LL ans=1;
    while(b!=0)
    {
        if(b%2==1) ans=ans*a%Mod;
        a=a*a%Mod;b/=2;
    }
    return ans;
}
int miu[1100000],prime[1100000],v[1100000];
LL f[1100000],nyf[1100000],nyg[1100000],sg[1100000];
LL sf[1100000],g[1100000];
void pre(int n)
{
    miu[1]=1;int m=0;
    for(int i=2;i<=n;i++)
    {
        if(v[i]==0)
        {
            v[i]=i;
            prime[++m]=i;
            miu[i]=-1;
        }
        for(int j=1;j<=m;j++)
        {
            if(prime[j]>n/i||prime[j]>v[i]) break;
            v[i*prime[j]]=prime[j];
            if(i%prime[j]==0) miu[i*prime[j]]=0;
            else miu[i*prime[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;
    sf[0]=1;
    for(int i=1;i<=n;i++) sf[i]=sf[i-1]*f[i]%Mod;
    nyf[0]=nyf[1]=1;
    LL d=p_mod(sf[n],Mod-2);
    for(int i=n;i>=1;i--) nyf[i]=d*sf[i-1]%Mod,d=d*f[i]%Mod;
    for(int i=1;i<=n;i++) g[i]=1;
    for(int i=2;i<=n;i++)
    {
        for(int j=1;i*j<=n;j++)
        {
            LL d;
            if(miu[j]==1) d=f[i];
            else if(miu[j]==0) d=1;
            else d=nyf[i];
            g[i*j]=g[i*j]*d%Mod;
        }
    }
    sg[0]=1;
    for(int i=1;i<=n;i++) sg[i]=sg[i-1]*g[i]%Mod;
    nyg[n]=p_mod(sg[n],Mod-2);
    for(int i=n;i>=1;i--) nyg[i-1]=nyg[i]*g[i]%Mod;
}
int main()
{
    pre(1000000);
    int T,n,m;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d",&n,&m);
        if(n>m) swap(n,m);
        LL ans=1;
        for(int i=1,j;i<=n;i=j+1)
        {
            j=min(n/(n/i),m/(m/i));
            ans=ans*p_mod(sg[j]*nyg[i-1]%Mod,LL(n/i)*(m/i)%(Mod-1))%Mod;
        }
        printf("%lld\n",ans);
    }
    return 0;
}