洛谷 P5518 [MtOI2019]幽灵乐团

东风谷 早苗(Kochiya Sanae)非常喜欢幽灵乐团的演奏，她想对她们的演奏评分。

因为幽灵乐团有 \(3\) 个人，所以我们可以用 \(3\) 个正整数 \(A,B,C\) 来表示出乐团演奏的分数，她们的演奏分数可以表示为

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

因为音乐在不同的部分会有不同的听觉感受，所以 \(type\) 会在 \(\{0,1,2\}\) 中发生变化，其中：

\[f(0)=1 \]

\[f(1)=i \times j \times k \]

\[f(2)=\gcd(i,j,k) \]

因为乐团的歌实在太好听了，导致分数特别高，所以她们的分数要对给定的正整数 \(p\) 取模。

因为有很多歌曲要演奏，所以早苗给出了 \(T\) 组询问。

\[1\leq A,B,C\leq 10^5 \ \ \ \ 10^7 \leq p \leq 1.05\times 10^9\ \ \ \ p\in \{ prime\} \ \ \ \ T =70 \]

---

刚开始写了个 \(n^\frac{7}{8}\log n\) 的，直接暴毙/kk

\[\begin{aligned} &\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c}\left(\frac{\text{lcm}(i,j)}{\gcd(i,k)}\right)^{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} \]

对于下面：

\[\begin{aligned}&\prod_{i=1}^a\prod_{j=1}^b\prod_{k=1}^cgcd(i,j)\\=&\prod_{i=1}^a\prod_{j=1}^bgcd(i,j)^k\end{aligned} \]

直接改为枚举 \(gcd(i,j)\) ，就有：

\[\prod_{i=1} i^{kg(a,b,i)} \]

其中 \(g(a,b,x)\) 表示 \(\sum_{i=1}^a\sum_{j=1}^b[gcd(i,j)==x]\) ，直接对这个反演：

\[\begin{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}\lfloor\frac{a}{d}\rfloor\lfloor\frac{b}{d}\rfloor\mu(d)\end{aligned} \]

然后对于上式就有：

\[\prod_{i=1}i^{kg(\frac{a}{i},\frac{b}{i},1)} \]

这样就可以整除分块套整除分块了，单次 \(O(n^\frac{3}{4}\log n)\) 。

type = 2

对于上面：

\[\begin{aligned}&\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c}(ij)^{ijk}\\=&\prod_{i=1}^{a}\prod_{j=1}^{b}(ij)^{ij\sum k} \ (\sum k=\sum_{x=1}^{k}x)\\=&\prod_{i=1}^ai^{i\sum j\sum k}\prod_{j=1}j^{j\sum i\sum k}\end{aligned} \]

\(O(n\log n)\) 预处理 \(f_n=\prod_{i=1}^ni^i\) 就可以 \(O(n)\) 计算了。

对于下面：

\[\begin{aligned}&\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c}gcd(i,j)^{ijk}\\=&\prod_{i=1}^a\prod_{j=1}^bgcd(i,j)^{ij\sum k}\end{aligned} \]

又可以直接枚举 \(gcd\) 了。

