「MtOI2019」幽灵乐团

Description

求$$\prod _{i=1}^A\prod _{j=1}^B\prod _{k=1}^C(\frac{[i,j]}{(i,k)}{)}^{f(type)}$$
$f(0)=1,f(1)=ijk,f(2)=(i,j,k)$
设T为数据组数
A,B,C<=10^5,T<=70

Solution

被车万骗进来的我做了一天的数学题
三合一差评
公式太多了会跳步，建议手推
所有的除法都是下取整
M表示min(a,b,c)或min(a,b)

显然分成两部分，一部分是${\prod }_{i=1}^A{\prod }_{j=1}^B{\prod }_{k=1}^C{i}^{f(type)}$
另一部分是${\prod }_{i=1}^A{\prod }_{j=1}^B{\prod }_{k=1}^C(i,j{)}^{f(type)}$
考虑第一部分
type=1直接预处理阶乘快速幂
type=2也是预处理${\prod }_{i=1}^n{i}^i$快速幂
type=3推式子：
${\prod }_{i=1}^A{\prod }_{j=1}^B{\prod }_{k=1}^C{i}^{(i,j,k)}$
${\prod }_{d=1}^M{\prod }_{i=1}^{A/d}(id{)}^{d{\sum }_{j=1}^{B/d}{\sum }_{k=1}^{C/d}[(i,j,k)=1]}$
${\prod }_{x=1}^M{\prod }_{d=1}^{M/x}{\prod }_{i=1}^{A/xd}(idx{)}^{d\mu (x)(B/xd)(C/xd)}$
${\prod }_{T=1}^M{\prod }_{i=1}^{A/T}(iT{)}^{{\sum }_{d|T}d\mu (T/x)(B/T)(C/T)}$
${\prod }_{T=1}^M{\prod }_{i=1}^{A/T}(iT{)}^{\phi (T)(B/T)(C/T)}$
预处理${\prod }_{i=1}^n{i}^{\phi (i)}$+分块即可

考虑第二部分
type=1求${\prod }_{i=1}^A{\prod }_{j=1}^B(i,j)$再快速幂
${\prod }_{d=1}^M{d}^{{\sum }_{i=1}^{A/d}{\sum }_{j=1}^{B/d}[(i,j)=1]}$
${\prod }_{x=1}^M{\prod }_{d=1}^{M/x}{d}^{(A/xd)(B/xd)\mu (x)}$
${\prod }_{T=1}^M{\prod }_{d|T}{d}^{\mu (T/d)(A/T)(B/T)}$
n log n预处理+分块
type=2同理，不再过多阐述
type=3
${\prod }_{i=1}^A{\prod }_{j=1}^B{\prod }_{k=1}^C(i,j{)}^{(i,j,k)}$
${\prod }_{d=1}^{min(A,B)}{d}^{({\sum }_{k=1}^C(k,d))({\sum }_{i=1}^{A/d}{\sum }_{j=1}^{B/d}[(i,j)=1)]}$
${\prod }_{x=1}^{min(A,B)}{\prod }_{d=1}^{min(A,B)/x}{d}^{\mu (x)({\sum }_{k=1}^C(k,d))(A/xd)(B/xd)}$
考虑${\sum }_{k=1}^C(k,d)$
${\sum }_{T|d}T{\sum }_{k=1}^{C/T}[(k,d/T)=1]$
${\sum }_{x|d}\mu (x){\sum }_{T|\frac{d}{x}}T(C/Tx)$
${\sum }_{i|d}(C/i){\sum }_{x|i}\mu (x)\frac{i}{x}$
${\sum }_{i|d}(C/i)\phi (i)$
原式=${\prod }_{i=1}^M{\prod }_{x=1}^{min(A,B)/i}{\prod }_{d=1}^{min(A,B)/ix}(id{)}^{(C/i)\phi (i)\mu (x)(A/idx)(B/idx)}$
分成两部分
${\prod }_{i=1}^M{i}^{(C/i)\phi (i){\sum }_{x=1}^{min(A,B)/i}{\sum }_{d=1}^{min(A,B)/ix}\mu (x)(A/idx)(B/idx)}$
考虑${\sum }_{x=1}^{min(A,B)/i}{\sum }_{d=1}^{min(A,B)/ix}\mu (x)(A/idx)(B/idx)$
${\sum }_{T=1}^{min(A,B)/i}(A/Ti)(B/Ti){\sum }_{d|T}\mu (d)$
后面那个显然只有T=1为1
所以这一部分就是${\prod }_{i=1}^M{i}^{(A/i)(B/i)(C/i)\phi (i)}$，直接分块
${\prod }_{i=1}^M{\prod }_{x=1}^{min(A,B)/i}{\prod }_{d=1}^{min(A,B)/ix}{d}^{(C/i)\phi (i)\mu (x)(A/idx)(B/idx)}$
${\prod }_{i=1}^M{\prod }_{T=1}^{min(A,B)/i}{\sum }_{d|T}{d}^{\mu (T/d)(A/iT)(B/iT)(C/i)\phi (i)}$，分块套分块即可

