「NOI 2016」循环之美

Problem

Description

给定三个十进制数 n, m, k ,求在 k 进制下,有多少个数值互不相等的纯循环小数可以用分数 \displaystyle\frac{x}{y} \left(x, y \in \mathbb{N^{+}}, x \le n, y \le m\right) 表示。
一个数是纯循环的,当且仅当其可以写成以下形式:

a.\dot{c_1} c_2 c_3 \ldots c_{p - 1} \dot{c_p}

其中, a 是一个整数, p \ge 1 ;对于 1 \le i \le p c_i k 进制下的一位数字。

需要特别注意的是,我们认为一个整数是纯循环的,因为它的小数部分可以表示成 0 的循环或是 k-1 的循环。

Constraints

1 \le n, m \le 10^{9}, 2 \le k \le 2000

Solution

Analysis

对于 x \perp y ,观察 \displaystyle\frac{x}{y} 进制转换的过程,若其为纯循环小数,设其循环节长度为 l ,则有:

\begin{aligned} z &= \frac{x}{y} k^{l} - \frac{x}{y} \in \mathbb{Z} \\ xk^{l} - x &= yz \\ xk^{l} &\equiv x \pmod{y} \\ k^{l} &\equiv 1 \pmod{y} \end{aligned}

那么 k^{l - 1} 即为 k 在模 y 意义下的逆元,而这是 k \perp y 的充要条件。

那么答案即为:

\begin{aligned} \sum_{i = 1}^{n} \sum_{j = 1}^{m} [i \perp j] [j \perp k] &= \sum_{i = 1}^{m} [i \perp k] \sum_{j = 1}^{n} [i \perp j] \\ &= \sum_{d = 1}^{\min(n, m)} [d \perp k] \mu(d) (n / d) \sum_{id \le m} [i \perp k] \end{aligned}

现在问题转化为求 [n \perp k] \mu(n) [n \perp k] 的前缀和。设:

\begin{aligned} f(n) &= [n \perp k] \\ F(n) &= \sum_{i = 1}^{n} f(i) \\ g(n) &= [n \perp k] \mu(n) \\ G(n) &= \sum_{i = 1}^{n} g(i) \end{aligned}

对于 F(n) ,由于 \gcd(x, y) = \gcd(y, x \bmod{y}) ,有:

\begin{aligned} F(n) &= \sum_{i = 1}^{n} [i \perp k] \\ &= \sum_{i = 1}^{n} [(i \bmod{k}) \perp k] \\ &= (n / k) F(k) + F(n \bmod{k}) \end{aligned}

预处理 F(1), F(2), \dots, F(k) 即可 O(1) 回答,时间复杂度 O(k) - O(1)

对于 G(n) ,发现 [n \perp k] 一项难以处理,考虑将其消掉。
观察到 [x \perp z] [y \perp z] = [xy \perp z] ,故可以将 g 卷上 f 来杜教筛:

\begin{aligned} (g \ast f)(n) &= \sum_{d | n} [d \perp k] \mu(d) \left[\frac{n}{d} \perp k\right] \\ &= \sum_{d | n} [n \perp k] \mu(d) \\ &= [n \perp k] \epsilon(n) \\ &= \epsilon(n) \\ f(1)G(n) &= \sum_{i = 1}^{n} (g \ast f)(i) - \sum_{i = 2}^{n} f(i) G(n / i) \\ G(n) &= 1 - \sum_{i = 2}^{n} f(i) G(n / i) \end{aligned}

时间复杂度 O\left(n^{\frac{2}{3}} + k\right)

Code

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
#include <cstdio>

using i64 = long long;

constexpr int maxk = 2000;
constexpr int maxs = 1000000; // n^{2 / 3}

template <typename _Tp>
inline _Tp Min(const _Tp &x, const _Tp &y) {
return x < y ? x : y;
}

inline int io() {
static int v;
return scanf("%d", &v), v;
}

int pd[6], dcnt;

inline void divide(int k) {
for (int i = 2; i <= k; ++i)
if (k % i == 0) {
pd[++dcnt] = i;
do k /= i; while (k % i == 0);
}
}

inline bool perp(const int &n) {
for (int i = 1; i <= dcnt; ++i)
if (n % pd[i] == 0)
return false;
return true;
}

int pri[maxs / 7], pcnt,
lpf[maxs + 1],
f[maxk + 1],
g[maxs + 1];

inline void precalcF(const int &k) {
f[1] = 1;
for (int i = 2; i <= k; ++i)
f[i] = f[i - 1] + perp(i);
}

inline void sieve(const int &n) {
g[1] = 1; i64 v;
for (int i = 2; i <= n; ++i) {
if (lpf[i] == 0)
pri[lpf[i] = ++pcnt] = i, g[i] = -1;
for (int j = 1; j < lpf[i] && (v = (i64)i * pri[j]) <= n; ++j)
lpf[v] = j, g[v] = -g[i];
if ((v = (i64)i * pri[lpf[i]]) <= n)
lpf[v] = lpf[i];
perp(i) ? g[i] += g[i - 1] : (g[i] = g[i - 1]);
}
}

const int global_n = io(), m = io(), k = io();

inline int F(const int &n) {
return (n / k) * f[k] + f[n % k];
}

namespace hash {

constexpr int hashsz = 1048576;

struct hashlnk {
int key, val; hashlnk *las;
inline hashlnk* init(const int &ik, const int &iv, hashlnk *const &ls) {
return key = ik, val = iv, las = ls, this;
}
} *las[hashsz];

inline int& find(const int &key) {
static hashlnk pool[maxs], *alloc_p = pool - 1;
const int pos = (key & 1048575);
for (hashlnk *o = las[pos]; o; o = o->las)
if (o->key == key)
return o->val;
las[pos] = (++alloc_p)->init(key, 0, las[pos]);
return las[pos]->val;
}

}

inline int G(const int &n) {
if (n <= maxs) return g[n];
int &mem = hash::find(n);
if (mem) return mem;
int ans = 1;
for (int i = 2, j; i <= n; i = j + 1) {
j = n / (n / i);
ans -= (F(j) - F(i - 1)) * G(n / i);
} return mem = ans;
}

int main() {
const int lim = Min(global_n, m);

divide(k);
precalcF(k);
sieve(Min(lim, maxs));

i64 ans = 0;
for (int i = 1, j; i <= lim; i = j + 1) {
j = Min(global_n / (global_n / i), m / (m / i));
ans += (i64)(G(j) - G(i - 1)) * (global_n / i) * F(m / i);
} printf("%lld", ans);

return 0;
}