[MtOI2019]幽灵乐团

题目

一个很暴力的辣鸡做法

考虑到两个数的(gcd)是所有质数次幂取(min)，两个数的( m lcm)是所有质数次幂取(max)，于是最后的答案一定是(prod p_i^{c_i})，而且这里最大的质数不会超过(n)，于是我们考虑算出每一个质数的次幂是多少

于是我们成功的将(prod)转换成了(sum)，指数上得对( m mod-1)取模

对于( m type=0)，我们枚举一个质数(p)，计算它的次幂

即为

[sum_{i=0}sum_{j=0}sum_{k=0}sum_{t_1=1}^{left lfloor frac{A}{p^i} ight floor}sum_{t_2=1}^{left lfloor frac{B}{p^j} ight floor}sum_{t_3=1}^{left lfloor frac{C}{p^k} ight floor}max(i,j)-min(i,k)[t_1perp p][t_2perp p][t_3perp p] ]

这个很显然啊，就是我们枚举当前这个质数的最大次幂，之后枚举这个拥有这个最大次幂的数，即为(p^i imes t_1)，但是如果(t_1)和(p)不互质，那么(i)就不是(p^i imes t)的最大次幂，于是还必须得保证(t_1perp p)

我们改变一下求和顺序即为

[sum_{i=0}sum_{j=0}sum_{k=0}max(i,j)-min(i,k)sum_{t_1=1}^{left lfloor frac{A}{p^i} ight floor}[t_1perp p]sum_{t_2=1}^{left lfloor frac{B}{p^j} ight floor}[t_2perp p]sum_{t_3=1}^{left lfloor frac{C}{p^k} ight floor}[t_3perp p] ]

又发现求(max(i,j))和(min(i,k))是相当独立的，于是我们可以强行拆成

[C imes sum_{i=0}sum_{j=0}max(i,j)sum_{t_1=1}^{left lfloor frac{A}{p^i} ight floor}[t_1perp p]sum_{t_2=1}^{left lfloor frac{B}{p^j} ight floor}[t_2perp p] ]

减去

[B imes sum_{i=0}sum_{k=0}min(i,k)sum_{t_1=1}^{left lfloor frac{A}{p^i} ight floor}[t_1perp p]sum_{t_3=1}^{left lfloor frac{C}{p^k} ight floor}[t_3perp p] ]

现在思考一下这个(sum_{t_1=1}^{left lfloor frac{A}{p^i} ight floor}[t_1perp p])到底怎么求，和一个质数互质只要不是这个质数的倍数就行了，于是我们简单的减去(1)到(left lfloor frac{A}{p^i} ight floor)里(p)的倍数即可，也就是(left lfloor frac{A}{p^i} ight floor-left lfloor frac{A}{p^{i+1}} ight floor)

也就是我们现在只需要求

[B imes sum_{i=0}sum_{j=0}max(i,j)(left lfloor frac{A}{p^i} ight floor-left lfloor frac{A}{p^{i+1}} ight floor)(left lfloor frac{B}{p^j} ight floor-left lfloor frac{B}{p^{j+1}} ight floor) ]

我们可以直接(O(log_p^2 A))来计算这个柿子，但是我们处理一下(left lfloor frac{A}{p^i} ight floor-left lfloor frac{A}{p^{i+1}} ight floor)的前缀和强行讨论一下也可以做到(O(log_p A))

一到(n)内的质数大概有(frac{n}{ln n})个，每一个质数的计算需要(O(log_p n))的复杂度，于是做一遍复杂度大概就是(O(n))的

对于( m type=1)，柿子变成了

[sum_{i=0}sum_{j=0}sum_{k=0}sum_{t_1=1}^{left lfloor frac{A}{p^i} ight floor}sum_{t_2=1}^{left lfloor frac{B}{p^j} ight floor}sum_{t_3=1}^{left lfloor frac{C}{p^k} ight floor}(max(i,j)-min(i,k)) imes t_1 imes t_2 imes t_3 imes p^{i+j+k}[t_1perp p][t_2perp p][t_3perp p] ]

强行一搞也能变成上面的形式，同样可以在(O(n))时间求出来

对于( m type=2)，我们考虑枚举上面的(gcd)

设(f(A,B,C))为第一问的答案，也就是三个上标分别取到(A,B,C)时的贡献

套路枚举(gcd)为(d)即(d)的倍数时的贡献

[sum_{d=1}^{min(A,B,C)}dsum_{d|i}mu(frac{i}{d})f(left lfloor frac{A}{d} ight floor,left lfloor frac{B}{d} ight floor,left lfloor frac{C}{d} ight floor) ]

就是(left lfloor frac{A}{d} ight floor,left lfloor frac{B}{d} ight floor,left lfloor frac{C}{d} ight floor)内的数乘上(d)分别不会超过(A,B,C)，且这些数的(gcd)都是(d)或(d)的倍数

再套路交换一下求和顺序

[sum_{i=1}^{min(A,B,C)}f(left lfloor frac{A}{i} ight floor,left lfloor frac{B}{i} ight floor,left lfloor frac{C}{i} ight floor)sum_{d|i}mu(frac{i}{d})d=sum_{i=1}^{min(A,B,C)}f(left lfloor frac{A}{i} ight floor,left lfloor frac{B}{i} ight floor,left lfloor frac{C}{i} ight floor)varphi(i) ]

于是直接大力整除分块，复杂度应该不会超过(O(nln n))

吧。

请大力卡常

代码

