LibreOJ 2085 「NOI2016」循环之美

首先,为了去重,我们可以增加一个限制:必须是最简分数。
然后我们看看什么情况会产生纯循环小数。

\(10\) 进制下,当 \(\gcd(x,y) = 1,\gcd(x,10) = 1\) 时,\(\dfrac y x\) 是纯循环小数。
于是我们猜测当 \(\gcd(x,y) = 1,\gcd(x,k) = 1\) 时,\(\dfrac y x\)\(k\) 进制下是纯循环小数。

证明:
\(\dfrac y x\)\(k\) 进制下是纯循环小数,当且仅当 \(\exists l \neq 0,\dfrac y x - \left\lfloor\dfrac y x\right\rfloor = \dfrac {yk^l} x - \left\lfloor \dfrac {yk^l} x \right\rfloor\)
这玩意看着很像取模的形式,所以两边同乘 \(x\),变成 \(y \equiv yk^l \pmod x\)

由于前面说了 \(\gcd(x,y) = 1\),所以 \(k^l \equiv 1 \pmod x\)
于是由欧拉定理得 \(\gcd(k,x) = 1\),一个特解为 \(l = \varphi(x)\)

证毕。

问题变成了 \(\sum\limits_{i = 1}^n\sum\limits_{j = 1}^m [\gcd(i,j) = 1][\gcd(j,k) = 1]\)

\[\begin{align*} & \sum\limits_{i = 1}^n\sum\limits_{j = 1}^m [\gcd(i,j) = 1][\gcd(j,k) = 1] \\ = & \sum\limits_{j = 1}^m [\gcd(j,k) = 1] \sum\limits_{i = 1}^n [\gcd(i,j) = 1] \\ = & \sum\limits_{j = 1}^m [\gcd(j,k) = 1] \sum\limits_{i = 1}^n \sum\limits_{d | i,d | j} \mu(d) \\ = & \sum\limits_{j = 1}^m [\gcd(j,k) = 1] \sum\limits_{d | j} \lfloor \dfrac n d \rfloor \mu(d) \\ = & \sum\limits_{d = 1}^{\min(n,m)} \lfloor \dfrac n d \rfloor \mu(d) \sum\limits_{jd \le m} [\gcd(jd,k) = 1] \\ = & \sum\limits_{d = 1}^{\min(n,m)} \lfloor \dfrac n d \rfloor [\gcd(d,k) = 1] \mu(d) \sum\limits_{jd \le m} [\gcd(j,k) = 1] \end{align*}\]

\(f(n) = \sum\limits_{i = 1}^n [\gcd(i,k) = 1]\)
由于 \(\gcd(i,k) = \gcd(i + k,k)\),有 \(f(n) = \lfloor \dfrac n k \rfloor \varphi(k) + f(n \bmod k)\)
发现 \(k\) 非常小,暴力求出 \(x \in [0,k]\)\(f(x)\) 值即可。

于是有 \[\begin{align*} & \sum\limits_{d = 1}^{\min(n,m)} \lfloor \dfrac n d \rfloor [\gcd(d,k) = 1] \mu(d) \sum\limits_{jd \le m} [\gcd(j,k) = 1] \\ = & \sum\limits_{d = 1}^{\min(n,m)} \lfloor \dfrac n d \rfloor f(\lfloor \dfrac m d \rfloor) [\gcd(d,k) = 1] \mu(d) \end{align*}\]

如果数论分块的话,就剩下后面的两项。
考虑怎么算它们的前缀和。

\(S(n,k) = \sum\limits_{i = 1}^n [\gcd(i,k) = 1] \mu(i)\)
\[\begin{align*} S(n,k) = & \sum\limits_{i = 1}^n [\gcd(i,k) = 1] \mu(i) \\ = & \sum\limits_{i = 1}^n \mu(i) \sum\limits_{d | i,d | k} \mu(d) \\ = & \sum\limits_{d | k} \mu(d) \sum\limits_{id \le n} \mu(id) \end{align*}\]

