題解-MtOI2019 幽靈樂團
題面
給定 \(p\)，\(Cnt\) 組測試資料，每次給 \(a,b,c\)，求
\[\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^c\left(\frac{{\rm lcm}(i,j)}{\gcd(i,k)}\right)^{f(t)}\bmod p \]
\(t\in\{0,1,2\}\)，\(f(0)=1,f(1)=ijk,f(2)=\gcd(i,j,k)\)。
資料範圍：\(1\le a,b,c\le 10^5\)，\(10^7\le p\le 105\cdot 10^7\)，\(p\in\mathbb{P}\)，\(1\le Cnt\le 70\)
。
蒟蒻語
初學莫比烏斯反演是很久之前了，最近集訓老師講了莫比烏斯反演。
那天蒟蒻沒去，神 \(\tt zhoukangyang\) 說：
老師抓了隔壁機房的 \(\texttt{x義x}\)（\(\tt orz\)）來現場做「MtOI2019 幽靈樂團」，結果他上來就開了題解……
於是神 \(\tt zky\) 和神【資料刪除】就幹了這題一個晚上。第二天蒟蒻去推了一個下午，終於推出來了，推了一面，細節很多，但不難。最後蒟蒻比他們兩個先做出來 \(\tt /yx\)。
蒟蒻解
\[\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^c\left(\frac{{\rm lcm}(i,j)}{\gcd(i,k)}\right)^{f(t)} =\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^c\left(\frac{ij}{\gcd(i,j)\gcd(i,k)}\right)^{f(t)} \]
所以只需算出:
\[\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^ci^{f(t)} \]
和
\[\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^c\gcd(i,j)^{f(t)} \]
即可，共 \(2\times 3=6\) 種式子。
t=0
\[\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^ci=(a!)^{bc} \]
預處理 \(n!\) 時間複雜度 \(\Theta(n)\)，計算 \(\Theta(1)\)。
\[\begin{split} &\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^c\gcd(i,j)\\ =&\left(\prod_{d=1}^{\min(a,b)}d^{\sum\limits_{i=1}^a\sum\limits_{j=1}^b[\gcd(i,j)=d]}\right)^c\\ =&\left(\prod_{d=1}^{\min(a,b)}d^{\sum\limits_{k=1}^{\lfloor\frac{a}{d}\rfloor}\mu(k)\lfloor\frac{a}{dk}\rfloor\lfloor\frac{b}{dk}\rfloor}\right)^c\\ =&\left(\prod_{T=1}^{\min(a,b)}\left(\prod_{d|T}d^{\mu(\frac{T}{d})}\right)^{\lfloor\frac{a}{T}\rfloor\lfloor\frac{b}{T}\rfloor}\right)^c \end{split} \]
預處理 \(\prod_{d|T}d^{\mu(\frac{T}{d})}\) 字首積時間複雜度 \(\Theta(n\log\log n)\)，分塊計算 \(\Theta(\sqrt{n})\)。
t=1
令 \(S(n)=\frac{n(n+1)}{2}\)。
\[\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^ci^{ijk}=\left(\prod_{i=1}^a i^i\right)^{S(b)S(c)} \]
預處理 \(\prod_{i=1}^a i^i\) 時間複雜度 \(\Theta(n)\)，計算 \(\Theta(1)\)。
未完待續
程式碼
#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 PredGcd(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<<"PredGcd("<<a<<','<<b<<"):"<<res<<'\n';
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<<'\n';
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<<'\n';
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<<'\n';
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<<'\n';
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<<'\n';
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(PredGcd(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<<'\n';
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]<<'\n';
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<<'\n';
}
// clock_t end=clock();
// cout<<db(end-start)/CLOCKS_PER_SEC<<'\n';
return 0;
}