洛谷 P5518 [MtOI2019]幽灵乐团

东风谷 早苗(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;
}