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

点此看题面

大致题意:\(\prod_{i=1}^n\prod_{j=1}^mf(gcd(i,j))\)

推式子

首先,按照套路我们枚举\(gcd\),得到:

\[\prod_{d=1}^{min(n,m)}f(d)^{\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}[gcd(i,j)=1]}\]

根据\(\sum_{d|n}\mu(d)=[n=1]\),我们可以再转化得到:

\[\prod_{d=1}^{min(n,m)}f(d)^{\sum_{p=1}^{\frac{min(n,m)}d}\mu(p)\lfloor\frac n{dp}\rfloor\lfloor\frac m{dp}\rfloor}\]

然后我们把对\(p\)的枚举移出来,得到:

\[\prod_{d=1}^{min(n,m)}\sum_{p=1}^{\lfloor\frac{min(n,m)}d\rfloor}f(d)^{\mu(p)\lfloor\frac n{dp}\rfloor\lfloor\frac m{dp}\rfloor}\]

这个式子如果直接用除法分块套除法分块,是\(O(n)\)的,但多组数据就\(TLE\)了。

所以我们要进一步优化,枚举\(D=dp\),就可以得到这样一个式子:

\[\prod_{D=1}^{min(n,m)}(\prod_{d|D}f(d)^{\mu({\frac Dd})})^{\lfloor\frac nD\rfloor\lfloor\frac mD\rfloor}\]

其中\(\prod_{d|D}f(d)^{\mu({\frac Dd})}\),我们可以枚举\(d\)以及其倍数\(D\)预处理。

再用除法分块求就是\(O(\sqrt n)\)的了。

代码

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 1000000
#define X 1000000007
#define XSum(x,y) ((x)+(y)>=X?(x)+(y)-X:(x)+(y))
#define Qinv(x) Qpow(x,X-2)
using namespace std;
int n,m,f[N+5],g[N+5],ig[N+5];
I int Qpow(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}//快速幂
template<int SZ> class LinearSiever//线性筛
{
    private:
        int Pt,P[SZ+5];
    public:
        int mu[SZ+5];
        I LinearSiever()
        {
            RI i,j;for(mu[1]=1,i=2;i<=SZ;++i)
            {
                !P[i]&&(mu[P[++Pt]=i]=-1);
                for(j=1;j<=Pt&&1LL*i*P[j]<=SZ;++j)
                    if(P[i*P[j]]=1,i%P[j]) mu[i*P[j]]=-mu[i];else break;
            }
        }
};LinearSiever<N> L;
int main()
{
    RI i,j,p;for(f[1]=g[1]=1,i=2;i<=N;++i) f[i]=XSum(f[i-2],f[i-1]),g[i]=1;//预处理斐波那契数
    for(i=2;i<=N;++i) for(p=Qinv(f[i]),j=1;1LL*i*j<=N;++j)//预处理
        g[i*j]=1LL*g[i*j]*(L.mu[j]?(~L.mu[j]?f[i]:p):1)%X;
    for(g[0]=ig[0]=i=1;i<=N;++i) g[i]=1LL*g[i]*g[i-1]%X,ig[i]=Qinv(g[i]);//求前缀积及逆元
    RI Tt,t,l,r,ans;scanf("%d",&Tt);W(Tt--)
    {
        for(scanf("%d%d",&n,&m),ans=l=1,t=min(n,m);l<=t;l=r+1)//除法分块
        {
            r=min(n/(n/l),m/(m/l)),
            ans=1LL*ans*Qpow(1LL*g[r]*ig[l-1]%X,1LL*(n/l)*(m/l)%(X-1))%X;//计算答案
        }printf("%d\n",ans);
    }return 0;
}
posted @ 2019-08-06 17:42  TheLostWeak  阅读(...)  评论(...编辑  收藏