题解 P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题

· · 题解

upd on 2021.8.7：修了个公式的问题

洛谷题面传送门

如果 Latex 炸了就到这里查看

一道究极恶心的毒瘤六合一题，式子推了我满满两面 A4 纸……

首先我们可以将式子拆成：

ans=\prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^C(\dfrac{ij}{\gcd(i,j)\gcd(i,j)})^{f(type)}

也就是说我们需要算出以下四项式子的值：

\prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^Ci^{f(type)} \prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^Cj^{f(type)} \prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^C(\dfrac{1}{\gcd(i,j)})^{f(type)} \prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^C(\dfrac{1}{\gcd(i,k)})^{f(type)}

显然前两项与后两项是等价的，因此我们只需算出：

f_1(type)=\prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^Ci^{f(type)}

与

f_2(type)=\prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^C(\gcd(i,j))^{f(type)}

即可求出答案。

考虑对 type 进行分类讨论，首先是 type=0，那么

f_1(0)=\prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^Ci

考虑每个 i 的贡献，稍微想想即可得到：

f_1(0)=(A!)^{BC}

然后是 f_2(0)，套路地枚举 \gcd(i,j) 有

\begin{aligned} f_2(0)&=\prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^C\gcd(i,j)\\ &=(\prod\limits_{d=1}^{\min(A,B)}d^{\sum\limits_{i=1}^A\sum\limits_{j=1}^B[\gcd(i,j)=d]})^C\\ &=(\prod\limits_{d=1}^{\min(A,B)}d^{\sum\limits_{i=1}^{A/d}\sum\limits_{j=1}^{B/d}[\gcd(i,j)=1]})^C\\ &=(\prod\limits_{d=1}^{\min(A,B)}d^{\sum\limits_{p}\mu(p)\lfloor\frac{A}{dp}\rfloor\lfloor\frac{B}{dp}\rfloor})^C\\ \end{aligned}

最右边那个 k 次方显然可以忽略掉不管它，最后求个快速幂即可。考虑对里面的 dp 进行二维的整除分块，那么答案的式子可以化为：

\begin{aligned} f_2(0)=(\prod\limits_{dp}g_1(dp)^{\lfloor\frac{A}{dp}\rfloor\lfloor\frac{B}{dp}\rfloor})^C \end{aligned}

其中

g_1(x)=\prod\limits_{d·p=x}d^{\mu(p)}

考虑整除分块，对于一段区间 [L,R]，满足 \forall x\in[L,R] 均有 \lfloor\dfrac{A}{x}\rfloor=\lfloor\dfrac{A}{L}\rfloor,\lfloor\dfrac{B}{x}\rfloor=\lfloor\dfrac{B}{L}\rfloor，我们这样计算它们的贡献：

\prod\limits_{i=L}^R(g_1(i)^{\lfloor\frac{A}{i}\rfloor\lfloor\frac{B}{i}\rfloor})^C=(\prod\limits_{i=L}^Rg_1(i))^{\lfloor\frac{A}{i}\rfloor\lfloor\frac{B}{i}\rfloor·C}

预处理前缀积即可 \mathcal O(1) 计算，时间复杂度 \mathcal O(\log n\sqrt{n})。

接下来是 type=1 的情况，个人感觉与 type=0 的情况大差不差，毕竟指数上都只与 i,j,k 本身而不涉及到它们的 \gcd 之类，只不过指数上枚举变量的次数稍微高了一点点，导致情况较于 type=0 略有一点点繁琐。

首先是 f_1(1)

f_1(1)=\prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^Ci^{ijk}

还是考虑每个 i 的贡献被累计了多少次：

\begin{aligned} f_1(0)&=(A!)^{\sum\limits_{j=1}^B\sum\limits_{k=1}^Cjk}\\ &=(A!)^{\frac{B(B+1)}{2}·\frac{C(C+1)}{2}} \end{aligned}

一波快速幂带走。

