「NOI 2016」循环之美-莫比乌斯反演 + 杜教筛

链接

UOJ 221
BZOJ 4652
LOJ 2085

题解

令 $x \perp y$ 表示 $x, y$ 互质,若 $\frac x y, x \perp y$ 是纯循环小数,则一定存在一个正整数 $l$ 使得
$$\begin{aligned} x \cdot k ^ l &\equiv x \ (\text{mod } y)\\ k ^ l &\equiv 1 \ (\text{mod } y) \end{aligned}$$
根据欧拉定理,若 $k \perp y$,则 $k ^ {\varphi(y)} \equiv 1 \ (\text{mod } y)$,由于逆元是唯一的,所以存在 $l$ 的充要条件是 $k \perp y$,故问题转化为求
$$\begin{aligned}&\sum_{x=1}^{n} \sum_{y=1}^{m}[x\perp y][y \perp k]\\=&\sum_{y=1}^{m}[y\perp k]\sum_{x=1}^{n}[x\perp y]\\ =& \sum_{y = 1} ^ m[y \perp k]\sum_{x = 1} ^ n[\gcd(x, y) = 1]\end{aligned}$$
后半部分是一个经典的莫比乌斯反演,我们可以得到
$$\begin{aligned}\sum_{y = 1} ^ m[y \perp k]\sum_{d | y} ^ n\mu(d)\lfloor \frac n d \rfloor \end{aligned}$$
交换一下枚举顺序,得
$$\begin{aligned}&\sum_{d=1}^{n}\mu(d)\lfloor \frac{n}{d} \rfloor \sum_{y=1}^{m}[d\ |\ y][y\perp k]\\=&\sum_{d=1}^{n}\mu(d)\lfloor \frac{n}{d}\rfloor \sum_{i=1}^{\lfloor \frac{m}{d} \rfloor}[id\perp k]\\=& \sum_{d=1}^{n}[d\perp k]\mu(d)\lfloor \frac{n}{d}\rfloor \sum_{i=1}^{\lfloor \frac{m}{d} \rfloor}[i\perp k]\end{aligned}$$

$$f(n, k) = \sum_{i = 1} ^ n[i \perp k]$$
由于
$$[a \perp b] = [a \text{ mod }b \perp b]$$

$$f(n, k) = \lfloor \frac n k \rfloor f(k, k) + f(n \text{ mod } k, k)$$
此时我们已经解决了原式的后半部分,现在考虑如何快速求出前半部分,令
$$g(n, k) = \sum_{i = 1} ^ n [i \perp k]\mu(i)$$
考虑 $k$ 的一个质因数 $p$,那么 $k$ 可以表示为 $p ^ cq$,由于 $[1, n]$ 只有与 $k$ 互质的才是有效值,那么若 $x \perp k$,则 $x \perp p$ 且 $x \perp q$,我们可以考虑从 $x \perp q$ 的取值中减去 $x$ 不与 $p$ 互质的部分,就可以得到 $x \perp k$ 的部分,故
$$\begin{aligned}g(n,k)&=\sum_{i=1}^n[i\perp q]\mu(i)-\sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[py\perp q]\mu(py) \\&=g(n,q)- \sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[y\perp q]\mu(py)\end{aligned}$$
当 $p \perp y$ 时,$\mu(py) = \mu(p)\mu(y)$,否则 $\mu(py) = 0$,则
$$\begin{aligned}g(n,k)&=g(n,q)- \sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[y\perp p][y\perp q]\mu(p)\mu(y)\\&=g(n,q)-\mu(p)\sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[y\perp k]\mu(y)\\&=g(n,q)+g(\lfloor\frac{n}{p}\rfloor,k) \end{aligned}$$
此时我们就可以递归求解了,考虑边界,当 $n = 0$ 或 $k = 1, n = 0$ 时,直接返回 $0$ 就可以了,$k = 1$ 时,就是莫比乌斯函数的前缀和,这时我们用杜教筛求就好了。
时间复杂度 $O(\omega(k)\sqrt{n} + n ^ {\frac 2 3})$

代码

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
/**
* Copyright (c) 2017, xehoth
* All rights reserved.
* 「NOI 2016」循环之美 22-07-2017
* 莫比乌斯反演 + 杜教筛
* @author xehoth
*/
#include <bits/stdc++.h>
#include <tr1/unordered_map>
namespace Task {
#define long long long
const int MAXN = 300010;
const int MAXK = 2010;
int prime[MAXN], cnt, n, m, k;
long mu[MAXN];
bool vis[MAXN];
inline void fastLinearSieveMethod(const int n) {
mu[1] = 1;
for (register int i = 2; i <= n; i++) {
if (!vis[i]) prime[cnt++] = i, mu[i] = -1;
for (register int j = 0, tmp; j < cnt && (tmp = i * prime[j]) <= n;
j++) {
vis[tmp] = true;
if (i % prime[j] == 0) {
mu[tmp] = 0;
break;
} else {
mu[tmp] = -mu[i];
}
}
}
for (register int i = 1; i <= n; i++) mu[i] += mu[i - 1];
}
typedef std::tr1::unordered_map<int, long> HashMap;
HashMap map[MAXK];
HashMap::iterator it;
const int sieveBlockSize = 300000;
HashMap buc;
inline long sieveMuMain(int x) {
if (x <= sieveBlockSize) return mu[x];
register long &cur = buc[x];
if (cur) return cur;
register long ret = 1;
for (register int i = 2, pos; i <= x; i = pos + 1)
pos = x / (x / i), ret -= (pos - i + 1) * sieveMuMain(x / i);
return cur = ret;
}
int f[MAXK], fac[MAXK], facCnt;
inline void initFunctionF() {
for (register int i = 1; i <= k; i++)
f[i] = f[i - 1] + (std::__gcd(i, k) == 1);
}
inline long getF(int x) { return (long)(x / k) * f[k] + f[x % k]; }
inline long getG(int x, int y) {
if (x <= 1) return x;
if (y == 0) return sieveMuMain(x);
register long &cur = map[y][x];
return cur ? cur : cur = getG(x, y - 1) + getG(x / fac[y], y);
}
inline void initFactor() {
register int now = k;
for (register int i = 2; i <= k; i++) {
if (now % i != 0) continue;
fac[++facCnt] = i;
while (now % i == 0) now /= i;
if (now == 1) break;
}
}
inline void solve() {
std::ios::sync_with_stdio(false);
std::cin.tie(NULL), std::cout.tie(NULL);
std::cin >> n >> m >> k;
fastLinearSieveMethod(sieveBlockSize);
initFunctionF(), initFactor();
register int min = std::min(n, m);
register long ans = 0;
for (register int i = 1, pos; i <= min; i = pos + 1) {
pos = std::min(n / (n / i), m / (m / i));
ans += (long)getF(m / i) * (n / i) *
(getG(pos, facCnt) - getG(i - 1, facCnt));
}
std::cout << ans;
}
}
int main() {
Task::solve();
return 0;
}
分享到