Dyd's Blog

He who has a strong enough why can bear almost any how.

luoguP5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题

Posted on Edited on
Symbols count in article: 8k Reading time ≈ 7 mins.

神仙 东方 莫反题

幽灵乐团 / 莫比乌斯反演基础练习题

以下所有除法/分数默认下取整，为方便记模数为 $P$

一眼望过去分讨题，然后发现每种都好难……
$$
\begin{aligned}
Ans
= & \prod_{i = 1}^{A} \prod_{j = 1}^{B} \prod_{k = 1}^{C} (\frac{lcm(i, j)}{\gcd(i, k)})^{f(type)} \\
= & \prod_{i = 1}^{A} \prod_{j = 1}^{B} \prod_{k = 1}^{C} (\frac{ij}{\gcd(i, k) \gcd(i, j)})^{f(type)} \\
\end{aligned}
$$
设：
$$
\begin{aligned}
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)}
\end{aligned}
$$
原式为
$$
Ans = \frac{g(A, B, C)g(B, A, C)}{h(A, B, C) h(A, C, B)}
$$
然后分讨

part 1

$$
\begin{aligned}
g(a, b, c)
= & \prod_{i = 1}^{a} \prod_{j = 1}^{b} \prod_{k = 1}^{c} i \\
= & \prod_{i = 1}^A i^{bc} \\
= & (a!)^{bc}
\end{aligned}
$$

$O(n)$ 预处理阶乘，快速幂可以单次 $O(\log P)$ 求出 $g$
$$
\begin{aligned}
h(a, b, c)
= & \prod_{i = 1}^{a} \prod_{j = 1}^{b} \prod_{k = 1}^{c} \gcd(i, j) \\
= & \prod_{i = 1}^{a} \prod_{j = 1}^{b} \gcd(i, j)^{c} \\
\end{aligned}
$$
看到连乘尝试取对数得：
$$
\begin{aligned}
\ln(h(a, b, c))
= & c \sum_{i = 1}^a \sum_{j = 1}^b \ln(\gcd(i, j)) \\
= & c \sum_{d = 1}^{t = \min(a, b)} \ln(d) \sum_{i = 1}^{\frac{a}{d}} \sum_{j = 1}^\frac{b}{d} [\gcd(i, j) == 1] \\
= & c \sum_{d = 1}^{t} \ln(d) \sum_{i = 1}^{\frac{a}{d}} \sum_{j = 1}^\frac{b}{d} \sum_{x | i \wedge x|j} \mu(x) \\
= & c \sum_{d = 1}^{t} \ln(d) \sum_{x}^{\frac{t}{d}} \sum_{i = 1}^{\frac{a}{dx}} \sum_{j = 1}^{\frac{b}{dx}} \mu(x) \\
= & c \sum_{d = 1}^{t} \ln(d) \sum_{x}^{\frac{t}{d}} \mu(x) \frac{a}{dx} \frac{b}{dx}\\
= & c \sum_{d = 1}^{t} \sum_{x}^{\frac{t}{d}} \ln(d) \mu(x) \frac{a}{dx} \frac{b}{dx}\\
= & c \sum_{p = 1}^{t} \sum_{d | p} \ln(d) \mu(\frac{p}{d}) \frac{a}{p} \frac{b}{p}\\
= & c \sum_{p = 1}^{t} \frac{a}{p} \frac{b}{p} \sum_{d | p} \ln(d) \mu(\frac{p}{d}) \\
\end{aligned}
$$
再取 $\exp$ ：
$$
\begin{aligned}
h(a, b, c)
= & (\prod_{p = 1}^{t} (\prod_{d | p} d^{\mu(\frac{p}{d})})^{\frac{a}{p} \frac{b}{p}})^c
\end{aligned}
$$
这个跟 shit 一样的东西竟然是可以算的，考虑到 $\mu$ 只有三个取值，只要 $O(n(\log P + \log n))$ 预处理中间的东西，再整除分块可以做到 $O(\sqrt{n} \log P)$

part 2