#include<bits/stdc++.h>
#define re register
#define LL long long
#define max(a,b) (a>b?a:b)
#define min(a,b) (a<b?a:b)
#pragma GCC optimize(3)
#pragma GCC optimize("-fcse-skip-blocks")
const int maxn=1e5+5;
int mod,P,T,R;
int is[maxn],p[maxn>>1],pre[maxn],phi[maxn],pw[maxn>>1];
int r[2],tax[2][20],pr[2][20];
inline int ksm(int a,int b) {
	int S=1;
	for(;b;b>>=1,a=1ll*a*a%mod) if(b&1) S=1ll*S*a%mod;
	return S;
}
inline int qm(int x) {return x>=P?x-P:x;}
inline int dqm(int x) {return x+=x>>31&P;}
inline int calc(int n,int p) {
	return dqm(pre[n]-1ll*pre[n/p]*p%P);
}
inline int out(int A) {
	int ans=1;
	for(re int i=1;i<=p[0]&&p[i]<=A;i++) 
		ans=1ll*ans*ksm(p[i],pw[i])%mod;
	return ans;
}
inline void clear(int A) {
	for(re int i=1;i<=p[0]&&p[i]<=A;i++) pw[i]=0;
}
inline void Solve(int A,int B,int C,int v) {
	int ans1=1;
	for(re int t=1;t<=p[0];t++) {
		if(p[t]>A&&p[t]>B&&p[t]>C) break;
		int now=0,tot=0;
		for(re int i=0,aa=A;aa;i++,aa/=p[t]) pr[0][i]=tax[0][i]=aa-aa/p[t],r[0]=i;
		for(re int j=0,bb=B;bb;j++,bb/=p[t]) pr[1][j]=tax[1][j]=bb-bb/p[t],r[1]=j;
		R=max(r[0],r[1]);
		for(re int i=1;i<=R;i++) pr[0][i]=qm(pr[0][i]+pr[0][i-1]);
		for(re int j=1;j<=R;j++) pr[1][j]=qm(pr[1][j]+pr[1][j-1]);
		for(re int i=0;i<=r[0];i++) 
			now=qm(now+1ll*tax[0][i]*pr[1][i]*i%P);
		for(re int j=1;j<=r[1];j++) 
			now=qm(now+1ll*tax[1][j]*pr[0][j-1]*j%P);
		for(re int i=1;i<=R;i++) pr[1][i]=0;
		for(re int i=r[0]+1;i<=R;i++) pr[0][i]=0;	
		now=1ll*now*C%P;
		for(re int k=0,cc=C;cc;k++,cc/=p[t]) pr[1][k]=tax[1][k]=cc-cc/p[t],r[1]=k;
		R=max(r[0],r[1]);
		for(re int i=r[0]+1;i<=R;i++) pr[0][i]=qm(pr[0][i-1]+pr[0][i]);
		for(re int k=1;k<=R;k++) pr[1][k]=qm(pr[1][k]+pr[1][k-1]);
		for(re int i=1;i<=r[0];i++) 
			tot=qm(tot+1ll*i*tax[0][i]*dqm(pr[1][r[1]]-pr[1][i-1])%P);
		for(re int k=1;k<=r[1];k++)
			tot=qm(tot+1ll*k*tax[1][k]*dqm(pr[0][r[0]]-pr[0][k])%P);
		for(re int i=0;i<=R;i++) pr[1][i]=pr[0][i]=0;
		now=dqm(now-1ll*B*tot%P);
		pw[t]=qm(pw[t]+1ll*now*v%P);
	}
}
inline int Calc(int n,int m,int h) {
	int U=min(min(n,m),h);
	for(re int l=2,r;l<=U;l=r+1) {
		r=min(n/(n/l),m/(m/l));
		r=min(r,h/(h/l));
		Solve(n/l,m/l,h/l,dqm(phi[r]-phi[l-1]));
	}
	return out(max(max(n,m),h));
}
int A,B,C;
int main() {
	scanf("%d%d",&T,&mod);P=mod-1;phi[1]=1;
	for(re int i=2;i<maxn;i++) {
		if(!is[i]) p[++p[0]]=i,phi[i]=i-1;
		for(re int j=1;j<=p[0]&&p[j]*i<maxn;j++) {
			is[p[j]*i]=1;if(i%p[j]==0) {
				phi[p[j]*i]=phi[i]*p[j];break;
			}
			phi[p[j]*i]=phi[i]*phi[p[j]];
		}
	}
	for(re int i=1;i<maxn;i++) 
		pre[i]=qm(pre[i-1]+i),phi[i]=qm(phi[i-1]+phi[i]);
	while(T--) {
		scanf("%d%d%d",&A,&B,&C);
		int ans2=1;
		for(re int t=1;t<=p[0];t++) {
			if(p[t]>A&&p[t]>B&&p[t]>C) break;
			int now=0,tot=0;
			for(re LL i=0,a=1;a<=A;i++,a*=p[t])
				for(re LL j=0,b=1;b<=B;j++,b*=p[t]) {
					int v=1ll*max(i,j)*a*b%P;
					tot=qm(tot+1ll*v*calc(A/a,p[t])*calc(B/b,p[t])%P*pre[C]%P);
			}
			for(re LL i=1,a=p[t];a<=A;i++,a*=p[t])
				for(re LL k=1,c=p[t];c<=C;k++,c*=p[t]) {
					int v=1ll*min(i,k)*a*c%P;
				    now=qm(now+1ll*v*calc(A/a,p[t])*calc(C/c,p[t])%P*pre[B]%P);
			}
			ans2=1ll*ans2*ksm(p[t],dqm(tot-now))%mod;
		}
		Solve(A,B,C,1);printf("%d ",out(max(max(A,B),C)));
		printf("%d %d
",ans2,Calc(A,B,C));clear(max(max(A,B),C));
	}
}