zoukankan      html  css  js  c++  java
  • P5518[MtOI2019]幽灵乐团【莫比乌斯反演,欧拉反演】

    正题

    题目链接:https://www.luogu.com.cn/problem/P5518


    题目大意

    \(T\)次给出\(A,B,C\)求以下三个式子

    \[\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^{C}\frac{lcm(i,j)}{gcd(i,k)} \]

    \[\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^{C}\left(\frac{lcm(i,j)}{gcd(i,k)}\right)^{i\times j\times k} \]

    \[\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^{C}\left(\frac{lcm(i,j)}{gcd(i,k)}\right)^{gcd(i,j,k)} \]

    \(1\leq T\leq 70,1\leq A,B,C\leq 10^5\)


    解题思路

    开始写了个\(O(Tn\log n)\)结果发现不能过,然后就多浪费了三个多小时
    在这里插入图片描述
    只需要用到两个反演的式子

    \[F(i)=\sum_{i|d}f(d)\Rightarrow f(i)=\sum_{i|d}\mu(\frac{d}{i})F(d) \]

    \[gcd(S)=\sum_{\forall x\in S,d|x}\varphi(d) \]

    然后因为推导过程出来冗长以外没有太多难的部分,所以推荐自己手推到不会的再翻题解。

    然后就开始吧,因为\(lcm(i,j)=\frac{ij}{gcd(i,j)}\)所以问题可以化为两个部分,求

    \[\left(\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^Cij\right)^{f(type)},\left(\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^C\frac{1}{gcd(i,j)}\right)^{f(type)} \]

    首先是第一个式子\(f(type)=1\)
    第一部分就是

    \[\prod_{i=1}^Ai^{B\times C}\times \prod_{j=1}^Bj^{A\times C} \]

    这个十分简单,我们预处理阶乘就可以做到\(O(\log P)\)
    然后第二部分考虑枚举约数

    \[\prod_{d=1}\frac{1}{d}^{C\times \sum_{i=1}^A\sum_{j=1}^B[gcd(i,j)=d]} \]

    然后莫反

    \[\prod_{d=1}\frac{1}{d}^{C\sum_{k=1}\mu(k)\lfloor\frac{A}{kd}\rfloor\lfloor\frac{B}{kd}\rfloor} \]

    显然的我们可以\(d,k\)都可以整除分块,预处理一下逆元的前缀乘积就可以快速计算区间逆元乘积了。

    第一个式子时间复杂度\(O(n^\frac{3}{4})\)

    然后第二个式子类似的,记\(S(n)=\frac{n\times (n+1)}{2}\)那么有

    \[\prod_{i=1}^Ai^{iS(B)S(C)}\times \prod_{j=1}^Bj^{jS(A)(C)} \]

    维护一个\(i^i\)的前缀积即可。
    第二部分

    \[\prod_{d=1}\frac{1}{d}^{S(C)\times \sum_{i=1}^A\sum_{j=1}^Bij[gcd(i,j)=d]}\Rightarrow \prod_{d=1}\frac{1}{d}^{S(C)\sum_{k=1}\mu(k)S(\lfloor\frac{A}{kd}\rfloor)S(\lfloor\frac{B}{kd}\rfloor)} \]

    需要提前处理\(\mu(i)\times i^2\)的前缀和,\(i^2,\frac{1}{i^2}\)的前缀积就好了

    第二个式子时间复杂度\(O(n^\frac{3}{4})\)

    主要的难点在第三个式子
    首先第一部分先只考虑\(i\)

    \[\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^Ci^{gcd(i,j)} \]

    考虑欧拉反演,枚举约数

    \[\prod_{d=1}\left(\prod_{i}^{\lfloor\frac{A}{d}\rfloor}i\times d\right)^{\varphi(d)\lfloor\frac{B}{d}\rfloor\lfloor\frac{C}{d}\rfloor} \]

    \(f_i=\prod_{j=1}^ii\)那么有

    \[\prod_{d=1}\left(f_{\lfloor\frac{A}{d}\rfloor}d^{\lfloor\frac{A}{d}\rfloor}\right)^{\varphi(d)\lfloor\frac{B}{d}\rfloor\lfloor\frac{C}{d}\rfloor} \]

    拆开来\(f\)\(d\)的部分

    \[\prod_{d=1}(f_{\lfloor\frac{A}{d}\rfloor})^{\varphi(d)\lfloor\frac{B}{d}\rfloor\lfloor\frac{C}{d}\rfloor}\times d^{\varphi(d)\lfloor\frac{A}{d}\rfloor\lfloor\frac{B}{d}\rfloor\lfloor\frac{C}{d}\rfloor} \]

    预处理出\(f\)数组,和\(g\)数组\(g_i=\prod_{j=1}^i j^{\varphi(j)}\)\(\varphi\)的前缀和就可以整除分块搞了

    然后是第二部分

    \[\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^C\frac{1}{gcd(i,j)}^{gcd(i,j,k)} \]

    你可以试一下枚举\(gcd(i,j)\),反正我推了好久都没有对出来/kk
    所以考虑枚举\(gcd(i,j,k)\)的约数然后欧拉反演里面再莫反

    \[\prod_{d=1}\left(\prod_{k=1}\frac{1}{dk}^{\sum_{z=1}\mu(z)\lfloor\frac{A}{dkz}\rfloor\lfloor\frac{B}{dkz}\rfloor}\right)^{\varphi(d)\lfloor\frac{C}{d}\rfloor} \]

    然后把\(\frac{1}{d}\)\(\frac{1}{k}\)分开处理

    \[\prod_{d=1}\left(\prod_{k=1}\frac{1}{k}^{\sum_{z=1}\mu(z)\lfloor\frac{A}{dkz}\rfloor\lfloor\frac{B}{dkz}\rfloor}\right)^{\varphi(d)\lfloor\frac{C}{d}\rfloor} \]

    \[\times \prod_{d=1}\left(\frac{1}{d}\right)^{\varphi(d)\lfloor\frac{C}{d}\rfloor\sum_{k=1}\sum_{z=1}\mu(z)\lfloor\frac{A}{dkz}\rfloor\lfloor\frac{B}{dkz}\rfloor} \]

    (至于为什么要这么分开,会注意到第一个式子如果我们进行两次整除分块,那么对于每个\(\frac{1}{dk}\)就会有两个影响它的指数\(\varphi(d)\)\(\sum_{z=1}\mu(z)\lfloor\frac{A}{dkz}\rfloor\lfloor\frac{B}{dkz}\rfloor\)所以很难处理,此时我们分开来就可以直接用区间的值来做了)

    注意到这样要做三次整除分块,十分地慢,但是考虑后面那个式子\(\sum_{z=1}\mu(z)\lfloor\frac{A}{dkz}\rfloor\lfloor\frac{B}{dkz}\rfloor\)对于一个固定的\(dk\)是一个确定的值的,并且因为是整除分块,所以不同的值不多我们可以考虑预处理这个东西,整除分块一个\(dk\)然后里面再整除分块计算就好了。

    还有拆开后处理\(\frac{1}{d}\)的式子需要预处理\(\frac{1}{d}^{\varphi(d)}\)的前缀积。

    然后就做完了,第三部分时间复杂度\(O(n+n^{\frac{3}{4}}\log n)\)

    实际上其实最后的式子两个部分可以约掉一些东西省一些常数,但是这样做也能过但是得开int。

    或者你可以用别的方法卡卡常反正我开int了。


    code

    因为中途long long改int所以代码巨丑

    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    //#define int long long
    using namespace std;
    const int N=1e5+10;
    int T,P,Phi,A,B,C,ans,cnt;bool v[N];
    int pri[N/10],mu[N],su[N],phi[N],f[N],g[N],h[N];
    int inv[N],inc[N],inw[N],fac[N],gac[N],hac[N];
    int power(int x,int b){
    	int ans=1;b=(b%Phi+Phi)%Phi;
    	while(b){
    		if(b&1)ans=ans*1ll*x%P;
    		x=x*1ll*x%P;b>>=1;
    	}
    	return ans; 
    }
    int Sum(int n)
    {return n*1ll*(n+1)/2%Phi;}
    int Tum(int n)
    {return n*1ll*(n+1)*1ll*(2ll*1ll*n+1)/6ll%Phi;} 
    void Prime(){
    	mu[1]=phi[1]=1;
    	for(int i=2;i<N;i++){
    		if(!v[i])pri[++cnt]=i,mu[i]=-1,phi[i]=i-1;
    		for(int j=1;j<=cnt&&i*1ll*pri[j]<N;j++){
    			v[i*1ll*pri[j]]=1;
    			if(i%pri[j]==0){
    				phi[i*1ll*pri[j]]=phi[i]*1ll*pri[j];
    				break;
    			}
    			mu[i*1ll*pri[j]]=-mu[i];
    			phi[i*1ll*pri[j]]=phi[i]*1ll*(pri[j]-1);
    		}
    	}
    	inc[1]=fac[0]=gac[0]=inw[0]=inv[0]=hac[0]=1;
    	for(int i=2;i<N;i++)inc[i]=P-inc[P%i]*1ll*(P/i)%P;
    	for(int i=1;i<N;i++)fac[i]=fac[i-1]*1ll*i%P;
    	for(int i=1;i<N;i++)gac[i]=gac[i-1]*1ll*power(i,i)%P;
    	for(int i=1;i<N;i++)hac[i]=hac[i-1]*1ll*power(i,i*1ll*i%Phi)%P;
    	for(int i=1;i<N;i++)inw[i]=inw[i-1]*1ll*power(inc[i],i*1ll*i%Phi)%P;
    	for(int i=1;i<N;i++)inv[i]=inv[i-1]*1ll*inc[i]%P;
    	for(int i=1;i<N;i++)su[i]=su[i-1]+mu[i];
    	g[0]=h[0]=1;
    	for(int i=1;i<N;i++)g[i]=g[i-1]*1ll*power(inc[i],phi[i])%P;
    	for(int i=1;i<N;i++)h[i]=h[i-1]*1ll*power(i,phi[i])%P;
    	for(int i=1;i<N;i++)(phi[i]+=phi[i-1])%=Phi;
    	return;
    }
    void qart1(int A,int B,int C){
    	int n=min(A,B);
    	for(int L=1,R;L<=n;L=R+1){
    		R=min(A/(A/L),B/(B/L));
    		int a=A/L,b=B/L,m=min(a,b),w=0;
    		for(int l=1,r;l<=m;l=r+1){
    			r=min(a/(a/l),b/(b/l));
    			(w+=(su[r]-su[l-1])*(a/l)*1ll*(b/l)%Phi)%=Phi;
    		}
    		ans=ans*1ll*power(inv[R]*1ll*fac[L-1]%P,w*1ll*C%Phi)%P;
    	}
    	return;
    }
    void part1(){
    	ans=power(fac[A],B*1ll*C%Phi)*1ll*power(fac[B],A*1ll*C%Phi)%P;
    	qart1(A,B,C);qart1(A,C,B);
    	printf("%d ",ans);return;
    }
    void qart2(int A,int B,int C){
    	int n=min(A,B);
    	for(int i=1;i<=n;i++)f[i]=(f[i-1]+mu[i]*1ll*i*1ll*i%Phi)%Phi;
    	for(int L=1,R;L<=n;L=R+1){
    		R=min(A/(A/L),B/(B/L));
    		int a=A/L,b=B/L,m=min(a,b),w=0;
    		for(int l=1,r;l<=m;l=r+1){
    			r=min(a/(a/l),b/(b/l));
    			(w+=Sum(a/l)*1ll*Sum(b/l)%Phi*1ll*(f[r]-f[l-1])%Phi)%=Phi;
    		}
    		w=w*1ll*Sum(C)%Phi;
    		ans=ans*1ll*power(inw[R]*1ll*hac[L-1]%P,w)%P;
    	}
    	return;
    }
    void part2(){
    	ans=power(gac[A],Sum(B)*1ll*Sum(C)%Phi)*1ll*power(gac[B],Sum(A)*1ll*Sum(C)%Phi)%P;
    	qart2(A,B,C);qart2(A,C,B);
    	printf("%d ",ans);return;
    }
    void qart3(int A,int B,int C){
    	int n=min(A,B);
    	for(int L=1,R;L<=n;L=R+1){
    		R=min(A/(A/L),B/(B/L));
    		int a=A/L,b=B/L,m=min(a,b),w=0;
    		for(int l=1,r;l<=m;l=r+1){
    			r=min(a/(a/l),b/(b/l));
    			(w+=1ll*(a/l)*(b/l)*(su[r]-su[l-1])%Phi)%=Phi;
    		}
    		for(int i=L;i<=R;i++)f[i]=w;
    	}
    	n=min(n,C);
    	for(int L=1,R;L<=n;L=R+1){
    		R=min(A/(A/L),min(B/(B/L),C/(C/L)));
    		int a=A/L,b=B/L,c=C/L,m=min(a,b),w=0,s=1;
    		for(int l=1,r;l<=m;l=r+1){
    			r=min(a/(a/l),b/(b/l));
    			(w+=1ll*f[l*1ll*L]*(r-l+1)%Phi)%=Phi; 
    			s=1ll*s*power(1ll*inv[r]*fac[l-1]%P,f[l*1ll*L])%P;
    		}
    		ans=ans*1ll*power(1ll*g[R]*power(g[L-1],P-2)%P,1ll*w*(C/L)%Phi)%P;
    		ans=ans*1ll*power(s,1ll*(phi[R]-phi[L-1])*(C/L)%Phi)%P;
    	}
    //	for(int i=1;i<=A;i++)
    //		for(int j=1;j<=B;j++)
    //			for(int k=1;k<=C;k++)
    //				ans=ans*1ll*power(inc[__gcd(i,j)],__gcd(__gcd(i,j),k))%P;
    	return;
    }
    void part3(){
    	ans=1;
    	int n=min(A,min(B,C));
    	for(int L=1,R;L<=n;L=R+1){
    		R=min(A/(A/L),min(B/(B/L),C/(C/L)));
    		int a=A/L,b=B/L,c=C/L;
    //		b=fac[b]*1ll*power(Sum(R)-Sum(L-1),b)%P;
    		ans=1ll*ans*power(fac[a],1ll*(phi[R]-phi[L-1])*b*c%Phi)%P;
    		ans=1ll*ans*power(fac[b],1ll*(phi[R]-phi[L-1])*a*c%Phi)%P;
    		ans=1ll*ans*power(h[R]*1ll*power(h[L-1],P-2)%P,2ll*a*b*c%Phi)%P;
    	}
    //	int k=ans;
    	qart3(A,B,C);qart3(A,C,B);
    //	else ans=1ll*ans*ans%P*power(k,P-2)%P;
    	printf("%d ",ans);return;
    }
    signed main()
    {
    	scanf("%d%d",&T,&P);
    	Phi=P-1;Prime();
    	while(T--){
    		scanf("%d%d%d",&A,&B,&C);
    //		A=B=C=1e5;
    		part1();
    		part2();
    		part3();
    		putchar('\n');
    	}
    	return 0;
    }
    
  • 相关阅读:
    怎么点击div之外的区域就隐藏这个div啊 找了很久,都没有很好解决
    ibatis 到 MyBatis区别
    MyBatis学习(一)一个简单的例子
    iBatis简单入门教程
    strut2的标签
    spring事务传播机制实例讲解
    ORACLE中Drop table cascade constraints之后果.
    oracle表的操作sql语句
    webService
    Oracle临时表
  • 原文地址:https://www.cnblogs.com/QuantAsk/p/15594319.html
Copyright © 2011-2022 走看看