Luogu2257 YY的GCD

Link

Description

给定N和M,设 p 为任一质数,求

\sum_{i=1}^{N} \sum_{j=1}^{M} [gcd(i,j) = p]

Sol

常规套路还是先把布尔表达式除一个 p

\sum_{i=1}^{N} \sum_{j=1}^{M} [gcd(\frac{i}{p},\frac{j}{p}) = 1]

然后改成枚举倍数

\sum_{i=1}^{\lfloor \frac{N}{p} \rfloor} \sum_{j=1}^{\lfloor \frac{M}{p} \rfloor} [gcd(i,j) = 1]

[gcd(i,j)=1] 这个东西你还是很熟悉，把他莫比乌斯反演一下

\sum_{i=1}^{\lfloor \frac{N}{p} \rfloor} \sum_{j=1}^{\lfloor \frac{M}{p} \rfloor} \sum_{k|gcd(i,j)} \mu(k)

可以直接枚举质数 p

\sum_{p}^{max(N,M)} \sum_{k=1}^{\lfloor \frac{N}{p} \rfloor} \mu(k) \times \lfloor \frac{N}{kp} \rfloor \times \lfloor \frac{M}{kp} \rfloor

但是你发现在 10^7 的范围内枚举质数还 10^4 组回答不现实。再优化，考虑直接枚举 k\times p

设 T = k \times p.那么

\sum_{T=1}^{N} \lfloor \frac{N}{T} \rfloor \lfloor \frac{M}{T} \rfloor \sum_{p|T \& p \in P} \mu(\frac{T}{p})

其中 P 表示质数集

现在就做完了

Code

#include <bits/stdc++.h>
using namespace std;

#define int long long

const int SIZE = 1e7 + 5;

int T, n, m, tot;
int isPrime[SIZE], prime[SIZE], mu[SIZE], ans[SIZE], sum[SIZE];

namespace ae86 {
    const int bufl = 1 << 15;
    char buf[bufl], *s = buf, *t = buf;
    inline int fetch() {
        if (s == t) { t = (s = buf) + fread(buf, 1, bufl, stdin); if (s == t) return EOF; }
        return *s++;
    }
    inline int read() {
        int a = 0, b = 1, c = fetch();
        while (!isdigit(c))b ^= c == '-', c = fetch();
        while (isdigit(c)) a = a * 10 + c - 48, c = fetch();
        return b ? a : -a;
    }
}
using ae86::read;

inline void getMu()
{
    int siz = SIZE - 5;
    isPrime[1] = mu[1] = 1;
    for (register int i = 2; i <= siz; ++ i) {
        if (!isPrime[i]) {
            prime[++ tot] = i, mu[i] = -1;
        }
        for (register int j = 1; j <= tot && i * prime[j] <= siz; ++ j) {
            isPrime[i * prime[j]] = 1;
            if (i % prime[j]) mu[i * prime[j]] = -mu[i];
            else {
                mu[i * prime[j]] = 0;
                break;
            }
        }
    }
}

signed main()
{
    // freopen("code.in", "r", stdin);
    // freopen("code.out", "w", stdout);
    T = read(); int siz = SIZE - 5;
    getMu();
    for (register int i = 1; i <= tot; ++ i) {
        for (int j = prime[i]; j <= siz; j += prime[i]) ans[j] += mu[j / prime[i]];
    }
    for (register int i = 1; i <= siz; ++ i) sum[i] = sum[i - 1] + ans[i];
    while (T --) {
        n = read(), m = read();
        int range = std::min(n, m), res = 0;
        for (register int l = 1, r = 0; l <= range; l = r + 1) {
            r = std::min(range, std::min(n / (n / l), m / (m / l)));
            res += (n / l) * (m / l) * (sum[r] - sum[l - 1]);
        }
        printf("%lld\n", res);
    }
    return 0;
}