题解-MtOI2019 幽灵乐团

题面

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(\log n)\)

\[\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}\log 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(\log n)\)

\[\begin{split}
&\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^c \gcd(i,j)^{ij}\\
=&\left(\prod_{d=1}^{\min(a,b)}d^{\sum\limits_{i=1}^a\sum\limits_{j=1}^b[\gcd(i,j)=d]ij}\right)^{S(c)}\\
=&\left(\prod_{d=1}^{\min(a,b)}d^{d^2\sum\limits_{k=1}^{\lfloor\frac{\min(a,b)}{d}\rfloor}\mu(k)k^2S(\lfloor\frac{a}{dk}\rfloor)S(\lfloor\frac{b}{dk}\rfloor)}\right)^{S(c)}\\
=&\left(\prod_{T=1}^{\min(a,b)}\left(\prod_{d|T}d^{\mu(\frac{T}{d})}\right)^{T^2S(\lfloor\frac{a}{T})\rfloor S(\lfloor\frac{b}{T}\rfloor)}\right)^{S(c)}\\
\end{split}
\]

\(\prod_{d|T}d^{\mu(\frac{T}{d})}\) 的基础上预处理 \(\left(\prod_{d|T}d^{\mu(\frac{T}{d})}\right)^{T^2}\) 时间复杂度 \(\Theta(n\log n)\),分块计算 \(\Theta(\sqrt{n}\log n)\)


t=2

\[\begin{split}
&\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^ci^{\gcd(i,j,k)}\\
=&\prod_{d=1}^{\min(a,b,c)}\prod_{i=1}^a i^{d\sum\limits_{j=1}^b\sum\limits_{k=1}^c[(i,j,k)=d]}\\
=&\prod_{d=1}^{\min(a,b,c)}\prod_{h=1}^{\lfloor\frac{\min(a,b,c)}{d}\rfloor}\left(\prod_{i=1}^{\lfloor\frac{a}{dh}\rfloor}(idh)^{d\mu(h)}\right)^{\lfloor\frac{a}{dh}\rfloor\lfloor\frac{b}{dh}\rfloor}\\
=&\prod_{T=1}^{\min(a,b,c)}\left(\lfloor\frac{a}{T}\rfloor!T^{\lfloor\frac{a}{T}\rfloor}\right)^{\varphi(T)\lfloor\frac{b}{T}\rfloor\lfloor\frac{c}{T}\rfloor}\\
=&\prod_{T=1}^{\min(a,b,c)}\left(\lfloor\frac{a}{T}\rfloor!\right)^{\varphi(T)\lfloor\frac{b}{T}\rfloor\lfloor\frac{c}{T}\rfloor}\times\prod_{T=1}^{\min(a,b,c)}\left(T^{\varphi(T)}\right)^{\lfloor\frac{a}{T}\rfloor\lfloor\frac{b}{T}\rfloor\lfloor\frac{c}{T}\rfloor}
\end{split}
\]

左边用已经预处理过的阶乘,右边预处理 \(T^{\varphi(T)}\) 前缀积时间复杂度 \(\Theta(n\log n)\),分块计算 \(\Theta(\sqrt{n}\log n)\)