其次是 f_2(1)

\begin{aligned} f_2(1)&=\prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^C\gcd(i,j)^{ijk}\\ &=(\prod\limits_{d=1}^{\min(A,B)}d^{\sum\limits_{i=1}^A\sum\limits_{j=1}^Bij[\gcd(i,j)=d]})^{\sum\limits_{k=1}^Ck}\\ &=(\prod\limits_{d=1}^{\min(A,B)}d^{\sum\limits_{i=1}^{A/d}\sum\limits_{j=1}^{B/d}ijd^2[\gcd(i,j)=1]})^{\frac{C(C+1)}{2}}\\ &=(\prod\limits_{d=1}^{\min(A,B)}d^{\sum\limits_{p}\mu(p)d^2p^2·s(\lfloor\frac{A}{dp}\rfloor)s(\lfloor\frac{B}{dp}\rfloor)})^{\frac{C(C+1)}{2}}\\ \end{aligned}

其中 s(i)=\dfrac{i(i+1)}{2}

那么我们还是枚举 dp，按照 f_2(0) 的套路设一个 g_2(x)，定义如下：

g_2(x)=\prod\limits_{d·p=x}d^{d^2p^2\mu(p)}

那么考虑对 dp 进行整除分块，那么

f_2(1)=(\prod\limits_{dp=1}^{\min(A,B)}g_2(dp)^{s(\lfloor\frac{A}{dp}\rfloor)s(\lfloor\frac{B}{dp}\rfloor)})^{\frac{C(C+1)}{2}}

预处理 g_2(dp) 的前缀积然后对 dp 整除分块即可。

最后是 type=2 的情况

先考虑 f_1(2)

f_1(2)=\prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^Ci^{\gcd(i,j,k)}

上来先把 \gcd 反演掉：

\begin{aligned} f_1(2)&=\prod\limits_{d=1}^{\min(A,B,C)}\prod\limits_{i=1}^{\lfloor\frac{A}{d}\rfloor}(id)^{d·\sum\limits_{j=1}^{\lfloor\frac{B}{d}\rfloor}\sum\limits_{k=1}^{\lfloor\frac{C}{d}\rfloor}[\gcd(i,j,k)=1]}\\ &=\prod\limits_{d=1}^{\min(A,B,C)}\prod\limits_{i=1}^{\lfloor\frac{A}{d}\rfloor}(id)^{d·\sum\limits_{p\mid i}\mu(p)\lfloor\frac{B}{dp}\rfloor\lfloor\frac{C}{dp}\rfloor} \end{aligned}

把 p 提到外面来

\begin{aligned} f_1(2)&=\prod\limits_{d=1}^{\min(A,B,C)}\prod\limits_{p}(\prod\limits_{i\in[1,\lfloor\frac{A}{d}\rfloor]\&p\mid i}(id)^{d·\mu(p)})^{\lfloor\frac{B}{dp}\rfloor\lfloor\frac{C}{dp}\rfloor} \end{aligned}

然后按照套路枚举 \dfrac{i}{p}

\begin{aligned} f_1(2)&=\prod\limits_{d=1}^{\min(A,B,C)}\prod\limits_{p}(\prod\limits_{i=1}^{\lfloor\frac{A}{dp}\rfloor}(idp)^{d·\mu(p)})^{\lfloor\frac{B}{dp}\rfloor\lfloor\frac{C}{dp}\rfloor} \end{aligned}

然后枚举 dp，根据 \mu*i=\varphi 可知 \sum\limits_{d·p=x}d·\mu(p)=\varphi(x)，于是

f_1(2)=\prod\limits_{x=1}^{\min(A,B,C)}(\prod\limits_{i=1}^{\lfloor\frac{A}{x}\rfloor}(ix)^{\varphi(x)})^{\lfloor\frac{B}{x}\rfloor\lfloor\frac{C}{x}\rfloor}

外面的东西显然整除分块一下就好了，里面的东西