乍一看推不下去了。
但是我们反观 \(\mu\) 函数,根据定义容易发现当 \(\gcd(x,y) \ne 1\) 时,\(\mu(xy) = 0\)

所以有 \[\begin{align*} S(n,k) = & \sum\limits_{d | k} \mu(d) \sum\limits_{id \le n} \mu(id) \\ = & \sum\limits_{d | k} \mu(d) \sum\limits_{id \le n} [\gcd(i,d) = 1] \mu(id) \\ = & \sum\limits_{d | k} \mu^2(d) \sum\limits_{id \le n} [\gcd(i,d) = 1] \mu(i) \\ = & \sum\limits_{d | k} \mu^2(d) S(\lfloor \dfrac n d \rfloor,d) \end{align*}\]

然后考虑递归求解 \(S\),容易发现 \(S(n,1) = \sum\limits_{i = 1}^n \mu(i)\),杜教筛即可。

于是最后对于 \(\sum\limits_{d = 1}^{\min(n,m)} \lfloor \dfrac n d \rfloor f(\lfloor \dfrac m d \rfloor) [\gcd(d,k) = 1] \mu(d)\) 数论分块即可。

略微卡常,少开 long long,算 \(S\) 的时候如果系数是 \(0\) 就不用往下算。
unordered_map 居然不能用 pair,偷懒换成了 map(

代码:

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
#include <cstdio>
#include <utility>
#include <algorithm>
#include <map>
#include <tr1/unordered_map>
using namespace std;
using namespace tr1;
const int N = 1e9;
const int K = 2e3;
const int CNT = 1e7;
int n,m,k;
int vis[CNT + 5],prime[CNT + 5],cnt,mu[CNT + 5];
int f[K + 5];
unordered_map<int,int> w;
map<pair<int,int>,long long> mem;
long long ans;
int calc(int n)
{
if(n <= CNT)
return mu[n];
if(w.count(n))
return w[n];
long long ret = 1;
for(register int l = 2,r;l <= n;l = r + 1)
{
r = n / (n / l);
ret -= (long long)(r - l + 1) * calc(n / l);
}
return w[n] = ret;
}
int F(int n)
{
return f[k] * (n / k) + f[n % k];
}
long long S(int n,int k)
{
if(n == 0 || n == 1)
return n;
if(k == 1)
return calc(n);
if(mem.count(make_pair(n,k)))
return mem[make_pair(n,k)];
long long ret = 0;
for(int i = 1;i * i <= k;++i)
if(!(k % i))
{
if(mu[i] - mu[i - 1])
ret += (mu[i] - mu[i - 1]) * (mu[i] - mu[i - 1]) * S(n / i,i);
if(i ^ (k / i) && mu[k / i] - mu[k / i - 1])
ret += (mu[k / i] - mu[k / i - 1]) * (mu[k / i] - mu[k / i - 1]) * S(n / (k / i),k / i);
}
return mem[make_pair(n,k)] = ret;
}
int main()
{
mu[1] = 1;
for(register int i = 2;i <= CNT;++i)
{
if(!vis[i])
mu[prime[++cnt] = i] = -1;
for(register int j = 1;j <= cnt && i * prime[j] <= CNT;++j)
{
vis[i * prime[j]] = 1;
if(!(i % prime[j]))
break;
else
mu[i * prime[j]] = -mu[i];
}
mu[i] += mu[i - 1];
}
scanf("%d%d%d",&n,&m,&k);
for(register int i = 1;i <= k;++i)
f[i] = f[i - 1] + (__gcd(i,k) == 1);
for(register int l = 1,r;l <= min(n,m);l = r + 1)
{
r = min(n / (n / l),m / (m / l));
ans += (S(r,k) - S(l - 1,k)) * (n / l) * F(m / l);
}
printf("%lld\n",ans);
}