记 $S(x) = \sum_{i = 1}^x i = \frac{x(x + 1)}{2}$
$$
\begin{aligned}
g(a, b, c)
= & \prod_{i = 1}^{a} \prod_{j = 1}^{b} \prod_{k = 1}^{c} i^{ijk} \\
= & \prod_{i = 1}^{a} (i^i)^{S(b) S(c)} \\
= & (\prod_{i = 1}^{a} i^i)^{S(b) S(c)} \\
\end{aligned}
$$
可以 $O(n \log P)$ 预处理， $O(\log P)$ 求
$$
\begin{aligned}
h(a, b, c)
= & \prod_{i = 1}^{a} \prod_{j = 1}^{b} \prod_{k = 1}^{c} \gcd(i, j)^{ijk} \\
= & \prod_{i = 1}^{a} \prod_{j = 1}^{b} (\gcd(i, j)^{ij})^{S(c)} \\
\end{aligned}
$$
还是取对数：
$$
\begin{aligned}
\ln(h(a, b, c))
= & S(c) \sum_{i = 1}^a \sum_{j = 1}^b ij \ln(\gcd(i, j)) \\
= & S(c) \sum_{d = 1}^{t = \min(a, b)} \sum_{i = 1}^{\frac{a}{d}} \sum_{j = 1}^{\frac{b}{d}} [\gcd(i, j) == 1] i j d^2 \ln(d) \\
= & S(c) \sum_{d = 1}^{t} d^2 \ln(d) \sum_{i = 1}^{\frac{a}{d}} \sum_{j = 1}^{\frac{b}{d}} i j \sum_{x | i \wedge x | j} \mu(x) \\
= & S(c) \sum_{d = 1}^{t} d^2 \ln(d) \sum_{x = 1}^{\frac{t}{d}} \mu(x) \sum_{i = 1}^{\frac{a}{dx}} \sum_{j = 1}^{\frac{b}{dx}} i j x^2\\
= & S(c) \sum_{d = 1}^{t} d^2 \ln(d) \sum_{x = 1}^{\frac{t}{d}} x^2 \mu(x) \sum_{i = 1}^{\frac{a}{dx}} \sum_{j = 1}^{\frac{b}{dx}} i j \\
= & S(c) \sum_{d = 1}^{t} \sum_{x = 1}^{\frac{t}{d}} d^2 \ln(d) x^2 \mu(x) S(\frac{a}{dx}) S(\frac{b}{dx}) \\
= & S(c) \sum_{d = 1}^{t} \sum_{x = 1}^{\frac{t}{d}} (dx)^2 \ln(d) \mu(x) S(\frac{a}{dx}) S(\frac{b}{dx}) \\
= & S(c) \sum_{p = 1}^{t} \sum_{d | p} p^2 \ln(d) \mu(\frac{p}{d}) S(\frac{a}{p}) S(\frac{b}{p}) \\
= & S(c) \sum_{p = 1}^{t} S(\frac{a}{p}) S(\frac{b}{p}) \sum_{d | p} \ln(d) \mu(\frac{p}{d}) p^2 \\
\end{aligned}
$$
然后 $\exp$
$$
\begin{aligned}
h(a, b, c)
= & (\prod_{p = 1}^t (\prod_{d | p} d^{\mu(\frac{p}{d})p^2} )^{S(\frac{a}{p}) S(\frac{b}{p})})^{S(c)}
\end{aligned}
$$
最里面的部分和 part 1 类似，但由于有 $p^2$ 导致复杂度变为 $O(n \log n \log P)$ ，每次询问还是整除分块 $O(\sqrt{n} \log P)$

part 3