如果写错了什么那也是无可奈何啊(摊手

Code

#include <cstdio>
#include <cstring>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
using namespace std;

typedef long long ll;

const int N=1e5+5;

int q,Mo,p[N],mu[N],phi[N];
bool vis[N];
ll f1[N],t1[N],f2[N],t2[N],s[N],f3[N],sp[N],Inv_t1[N],Inv_t2[N],a,b,c;

ll pwr(ll x,ll y) {
	y%=Mo-1;if (y<0) y+=Mo-1;
	ll z=1;
	for(;y;y>>=1,x=x*x%Mo)
		if (y&1) z=z*x%Mo;
	return z;
}

ll Inv(ll x) {return pwr(x,Mo-2);}

ll S(ll x) {return x*(x+1)/2%(Mo-1);}

void pre(int N) {
	mu[1]=phi[1]=1;
	fo(i,2,N) {
		if (!vis[i]) p[++p[0]]=i,mu[i]=-1,phi[i]=i-1;
		fo(j,1,p[0]) {
			int k=i*p[j];if (k>N) break;
			vis[k]=1;if (!(i%p[j])) {phi[k]=phi[i]*p[j];break;}
			mu[k]=-mu[i];phi[k]=phi[i]*(p[j]-1);
		}
	}

	f1[0]=f2[0]=f3[0]=1;
	fo(i,1,N) {
		f1[i]=f1[i-1]*i%Mo;
		f2[i]=f2[i-1]*pwr(i,i)%Mo;
		f3[i]=f3[i-1]*pwr(i,phi[i])%Mo;
	}

	fo(i,0,N) t1[i]=t2[i]=1;
	fo(d,1,N) fo(x,1,N/d) {
		int now=pwr(d,mu[x]);
		t1[d*x]=t1[d*x]*now%Mo;
		t2[d*x]=t2[d*x]*now%Mo;
	}
	Inv_t1[0]=Inv_t2[0]=1;
	fo(i,1,N) {
		t1[i]=t1[i-1]*t1[i]%Mo;
		t2[i]=t2[i-1]*pwr(t2[i],(ll)i*i)%Mo;
		Inv_t1[i]=Inv(t1[i]);Inv_t2[i]=Inv(t2[i]);
		sp[i]=(sp[i-1]+phi[i])%(Mo-1);
	}
}

namespace Type_1{
	ll calc(ll n,ll m) {
		if (n>m) swap(n,m);
		ll ans=1;
		for(ll l=1,r=0;l<=n;l=r+1) {
			r=min(n/(n/l),m/(m/l));
			ans=ans*pwr(t1[r]*Inv_t1[l-1]%Mo,(n/l)*(m/r))%Mo;
		}
		return ans;
	}
}

namespace Type_2{
	ll calc(ll n,ll m) {
		if (n>m) swap(n,m);
		ll ans=1;
		for(ll l=1,r=0;l<=n;l=r+1) {
			r=min(n/(n/l),m/(m/l));
			ans=ans*pwr(t2[r]*Inv_t2[l-1]%Mo,S(n/l)*S(m/l))%Mo;
		}
		return ans;
	}
}

namespace Type_3{
	ll calc(ll a,ll b) {
		if (a>b) swap(a,b);
		ll ans=1;
		for(ll l=1,r=0;l<=a;l=r+1) {
			r=min(a/(a/l),b/(b/l));
			ans=ans*pwr(t1[r]*Inv_t1[l-1]%Mo,(a/l)*(b/l))%Mo;
		}
		return ans;
	}
	ll calc_1(ll a,ll b,ll c) {
		ll ans=1;int m=min(min(a,b),c);
		for(ll l=1,r=0;l<=m;l=r+1) {
			r=min(min(a/(a/l),b/(b/l)),c/(c/l));
			ll ret=(sp[r]-sp[l-1])%(Mo-1);
			ans=ans*pwr(f1[a/l],(b/l)*(c/l)%(Mo-1)*ret)%Mo;
		}
		return ans;
	}
	ll calc_2(ll a,ll b,ll c) {
		ll ans=1;int m=min(min(a,b),c);
		for(ll l=1,r=0;l<=m;l=r+1) {
			r=min(min(a/(a/l),b/(b/l)),c/(c/l));
			ans=ans*pwr(calc(a/l,b/l),(sp[r]-sp[l-1])*(c/l))%Mo;
		}
		return ans;
	}
}

int main() {
	scanf("%d%d",&q,&Mo);
	pre(1e5);
	for(;q;q--) {
		scanf("%lld%lld%lld",&a,&b,&c);
		ll ans=pwr(f1[a],b*c)%Mo*pwr(f1[b],a*c)%Mo;
		ans=ans*Inv(pwr(Type_1::calc(a,b),c))%Mo;
		ans=ans*Inv(pwr(Type_1::calc(a,c),b))%Mo;
		printf("%lld ",ans);
		ans=pwr(f2[a],S(b)*S(c))*pwr(f2[b],S(a)*S(c))%Mo;
		ans=ans*Inv(pwr(Type_2::calc(a,b),S(c)))%Mo;
		ans=ans*Inv(pwr(Type_2::calc(a,c),S(b)))%Mo;
		printf("%lld ",ans);
		ans=Type_3::calc_1(a,b,c)*Type_3::calc_1(b,a,c)%Mo;
		ans=ans*Inv(Type_3::calc_2(a,b,c))%Mo;
		ans=ans*Inv(Type_3::calc_2(a,c,b))%Mo;
		printf("%lld\n",ans);
	}
	return 0;
}