\[\begin{aligned}\prod_{i=1}i^{\sum kg'(a,b,i)}\end{aligned} \]

\(g'(a,b,x)\) 表示 \(\sum_{i=1}^a\sum_{j=1}^bij[gcd(i,j)==x]\)

其实和上面的还是基本一样的套路，注意到 \(g'(a,b,x)= x^2g'(\frac{a}{x},\frac{b}{x},1)\) ，所以可以写成这样：

\[\begin{aligned}&x^2\sum_{d=1}d^2\sum_{i=1}^{\frac{a}{d}}\sum_{j=1}^{\frac{b}{d}}ij\mu(d)\\=&x^2\sum_{d=1}d^2\mu(d)\sum \frac{a}{d}\sum \frac{b}{d}\end{aligned} \]

然后就可以整除分块套整除分块了，预处理出来 \(\sum_{i}i^2\mu(d)\) 就可以单次 \(O(n^\frac{3}{4}\log n)\) 计算了。

type = 3

对于上面：

\[\begin{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}^b\sum_{k=1}^cgcd(i,j,k)}\prod_{j=1}^bj^{\sum_{i=1}^a\sum_{k=1}^cgcd(i,j,k)}\end{aligned} \]

发现左右两边是一样的，所以只讨论左边的值，直接对幂反演：

\[\begin{aligned}G(x)&=\sum_{j=1}^b\sum_{k=1}^cgcd(x,j,k)\\&=\sum_{j=1}^b\sum_{k=1}^c\sum_{d|x,d|j,d|k}\varphi(x)\\&=\sum_{d|x}\varphi(d)\lfloor\frac{b}{d}\rfloor\lfloor\frac{c}{d}\rfloor\end{aligned} \]

带回去原式（只有左半边）：

\[\begin{aligned}=&\prod_{i=1}^ai^{\sum_{d|x}\varphi(d)\lfloor\frac{b}{d}\rfloor\lfloor\frac{c}{d}\rfloor}\\=&\prod_{d=1}^a\prod_{i=1}^\frac{a}{d}(id)^{\varphi(d)\lfloor\frac{c}{d}\rfloor\lfloor\frac{b}{d}\rfloor}\\=&\prod_{d=1}^{a}(\lfloor\frac{a}{d}\rfloor!)^{\varphi(d)\lfloor\frac{c}{d}\rfloor\lfloor\frac{b}{d}\rfloor}\times d^{{\varphi(d)\lfloor\frac{b}{d}\rfloor\lfloor\frac{c}{d}\rfloor\lfloor\frac{b}{d}\rfloor}}\end{aligned} \]

这样子两边就可以分别整除分块了， \(O(n\log n)\) 预处理出来 \(n!\) 和\(g_n=\prod_{i=1}^ni^{\varphi(i)}\) 就可以单次 \(\sqrt{n}\log n\) 计算了。

然后对于下面：

\[\begin{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)\) ，就有：

\[\begin{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))^{d\sum_{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))^{d\mu(t)\lfloor\frac{c}{td}\rfloor}\end{aligned} \]

把 \(td\) 和 \(gcd(i,j)\) 分开考虑，对于 \(td\) 部分，就有：

\[\begin{aligned}&\prod_{d=1}\prod_{t=1}\prod_{i=1}^\frac{a}{td}\prod_{j=1}^\frac{b}{td}(td)^{d\mu(t)\lfloor\frac{c}{td}\rfloor}\\=&\prod_{T=1}\prod_{i=1}^{\frac{a}{T}}\prod_{i=1}^{\frac{b}{T}}T^{\sum_{t|T}t\mu(\frac{T}{t})\lfloor\frac{c}{T}\rfloor}\end{aligned} \]

注意到 \(\mu*id=\varphi\) ，所以就有：

\[\prod_{T=1}T^{\varphi(T){\lfloor\frac{a}{T}\rfloor}\lfloor\frac{b}{T}\rfloor} \]

根据之前预处理出来的东西可以单次 \(O(\sqrt{n}\log n)\) 计算。

然后再考虑 \(gcd(i,j)\) 部分，就有：

\[\begin{aligned}&\prod_{d=1}\prod_{t=1}\prod_{i=1}^\frac{a}{td}\prod_{j=1}^\frac{b}{td}gcd(i,j)^{d\mu(t)\lfloor\frac{c}{td}\rfloor}\\=&\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)d\mu(t)\lfloor\frac{c}{td}\rfloor}\\=&\prod_{d=1}\prod_{t=1}\prod_{h=1}\prod_{g=1}gh^{d\mu(h)\lfloor\frac{a}{tdgh}\rfloor\lfloor\frac{b}{tdgh}\rfloor\mu(t)\lfloor\frac{c}{td}\rfloor}\\=&\prod_{T=1(td)}(\prod_{t=1(gh)}(\prod_{d|t}d^{\mu(\frac{t}{d})})^{\lfloor\frac{a}{tT}\rfloor\lfloor\frac{b}{tT}\rfloor})^{\sum_{k|T}\mu(k)\frac{T}{k}\lfloor\frac{c}{T}\rfloor}\\=&\prod_{T=1}(\prod_{t=1}(\prod_{d|t}d^{\mu(\frac{t}{d})})^{\lfloor\frac{a}{tT}\rfloor\lfloor\frac{b}{tT}\rfloor})^{\varphi(T)\lfloor\frac{c}{T}\rfloor}\end{aligned} \]

设 \(G'_t=\prod_{d|t}d^\mu(\frac{t}{d})\) ，这个可以 \(O(n\log n)\) 时间预处理出来，然后式子就变成了：

\[\prod_{T=1}(\prod_{t=1}(G'_t)^{\lfloor\frac{a}{tT}\rfloor\lfloor\frac{b}{tT}\rfloor})^{\varphi(T)\lfloor\frac{c}{T}\rfloor} \]

这个东西就又是一个整除分块套整除分块了，可以单词 \(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;
}