[NOI2016] 循环之美

填坑。

先考虑题面，若 \(i \over j\) 最简，必定有 \(i \perp j\)。 然后要在 \(k\) 进制下循环，设循环节长为 \(x\)，有小学奥数经典套路：

\[{i\over j} - \lfloor{i \over j}\rfloor = {ik^x \over j} - \lfloor{ik^x\over j}\rfloor \]

两边都乘 \(j\), 再对 \(j\) 取模，约分掉 \(i\) 得到：

\[k^x \equiv 1 \pmod j \]

由欧拉定理知，满足 \(j\perp k\) 一定有解。而不满足的话，\(\gcd(k,j)\neq 1\), \(\gcd(k^x,j)\neq 1\)，而 \(k^x\bmod j = 1\), \(\gcd(k^x,j)=\gcd(k^x,k^x\bmod j)=\gcd(k^x,1)=1\)，前后矛盾，无解。

因此我们可以得到式子：

\[\sum_{i=1}^n\sum_{j=1}^m [i ⊥ j][j ⊥ k] \]

注意 \(i\perp j\) 的含义是 \(\gcd(i,j)=1\)，我们得到：

\[\sum_{i=1}^n\sum_{j=1}^m \sum_{d|\gcd(i,j)} \mu(d) [j ⊥ k] \]

所以有：

\[\sum_{d}\mu(d)\sum_{i=1}^{n\over d} \sum_{j = 1}^{m\over d}[jd⊥ k] \]

\([jd⊥k]\) 可以拆成 \([j⊥k]\) 和 \([d⊥k]\)。后面的式子和 \(i\) 没什么关系，直接乘上 \(n\over d\):

\[\sum_{d}\mu(d){n\over d}[d⊥k]\sum_{j=1}^{m\over d} [j⊥k] \]

我们令 \(g(n)\) 表示后面那个和式，即：

\[g(n)=\sum_{i=1}^n [j⊥k] \]

由于 \(\gcd(i,k)=\gcd(i + k, k)\)，所以：

\[\sum_{i=ak+1}^{(a+1)k}[i⊥k]=\sum_{i=1}^k [i⊥k]=\varphi(k) \]

我们轻松得到：\(g(n)=g(n \bmod k) + {n\over k}\varphi(k)\)，\(\Theta(k\log k)\) 预处理出 \(g(1)~g(k)\) 后可以 \(\Theta(1)\) 求。

带回原式得到：

\[\sum_{d}{n\over d}g({m\over d})\mu(d)[d⊥k] \]

前面两项可以整除分块搞，考虑后面的怎么求前缀和，即求。

\[f(n)=\sum_{i=1}^n \mu(i)[i⊥k] \]

直接推有点难。考虑杜教筛。难点在于求和 \([i ⊥ k]\) 卷起来后的前缀和，简单推一下：

\[\sum_{i=1}^n \sum_{d|i} \mu(d)[d⊥k][{i\over d} ⊥ k]=\sum_{i=1}^n [i⊥k][i=1]=1 \]

于是我们就愉快的 AC 了。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <map>

using namespace std;

const int N = 2e3 + 5, Lim = 1e6 + 5, lim = 1e6;

int mu[Lim], g[N];

inline int gcd(int a, int b) {
	return !b ? a : gcd(b, a % b);
}

int prime[Lim], top = 0, notprime[Lim], phk = 0, n, m, k, f[Lim];

map<int, int> mp;

inline void pre() {
	mu[1] = 1; notprime[1] = 1;
	for (int i = 2; i <= lim; ++i) {
		if (!notprime[i]) prime[++top] = i, mu[i] = -1;
		for (int j = 1; j <= top && i * prime[j] <= Lim; ++j) {
			notprime[i * prime[j]] = 1;
			if (i % prime[j]) mu[i * prime[j]] = -mu[i];
			else { mu[i * prime[j]] = 0; break; }
		}
	} for (int i = 1; i <= k; ++i) g[i] = g[i - 1] + (gcd(i, k) == 1);
	phk = g[k];
	for (int i = 1; i <= lim; ++i) f[i] = f[i - 1] + (gcd(i, k) == 1) * mu[i];
}

typedef long long ll;

inline int G(int x) {
	return (x / k) * phk + g[x % k]; 
}

inline int min_(int a, int b) {
	return a < b ? a : b;
}

inline int F(int n) {
	if (n <= lim) return f[n];
	if (mp[n]) return mp[n];
	int ret = 1;
	for (int l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ret -= F(n / l) * (G(r) - G(l - 1));
	} return mp[n] = ret;
}

int main() {
	scanf("%d%d%d", &n, &m, &k);
	pre(); ll ans = 0;
	for (int l = 1, r; l <= min_(n, m); l = r + 1) {
		r = min_(n / (n / l), m / (m / l));
		ans += 1ll * (n / l) * G(m / l)  * (F(r) - F(l - 1));
	} printf("%lld\n", ans);
	return 0;
}