【题目描述】

给定 $n. m, k$,求: $$ \sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) = k] $$ $n, m, k \leq 5 \times 10^4$

有 $T$ 组询问,$T \leq 5 \times 10^4$。

【题目链接】

BZOJ 1011 Zap

【解题思路】

考虑到 $\gcd(i, j) = k$ 等价于 $\gcd(\frac i k, \frac j k) = 1$,那么我们先将 $n, m$ 分别除以 $k$,为方便,将除后的结果仍记作 $n, m$ ,那么我们要求的就是:

$$ \text{ans}(n, m) = \sum_{i = 1}^n\sum_{j = 1}^m[\gcd(i, j) = 1] $$ 通用思路,将$[\gcd(i, j) = 1]$ 通过 $\sum\limits_{d \mid x}^x \mu(d) = [x = 1]$ 即$\mu \times u = e$替换: $$ \text{ans}(n, m) =\sum_{i = 1}^n\sum_{j = 1}^m\sum_{d | \gcd(i, j)} \mu(d) $$ 接下来更改求和指标,改为外层枚举 $d$,并只有在 $d \mid \gcd(i, j)$ 时将$\mu(d)$计入答案,即: $$ ans(n, m) = \sum_{d = 1}^{min(n, m)}\mu(d) \sum_{i = 1}^n \sum_{j = 1}^m [d \mid \gcd(i, j)] $$ 考虑后面那两个 $\sum$,其实它们计算的是 $d$ 成为 $\gcd(i, j)$ 的次数,这个其实是可以快速计算的。

我们知道 $d \mid \gcd(i, j)$,当且仅当 $d \mid i \land d \mid j$ ,而这两个条件是独立的,我们可以分别计算它们各自成立的次数,然后将其相乘,即: $$ \begin{align*} \sum_{i = 1}^n \sum_{j = 1}^m[d \mid \gcd(i, j)] &= \sum_{i = 1}^n \sum_{j = 1}^m[d \mid i \land d \mid j] \\ &= \sum_{i = 1}^n [d \mid i] \cdot \sum_{j = 1}^m [d \mid j] \end{align*} $$ 而这样一个式子 $\sum\limits_{i = 1}^n[d \mid i]$ ,其含义是 $1 \sim n$ 中,$d$ 的倍数的个数,那么其值也就很明了了: $$ \sum_{i = 1}^n[d \mid i] = \left \lfloor \frac n d \right \rfloor $$ 至此,我们的式子成了这个样子: $$ \text {ans}(n, m) = \sum_{d}^{\min(n, m)}\mu(d) \cdot \left \lfloor \frac n d \right \rfloor \cdot \left \lfloor \frac m d \right \rfloor $$ 由于 $\left \lfloor \frac n d \right \rfloor$ 至多有 $2\sqrt n$ 种取值,$\left \lfloor \frac m d \right \rfloor$同理,所以我们可以枚举使其各自取值相同的 $d$ 的区间,然后并在一起计算,一起计算过程的进行需要 $\mu$ 的区间和,可以通过预处理 $\mu$ 的前缀和来计算 。

总结:

利用 $\mu \times u = e$ 对形如 $[\gcd(i, j) = 1]$ 的式子进行代换。

更改求和指标,枚举约数。

考虑每个约数满足条件的次数。

最后一般是利用下取整的性质,成块统计答案。

【AC 代码】

#include <cstdio>
#include <algorithm>

typedef long long int64;

const int MAX_N = 50000;

bool isNotPrime[MAX_N + 1];
int primes[MAX_N], m;
int64 mu[MAX_N + 1];

inline void seive(int n = MAX_N){
    isNotPrime[1] = true, mu[1] = 1;
    for(int i = 2; i <= n; i++){
        if(!isNotPrime[i]){
            primes[m++] = i;
            mu[i] = -1;
        }

        for(int j = 0, p, x; j < m; j++){
            if((x = i * (p = primes[j])) > n) break;
            isNotPrime[x] = true;
            if(i % p == 0){
                mu[x] = 0;
                break;
            } else{
                mu[x] = -mu[i];
            }
        }
    }

    for(int i = 1; i <= n; i++) mu[i] += mu[i - 1];
}

inline int64 solve(int n, int m, int k){
    n /= k, m /= k;
    if(n > m) std::swap(n, m);

    int64 ans = 0;
    for(int l = 1, r; l <= n; l = r + 1){
        r = std::min(n / (n / l), m / (m / l));
        ans += (mu[r] - mu[l - 1]) * (int64)(n / l) * (m / l);
    }

    return ans;
}

int main(){
    seive();

    int T_T;
    scanf("%d", &T_T);
    while(T_T--){
        int n, m, k;
        scanf("%d%d%d", &n, &m, &k);
        printf("%lld\n", solve(n, m, k));
    }

    return 0;
}