首先，一个分数 $\frac{p}{q}$ 在 $k$ 进制下是纯循环小数，就是 $\exists t$，$p \equiv p \times k^t \pmod q$，其中 $p \bot q$，那么即为 $k^t \equiv 1 \pmod q$，要存在解，就得 $k \bot q$。
所以答案就是 $\sum_{i=1}^n\sum_{j=1}^m [i\bot j][k\bot q]$
$$\sum_{i=1}^n\sum_{j=1}^m [i\bot j][k\bot j]$$
$$=\sum_{i=1}^n\sum_{j=1}^m[k\bot j]\sum_{d|(i,j)}\mu(d)$$
$$=\sum_{d=1}^{\min\{n,m\}}[d\bot k]\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[j\bot k]$$
$$=\sum_{d=1}^{\min\{n,m\}}[d\bot k]\mu(d)\lfloor\frac{n}{d}\rfloor s(\lfloor\frac{m}{d}\rfloor)$$
其中 $s(n)=\sum_{i=1}^{n}[i\bot k]=\lfloor\frac{n}{k}\rfloor s(k) + s(n$ $\text{mod}$ $k)$，这个暴力预处理即可。
要求 $f(n)=\mu(n)[n\bot k]$ 的前缀和，设 $g(n)=[n\bot k]$。
$$(f * g)(n)$$
$$=\sum_{d\mid n}\mu(d)[d\bot k][\frac{n}{d}\bot k]$$
$$=[n\bot k]\sum_{d\mid n}\mu(d)$$
$$=[n \bot k][n=1]$$
所以 $\sum_{i=1}^n (f * g)(i)=1$
所以 $f$ 的前缀和 $S(n)=1-\sum_{i=2}^{n}[i\bot k]S(\lfloor\frac{n}{i}\rfloor)$
至此解决。

#include <bits/stdc++.h>

const int N = 1e6, XN = N + 7;
int prin, prime[XN], mu[XN];
int fs[XN], gs[XN];

int gcd(int a, int b) {
    while (b) {
        a %= b;
        std::swap(a, b);
    }
    return a;
}

void init(int k) {
    static bool vis[XN];
    mu[1] = 1;
    for (int i = 2; i <= N; i++) {
        if (!vis[i]) {
            prime[++prin] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= prin && i * prime[j] <= N; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= N; i++) {
        fs[i] = fs[i - 1] + mu[i] * (gcd(i, k) == 1);
        gs[i] = gs[i - 1] + (gcd(i, k) == 1);
    }
}

int k, n, m;

int GS(int n) {
    return (n / k) * gs[k] + gs[n % k];
}
#define ll long long
std::unordered_map<int, ll> mp;

ll solve(int n) {
    if (n <= N) return fs[n];
    if (mp.count(n)) return mp[n];
    ll ans = 1;
    for (int i = 2, j; i <= n; i = j + 1) {
        j = n / (n / i);
        ans -= 1LL * (GS(j) - GS(i - 1)) * solve(n / i);
    }
    return mp[n] = ans;
}

signed main() {
    scanf("%d%d%d", &n, &m, &k);
    init(k);
    ll ans = 0;
    for (int i = 1, j; i <= std::min(n, m); i = j + 1) {
        j = std::min(n / (n / i), m / (m / i));
        ans += 1LL * (solve(j) - solve(i - 1)) * (n / i) * GS(m / i);
    }
    printf("%lld\n", ans);
    return 0;
}

View Code