「BZOJ 2820」YY的GCD - 莫比乌斯反演

1 \leq x \leq N 1 \leq y \leq M \gcd(x, y) 为质数的 (x, y) 数量。

链接

BZOJ 2820

题解

首先,我们只需要处理 N \leq M 的情况,当 M \lt N 的情况只需要交换 N M 即可。

设小于 N 质数为 p_1, p_2, …, p_n ,则答案为

\sum\limits_{k = 1} ^ {n} \sum\limits_{i = 1} ^ {N} \sum\limits_{j = 1} ^ {M} [\gcd(i, j) = p_k]

根据莫比乌斯反演,我们有

\sum\limits_{k = 1} ^ {n} \sum\limits_{i = 1} ^ {N} \sum\limits_{j = 1} ^ {M} [\gcd(i, j) = p_k] = \sum\limits_{k = 1} ^ {n} \sum\limits_{d = 1} ^ {\lfloor \frac{N}{p_k} \rfloor} \mu(d) \lfloor \frac{N}{p_k \times d} \rfloor \lfloor \frac{M}{p_k \times d} \rfloor

T = p_k \times d ,我们在外层枚举 T 然后对每个质因子计算 \mu

\sum\limits_{T = 1} ^ {N} \lfloor \frac{N}{T} \rfloor \lfloor \frac{M}{T} \rfloor \sum\limits_{k \mid T} \mu(\frac{T}{k})

f(T) = \sum\limits_{k \mid T} \mu(\frac{T}{k})

T = p_1 ^ {x_1} \times p_2 ^ {x_2} \times … p_n ^ {x_n}

T' = p_1 ^ {x_1 - 1} \times p_2 ^ {x_2} \times … p_n ^ {x_n}

考虑线性筛求 \mu 的过程,当 T' \ {\rm mod} \ p_1 = 0

\begin{align} f(T') &= \sum\limits_{i = 2} ^ {k} \mu(\frac{T'}{p_i}) \\ f(T) &= \sum\limits_{i = 1} ^ {k} \mu(\frac{T}{p_i}) \\ &= \mu(\frac{T}{p_i}) + \sum\limits_{i = 2} ^ {k} \mu(\frac{T}{p_i}) \\ &= \mu(T') + \sum\limits_{i = 2} ^ {k} \mu(\frac{T'}{p_i} \times p_1) \\ &= \mu(T') + \sum\limits_{i = 2} ^ {k} \mu(\frac{T'}{p_i}) \times \mu(p_1) \\ \mu(p_i) &= 1 \\ f(T) &= \mu(T') - f(T') \end{align}

x_1 \gt 1

\begin{align} \mu(p_1 ^ {x_1}) &= 0 \\ f(T) &= \sum\limits_{i = 1} ^ {k} \mu(\frac{T}{p_i}) \\ &= \mu(\frac{T}{p_1}) + \sum\limits_{i = 2} ^ {k} \mu(\frac{T}{p_i}) \\ &= \mu(T') + \sum\limits_{i = 2} ^ {k} \mu(\frac{T}{p_i \times p_1 ^ {x_1}} \times p_1 ^ {x_1}) \\ &= \mu(T') + \sum\limits_{i = 2} ^ {k} \mu(\frac{T}{p_i \times p_1 ^ {x_1}}) \times \mu(p_1 ^ {x_1}) \\ &= \mu(T') \end{align}

线性筛后分块处理即可。

代码

#include <cstdio>
#include <cmath>
#include <vector>

const int MAXT = 10000;
const int MAXN = 10000000;
const int MAXM = 10000000;

bool isNotPrime[MAXN + 1];
int miu[MAXN + 1], f[MAXN + 1], s[MAXN + 1];
std::vector<int> primes;
inline void getPrimes() {
    primes.reserve(MAXN);
    isNotPrime[0] = isNotPrime[1] = true;
    for (int i = 2; i <= MAXN; i++) {
        if (!isNotPrime[i]) {
            primes.push_back(i);
            miu[i] = -1;
            f[i] = 1;
        }
        for (std::vector<int>::const_iterator p = primes.begin(); p != primes.end() && i * *p <= MAXN; p++) {
            isNotPrime[i * *p] = true;
            if (i % *p == 0) {
                miu[i * *p] = 0;
                f[i * *p] = miu[i];
                break;
            } else {
                miu[i * *p] = -miu[i];
                f[i * *p] = miu[i] - f[i];
            }
        }
    }

    for (int i = 1; i <= MAXN; i++) s[i] = s[i - 1] + f[i];
}

inline long long solve(const int n, const int m) {
    long long ans = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = std::min(n / (n / l), m / (m / l));
        ans += (s[r] - s[l - 1]) * static_cast<long long>(n / l) * static_cast<long long>(m / l);
    }

    return ans;
}

int main() {
    getPrimes();
    int t;
    scanf("%d", &t);
    while (t--) {
        int n, m;
        scanf("%d %d", &n, &m);
        if (n > m) std::swap(n, m);
        printf("%lld\n", solve(n, m));
    }

    return 0;
}