\prod\limits_{i=1}^{\lfloor\frac{A}{x}\rfloor}(ix)^{\varphi(x)}

显然等于

\begin{aligned} &\prod\limits_{i=1}^{\lfloor\frac{A}{x}\rfloor}i^{\varphi(x)}·x^{\varphi(x)}\\ =&((\lfloor\dfrac{A}{x}\rfloor)!)^{\varphi(x)}·(x^{\varphi(x)})^{\lfloor\frac{A}{x}\rfloor} \end{aligned}

然后套路地预处理 g_3=x^{\varphi(x)} 的前缀积，以及 \varphi(x) 的前缀和即可在整除分块的过程中 \mathcal O(1) 求出式子的值，注意 \varphi(x) 的前缀和应 \bmod(MOD-1) 而不是 \bmod MOD，因为 \varphi(x) 的前缀和作用在指数上。

时间复杂度 \sqrt{n}\log n

然后是最精神污染的一个式子：

f_2(2)=\prod\limits_{i=1}^A\prod\limits_{j=1}^B\prod\limits_{k=1}^C(\gcd(i,j))^{\gcd(i,j,k)}

按照 P1587 的套路，碰到两个 \gcd 咱们最好不要莽，要一个个反演，因此考虑先反演下面这个 \gcd：

\begin{aligned} f_2(2)&=\prod\limits_{d=1}^{\min(A,B)}d^{\sum\limits_{i=1}^{\lfloor\frac{A}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{B}{d}\rfloor}[\gcd(i,j)=1]\sum\limits_{k=1}^C\gcd(d,k)}\\ &=\prod\limits_{d=1}^{\min(A,B)}d^{\sum\limits_{p}^{\lfloor\frac{\min(A,B)}{d}\rfloor}\mu(p)\lfloor\frac{A}{dp}\rfloor\lfloor\frac{B}{dp}\rfloor\sum\limits_{k=1}^C\gcd(d,k)} \end{aligned}

发现后面有个 \sum\limits_{k=1}^C\gcd(d,k)，考虑对这个东西单独推个式子：

\begin{aligned} g_4(d,C)&=\sum\limits_{k=1}^C\gcd(d,k)\\ &=\sum\limits_{x\mid d}x·\sum\limits_{i=1}^{\lfloor\frac{C}{x}\rfloor}[i\perp\dfrac{d}{x}]\\ &=\sum\limits_{x\mid d}x·\sum\limits_{y\mid\frac{d}{x}}\mu(y)\lfloor\dfrac{C}{xy}\rfloor \end{aligned}

套路地枚举 xy=T 可得：

g_4(d,C)=\sum\limits_{T\mid d}\lfloor\dfrac{C}{T}\rfloor\sum\limits_{x\mid T}x\mu(\dfrac{T}{x})

这已经是本题中第二次看到这个式子了：

\sum\limits_{x\mid T}x\mu(\dfrac{T}{x})=\varphi(T)

于是

g_4(d,C)=\sum\limits_{T\mid d}\lfloor\dfrac{C}{T}\rfloor\varphi(T)

带回去

f_2(2)=\prod\limits_{d=1}^{\min(A,B)}d^{\sum\limits_{p}^{\lfloor\frac{\min(A,B)}{d}\rfloor}\mu(p)\lfloor\frac{A}{dp}\rfloor\lfloor\frac{B}{dp}\rfloor\sum\limits_{T\mid d}\lfloor\frac{C}{T}\rfloor\varphi(T)}

按照这里总结出来的套路，看到先枚举 i 再枚举 j\mid i 的求和/积式我们可以考虑交换求和/积的顺序，先枚举 j 再枚举 i，这样会出现下取整，就可以整除分块了。

因此考虑先枚举 T 再枚举 d，有：

f_2(2)=\prod\limits_{T=1}^{\min(A,B)}\prod\limits_{d=1}^{\lfloor\frac{\min(A,B)}{T}\rfloor}(dT)^{\sum\limits_{p}\mu(p)\lfloor\frac{A}{dTp}\rfloor\lfloor\frac{B}{dTp}\rfloor\lfloor\frac{C}{T}\rfloor\varphi(T)}

