洛谷 P5518 - [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题（莫比乌斯反演+整除分块）
一道究极恶心的毒瘤六合一题，式子推了我满满两面 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;
}
