循环之美

循环之美

首先我们得知道一个结论:当分母与进制数互质时,它是一个纯循环小数.知道这个结论后,题目就是要我们求下式: \[ \begin{aligned}ans=&\sum_{x = 1}^n\sum_{y=1}^m\left[x\perp y\right]\left[y \perp k\right]\\=&\sum_{x=1}^n\sum_{y=1,k\perp y}^m\sum_{d|x,d|y}\mu(d)\\=&\sum_{k\perp d}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\left[k\perp y\right]\\=&\sum_{k\perp d}\mu(d)\left\lfloor\frac{n}{d}\right\rfloor\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\left[k\perp y\right]\end{aligned} \]\[ \begin{aligned}f(n)=&\sum_{i=1}^n\left[k\perp i\right]\end{aligned} \]

因为\(gcd(n, k) =gcd(n\bmod k, k)\), 所以有: \[ f(n) = \left\lfloor\frac{n}{k}\right\rfloor f(k)+f(n \bmod k) \] 所以我们直接预处理\(f\)的前\(k\)项就好了, 现在问题在于求: \[ \begin{aligned}S(n,k)=&\sum_{d = 1}^n\left[k\perp d\right]\mu(d)\\=&\sum_{d=1}^n\mu(d)\sum_{i|d, i|k}\mu(i)\\=&\sum_{i|k}\mu(i)\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}\mu(d\times i)\\=&\sum_{i|k}\mu^2(i)\sum_{d=1}^{\lfloor\frac{n}{i}\rfloor}\left[i\perp d\right]\mu(d)\\=&\sum_{i|k}\mu^2(i)S(\left\lfloor\frac{n}{i}\right\rfloor, i)\end{aligned} \] 又可以递归了.倒数第四个转化妙妙啊.

边界就是当\(k=1\)时,直接返回\(\mu\)的前缀和就好了,这个直接杜教筛一下就\(ok\)了.

复杂度不太会证.所以\(O(能过)\)

代码:

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
/*
* 4652.cpp
* This file is part of BZOJ_4652
*
* Copyright (C) 2019 - ViXbob
*
* BZOJ_4652 is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* BZOJ_4652 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with BZOJ_4652. If not, see <http://www.gnu.org/licenses/>.
*/

// Begin at 2019-04-25 18:22:47 - ViXbob

/**
* There is no end though there is a start in space. ---Infinity.
* It has own power, it ruins, and it goes though there is a start also in the star. ---Finite.
* Only the person who was wisdom can read the most foolish one from the history.
* The fish that lives in the sea doesn't know the world in the land.
* It also ruins and goes if they have wisdom.
* It is funnier that man exceeds the speed of light than fish start living in the land.
* It can be said that this is an final ultimatum from the god to the people who can fight.
*
* Steins;Gate
*/
#include <bits/stdc++.h>
#define rep(i, j, k) for(int i = j; i <= k; ++i)
#define dep(i, j, k) for(int i = j; i >= k; --i)

using namespace std;

typedef long long ll;

const int maxn = 2e3 + 5;
const int maxm = 1e6 + 5;

int n, m, k, f[maxn], prime[maxm], mu[maxm];
bool vis[maxm];
unordered_map<int, int> Mu;
vector<int> fac[maxn];

inline int read() {
char ch = getchar(); int u = 0, f = 1;
while(!isdigit(ch)) { if(ch == '-') f = -1; ch = getchar(); }
while(isdigit(ch)) { u = u * 10 + ch - 48; ch = getchar(); } return u * f;
}

inline void PreSolve(int maxk) {
mu[1] = vis[1] = 1;
rep(i, 2, maxm - 1) {
if(!vis[i]) prime[++prime[0]] = i, mu[i] = -1;
for(int j = 1; j <= prime[0] && prime[j] * i < maxm; j++) {
vis[prime[j] * i] = 1;
if(i % prime[j]) mu[prime[j] * i] = -mu[i]; else break;
}
}
rep(i, 1, maxm - 1) mu[i] += mu[i - 1];
f[1] = 1; rep(i, 2, maxk) f[i] = f[i - 1] + (__gcd(i, k) == 1);
rep(i, 1, maxk) for(int j = i; j <= maxk; j += i) fac[j].push_back(i);
}

inline int get_mu(int n, int rnt = 1) {
if(n < maxm) return mu[n];
if(Mu.count(n)) return Mu[n];
for(int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
rnt -= get_mu(n / l) * (r - l + 1);
} return (Mu[n] = rnt);
}

inline ll F(int n) { return (n / k) * f[k - 1] + f[n % k]; }

inline int g(int n, int k) {
// cerr << n << " " << k << endl;
if(k == 1) return get_mu(n);
int rnt = 0;
for(auto v : fac[k]) {
// cerr << n << " / " << v << " = " << n / v << endl;
if(n / v <= 0) break;
if(mu[v] - mu[v - 1]) rnt += g(n / v, v);
} return rnt;
}

int main() {
n = read(); m = read(); k = read();
PreSolve(k); ll ans = 0;
for(int l = 1, r, lst = 0; l <= min(n, m); l = r + 1) {
r = min(n / (n / l), m / (m / l));
// cerr << "l = " << l << ", r = " << r << endl;
int now = g(r, k);
// cerr << "now = " << now << endl;
ans += (now - lst) * (n / l) * F(m / l);
lst = now;
} cout << ans << endl;
return 0;
}
/*
1 1 1.00000000
1 2 0.50000000
1 3 0.33333333
1 4 0.25000000
1 5 0.20000000
1 6 0.16666666
2 1 2.00000000
2 3 0.66666666
2 5 0.40000000
*/