题面
给定 (p)，(Cnt) 组测试数据，每次给 (a,b,c)，求
[prod_{i=1}^aprod_{j=1}^bprod_{k=1}^cleft(frac{{ m lcm}(i,j)}{gcd(i,k)} ight)^{f(t)}mod p ](tin{0,1,2})，(f(0)=1,f(1)=ijk,f(2)=gcd(i,j,k))。
数据范围：(1le a,b,cle 10^5)，(10^7le ple 105cdot 10^7)，(pinmathbb{P})，(1le Cntle 70)。
蒟蒻语
初学莫比乌斯反演是很久之前了，最近集训老师讲了莫比乌斯反演。
那天蒟蒻没去，神 ( t zhoukangyang) 说：
老师抓了隔壁机房的 ( exttt{x义x})（( t orz)）来现场做「MtOI2019 幽灵乐团」，结果他上来就开了题解……
于是神 ( t zky) 和神【数据删除】就干了这题一个晚上。第二天蒟蒻去推了一个下午，终于推出来了，推了一面，细节很多，但不难。最后蒟蒻比他们两个先做出来 ( t /yx)。
蒟蒻解
[prod_{i=1}^aprod_{j=1}^bprod_{k=1}^cleft(frac{{
m lcm}(i,j)}{gcd(i,k)}
ight)^{f(t)}
=prod_{i=1}^aprod_{j=1}^bprod_{k=1}^cleft(frac{ij}{gcd(i,j)gcd(i,k)}
ight)^{f(t)}
]
所以只需算出:
[prod_{i=1}^aprod_{j=1}^bprod_{k=1}^ci^{f(t)}
]
和
[prod_{i=1}^aprod_{j=1}^bprod_{k=1}^cgcd(i,j)^{f(t)}
]
即可，共 (2 imes 3=6) 种式子。
t=0
[prod_{i=1}^aprod_{j=1}^bprod_{k=1}^ci=(a!)^{bc}
]
预处理 (n!) 时间复杂度 (Theta(n))，计算 (Theta(log n))。
[egin{split}
&prod_{i=1}^aprod_{j=1}^bprod_{k=1}^cgcd(i,j)\
=&left(prod_{d=1}^{min(a,b)}d^{sumlimits_{i=1}^asumlimits_{j=1}^b[gcd(i,j)=d]}
ight)^c\
=&left(prod_{d=1}^{min(a,b)}d^{sumlimits_{k=1}^{lfloorfrac{a}{d}
floor}mu(k)lfloorfrac{a}{dk}
floorlfloorfrac{b}{dk}
floor}
ight)^c\
=&left(prod_{T=1}^{min(a,b)}left(prod_{d|T}d^{mu(frac{T}{d})}
ight)^{lfloorfrac{a}{T}
floorlfloorfrac{b}{T}
floor}
ight)^c
end{split}
]
预处理 (prod_{d|T}d^{mu(frac{T}{d})}) 前缀积时间复杂度 (Theta(nloglog n))，分块计算 (Theta(sqrt{n}log n))。
t=1
令 (S(n)=frac{n(n+1)}{2})。
[prod_{i=1}^aprod_{j=1}^bprod_{k=1}^ci^{ijk}=left(prod_{i=1}^a i^i
ight)^{S(b)S(c)}
]
预处理 (prod_{i=1}^a i^i) 时间复杂度 (Theta(n))，计算 (Theta(log n))。
[egin{split}
&prod_{i=1}^aprod_{j=1}^bprod_{k=1}^c gcd(i,j)^{ij}\
=&left(prod_{d=1}^{min(a,b)}d^{sumlimits_{i=1}^asumlimits_{j=1}^b[gcd(i,j)=d]ij}
ight)^{S(c)}\
=&left(prod_{d=1}^{min(a,b)}d^{d^2sumlimits_{k=1}^{lfloorfrac{min(a,b)}{d}
floor}mu(k)k^2S(lfloorfrac{a}{dk}
floor)S(lfloorfrac{b}{dk}
floor)}
ight)^{S(c)}\
=&left(prod_{T=1}^{min(a,b)}left(prod_{d|T}d^{mu(frac{T}{d})}
ight)^{T^2S(lfloorfrac{a}{T})
floor S(lfloorfrac{b}{T}
floor)}
ight)^{S(c)}\
end{split}
]
在 (prod_{d|T}d^{mu(frac{T}{d})}) 的基础上预处理 (left(prod_{d|T}d^{mu(frac{T}{d})} ight)^{T^2}) 时间复杂度 (Theta(nlog n))，分块计算 (Theta(sqrt{n}log n))。
t=2
[egin{split}
&prod_{i=1}^aprod_{j=1}^bprod_{k=1}^ci^{gcd(i,j,k)}\
=&prod_{d=1}^{min(a,b,c)}prod_{i=1}^a i^{dsumlimits_{j=1}^bsumlimits_{k=1}^c[(i,j,k)=d]}\
=&prod_{d=1}^{min(a,b,c)}prod_{h=1}^{lfloorfrac{min(a,b,c)}{d}
floor}left(prod_{i=1}^{lfloorfrac{a}{dh}
floor}(idh)^{dmu(h)}
ight)^{lfloorfrac{a}{dh}
floorlfloorfrac{b}{dh}
floor}\
=&prod_{T=1}^{min(a,b,c)}left(lfloorfrac{a}{T}
floor!T^{lfloorfrac{a}{T}
floor}
ight)^{varphi(T)lfloorfrac{b}{T}
floorlfloorfrac{c}{T}
floor}\
=&prod_{T=1}^{min(a,b,c)}left(lfloorfrac{a}{T}
floor!
ight)^{varphi(T)lfloorfrac{b}{T}
floorlfloorfrac{c}{T}
floor} imesprod_{T=1}^{min(a,b,c)}left(T^{varphi(T)}
ight)^{lfloorfrac{a}{T}
floorlfloorfrac{b}{T}
floorlfloorfrac{c}{T}
floor}
end{split}
]
左边用已经预处理过的阶乘，右边预处理 (T^{varphi(T)}) 前缀积时间复杂度 (Theta(nlog n))，分块计算 (Theta(sqrt{n}log n))。
[egin{split}
&prod_{i=1}^aprod_{j=1}^bprod_{k=1}^cgcd(i,j)^{gcd(i,j,k)}\
=&prod_{d=1}^{min(a,b,c)}prod_{i=1}^aprod_{j=1}^bgcd(i,j)^{dsumlimits_{k=1}^c[gcd(i,j,k)=d]}\
=&prod_{d=1}^{min(a,b,c)}prod_{h=1}^{lfloorfrac{min(a,b,c)}{d}
floor}prod_{i=1}^{lfloorfrac{a}{dh}
floor}prod_{j=1}^{lfloorfrac{b}{dh}
floor}left(gcd(i,j)dh
ight)^{dmu(h)lfloorfrac{c}{dh}
floor}\
=&prod_{T=1}^{min(a,b,c)}T^{sum_{d|T}dmu(frac{T}{d})lfloorfrac{a}{T}
floorlfloorfrac{b}{T}
floorlfloorfrac{c}{T}
floor} imes prod_{T=1}^{min(a,b,c)}left(prod_{i=1}^{lfloorfrac{a}{T}
floor}prod_{j=1}^{lfloorfrac{b}{T}
floor}gcd(i,j)^{sum_{d|T}dmu(frac{T}{d})}
ight)^{lfloorfrac{c}{T}
floor}\
=&prod_{T=1}^{min(a,b,c)}T^{varphi(T)lfloorfrac{a}{T}
floorlfloorfrac{b}{T}
floorlfloorfrac{c}{T}
floor} imes prod_{T=1}^{min(a,b,c)}left(prod_{i=1}^{lfloorfrac{a}{T}
floor}prod_{j=1}^{lfloorfrac{b}{T}
floor}gcd(i,j)^{varphi(T)}
ight)^{lfloorfrac{c}{T}
floor}\
=&prod_{T=1}^{min(a,b,c)}T^{varphi(T)lfloorfrac{a}{T}
floorlfloorfrac{b}{T}
floorlfloorfrac{c}{T}
floor} imes prod_{T=1}^{min(a,b,c)}left(prod_{G=1}^{lfloorfrac{min(a,b,c)}{T}
floor}left(prod_{e|G}e^{mu(frac{G}{e})}
ight)^{lfloorfrac{a}{TG}
floorlfloorfrac{b}{TG}
floor}
ight)^{lfloorfrac{c}{T}
floorvarphi(T)}\
end{split}
]
有点跳步，但是有了这点拨并不难推，毕竟前人之述备矣，然则有些细节，确实该细讲。
(T^{varphi(T)}) 和 (prod_{e|G}e^{mu(frac{G}{e})}) 已经预处理过了，分块套分块 (Theta(nlog n))。
细节&优化
注意要模 (p)，但是根据欧拉定理，指数要模 (varphi(p)=p-1)。
看到一个式子怎么知道该筛什么？看看哪些东西不可以分块统一算。
怎么 (Theta(nlog log n)) 预处理 (prod_{d|T}d^{mu(frac{T}{d})})？
利用积性函数的性质，同时计算这东西和这东西的逆元（注意还没有求前缀积）：
for(int i=1;i<=N;i++) f[i]=i,ivf[i]=in[i];
for(int j=0;j<cp;j++)for(int i=N/p[j];i>=1;i--)
f[i*p[j]]=(ll)f[i*p[j]]*ivf[i]%P,ivf[i*p[j]]=(ll)ivf[i*p[j]]*f[i]%P;
现在的时间复杂度是 (Theta(Tnlog n)) 我卡过去了，但是是可以优化到 (Theta(n^{frac{4}{3}}+Tn^{frac{2}{3}}log n))。
注意这个式子里面有什么：
[prod_{T=1}^{min(a,b,c)}left(prod_{G=1}^{lfloorfrac{min(a,b,c)}{T}
floor}left(prod_{e|G}e^{mu(frac{G}{e})}
ight)^{lfloorfrac{a}{TG}
floorlfloorfrac{b}{TG}
floor}
ight)^{lfloorfrac{c}{T}
floorvarphi(T)}
]
[egin{split}
g(a,b)=&prod_{T=1}^{min(a,b)}left(prod_{d|T}d^{mu(frac{T}{d})}
ight)^{lfloorfrac{a}{T}
floorlfloorfrac{b}{T}
floor}\
=&prod_{i=1}^aprod_{j=1}^b gcd(i,j)\
end{split}
]
所以可以 (Theta(n^{frac{4}{3}})) 用类似辗转相除的方法预处理 (a,ble n^{frac{2}{3}}) 的 (g(a,b)) 以及前缀积，然后总时间复杂度就 (Theta(n^{frac{4}{3}}+Tn^{frac{2}{3}}log n)) 了。
for(int i=1;i<=K;i++){
for(int j=1;j<i;j++){
g[i][j]=g[max(j,i-j)][min(j,i-j)];
sg[i][j]=(ll)sg[i-1][j]*sg[i][j-1]%P*isg[i-1][j-1]%P*g[i][j]%P;
isg[i][j]=(ll)isg[i-1][j]*isg[i][j-1]%P*sg[i-1][j-1]%P*in[g[i][j]]%P;
}
g[i][i]=i;
sg[i][i]=(ll)sg[i][i-1]*sg[i][i-1]%P*isg[i-1][i-1]%P*i%P;
isg[i][i]=(ll)isg[i][i-1]*isg[i][i-1]%P*sg[i-1][i-1]%P*in[i]%P;
}
代码
#include <bits/stdc++.h>
using namespace std;
//Start
typedef long long ll;
typedef double db;
#define mp(a,b) make_pair(a,b)
#define x first
#define y second
#define be(a) a.begin()
#define en(a) a.end()
#define sz(a) int((a).size())
#define pb(a) push_back(a)
const int inf=0x3f3f3f3f;
const ll INF=0x3f3f3f3f3f3f3f3f;
//Data
const int N=1e5;
int P,fac[N+1],faci[N+1],in[N+1];
int sum(int a){return (ll)a*(a+1)/2%(P-1);}
int Pow(int a,int x){
if(!a) return 0; int res=1;
for(;x;a=(ll)a*a%P,x>>=1)if(x&1) res=(ll)res*a%P;
return res;
}
//Sieve
const int K=2500;
bool np[N+1];
int mu[N+1],phi[N+1],cp,p[N];
int f[N+1],ivf[N+1],fii[N+1],ifii[N+1],dp[N+1],idp[N+1];
int g[K+1][K+1],sg[K+1][K+1],isg[K+1][K+1];
void Sieve(){
np[1]=true,mu[1]=phi[1]=fii[0]=ifii[0]=dp[0]=idp[0]=f[0]=ivf[0]=in[1]=fac[0]=faci[0]=1;
for(int i=2;i<=N;i++) in[i]=(ll)(P-P/i)*in[P%i]%P;
for(int i=1;i<=N;i++) fac[i]=(ll)fac[i-1]*i%P,faci[i]=(ll)faci[i-1]*Pow(i,i)%P;
for(int i=2;i<=N;i++){
if(!np[i]) mu[i]=-1,phi[i]=i-1,p[cp++]=i;
for(int j=0;j<cp&&i*p[j]<=N;j++){
np[i*p[j]]=1;
if(i%p[j]==0){mu[i*p[j]]=0,phi[i*p[j]]=phi[i]*p[j];break;}
mu[i*p[j]]=mu[i]*mu[p[j]],phi[i*p[j]]=phi[i]*phi[p[j]];
}
}
for(int i=1;i<=N;i++) f[i]=i,ivf[i]=in[i];
for(int j=0;j<cp;j++)for(int i=N/p[j];i>=1;i--)
f[i*p[j]]=(ll)f[i*p[j]]*ivf[i]%P,ivf[i*p[j]]=(ll)ivf[i*p[j]]*f[i]%P;
for(int i=1;i<=N;i++){
fii[i]=Pow(f[i],(ll)i*i%(P-1)),ifii[i]=Pow(ivf[i],(ll)i*i%(P-1));
dp[i]=Pow(i,phi[i]%(P-1)),idp[i]=Pow(in[i],phi[i]%(P-1));
}
for(int i=2;i<=N;i++){
f[i]=(ll)f[i]*f[i-1]%P,ivf[i]=(ll)ivf[i]*ivf[i-1]%P;
fii[i]=(ll)fii[i]*fii[i-1]%P,ifii[i]=(ll)ifii[i]*ifii[i-1]%P;
dp[i]=(ll)dp[i]*dp[i-1]%P,idp[i]=(ll)idp[i]*idp[i-1]%P;
(phi[i]+=phi[i-1])%=(P-1);
}
for(int i=0;i<=K;i++) g[i][0]=g[0][i]=i,sg[i][0]=sg[0][i]=isg[i][0]=isg[0][i]=1;
for(int i=1;i<=K;i++){
for(int j=1;j<i;j++){
g[i][j]=g[max(j,i-j)][min(j,i-j)];
sg[i][j]=(ll)sg[i-1][j]*sg[i][j-1]%P*isg[i-1][j-1]%P*g[i][j]%P;
isg[i][j]=(ll)isg[i-1][j]*isg[i][j-1]%P*sg[i-1][j-1]%P*in[g[i][j]]%P;
}
g[i][i]=i;
sg[i][i]=(ll)sg[i][i-1]*sg[i][i-1]%P*isg[i-1][i-1]%P*i%P;
isg[i][i]=(ll)isg[i][i-1]*isg[i][i-1]%P*sg[i-1][i-1]%P*in[i]%P;
}
}
//Mobius
int ProdGcd(int a,int b){
int res=1;
if(a<=K&&b<=K) res=sg[max(a,b)][min(a,b)];
else {
for(int d=1,r;d<=min(a,b);d=r+1){
r=min(a/(a/d),b/(b/d));
res=(ll)res*Pow((ll)f[r]*ivf[d-1]%P,(ll)(a/d)*(b/d)%(P-1))%P;
}
}
// cout<<"ProdGcd("<<a<<','<<b<<"):"<<res<<'
';
return res;
}
int Geti0(int a,int b,int c){
int res=Pow(fac[a],(ll)b*c%(P-1));
// cout<<"Geti0("<<a<<','<<b<<','<<c<<"):"<<res<<'
';
return res;
}
int Getg0(int a,int b,int c){
int res=1;
for(int d=1,r;d<=min(a,b);d=r+1){
r=min(a/(a/d),b/(b/d));
res=(ll)res*Pow((ll)f[r]*ivf[d-1]%P,(ll)(a/d)*(b/d)%(P-1))%P;
}
res=Pow(res,c%(P-1));
// cout<<"Getg0("<<a<<','<<b<<','<<c<<"):"<<res<<'
';
return res;
}
int Geti1(int a,int b,int c){
int res=Pow(faci[a],(ll)sum(b)*sum(c)%(P-1));
// cout<<"Geti1("<<a<<','<<b<<','<<c<<"):"<<res<<'
';
return res;
}
int Getg1(int a,int b,int c){
int res=1;
for(int d=1,r;d<=min(a,b);d=r+1){
r=min(a/(a/d),b/(b/d));
res=(ll)res*Pow((ll)fii[r]*ifii[d-1]%P,(ll)sum(a/d)*sum(b/d)%(P-1))%P;
}
res=Pow(res,sum(c)%(P-1));
// cout<<"Getg1("<<a<<','<<b<<','<<c<<"):"<<res<<'
';
return res;
}
int Geti2(int a,int b,int c){
int res=1;
for(int d=1,r;d<=min(a,min(b,c));d=r+1){
r=min(a/(a/d),min(b/(b/d),c/(c/d)));
res=(ll)res*Pow(fac[a/d],(ll)(b/d)*(c/d)%(P-1)*(phi[r]-phi[d-1]+P-1)%(P-1))%P;
res=(ll)res*Pow((ll)dp[r]*idp[d-1]%P,(ll)(a/d)*(b/d)%(P-1)*(c/d)%(P-1))%P;
}
// cout<<"Geti2("<<a<<','<<b<<','<<c<<"):"<<res<<'
';
return res;
}
int Getg2(int a,int b,int c){
int res=1;
for(int d=1,r;d<=min(a,min(b,c));d=r+1){
r=min(a/(a/d),min(b/(b/d),c/(c/d)));
res=(ll)res*Pow(ProdGcd(a/d,b/d),(ll)(c/d)*(phi[r]-phi[d-1]+P-1)%(P-1))%P;
res=(ll)res*Pow((ll)dp[r]*idp[d-1]%P,(ll)(a/d)*(b/d)%(P-1)*(c/d)%(P-1))%P;
}
// cout<<"Getg2("<<a<<','<<b<<','<<c<<"):"<<res<<'
';
return res;
}
//Main
int main(){
ios::sync_with_stdio(0);
cin.tie(0),cout.tie(0);
int cacnt; cin>>cacnt>>P;
// clock_t start=clock();
Sieve();// cout<<f[2]<<'
';
while(cacnt--){
int a,b,c,ans0,ans1,ans2;
cin>>a>>b>>c;
ans0=(ll)Geti0(a,b,c)*Geti0(b,a,c)%P*Pow(Getg0(a,b,c),P-2)%P*Pow(Getg0(a,c,b),P-2)%P;
ans1=(ll)Geti1(a,b,c)*Geti1(b,a,c)%P*Pow(Getg1(a,b,c),P-2)%P*Pow(Getg1(a,c,b),P-2)%P;
ans2=(ll)Geti2(a,b,c)*Geti2(b,a,c)%P*Pow(Getg2(a,b,c),P-2)%P*Pow(Getg2(a,c,b),P-2)%P;
cout<<ans0<<' '<<ans1<<' '<<ans2<<'
';
}
// clock_t end=clock();
// cout<<db(end-start)/CLOCKS_PER_SEC<<'
';
return 0;
}
祝大家学习愉快！