此时 $g, h$ 都不太好算，我们直接计算原式：
$$
\begin{aligned}
& \prod_{i = 1}^{A} \prod_{j = 1}^{B} \prod_{k = 1}^{C} (\frac{lcm(i, j)}{\gcd(i, k)})^{\gcd(i, j, k)} \\
\end{aligned}
$$
取对数：
$$
\begin{aligned}
& \sum_{i = 1}^{A} \sum_{j = 1}^{B} \sum_{k = 1}^{C} \gcd(i, j, k) \ln(\frac{lcm(i, j)}{\gcd(i, k)}) \\
= & \sum_{d = 1}^{t = \min(a, b, c)} \sum_{i = 1}^{\frac{A}{d}} \sum_{j = 1}^{\frac{B}{d}} \sum_{k = 1}^{\frac{C}{d}} [\gcd(i, j, k) == 1] d \ln(\frac{lcm(id, jd)}{\gcd(id, kd)}) \\
= & \sum_{d = 1}^{t} \sum_{i = 1}^{\frac{A}{d}} \sum_{j = 1}^{\frac{B}{d}} \sum_{k = 1}^{\frac{C}{d}} d \ln(\frac{lcm(i, j)}{\gcd(i, k)}) \sum_{x | i \wedge x | j \wedge x | k} \mu(x) \\
= & \sum_{d = 1}^{t} \sum_{x = 1}^{\frac{t}{d}} \sum_{i = 1}^{\frac{A}{dx}} \sum_{j = 1}^{\frac{B}{dx}} \sum_{k = 1}^{\frac{C}{dx}} d \mu(x) \ln(\frac{lcm(i, j)}{\gcd(i, k)}) \\
= & \sum_{p = 1}^{t} \sum_{d | p} \sum_{i = 1}^{\frac{A}{p}} \sum_{j = 1}^{\frac{B}{p}} \sum_{k = 1}^{\frac{C}{p}} d \mu(\frac{p}{d}) \ln(\frac{lcm(i, j)}{\gcd(i, k)}) \\
= & \sum_{p = 1}^{t} \sum_{d | p} d \mu(\frac{p}{d}) \sum_{i = 1}^{\frac{A}{p}} \sum_{j = 1}^{\frac{B}{p}} \sum_{k = 1}^{\frac{C}{p}} \ln(\frac{lcm(i, j)}{\gcd(i, k)}) \\
= & \sum_{p = 1}^{t} \varphi(p) \sum_{i = 1}^{\frac{A}{p}} \sum_{j = 1}^{\frac{B}{p}} \sum_{k = 1}^{\frac{C}{p}} \ln(\frac{lcm(i, j)}{\gcd(i, k)}) \\
\end{aligned}
$$
然后 $exp$
$$
\prod_{p = 1}^t (\prod_{i = 1}^{\frac{A}{p}} \prod_{j = 1}^{\frac{B}{p}} \prod_{k = 1}^{\frac{C}{p}} \frac{lcm(i, j)}{\gcd(i, k)})^{\varphi(p)}
$$
发现里面的东西就是 part 1，于是预处理出 $\varphi$ 的前缀和，整除分块是 $O(\sqrt{n} \log P)$ 的，里面又是个 $O(\sqrt{n} \log P)$ ，总复杂度 $O(n \log^2 P)$

代码