\[\begin{split}
&\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^c\gcd(i,j)^{\gcd(i,j,k)}\\
=&\prod_{d=1}^{\min(a,b,c)}\prod_{i=1}^a\prod_{j=1}^b\gcd(i,j)^{d\sum\limits_{k=1}^c[\gcd(i,j,k)=d]}\\
=&\prod_{d=1}^{\min(a,b,c)}\prod_{h=1}^{\lfloor\frac{\min(a,b,c)}{d}\rfloor}\prod_{i=1}^{\lfloor\frac{a}{dh}\rfloor}\prod_{j=1}^{\lfloor\frac{b}{dh}\rfloor}\left(\gcd(i,j)dh\right)^{d\mu(h)\lfloor\frac{c}{dh}\rfloor}\\
=&\prod_{T=1}^{\min(a,b,c)}T^{\sum_{d|T}d\mu(\frac{T}{d})\lfloor\frac{a}{T}\rfloor\lfloor\frac{b}{T}\rfloor\lfloor\frac{c}{T}\rfloor}\times \prod_{T=1}^{\min(a,b,c)}\left(\prod_{i=1}^{\lfloor\frac{a}{T}\rfloor}\prod_{j=1}^{\lfloor\frac{b}{T}\rfloor}\gcd(i,j)^{\sum_{d|T}d\mu(\frac{T}{d})}\right)^{\lfloor\frac{c}{T}\rfloor}\\
=&\prod_{T=1}^{\min(a,b,c)}T^{\varphi(T)\lfloor\frac{a}{T}\rfloor\lfloor\frac{b}{T}\rfloor\lfloor\frac{c}{T}\rfloor}\times \prod_{T=1}^{\min(a,b,c)}\left(\prod_{i=1}^{\lfloor\frac{a}{T}\rfloor}\prod_{j=1}^{\lfloor\frac{b}{T}\rfloor}\gcd(i,j)^{\varphi(T)}\right)^{\lfloor\frac{c}{T}\rfloor}\\
=&\prod_{T=1}^{\min(a,b,c)}T^{\varphi(T)\lfloor\frac{a}{T}\rfloor\lfloor\frac{b}{T}\rfloor\lfloor\frac{c}{T}\rfloor}\times \prod_{T=1}^{\min(a,b,c)}\left(\prod_{G=1}^{\lfloor\frac{\min(a,b,c)}{T}\rfloor}\left(\prod_{e|G}e^{\mu(\frac{G}{e})}\right)^{\lfloor\frac{a}{TG}\rfloor\lfloor\frac{b}{TG}\rfloor}\right)^{\lfloor\frac{c}{T}\rfloor\varphi(T)}\\
\end{split}
\]

有点跳步,但是有了这点拨并不难推,毕竟前人之述备矣,然则有些细节,确实该细讲。

\(T^{\varphi(T)}\)\(\prod_{e|G}e^{\mu(\frac{G}{e})}\) 已经预处理过了,分块套分块 \(\Theta(n\log n)\)


细节&优化

注意要模 \(p\),但是根据欧拉定理,指数要模 \(\varphi(p)=p-1\)

看到一个式子怎么知道该筛什么?看看哪些东西不可以分块统一算。

怎么 \(\Theta(n\log \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(Tn\log n)\) 我卡过去了,但是是可以优化到 \(\Theta(n^{\frac{4}{3}}+Tn^{\frac{2}{3}}\log n)\)

注意这个式子里面有什么:

\[\prod_{T=1}^{\min(a,b,c)}\left(\prod_{G=1}^{\lfloor\frac{\min(a,b,c)}{T}\rfloor}\left(\prod_{e|G}e^{\mu(\frac{G}{e})}\right)^{\lfloor\frac{a}{TG}\rfloor\lfloor\frac{b}{TG}\rfloor}\right)^{\lfloor\frac{c}{T}\rfloor\varphi(T)}
\]

\[\begin{split}
g(a,b)=&\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}\\
=&\prod_{i=1}^a\prod_{j=1}^b \gcd(i,j)\\
\end{split}
\]

所以可以 \(\Theta(n^{\frac{4}{3}})\) 用类似辗转相除的方法预处理 \(a,b\le 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<<'\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(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<<'\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;
}


祝大家学习愉快!

原文链接: https://www.cnblogs.com/George1123/p/13341419.html

欢迎关注

微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍;

也有高质量的技术群,里面有嵌入式、搜广推等BAT大佬

    题解-MtOI2019 幽灵乐团

原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/368648

非原创文章文中已经注明原地址,如有侵权,联系删除

关注公众号【高性能架构探索】,第一时间获取最新文章

转载文章受原作者版权保护。转载请注明原作者出处!

(0)
上一篇 2023年3月2日 下午6:31
下一篇 2023年3月2日 下午6:32

相关推荐