然后考虑拆开来：

f_2(2)=\prod\limits_{T=1}^{\min(A,B)}(T^{\varphi(T)\lfloor\frac{C}{T}\rfloor})^{\sum\limits_{d}\sum\limits_{p}\mu(p)\lfloor\frac{A}{dTp}\rfloor\lfloor\frac{B}{dTp}\rfloor}·\prod\limits_{d=1}^{\lfloor\frac{\min(A,B)}{T}\rfloor}(d^{\sum\limits_{p}\mu(p)\lfloor\frac{A}{dTp}\rfloor\lfloor\frac{B}{dTp}\rfloor})^{\varphi(T)\lfloor\frac{C}{T}\rfloor}

先考虑前面的式子：

(T^{\varphi(T)\lfloor\frac{C}{T}\rfloor})^{\sum\limits_{d}\sum\limits_{p}\mu(p)\lfloor\frac{A}{dTp}\rfloor\lfloor\frac{B}{dTp}\rfloor}

考虑枚举 dp=x，那么：

\begin{aligned} \text{原式}&=(T^{\varphi(T)\lfloor\frac{C}{T}\rfloor})^{\sum\limits_{x}\sum\limits_{p\mid x}\mu(p)\lfloor\frac{A}{Tx}\rfloor\lfloor\frac{B}{Tx}\rfloor}\\ &=(T^{\varphi(T)\lfloor\frac{C}{T}\rfloor})^{\epsilon(x)\lfloor\frac{A}{Tx}\rfloor\lfloor\frac{B}{Tx}\rfloor}\\ &=(T^{\varphi(T)\lfloor\frac{C}{T}\rfloor})^{\lfloor\frac{A}{T}\rfloor\lfloor\frac{B}{T}\rfloor} \end{aligned}

对 T 整除分块一下，预处理 \varphi(T) 的前缀和即可 \mathcal O(1) 算出。

然后是后面的式子（胜利就在眼前！）

\prod\limits_{T=1}^{\min(A,B)}\prod\limits_{d=1}^{\lfloor\frac{\min(A,B)}{T}\rfloor}(d^{\sum\limits_{p}\mu(p)\lfloor\frac{A}{dTp}\rfloor\lfloor\frac{B}{dTp}\rfloor})^{\varphi(T)\lfloor\frac{C}{T}\rfloor}

还是对 T 整除分块，然后枚举 dp=x，那么上面的式子可以写成：

\prod\limits_{T=1}^{\min(A,B)}\prod\limits_{x=1}^{\lfloor\frac{\min(A,B)}{T}\rfloor}((\prod\limits_{d\mid x}d^{\mu(p)})^{\lfloor\frac{A}{Tx}\rfloor\lfloor\frac{B}{Tx}\rfloor})^{\varphi(T)\lfloor\frac{C}{T}\rfloor}

发现最里面的括号的东西就是在求 f_2(0) 时用到的 g_1(x)，那么我们再套一个对 x 的整除分块即可。

根据整除分块里再套一个整除分块复杂度是 \sum\limits_{x,\exists k,s.t.\lfloor\frac{n}{k}\rfloor=x}\sqrt{x}=n^{0.75} 可知这一部分复杂度为 n^{0.75}\log n

于是这题就做完了，时间复杂度 n\log n+Tn^{0.75}\log n。

关于此题的常数问题，由于取模运算较多，可以使用快速取模优化常数，具体可见 chenxia25 神仙的这篇博客。