注意指数要模的是 $P - 1$ ，筛 $\mu$ 的时候 $-1$ 不要存为 $P - 1$ ，因为 $\mu$ 有可能在指数上，还要预处理的 $\varphi$ 也是模 $P - 1$ 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
#include <cstdio>
#include <algorithm>
using LL = long long;
const int N = 1e5 + 100;
int P;
int fac[N], mu[N], dm[N], idm[N], tt[N], ii[N], dmp[N], idmp[N], phi[N];
int ct = 0, pri[N];
void adj(int &x){ x += (x >> 31) & P; }
int qpow(int x, int y = P - 2)
{
	int res = 1;
	for (; y; y >>= 1, x = LL(x) * x % P) if (y & 1) res = LL(res) * x % P;
	return res;
}
void init()
{
	fac[0] = fac[1] = ii[0] = ii[1] = 1;
	for (int i = 2; i < N; ++i)
	{
		fac[i] = LL(i) * fac[i - 1] % P;
		ii[i] = qpow(i, i);
		ii[i] = LL(ii[i]) * ii[i - 1] % P;
	}
	mu[1] = phi[1] = 1;
	for (int i = 2; i < N; ++i)
	{
		if (!pri[i]) pri[++ct] = i, mu[i] = -1, phi[i] = i - 1;
		for (int j = 1; j <= ct && LL(i) * pri[j] < N; ++j)
		{
			pri[i * pri[j]] = -1;
			if (i % pri[j] == 0)
			{
				phi[i * pri[j]] = pri[j] * phi[i];
				mu[i * pri[j]] = 0;
				break;
			}
			mu[i * pri[j]] = -mu[i];
			phi[i * pri[j]] = (pri[j] - 1) * phi[i]; 
		}
	}
	for (int i = 2; i < N; ++i) phi[i] += phi[i - 1], phi[i] %= P - 1;
	for (int i = 0; i < N; ++i) dm[i] = idm[i] = dmp[i] = idmp[i] = 1;
	for (int i = 2, t, p; i < N; ++i)
	{
		t = qpow(i, P - 2);
		for (int j = 1; LL(i) * j < N; ++j)
		{
			p = i * j;
			if (mu[j] == -1)
			{
				dm[p] = LL(t) * dm[p] % P;
				dmp[p] = LL(qpow(t, LL(p) * p % (P - 1))) * dmp[p] % P;
			}
			else if (mu[j] == 1)
			{
				dm[p] = LL(i) * dm[p] % P;
				dmp[p] = LL(qpow(i, LL(p) * p % (P - 1))) * dmp[p] % P;
			}
		}	
	}
	for (int i = 3; i < N; ++i)
	{
		tt[i] = dm[i];
		dm[i] = LL(dm[i]) * dm[i - 1] % P;
	}
	idm[N - 1] = qpow(dm[N - 1]);
	for (int i = N - 2; i > 1; --i) idm[i] = LL(tt[i + 1]) * idm[i + 1] % P;
	for (int i = 3; i < N; ++i)
	{
		tt[i] = dmp[i];
		dmp[i] = LL(dmp[i]) * dmp[i - 1] % P;
	}
	idmp[N - 1] = qpow(dmp[N - 1]);
	for (int i = N - 2; i > 1; --i) idmp[i] = LL(tt[i + 1]) * idmp[i + 1] % P;
}
int g1(int a, int b, int c){ return qpow(fac[a], LL(b) * c % (P - 1)); }
int h1(int a, int b, int c)
{
	int res = 1, t = std::min(a, b), x;
	for (int l = 2, r; l <= t; l = r + 1)
	{
		r = std::min(a / (a / l), b / (b / l));
		x = LL(dm[r]) * idm[l - 1] % P;
		res = LL(res) * qpow(x, LL(a / l) * (b / l) % (P - 1)) % P;
	}
	return qpow(res, c);
}
int calc1(int a, int b, int c)
{
	int t1 = LL(g1(a, b, c)) * g1(b, a, c) % P, t2 = LL(h1(a, b, c)) * h1(a, c, b) % P;
	t2 = qpow(t2);
	return LL(t1) * t2 % P;
}
int S(int x){ return LL(x + 1) * x / 2 % (P - 1); }
int g2(int a, int b, int c){ return qpow(ii[a], LL(S(b)) * S(c) % (P - 1)); }
int h2(int a, int b, int c)
{
	int res = 1, t = std::min(a, b), x;
	for (int l = 2, r; l <= t; l = r + 1)
	{
		r = std::min(a / (a / l), b / (b / l));
		x = LL(dmp[r]) * idmp[l - 1] % P;
		res = LL(res) * qpow(x, LL(S(a / l)) * S(b / l) % (P - 1)) % P;
	}
	return qpow(res, S(c));
}
int calc2(int a, int b, int c)
{
	int t1 = LL(g2(a, b, c)) * g2(b, a, c) % P, t2 = LL(h2(a, b, c)) * h2(a, c, b) % P;
	t2 = qpow(t2);
	return LL(t1) * t2 % P; 
}
int calc3(int a, int b, int c)
{
	int res = 1, t = std::min(a, std::min(b, c)), x, y;
	for (int l = 1, r; l <= t; l = r + 1)
	{
		r = std::min(a / (a / l), std::min(b / (b / l), c / (c / l)));
		x = calc1(a / l, b / l, c / l);
		y = (phi[r] - phi[l - 1] + P - 1) % (P - 1);
		res = LL(res) * qpow(x, y) % P;
	}
	return res;
}
int main()
{
	int T, A, B, C;
	scanf("%d %d", &T, &P);
	init();
	while (T--)
	{
		scanf("%d %d %d", &A, &B, &C);
		printf("%d %d %d\n", calc1(A, B, C), calc2(A, B, C), calc3(A, B, C));
	}
	return 0;
}