「题解」Luogu P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题
P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题
Description
-
多测。
-
在一开始给定模数 \(p\)。
-
每次给定 \(3\) 个正整数 \(a, b, c\)，求
\[\prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \left(\dfrac{\operatorname{lcm}(i, j)}{\gcd(i, k)} \right)^{f(type)} \bmod p \]其中 \(type\in \{0, 1, 2\}\)：
\[f(0) = 1 \\ f(1) = ijk \\ f(2) = \gcd(i, j, k) \]每次分别回答 \(type = 0/ 1/ 2\) 时的答案。
-
\(1\le a, b, c\le 10^5, 10^7 \le p \le 1.05\times 10^9, p\in \mathbb{P}, T = 70\)。
Solution
以下令 \(S(n) = \sum_{i = 1}^n i\)。
参考了 peterwuyihong 大佬提供的思路，在 \(\prod\) 较多较复杂时，可以对原式去 \(\ln\)，将其变成 \(\sum\) 的形式，处理完之后再使用 \(\exp\) 转回 \(\prod\) 的形式。
\[ans = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \left(\dfrac{ij}{\gcd(i, j)\gcd(i, k)} \right)^{f(type)}
\]
发现每一项只与 \(2\) 个 \(\prod\) 有关，于是设：
\[g(a, b, c) = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c i^{f(type)} \\
h(a, b, c) = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \gcd(i, j)^{f(type)}
\]
则
\[ans = \dfrac{g(a, b, c) g(b, a, c)}{h(a, b, c) h(a, c, b)}
\]
分类讨论。
type = 0
\[\begin{aligned}
g(a, b, c)
& = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c i^1 \\
& = (a!)^{bc}
\end{aligned}
\]
预处理阶乘做到单次 \(\Omicron(\log p)\)。
\[h(a, b, c) = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \gcd(i, j)^1
\]
不妨设 \(a\le b\)。
取 \(\ln\) 得
\[\begin{aligned}
\sum_{i = 1}^a \sum_{j = 1}^b \sum_{k = 1}^c \ln(\gcd(i, j))
& = c \sum_{d = 1}^a \ln(d) \sum_{i = 1}^a \sum_{j = 1}^b [\gcd(i, j) = d] \\
& = c \sum_{d = 1}^a \ln(d) \sum_{k = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \mu(k) \left\lfloor\dfrac{a}{dk}\right\rfloor \left\lfloor\dfrac{b}{dk}\right\rfloor \\
& = c \sum_{T = 1}^a \left\lfloor\dfrac{a}{T}\right\rfloor \left\lfloor\dfrac{b}{T}\right\rfloor \sum_{d\mid T} \ln(d) \mu\left(\dfrac{T}{d}\right)
\end{aligned}
\]
取 \(\exp\) 得
\[\left[\prod_{T = 1}^a \left(\prod_{d \mid T} d^{\mu\left(\frac{T}{d}\right)} \right)^{\left\lfloor\frac{a}{T}\right\rfloor \left\lfloor\frac{b}{T}\right\rfloor} \right]^c
\]
中间那个可以 \(\Omicron(n\ln n\log p)\) 预处理，结合整除分块即可 \(\Omicron(\sqrt{n}\log p)\)。
type = 1
\[\begin{aligned}
g(a, b, c)
& = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c i^{ijk} \\
& = \left(\prod_{i = 1}^a i^i \right)^{S(b) S(c)}
\end{aligned}
\]
\(\Omicron(n\log p)\) 预处理，\(\Omicron(\log p)\) 回答。
\[h(a, b, c) = \prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \gcd(i, j)^{ijk}
\]
不妨设 \(a\le b\)。
取 \(\ln\) 得
\[\begin{aligned}
\sum_{i = 1}^a \sum_{j = 1}^b \sum_{k = 1}^c \ln(\gcd(i, j)) ijk
& = \sum_{d = 1}^a \ln(d) \sum_{i = 1}^a i \sum_{j = 1}^b j \sum_{k = 1}^c k [\gcd(i, j) = d] \\
& = S(c) \sum_{d = 1}^a \ln(d) \sum_{i = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} id \sum_{j = 1}^{\left\lfloor\frac{b}{d}\right\rfloor} jd [\gcd(i, j) = 1] \\
& = S(c) \sum_{d = 1}^a \ln(d) \sum_{k = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \mu(k) \sum_{i = 1}^{\left\lfloor\frac{a}{dk}\right\rfloor} idk \sum_{j = 1}^{\left\lfloor\frac{b}{dk}\right\rfloor} jdk \\
& = S(c) \sum_{d = 1}^a \ln(d) d^2 \sum_{k = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \mu(k) k^2 S\left(\left\lfloor\dfrac{a}{dk}\right\rfloor \right) S\left(\left\lfloor\dfrac{b}{dk}\right\rfloor \right) \\
& = S(c) \sum_{T = 1}^a S\left(\left\lfloor\dfrac{a}{T}\right\rfloor \right) S\left(\left\lfloor\dfrac{b}{T}\right\rfloor \right) \sum_{d\mid T} \ln(d) d^2 \mu\left(\dfrac{T}{d}\right) \left(\dfrac{T}{d}\right)^2 \\
& = S(c) \sum_{T = 1}^a T^2 S\left(\left\lfloor\dfrac{a}{T}\right\rfloor\right) S\left(\left\lfloor\dfrac{b}{T}\right\rfloor\right) \sum_{d\mid T} \ln(d) \mu\left(\dfrac{T}{d}\right)
\end{aligned}
\]
取 \(\exp\) 得
\[\left[\prod_{T = 1}^a \left(\prod_{d\mid T} d^{\mu\left(\frac{T}{d}\right) T^2} \right)^{S\left(\left\lfloor\frac{a}{T}\right\rfloor \right) S\left(\left\lfloor\frac{b}{T}\right\rfloor \right)} \right]^{S(c)}
\]
照样是 \(\Omicron(n\ln n\log p)\) 的预处理和 \(\Omicron(\sqrt{n} \log p)\) 的回答。
type = 2
直接计算原式（注意原式不可调换 \(a, b\)）
\[\prod_{i = 1}^a \prod_{j = 1}^b \prod_{k = 1}^c \left(\dfrac{\operatorname{lcm}(i, j)}{\gcd(i, k)} \right)^{\gcd(i, j, k)}
\]
取 \(\ln\) 得
\[\begin{aligned}
\sum_{i = 1}^a \sum_{j = 1}^b \sum_{k = 1}^c \ln\left(\dfrac{\operatorname{lcm}(i, j)}{\gcd(i, k)} \right) \gcd(i, j, k)
& = \sum_{d = 1}^{\min(a, b, c)} \varphi(d) \sum_{i = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \sum_{j = 1}^{\left\lfloor\frac{b}{d}\right\rfloor} \sum_{k = 1}^{\left\lfloor\frac{c}{d}\right\rfloor} \ln\left(\dfrac{\operatorname{lcm}(id, jd)}{\gcd(id, kd)} \right) \\
& = \sum_{d = 1}^{\min(a, b, c)} \varphi(d) \sum_{i = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \sum_{j = 1}^{\left\lfloor\frac{b}{d}\right\rfloor} \sum_{k = 1}^{\left\lfloor\frac{c}{d}\right\rfloor} \ln\left(\dfrac{\operatorname{lcm}(i, j)}{\gcd(i, k)} \right)
\end{aligned}
\]
取 \(\exp\) 得E
\[\prod_{d = 1}^{\min(a, b, c)} \left(\prod_{i = 1}^{\left\lfloor\frac{a}{d}\right\rfloor} \prod_{j = 1}^{\left\lfloor\frac{b}{d}\right\rfloor} \prod_{k = 1}^{\left\lfloor\frac{c}{d}\right\rfloor} \dfrac{\operatorname{lcm}(i, j)}{\gcd(i, k)} \right)^{\varphi(d)}
\]
发现没？中间那一块就是 \(type = 0\) 时的答案！
于是整除分块套整除分块实现 \(\Omicron(\sqrt{n} \log p \cdot \sqrt{n}\log p) = \Omicron(n\log^2 p)\)。
Code
tips : 注意哪些是指数，要 \(\bmod p - 1\)。
// 18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#include <cstring>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;
namespace IO
{
}
using IO::read;
int p, q;
int add1(int a, int b) {return (a + b) % p;}
int add2(int a, int b) {return (a + b) % q;}
int sub1(int a, int b) {return (a - b + p) % p;}
int sub2(int a, int b) {return (a - b + q) % q;}
int mul1(int a, int b) {return (ll)a * b % p;}
int mul2(int a, int b) {return (ll)a * b % q;}
int qpow(int a, int b) {if (b < 0) b += q; int base = a, ans = 1; while (b) {if (b & 1) ans = mul1(ans, base); base = mul1(base, base); b >>= 1;} return ans;}
int inv(int a) {return qpow(a, p - 2);}
const int MAXN = 1e5 + 5;
typedef int arr[MAXN];
const int N = 1e5;
arr prime, mu, phi, phiS, fac, dyy, dyyP, dyyPP, dyyPPI, dyyPI, indexP, xhj, xhjP, xhjPP, xhjPPI, xhjPI;
bool vis[MAXN];
void pre()
{
mu[1] = phi[1] = 1;
for (int i = 2; i <= N; i++)
{
if (!vis[i])
{
prime[++prime[0]] = i;
mu[i] = -1;
phi[i] = i - 1;
}
for (int j = 1; j <= prime[0] && i * prime[j] <= N; j++)
{
vis[i * prime[j]] = true;
if (i % prime[j] == 0)
{
mu[i * prime[j]] = 0;
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
mu[i * prime[j]] = mu[i] * mu[prime[j]];
phi[i * prime[j]] = phi[i] * phi[prime[j]];
}
}
for (int i = 1; i <= N; i++)
{
phiS[i] = add2(phiS[i - 1], phi[i]);
}
fac[0] = 1;
for (int i = 1; i <= N; i++)
{
fac[i] = mul1(fac[i - 1], i);
}
for (int i = 1; i <= N; i++)
{
dyy[i] = 1;
}
for (int i = 1; i <= N; i++)
{
for (int j = 1; i * j <= N; j++)
{
dyy[i * j] = mul1(dyy[i * j], qpow(i, mu[j]));
}
}
dyyP[0] = dyyPP[0] = 1;
for (int i = 1; i <= N; i++)
{
dyyP[i] = mul1(dyyP[i - 1], dyy[i]);
dyyPP[i] = mul1(dyyPP[i - 1], dyyP[i]);
}
dyyPPI[N] = inv(dyyPP[N]);
dyyPI[0] = 1;
for (int i = N - 1; i >= 0; i--)
{
dyyPPI[i] = mul1(dyyPPI[i + 1], dyyP[i + 1]);
dyyPI[i + 1] = mul1(dyyPPI[i + 1], dyyPP[i]);
}
indexP[0] = 1;
for (int i = 1; i <= N; i++)
{
indexP[i] = mul1(indexP[i - 1], qpow(i, i));
}
for (int i = 1; i <= N; i++)
{
xhj[i] = 1;
}
for (int i = 1; i <= N; i++)
{
for (int j = 1; i * j <= N; j++)
{
xhj[i * j] = mul1(xhj[i * j], qpow(i, mul2(mu[j], mul2(i * j, i * j))));
}
}
xhjP[0] = xhjPP[0] = 1;
for (int i = 1; i <= N; i++)
{
xhjP[i] = mul1(xhjP[i - 1], xhj[i]);
xhjPP[i] = mul1(xhjPP[i - 1], xhjP[i]);
}
xhjPPI[N] = inv(xhjPP[N]);
xhjPI[0] = 1;
for (int i = N - 1; i >= 0; i--)
{
xhjPPI[i] = mul1(xhjPPI[i + 1], xhjP[i + 1]);
xhjPI[i + 1] = mul1(xhjPPI[i + 1], xhjPP[i]);
}
}
int S(int n)
{
return (ll)n * (n + 1) / 2 % q;
}
int GetSum_phi(int l, int r)
{
return sub2(phiS[r], phiS[l - 1]);
}
int GetPro_dyy(int l, int r)
{
return mul1(dyyP[r], dyyPI[l - 1]);
}
int GetPro_xhj(int l, int r)
{
return mul1(xhjP[r], xhjPI[l - 1]);
}
struct Type0
{
int g(int a, int b, int c)
{
return qpow(fac[a], mul2(b, c));
}
int h(int a, int b, int c)
{
if (a > b)
{
swap(a, b);
}
int res = 1;
for (int l = 1, r; l <= a; l = r + 1)
{
int k1 = a / l, k2 = b / l;
r = min(a / k1, b / k2);
res = mul1(res, qpow(GetPro_dyy(l, r), mul2(k1, k2)));
}
return qpow(res, c);
}
int GetAns(int a, int b, int c)
{
int x = mul1(g(a, b, c), g(b, a, c)), y = mul1(h(a, b, c), h(a, c, b));
return mul1(x, inv(y));
}
}type0;
struct Type1
{
int g(int a, int b, int c)
{
return qpow(indexP[a], mul2(S(b), S(c)));
}
int h(int a, int b, int c)
{
if (a > b)
{
swap(a, b);
}
int res = 1;
for (int l = 1, r; l <= a; l = r + 1)
{
int k1 = a / l, k2 = b / l;
r = min(a / k1, b / k2);
res = mul1(res, qpow(GetPro_xhj(l, r), mul2(S(k1), S(k2))));
}
return qpow(res, S(c));
}
int GetAns(int a, int b, int c)
{
int x = mul1(g(a, b, c), g(b, a, c)), y = mul1(h(a, b, c), h(a, c, b));
return mul1(x, inv(y));
}
}type1;
struct Type2
{
int GetAns(int a, int b, int c)
{
int res = 1, upbound = min(a, min(b, c));
for (int l = 1, r; l <= upbound; l = r + 1)
{
int k1 = a / l, k2 = b / l, k3 = c / l;
r = min(a / k1, min(b / k2, c / k3));
res = mul1(res, qpow(type0.GetAns(k1, k2, k3), GetSum_phi(l, r)));
}
return res;
}
}type2;
signed main()
{
int t = read();
p = read(), q = p - 1;
pre();
while (t--)
{
int a = read(), b = read(), c = read();
printf("%d %d %d\n", type0.GetAns(a, b, c), type1.GetAns(a, b, c), type2.GetAns(a, b, c));
}
return 0;
}
