链接 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$ 时，就是莫比乌斯函数的前缀和，这时我们用杜教筛求就好了。
代码 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 #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 ; }
