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

链接

UOJ 221
BZOJ 4652
LOJ 2085

题解

xyx \perp y 表示 x,yx, y 互质,若 xy,xy\frac x y, x \perp y 是纯循环小数,则一定存在一个正整数 ll 使得

xklx (mod y)kl1 (mod y)\begin{aligned} x \cdot k ^ l &\equiv x \ (\text{mod } y)\\ k ^ l &\equiv 1 \ (\text{mod } y) \end{aligned}

根据欧拉定理,若 kyk \perp y,则 kφ(y)1 (mod y)k ^ {\varphi(y)} \equiv 1 \ (\text{mod } y),由于逆元是唯一的,所以存在 ll 的充要条件是 kyk \perp y,故问题转化为求

x=1ny=1m[xy][yk]=y=1m[yk]x=1n[xy]=y=1m[yk]x=1n[gcd(x,y)=1]\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}

后半部分是一个经典的莫比乌斯反演,我们可以得到

y=1m[yk]dynμ(d)nd\begin{aligned}\sum_{y = 1} ^ m[y \perp k]\sum_{d | y} ^ n\mu(d)\lfloor \frac n d \rfloor \end{aligned}

交换一下枚举顺序,得

d=1nμ(d)ndy=1m[d  y][yk]=d=1nμ(d)ndi=1md[idk]=d=1n[dk]μ(d)ndi=1md[ik]\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)=i=1n[ik]f(n, k) = \sum_{i = 1} ^ n[i \perp k]

由于

[ab]=[a mod bb][a \perp b] = [a \text{ mod }b \perp b]

f(n,k)=nkf(k,k)+f(n mod k,k)f(n, k) = \lfloor \frac n k \rfloor f(k, k) + f(n \text{ mod } k, k)

此时我们已经解决了原式的后半部分,现在考虑如何快速求出前半部分,令

g(n,k)=i=1n[ik]μ(i)g(n, k) = \sum_{i = 1} ^ n [i \perp k]\mu(i)

考虑 kk 的一个质因数 pp,那么 kk 可以表示为 pcqp ^ cq,由于 [1,n][1, n] 只有与 kk 互质的才是有效值,那么若 xkx \perp k,则 xpx \perp pxqx \perp q,我们可以考虑从 xqx \perp q 的取值中减去 xx 不与 pp 互质的部分,就可以得到 xkx \perp k 的部分,故

g(n,k)=i=1n[iq]μ(i)y=1np[pyq]μ(py)=g(n,q)y=1np[yq]μ(py)\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}

pyp \perp y 时,μ(py)=μ(p)μ(y)\mu(py) = \mu(p)\mu(y),否则 μ(py)=0\mu(py) = 0,则

g(n,k)=g(n,q)y=1np[yp][yq]μ(p)μ(y)=g(n,q)μ(p)y=1np[yk]μ(y)=g(n,q)+g(np,k)\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=0n = 0k=1,n=0k = 1, n = 0 时,直接返回 00 就可以了,k=1k = 1 时,就是莫比乌斯函数的前缀和,这时我们用杜教筛求就好了。
时间复杂度 O(ω(k)n+n23)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;
}
文章目录
  1. 1. 链接
  2. 2. 题解
  3. 3. 代码