const int MAXV=1e5;
ll mod;
int pr[MAXV/5+5],prcnt=0,vis[MAXV+5],mu[MAXV+5],smu[MAXV+5],phi[MAXV+5];
int fac[MAXV+5],prd[MAXV+5],prd_inv[MAXV+5],prd_sq[MAXV+5],prd_sq_inv[MAXV+5];
int pre_ii[MAXV+5],inv_pre_ii[MAXV+5],pre_mul[MAXV+5],inv_pre_mul[MAXV+5];
int inv[MAXV+5],prd_phi[MAXV+5],prd_phi_inv[MAXV+5];
int sphi[MAXV+5];
__int128_t _base1=1,_base2=1;
inline int mol1(__int128_t x){return x-mod*(_base1*x>>64);}
inline int mol2(__int128_t x){return x-(mod-1)*(_base2*x>>64);}
int qpow(int x,int e){
    if(e<0) e+=mod-1;int ret=1;
    for(;e;e>>=1,x=mol1(1ll*x*x)) if(e&1) ret=mol1(1ll*ret*x);
    return ret;
}
int work(int x,int y){return (!y)?1:((~y)?x:inv[x]);}
void sieve(int n){
    for(int i=(inv[0]=inv[1]=1)+1;i<=n;i++) inv[i]=mol1(1ll*inv[mod%i]*(mod-mod/i));
    mu[1]=phi[1]=1;
    for(int i=(fac[0]=1);i<=n;i++) fac[i]=mol1(1ll*fac[i-1]*i);
    for(int i=2;i<=n;i++){
        if(!vis[i]) mu[i]=-1,pr[++prcnt]=i,phi[i]=i-1;
        for(int j=1;j<=prcnt&&pr[j]*i<=n;j++){
            vis[pr[j]*i]=1;
            if(i%pr[j]==0){phi[i*pr[j]]=phi[i]*pr[j];break;}
            mu[i*pr[j]]=-mu[i];phi[i*pr[j]]=phi[i]*phi[pr[j]];
        }
    }
    for(int i=1;i<=n;i++) smu[i]=smu[i-1]+mu[i];
    for(int i=0;i<=n;i++) prd[i]=prd_sq[i]=1;
    for(int i=1;i<=n;i++) for(int j=i;j<=n;j+=i) prd[j]=mol1(1ll*prd[j]*work(i,mu[j/i]));
    for(int i=1;i<=n;i++) prd[i]=mol1(1ll*prd[i-1]*prd[i]);
    for(int i=0;i<=n;i++) prd_inv[i]=qpow(prd[i],-1);
    for(int i=1;i<=n;i++) for(int j=i;j<=n;j+=i)
        prd_sq[j]=mol1(1ll*prd_sq[j]*
        qpow(qpow(i,mol2(1ll*i*i)),mol2(1ll*mu[j/i]*(j/i)*(j/i))));
    for(int i=1;i<=n;i++) prd_sq[i]=mol1(1ll*prd_sq[i-1]*prd_sq[i]);
    for(int i=0;i<=n;i++) prd_sq_inv[i]=qpow(prd_sq[i],-1);
    pre_ii[0]=1;for(int i=1;i<=n;i++) pre_ii[i]=mol1(1ll*pre_ii[i-1]*qpow(i,i));
    for(int i=0;i<=n;i++) inv_pre_ii[i]=qpow(pre_ii[i],-1);
    pre_mul[0]=1;for(int i=1;i<=n;i++) pre_mul[i]=mol1(1ll*pre_mul[i-1]*work(i,mu[i]));
    for(int i=0;i<=n;i++) inv_pre_mul[i]=qpow(pre_mul[i],-1);
    prd_phi[0]=1;for(int i=1;i<=n;i++) prd_phi[i]=mol1(1ll*prd_phi[i-1]*qpow(i,phi[i]));
    for(int i=0;i<=n;i++) prd_phi_inv[i]=qpow(prd_phi[i],-1);
    for(int i=1;i<=n;i++) sphi[i]=mol2(sphi[i-1]+phi[i]);
}
int calc1(int x,int y){
    int res=1;
    for(int l=1,r;l<=min(x,y);l=r+1){
        r=min(x/(x/l),y/(y/l));
        res=1ll*res*qpow(1ll*prd[r]*prd_inv[l-1]%mod,1ll*(x/l)*(y/l)%(mod-1))%mod;
    }
    return res;
}
int solve1(int a,int b,int c){
    int res=1ll*qpow(calc1(a,b),-c)*qpow(calc1(a,c),-b)%mod;
    res=1ll*res*qpow(fac[a],1ll*b*c%(mod-1))%mod;
    res=1ll*res*qpow(fac[b],1ll*a*c%(mod-1))%mod;
    return res;
}
ll get(int x){return mol2(1ll*x*(x+1)/2);}
int calc2(int x,int y){
    int res=1;
    for(int l=1,r;l<=min(x,y);l=r+1){
        r=min(x/(x/l),y/(y/l));
        res=1ll*res*qpow(1ll*prd_sq[r]*prd_sq_inv[l-1]%mod,
        1ll*get(x/l)*get(y/l)%(mod-1))%mod;
    }
    return res;
}
int solve2(int a,int b,int c){
    int res=1;
    res=1ll*res*qpow(pre_ii[a],(1ll*b*(b+1)>>1)%(mod-1))%mod;
    res=1ll*res*qpow(pre_ii[b],(1ll*a*(a+1)>>1)%(mod-1))%mod;
    res=1ll*res*qpow(calc2(a,b),mod-2)%mod;
    res=qpow(res,(1ll*c*(c+1)>>1)%(mod-1));
    res=1ll*res*qpow(calc2(a,c),-(1ll*b*(b+1)>>1)%(mod-1))%mod;
    return res;
}
int calc3(int x,int y,int z){
    int res=1;
    for(int l=1,r;l<=x;l=r+1){
        r=1e9;
        if(x/l) chkmin(r,x/(x/l));
        if(y/l) chkmin(r,y/(y/l));
        if(z/l) chkmin(r,z/(z/l));
        int mul=qpow(1ll*prd_phi[r]*prd_phi_inv[l-1]%mod,x/l);
        mul=1ll*mul*qpow(fac[x/l],sphi[r]-sphi[l-1])%mod;
        res=1ll*res*qpow(mul,1ll*(y/l)*(z/l)%(mod-1))%mod;
    } return res;
}
int calc4(int x,int y,int z){
    int res=1;
    for(int l=1,r;l<=min(x,y);l=r+1){
        r=1e9;
        if(x/l) chkmin(r,x/(x/l));
        if(y/l) chkmin(r,y/(y/l));
        if(z/l) chkmin(r,z/(z/l));
        res=1ll*res*qpow(1ll*prd_phi[r]*prd_phi_inv[l-1]%mod,1ll*(x/l)*(y/l)*(z/l)%(mod-1))%mod;
        int X=x/l,Y=y/l,Z=z/l,sm=1ll*Z*(sphi[r]-sphi[l-1])%(mod-1);
        for(int L=1,R;L<=min(X,Y);L=R+1){
            R=1e9;
            if(X/L) chkmin(R,X/(X/L));
            if(Y/L) chkmin(R,Y/(Y/L));
            res=1ll*res*qpow(1ll*prd[R]*prd_inv[L-1]%mod,1ll*(X/L)*(Y/L)*sm%(mod-1))%mod;
        }
    }
    return res;
}
int solve3(int x,int y,int z){
    return 1ll*calc3(x,y,z)*calc3(y,x,z)%mod*
    qpow(calc4(x,y,z),-1)%mod*qpow(calc4(x,z,y),-1)%mod;
}
int main(){
    int qu;scanf("%d%lld",&qu,&mod);
    _base1=(_base1<<64)/mod;_base2=(_base2<<64)/(mod-1);
    sieve(MAXV);
    while(qu--){
        int a,b,c;scanf("%d%d%d",&a,&b,&c);
        printf("%d %d %d\n",solve1(a,b,c),solve2(a,b,c),solve3(a,b,c));
    }
    return 0;
}