东风谷 早苗(Kochiya Sanae)非常喜欢幽灵乐团的演奏，她想对她们的演奏评分。
因为幽灵乐团有 (3) 个人，所以我们可以用 (3) 个正整数 (A,B,C) 来表示出乐团演奏的分数，她们的演奏分数可以表示为
[prod_{i=1}^{A}prod_{j=1}^{B}prod_{k=1}^{C}left(frac{ ext{lcm}(i,j)}{gcd(i,k)}
ight)^{f(type)}
]
因为音乐在不同的部分会有不同的听觉感受，所以 (type) 会在 ({0,1,2}) 中发生变化，其中：
[f(0)=1
]
[f(1)=i imes j imes k
]
[f(2)=gcd(i,j,k)
]
因为乐团的歌实在太好听了，导致分数特别高，所以她们的分数要对给定的正整数 (p) 取模。
因为有很多歌曲要演奏，所以早苗给出了 (T) 组询问。
[1leq A,B,Cleq 10^5 10^7 leq p leq 1.05 imes 10^9 pin { prime} T =70
]
刚开始写了个 (n^frac{7}{8}log n) 的，直接暴毙/kk
[egin{aligned}
&prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c}left(frac{ ext{lcm}(i,j)}{gcd(i,k)}
ight)^{f(type)}\=&prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c}frac{(ij)^{f(type)}}{(gcd(i,j)gcd(i,k))^{f(type)}}end{aligned}]
上面和下面的就可以分开考虑，我们对 (f(type)) 分情况做。
注意到 (prod gcd(i,j)) 和 (prod gcd(i,k)) 本质上是等价的，所以我们后面只讨论 (gcd(i,j)) 。
type = 1
[prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c}ij=(a!)^{bc}(b!)^{ac}
]
对于下面：
[egin{aligned}&prod_{i=1}^aprod_{j=1}^bprod_{k=1}^cgcd(i,j)\=&prod_{i=1}^aprod_{j=1}^bgcd(i,j)^kend{aligned}
]
直接改为枚举 (gcd(i,j)) ，就有：
[prod_{i=1} i^{kg(a,b,i)}
]
其中 (g(a,b,x)) 表示 (sum_{i=1}^asum_{j=1}^b[gcd(i,j)==x]) ，直接对这个反演：
[egin{aligned}&=sum_{i=1}^{frac{a}{x}}sum_{j=1}^{frac{b}{x}}[gcd(i,j)==1]\&=sum_{i=1}^{frac{a}{x}}sum_{j=1}^{frac{b}{x}}sum_{d|i,d|j}mu(d)\&=sum_{d=1}lfloorfrac{a}{d}
floorlfloorfrac{b}{d}
floormu(d)end{aligned}
]
然后对于上式就有：
[prod_{i=1}i^{kg(frac{a}{i},frac{b}{i},1)}
]
这样就可以整除分块套整除分块了，单次 (O(n^frac{3}{4}log n)) 。
type = 2
对于上面：
[egin{aligned}&prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c}(ij)^{ijk}\=&prod_{i=1}^{a}prod_{j=1}^{b}(ij)^{ijsum k} (sum k=sum_{x=1}^{k}x)\=&prod_{i=1}^ai^{isum jsum k}prod_{j=1}j^{jsum isum k}end{aligned}
]
(O(nlog n)) 预处理 (f_n=prod_{i=1}^ni^i) 就可以 (O(n)) 计算了。
对于下面：
[egin{aligned}&prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c}gcd(i,j)^{ijk}\=&prod_{i=1}^aprod_{j=1}^bgcd(i,j)^{ijsum k}end{aligned}
]
又可以直接枚举 (gcd) 了。
[egin{aligned}prod_{i=1}i^{sum kg'(a,b,i)}end{aligned}
]
(g'(a,b,x)) 表示 (sum_{i=1}^asum_{j=1}^bij[gcd(i,j)==x])
其实和上面的还是基本一样的套路，注意到 (g'(a,b,x)= x^2g'(frac{a}{x},frac{b}{x},1)) ，所以可以写成这样：
[egin{aligned}&x^2sum_{d=1}d^2sum_{i=1}^{frac{a}{d}}sum_{j=1}^{frac{b}{d}}ijmu(d)\=&x^2sum_{d=1}d^2mu(d)sum frac{a}{d}sum frac{b}{d}end{aligned}
]
然后就可以整除分块套整除分块了，预处理出来 (sum_{i}i^2mu(d)) 就可以单次 (O(n^frac{3}{4}log n)) 计算了。
type = 3
对于上面：
[egin{aligned}&prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c}(ij)^{gcd(i,j,k)}\=&prod_{i=1}^{a}prod_{j=1}^{b}(ij)^{sum_{k=1}^cgcd(i,j,k)}\=&prod_{i=1}^ai^{sum_{j=1}^bsum_{k=1}^cgcd(i,j,k)}prod_{j=1}^bj^{sum_{i=1}^asum_{k=1}^cgcd(i,j,k)}end{aligned}
]
发现左右两边是一样的，所以只讨论左边的值，直接对幂反演：
[egin{aligned}G(x)&=sum_{j=1}^bsum_{k=1}^cgcd(x,j,k)\&=sum_{j=1}^bsum_{k=1}^csum_{d|x,d|j,d|k}varphi(x)\&=sum_{d|x}varphi(d)lfloorfrac{b}{d}
floorlfloorfrac{c}{d}
floorend{aligned}
]
带回去原式（只有左半边）：
[egin{aligned}=&prod_{i=1}^ai^{sum_{d|x}varphi(d)lfloorfrac{b}{d}
floorlfloorfrac{c}{d}
floor}\=&prod_{d=1}^aprod_{i=1}^frac{a}{d}(id)^{varphi(d)lfloorfrac{c}{d}
floorlfloorfrac{b}{d}
floor}\=&prod_{d=1}^{a}(lfloorfrac{a}{d}
floor!)^{varphi(d)lfloorfrac{c}{d}
floorlfloorfrac{b}{d}
floor} imes d^{{varphi(d)lfloorfrac{b}{d}
floorlfloorfrac{c}{d}
floorlfloorfrac{b}{d}
floor}}end{aligned}
]
这样子两边就可以分别整除分块了， (O(nlog n)) 预处理出来 (n!) 和(g_n=prod_{i=1}^ni^{varphi(i)}) 就可以单次 (sqrt{n}log n) 计算了。
然后对于下面：
[egin{aligned}&prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c}gcd(i,j)^{gcd(i,j,k)}end{aligned}
]
如果枚举 (gcd(i,j)) 就会变成跟我刚开始一样的极丑复杂度，所以我们考虑枚举 (gcd(i,j,k)) ，就有：
[egin{aligned}&prod_{i=1}^{a}prod_{j=1}^{b}prod_{k=1}^{c}gcd(i,j)^{sum_{d=1}d[gcd(i,j,k)==d]}\=&prod_{d=1}prod_{i=1}^{frac{a}{d}}prod_{j=1}^{frac{b}{d}}(dgcd(i,j))^{dsum_{k=1}^frac{c}{d}[gcd(i,j,k)==1]}\=&prod_{d=1}prod_{t=1}prod_{i=1}^{frac{a}{td}}prod_{j=1}^{frac{b}{td}}(tdgcd(i,j))^{dmu(t)lfloorfrac{c}{td}
floor}end{aligned}
]
把 (td) 和 (gcd(i,j)) 分开考虑，对于 (td) 部分，就有：
[egin{aligned}&prod_{d=1}prod_{t=1}prod_{i=1}^frac{a}{td}prod_{j=1}^frac{b}{td}(td)^{dmu(t)lfloorfrac{c}{td}
floor}\=&prod_{T=1}prod_{i=1}^{frac{a}{T}}prod_{i=1}^{frac{b}{T}}T^{sum_{t|T}tmu(frac{T}{t})lfloorfrac{c}{T}
floor}end{aligned}
]
注意到 (mu*id=varphi) ，所以就有：
[prod_{T=1}T^{varphi(T){lfloorfrac{a}{T}
floor}lfloorfrac{b}{T}
floor}
]
根据之前预处理出来的东西可以单次 (O(sqrt{n}log n)) 计算。
然后再考虑 (gcd(i,j)) 部分，就有：
[egin{aligned}&prod_{d=1}prod_{t=1}prod_{i=1}^frac{a}{td}prod_{j=1}^frac{b}{td}gcd(i,j)^{dmu(t)lfloorfrac{c}{td}
floor}\=&prod_{d=1}prod_{t=1}prod_{i=1}^frac{a}{td}prod_{j=1}^frac{b}{td}prod_{g=1}g^{sum_{h=1}mu(h)dmu(t)lfloorfrac{c}{td}
floor}\=&prod_{d=1}prod_{t=1}prod_{h=1}prod_{g=1}gh^{dmu(h)lfloorfrac{a}{tdgh}
floorlfloorfrac{b}{tdgh}
floormu(t)lfloorfrac{c}{td}
floor}\=&prod_{T=1(td)}(prod_{t=1(gh)}(prod_{d|t}d^{mu(frac{t}{d})})^{lfloorfrac{a}{tT}
floorlfloorfrac{b}{tT}
floor})^{sum_{k|T}mu(k)frac{T}{k}lfloorfrac{c}{T}
floor}\=&prod_{T=1}(prod_{t=1}(prod_{d|t}d^{mu(frac{t}{d})})^{lfloorfrac{a}{tT}
floorlfloorfrac{b}{tT}
floor})^{varphi(T)lfloorfrac{c}{T}
floor}end{aligned}
]
设 (G'_t=prod_{d|t}d^mu(frac{t}{d})) ，这个可以 (O(nlog n)) 时间预处理出来，然后式子就变成了：
[prod_{T=1}(prod_{t=1}(G'_t)^{lfloorfrac{a}{tT}
floorlfloorfrac{b}{tT}
floor})^{varphi(T)lfloorfrac{c}{T}
floor}
]
这个东西就又是一个整除分块套整除分块了，可以单词 (O(n^frac{3}{4}log n)) 计算。
Code
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std;
const int N = 1e5;
const int M = 400;
int T,p,a,b,c,fac[N + 5],prime[N + 5],pcnt,vis[N + 5],is[N + 5],mul[N + 5],fp,sm2[N + 5],ifm[N + 5],ipd[N + 5],inv[N + 5],fi[N + 5],fii[N + 5],phi[N + 5],sp[N + 5],f[N + 5],pd[N + 5],g[N + 5],fm[N + 5];
int mypow(int a,int x){int s = 1;for (;x;x & 1 ? s = 1ll * s * a % p : 0,a = 1ll * a * a % p,x >>= 1);return s;}
int gcd(int a,int b)
{
if (!b)
return a;
return gcd(b,a % b);
}
void prework()
{
fac[0] = 1;
for (int i = 1;i <= N;i++)
fac[i] = 1ll * fac[i - 1] * i % p;
inv[0] = inv[1] = 1;
for (int i = 2;i <= N;i++)
inv[i] = 1ll * inv[p % i] * (p - p / i) % p;
is[0] = 1;
for (int i = 1;i <= N;i++)
is[i] = 1ll * is[i - 1] * inv[i] % p;
vis[1] = mul[1] = phi[1] = 1;
for (int i = 2;i <= N;i++)
{
if (!vis[i])
{
prime[++pcnt] = i;
mul[i] = -1;
phi[i] = i - 1;
}
for (int j = 1;j <= pcnt && prime[j] * i <= N;j++)
{
vis[prime[j] * i] = 1;
mul[prime[j] * i] = -mul[i];
phi[prime[j] * i] = phi[i] * (prime[j] - 1);
if (i % prime[j] == 0)
{
mul[prime[j] * i] = 0;
phi[prime[j] * i] = phi[i] * prime[j];
break;
}
}
}
fp = p - 1;
for (int i = 0;i <= N;i++)
fm[i] = 1;
for (int i = 1;i <= N;i++)
for (int j = 1;j * i <= N;j++)
{
if (mul[j] == -1)
fm[i * j] = 1ll * fm[i * j] * inv[i] % p;
if (mul[j] == 1)
fm[i * j] = 1ll * fm[i * j] * i % p;
}
for (int i = 1;i <= N;i++)
fm[i] = 1ll * fm[i - 1] * fm[i] % p;
for (int i = 0;i <= N;i++)
ifm[i] = mypow(fm[i],p - 2);
for (int i = 1;i <= N;i++)
{
sm2[i] = (sm2[i - 1] + 1ll * i * i * mul[i] % fp) % fp;
mul[i] += mul[i - 1];
}
fi[0] = 1;
for (int i = 1;i <= N;i++)
fi[i] = 1ll * fi[i - 1] * mypow(i,i) % p;
fii[0] = 1;
for (int i = 1;i <= N;i++)
fii[i] = 1ll * fii[i - 1] * mypow(mypow(i,i),i) % p;
for (int i = 1;i <= N;i++)
sp[i] = (sp[i - 1] + phi[i]) % fp;
pd[0] = 1;
for (int i = 1;i <= N;i++)
pd[i] = 1ll * pd[i - 1] * mypow(i,phi[i]) % p;
for (int i = 0;i <= N;i++)
ipd[i] = mypow(pd[i],p - 2) % p;
}
int calc1(int n,int m)
{
int ans = 0;
for (int l = 1,r;l <= min(n,m);l = r + 1)
{
r = min(n / (n / l),m / (m / l));
ans += 1ll * (mul[r] - mul[l - 1]) * (n / l) % fp * (m / l) % fp;
ans %= fp;
}
ans = (ans + fp) % fp;
return ans;
}
void solve1(int a,int b,int c)
{
int ans = mypow(fac[a],1ll * b * c % fp),s = 1;
ans = 1ll * ans * mypow(fac[b],1ll * a * c % fp) % p;
for (int l = 1,r;l <= min(a,b);l = r + 1)
{
r = min(a / (a / l),b / (b / l));
s = 1ll * s * mypow(mypow(1ll * fac[r] * is[l - 1] % p,calc1(a / l,b / l)),c) % p;
}
for (int l = 1,r;l <= min(a,c);l = r + 1)
{
r = min(a / (a / l),c / (c / l));
s = 1ll * s * mypow(mypow(1ll * fac[r] * is[l - 1] % p,calc1(a / l,c / l)),b) % p;
}
ans = 1ll * ans * mypow(s,p - 2) % p;
ans = (ans + p) % p;
cout<<ans<<" ";
}
int sum1(int a){return 1ll * a * (a + 1) / 2 % fp;}
int sum2(int a){return 1ll * a * (a + 1) / 2 * (2 * a + 1) / 3 % fp;}
int calc2(int n,int m)
{
int ans = 0;
for (int l = 1,r;l <= min(n,m);l = r + 1)
{
r = min(n / (n / l),m / (m / l));
ans += 1ll * sum1(n / l) * sum1(m / l) % fp * (sm2[r] - sm2[l - 1]) % fp;
ans %= fp;
}
ans = (ans + fp) % fp;
return ans;
}
void solve2(int a,int b,int c)
{
int ans = 1,s = 1,fj = 1;
ans = 1ll * mypow(fi[a],sum1(b)) * mypow(fi[b],sum1(a)) % p;
ans = mypow(ans,sum1(c));
s = 1;
for (int l = 1,r;l <= min(a,b);l = r + 1)
{
r = min(a / (a / l),b / (b / l));
int now = calc2(a / l,b / l);
s = 1ll * s * mypow(1ll * fii[r] * mypow(fii[l - 1],p - 2) % p,1ll * now * sum1(c) % fp) % p;
}
for (int l = 1,r;l <= min(a,c);l = r + 1)
{
r = min(a / (a / l),c / (c / l));
int now = calc2(a / l,c / l);
s = 1ll * s * mypow(1ll * fii[r] * mypow(fii[l - 1],p - 2) % p,1ll * now * sum1(b) % fp) % p;
}
ans = 1ll * ans * mypow(s,p - 2) % p;
ans = (ans + p) % p;
cout<<ans<<" ";
}
int calc3(int a,int b,int c)
{
int ans = 1;
for (int l = 1,r;l <= min(b,min(a,c));l = r + 1)
{
r = min(a / (a / l),min(b / (b / l),c / (c / l)));
ans = 1ll * ans * mypow(fac[b / l],1ll * ((sp[r] - sp[l - 1]) % fp + fp) % fp * (a / l) % fp * (c / l) % fp) % p;
}
for (int l = 1,r;l <= min(a,min(b,c));l = r + 1)
{
r = min(a / (a / l),min(b / (b / l),c / (c / l)));
ans = 1ll * ans * mypow(1ll * pd[r] * ipd[l - 1] % p,1ll * (a / l) * (b / l) % fp * (c / l) % fp) % p;
}
return ans;
}
int calc5(int a,int b)
{
int ans = 1;
for (int l = 1,r;l <= min(a,b);l = r + 1)
{
r = min(a / (a / l),b / (b / l));
ans = 1ll * ans * mypow(1ll * fm[r] * ifm[l - 1] % p,1ll * (a / l) * (b / l) % fp) % p;
}
return ans;
}
int calc4(int a,int b,int c)
{
int ans = 1;
for (int l = 1,r;l <= min(a,min(b,c));l = r + 1)
{
r = min(a / (a / l),min(b / (b / l),c / (c / l)));
ans = 1ll * ans * mypow(1ll * pd[r] * ipd[l - 1] % p,1ll * (a / l) * (b / l) % fp * (c / l) % fp) % p;
}
for (int l = 1,r;l <= min(a,min(b,c));l = r + 1)
{
r = min(a / (a / l),min(b / (b / l),c / (c / l)));
ans = 1ll * ans * mypow(calc5(a / l,b / l),1ll * ((sp[r] - sp[l - 1]) % fp + fp) % fp * (c / l) % fp) % p;
}
return ans;
}
void solve3(int a,int b,int c)
{
int ans = 1,s = 1;
ans = 1ll * calc3(a,b,c) * calc3(b,a,c) % p;
s = 1ll * calc4(a,b,c) * calc4(a,c,b) % p;
ans = 1ll * ans * mypow(s,p - 2) % p;
cout<<ans<<endl;
}
int main()
{
scanf("%d%d",&T,&p);
prework();
while (T--)
{
scanf("%d%d%d",&a,&b,&c);
solve1(a,b,c);
solve2(a,b,c);
solve3(a,b,c);
}
return 0;
}