-
P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题
看完数据范围((T=70))大概是对于每个 (type) 做到 (O(n)) 以内。
type=0
[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}dfrac{operatorname{lcm(i,j)}}{gcd(i,k)}\ =prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}dfrac{i imes j}{gcd(i,j) imes gcd(i,k)}\ ]分母：
[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}gcd(i,j)\ =(prod_{d=1}^{A}d^{sum_{i=1}^{A}sum_{j=1}^{B}[gcd(i,j)==d]})^C\ =(prod_{d=1}^{A}d^{f(frac{n}{d},frac{m}{d})})^C ]其中
[f(n,m)=sum_{i=1}^{n} sum_{j=1}^{m}[gcd(i,j)==1]\ =sum_{x=1}^{n}mu(x)frac{n}{x}dfrac{m}{x} ](2) 层整除分块就可以 (O(n+sqrt{n}log n)) 了。
分子：
[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}ij\ =(prod_{i=1}^{A}prod_{j=1}^{B}ij)^C\ =((A!)^{B} imes(B!)^{A})^C ]预处理阶乘可 (O(log n)) 计算。
type=1
[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}(dfrac{i imes j}{gcd(i,j) imes gcd(i,k)})^{ijk}\ ]分子：
[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}(ij)^{ijk}\ =(prod_{i=1}^{A}(i^i)^{1+2+cdots A}prod_{j=1}^{B}(j^j)^{1+2+cdots+A})^{1+2+cdots+C}\ ]令 (jp_n=prod_{i=1}^{n}i^i) （快速幂预处理），则原式
[=((jp_A)^{1+2+cdots+B}(jp_B)^{1+2+cdots+A})^{1+2+cdots+C}\ =((jp_A)^{B imes(B+1)/2} imes (jp_B)^{A imes(A+1)/2})^{C imes(C+1)/2} ]分母：
[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}gcd(i,j)^{ijk}\ =(prod_{i=1}^{A}prod_{j=1}^{B}gcd(i,j)^{ij})^{C imes(C+1)/2}\ =(prod_{d=1}^{A}d^{f(A,B,d)})^{C imes(C+1)/2}\ ]其中
[f(n,m,d)=sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)==d]ij\ =d^2sum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{m}{d}}[gcd(i,j)==1]ij\ =d^2sum_{x=1}^{frac{n}{d}}mu(x)x^2sum_{i=1}^{frac{n}{dx}}sum_{j=1}^{frac{m}{dx}}ij\ ]带回原式换个写法
[ =(prod_{d=1}^{A}d^{f(A,B,d)})^{C imes(C+1)/2}\=(prod_{d=1}^{A}d^{g(frac{A}{d},frac{B}{d})*d^2})^{C imes(C+1)/2}\=(prod_{d=1}^{A}(d^{(d^2)})^{g(frac{A}{d},frac{B}{d})})^{C imes(C+1)/2}
]
其中
[ g(n,m)=sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)==1]ij
]
快速幂预处理 (d^{(d^2)}) 就可以整除分块套整除分块 (O(n)) 了
type=2
请确认您掌握狄利克雷卷积，这部分大量运用狄利克雷卷积来化简式子
[ prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}(dfrac{i imes j}{gcd(i,j) imes gcd(i,k)})^{gcd(i,j,k)}\
]
分子：
[ prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}i^{gcd(i,j,k)}\
=prod_{d=1}^{A}prod_{i=1}^{frac{A}{d}}(id)^{sum_{j=1}^{B}sum_{k=1}^{C}{[gcd(i,j,k)==d}]d}\
=prod_{d=1}^{A}prod_{i=1}^{frac{A}{d}}(id)^{sum_{j=1}^{frac{B}{d}}sum_{k=1}^{frac{C}{d}}{[gcd(i,j,k)==1}]d}\
=prod_{d=1}^{A}prod_{i=1}^{frac{A}{d}}(id)^{sum_{j=1}^{frac{B}{d}}sum_{k=1}^{frac{C}{d}}sum_{x|gcd(i,j,k)}mu(x)d}\
=prod_{d=1}^{A}prod_{x=1}^{frac{A}{d}}prod_{i=1}^{frac{A}{dx}}(idx)^{mu(x)dfrac{B}{dx}frac{C}{dx}}\
=prod_{T=1}^{A}prod_{d|T}^{}prod_{i=1}^{frac{A}{T}}(iT)^{mu(frac{T}{d})dfrac{B}{T}frac{C}{T}}\
]
注意到指数有一个 (sum_{d|T}mu(dfrac{T}{d})d=mu*I=varphi)
[ prod_{T=1}^{A}prod_{d|T}^{}prod_{i=1}^{frac{A}{T}}(iT)^{mu(frac{T}{d})dfrac{B}{T}frac{C}{T}}\
=prod_{T=1}^{A}prod_{i=1}^{frac{A}{T}}(iT)^{varphi(T)frac{B}{T}frac{C}{T}}\
=prod_{T=1}^{A}(jc_{frac{A}{T}}T^{frac{A}{T}})^{varphi(T)frac{B}{T}frac{C}{T}}\
]
预处理 (T^{varphi(T)}) 的前缀积和 (sumvarphi(i) \% (mod-1)) 整除分块即可。
分母：
[ prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}gcd(i,j)^{gcd(i,j,k)}\
=prod_{d=1}^{A}prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}[gcd(i,j)==d]d^{gcd(d,k)}\
=prod_{d=1}^{A}d^{sum_{i=1}^{A}sum_{j=1}^{B}[gcd(i,j)==d]sum_{k=1}^{C}gcd(d,k)}\
=prod_{d=1}^{A}d^{sum_{x=1}^{frac{A}{d}}mu(x)frac{A}{dx}frac{B}{dx}sum_{k=1}^{C}gcd(d,k)}\
=prod_{d=1}^{A}d^{sum_{T=1,d|T}^{A}mu(frac{T}{d})frac{A}{T}frac{B}{T}sum_{k=1}^{C}gcd(d,k)}\
]
单独提出后一部分来算。
[ sum_{k=1}^{C}gcd(d,k)\
=sum_{x|d}sum_{k=1}^{C}[gcd(d,k)==x]x\
=sum_{x|d}xsum_{k=1}^{frac{C}{x}}[gcd(dfrac{d}{x},k)==1]\
=sum_{x|d}xsum_{k=1}^{frac{C}{x}}sum_{t|gcd(frac{d}{x},k)}mu(t)\
=sum_{x|d}xsum_{t=1}^{frac{C}{x}}mu(t)dfrac{C}{tx}\
=sum_{x|d}xsum_{T=1,x|T}^{C}mu(dfrac{T}{x})dfrac{C}{T}\
=sum_{T=1,T|d}^{C}dfrac{C}{T}sum_{x|T}xmu(dfrac{T}{x})\
=sum_{T|d}varphi(T)dfrac{C}{T}
]
带回去
[ =prod_{d=1}^{A}d^{sum_{T=1,d|T}^{A}mu(frac{T}{d})frac{A}{T}frac{B}{T}sum_{k=1}^{C}gcd(d,k)}\
=prod_{d=1}^{A}d^{sum_{T=1,d|T}^{A}mu(frac{T}{d})frac{A}{T}frac{B}{T}sum_{x|d}varphi(x)frac{C}{x}}\
=prod_{T=1}^{A}(prod_{d=1,d|T}^{A}d^{mu(frac{T}{d})sum_{x|d}varphi(x)frac{C}{x}})^{frac{A}{T}frac{B}{T}}\
=prod_{T=1}^{A}(prod_{d|T}prod_{x|d}d^{mu(frac{T}{d})varphi(x)frac{C}{x}})^{frac{A}{T}frac{B}{T}}\
=prod_{T=1}^{A}(prod_{d|T}prod_{x|d}d^{mu(frac{T}{d})varphi(x)frac{C}{x}})^{frac{A}{T}frac{B}{T}}\
]
(x|d|T) 连着整除好奇怪。。。
考虑到 (d=xcdot dfrac{d}{x}) ，分成 (2) 部分算
部分一
[prod_{T=1}^{A}prod_{d|T}prod_{x|d}(dfrac{d}{x})^{mu(frac{T}{d})varphi(x)frac{C}{x}frac{A}{T}frac{B}{T}}\
=prod_{x=1}^{A}prod_{T=1}^{frac{A}{x}}prod_{x|d,d|Tx}(dfrac{d}{x})^{mu(frac{Tx}{d})varphi(x)frac{C}{x}frac{A}{Tx}frac{B}{Tx}}\
=prod_{x=1}^{A}prod_{T=1}^{frac{A}{x}}prod_{d|T}d^{mu(frac{T}{d})varphi(x)frac{C}{x}frac{A}{Tx}frac{B}{Tx}}\
=prod_{x=1}^{A}(prod_{T=1}^{frac{A}{x}}(prod_{d|T}d^{mu(frac{T}{d})})^{frac{A}{Tx}frac{B}{Tx}})^{varphi(x)frac{C}{x}}\
]
(O(nlog n)) 预处理 (d^{mu(frac{T}{d})}) 外边两层整除分块即可 (O(n))
小优化：(mu(x)) 非零的个数并不多，实测 (10000) 不到一点，虽然有个 (log) ，合起来近似 (O(n)) ，判一下 (mu) 的值可以快一点。
部分二
[ prod_{T=1}^{A}prod_{d|T}prod_{x|d}x^{mu(frac{T}{d})varphi(x)frac{C}{x}frac{A}{T}frac{B}{T}}\
=prod_{x=1}^{A}prod_{T=1}^{frac{A}{x}}prod_{x|d,d|Tx}x^{mu(frac{Tx}{d})varphi(x)frac{C}{x}frac{A}{Tx}frac{B}{Tx}}\
=prod_{x=1}^{A}prod_{T=1}^{frac{A}{x}}prod_{d|T}x^{mu(frac{T}{d})varphi(x)frac{C}{x}frac{A}{Tx}frac{B}{Tx}}\
=prod_{x=1}^{A}x^{sum_{T=1}^{frac{A}{x}}sum_{d|T}mu(frac{T}{d})varphi(x)frac{C}{x}frac{A}{Tx}frac{B}{Tx}}\
=prod_{x=1}^{A}x^{varphi(x)sum_{T=1}^{frac{A}{x}}frac{C}{x}frac{A}{Tx}frac{B}{Tx}sum_{d|T}mu(frac{T}{d})}\
]
(sum_{d|T}mu(dfrac{T}{d})=mu*I=epsilon)
所以只用计算 (T=1) 时的情形即可
[ =prod_{x=1}^{A}x^{varphi(x)sum_{T=1}^{frac{A}{x}}frac{C}{x}frac{A}{Tx}frac{B}{Tx}sum_{d|T}mu(frac{T}{d})}\
=sum_{x=1}^{A}x^{varphi(x)frac{C}{x}frac{A}{x}frac{B}{x}}
]
整除分块就够了！
注意事项
-
指数取模要模 (mod-1)
-
long long别少开，1ll别少乘 -
整除分块内层要快速幂的时候预处理逆元，否则单次询问变成 (O(nlog n))
心力憔悴啊，毒瘤。
提供一组大样例以供调试
Input
3 998244353
6 6 6
100000 100000 100000
99998 99999 100000
Output
570465346 682784945 524427235
862376103 371412326 358248208
103815203 127873959 745848036
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef double db;
#define pb(x) push_back(x)
#define mkp(x,y) make_pair(x,y)
//#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
//char buf[1<<21],*p1=buf,*p2=buf;
inline int read() {
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)) {if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch))x=x*10+(ch^48),ch=getchar();
return x*f;
}
const int N=100005;
int T,mod,A,B,C,iv6;
int mu[N],pri[N/9],cnt,jc[N],jp[N],sm[N],jk[N],st[N],phi[N],fp[N],vf[N];
bool vis[N];
int qpow(int n,int k,int res=1){
for(;k;k>>=1,n=1ll*n*n%mod)
if(k&1)res=1ll*n*res%mod;
return res;
}
void Sieve(const int N=100000){
mu[1]=1,phi[1]=1,jc[0]=1,jp[0]=1,jk[0]=1,st[0]=1,fp[0]=1,fp[1]=1,vf[0]=1;
for(int i=2;i<=N;++i){
fp[i]=1;
if(!vis[i])pri[++cnt]=i,mu[i]=-1,phi[i]=i-1;
for(int j=1;j<=cnt&&i*pri[j]<=N;++j){
vis[i*pri[j]]=1;
if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}
mu[i*pri[j]]=-mu[i],phi[i*pri[j]]=phi[i]*phi[pri[j]];
}
}
for(int i=1;i<=N;++i){
if(!mu[i])continue;
for(int j=1;i*j<=N;++j)
if(mu[i]==1)fp[i*j]=1ll*fp[i*j]*j%mod;
else fp[i*j]=1ll*fp[i*j]*qpow(j,mod-2)%mod;
}
for(int i=1;i<=N;++i)fp[i]=1ll*fp[i]*fp[i-1]%mod,vf[i]=qpow(fp[i],mod-2);
for(int i=1;i<=N;++i)
jc[i]=1ll*jc[i-1]*i%mod,
jp[i]=1ll*jp[i-1]*qpow(i,i)%mod,
sm[i]=(sm[i-1]+1ll*mu[i]*i%(mod-1)*i%(mod-1))%(mod-1),
jk[i]=1ll*jk[i-1]*qpow(i,1ll*i*i%(mod-1))%mod,
mu[i]+=mu[i-1],
st[i]=1ll*st[i-1]*qpow(i,phi[i])%mod,
phi[i]=(phi[i]+phi[i-1])%(mod-1);
/*
Attention:
mu:∑μ
jc :阶乘 mod mod
jp : i^i 的前缀积 mod mod
sm：∑μ(x)*x*x mod (mod-1)
jk:∏i^(i^2) mod (mod-1)
phi:∑φ mod (mod-1)
st:πT^phi(T) mod mod
fp:∑π(d|T)d^μ(T/d) mod mod
vf:qpow(fp,mod-2)
*/
}
namespace Task0{
int fz,fm;
int f(int n,int m){
if(n>m)n^=m^=n^=m;
int res=0;
for(int l=1,r;l<=n;l=r+1)
r=min(n/(n/l),m/(m/l)),res=(res+1ll*(n/l)*(m/l)%(mod-1)*(mu[r]-mu[l-1])%(mod-1))%(mod-1);
return (res+mod-1)%(mod-1);
}
int g(int n,int m){
if(n>m)n^=m^=n^=m;
int res=1;
for(int l=1,r;l<=n;l=r+1)
r=min(n/(n/l),m/(m/l)),res=1ll*res*qpow(1ll*jc[r]*qpow(jc[l-1],mod-2)%mod,f(n/l,m/l))%mod;
return (res+mod)%mod;
}
void main(){
fz=qpow(1ll*qpow(jc[A],B)*qpow(jc[B],A)%mod,C)%mod;
fm=qpow(1ll*qpow(g(A,B),C)*qpow(g(A,C),B)%mod,mod-2);
printf("%lld ",1ll*fz*fm%mod);
}
}
namespace Task1{
int fz,fm;
int s(int x,int y){
return (1ll*x*(x+1)/2%(mod-1))*(1ll*y*(y+1)/2%(mod-1))%(mod-1);
}
int s2(int x){
return 1ll*x*(x+1)%mod*(x+x+1)%mod*iv6%mod;
}
int f(int n,int m){
if(n>m)n^=m^=n^=m;
int res=0;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
res=(res+1ll*s(n/l,m/l)*(sm[r]-sm[l-1])%(mod-1))%(mod-1);
}
return (res+mod-1)%(mod-1);
}
int g(int n,int m){
if(n>m)n^=m^=n^=m;
int res=1;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
res=1ll*res*qpow(1ll*jk[r]*qpow(jk[l-1],mod-2)%mod,f(n/l,m/l))%mod;
}
return (res+mod)%mod;
}
void main(){
fz=qpow(1ll*qpow(jp[A],1ll*B*(B+1)/2%(mod-1))*qpow(jp[B],1ll*A*(A+1)/2%(mod-1))%mod,1ll*C*(C+1)/2%(mod-1));
fm=qpow(1ll*qpow(g(A,B),1ll*C*(C+1)/2%(mod-1))*qpow(g(A,C),1ll*B*(B+1)/2%(mod-1))%mod,mod-2);
printf("%lld ",1ll*fz*fm%mod);
}
}
namespace Task2{
int fz,fm;
int f(int A,int B,int C){
int res=1;
for(int l=1,r,mx=min(A,min(B,C));l<=mx;l=r+1){
r=min(A/(A/l),min(B/(B/l),C/(C/l)));
res=1ll*res*qpow(jc[A/l],1ll*(B/l)*(C/l)%(mod-1)*(phi[r]-phi[l-1]+mod-1)%(mod-1))%mod;
res=1ll*res*qpow(1ll*st[r]*qpow(st[l-1],mod-2)%mod,1ll*(A/l)*(B/l)*(C/l)%(mod-1))%mod;
}
return res;
}
int h(int n,int m){
if(n>m)n^=m^=n^=m;
int res=1;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
res=1ll*res*qpow(1ll*fp[r]*vf[l-1]%mod,1ll*(n/l)*(m/l)%(mod-1))%mod;
}
return res;
}
int g(int A,int B,int C){
int res=1;
for(int l=1,r,mx=min(A,min(B,C));l<=mx;l=r+1){
r=min(A/(A/l),min(B/(B/l),C/(C/l)));
res=1ll*res*qpow(1ll*st[r]*qpow(st[l-1],mod-2)%mod,1ll*(A/l)*(B/l)*(C/l)%(mod-1))%mod;
res=1ll*res*qpow(h(A/l,B/l),1ll*(phi[r]-phi[l-1]+mod-1)*(C/l)%(mod-1))%mod;
}
return res;
}
void main(){
fz=1ll*f(A,B,C)*f(B,A,C)%mod;
fm=qpow(1ll*g(A,B,C)*g(A,C,B)%mod,mod-2);
printf("%lld ",1ll*fz*fm%mod);
}
}
signed main(){
T=read(),mod=read(),iv6=qpow(6,mod-2),Sieve();
while(T--) {
A=read(),B=read(),C=read();
Task0::main(),Task1::main(),Task2::main(),puts("");
}
}