Luogu3455 [POI2007]ZAP-Queries

Link

Description

给定a,b,d, 求

\sum_{i=1}^{a} \sum_{j=1}^{b} [gcd(i,j) = d]

Sol

还是老套路,先直接把中括号的东西除一个 d

\sum_{i=1}^{a} \sum_{j=1}^{b} [gcd(\frac{i}{d},\frac{j}{d}) = 1]

然后改成枚举约数的

\sum_{i=1}^{\lfloor \frac{a}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{b}{d} \rfloor} [gcd(i,j) = 1]

你可能会想把后面那个 [gcd(i,j)=1] 改成 \varphi 函数,但事实上对这个题一点帮助都没有

考虑使用莫比乌斯反演最常用的性质: \sum_{d|n} \mu(d) = [n=1]

相当于上面这个式子的 n 就是 gcd(i,j) ,枚举 k | gcd(i,j) ,对 \mu(k) 求和

\sum_{i=1}^{\lfloor \frac{a}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{b}{d} \rfloor} \sum_{k|gcd(i,j)} \mu(k)

然后改成枚举 gcd(i,j) 的约数

\sum_{k=1}^{\lfloor \frac{a}{d} \rfloor} \sum_{i=1}^{\lfloor {\frac{a}{dk} \rfloor}} \sum_{j=1}^{\lfloor \frac{b}{dk} \rfloor} \mu(k)

事实上第一个和式 \sum_{k=1}^{\lfloor \frac{a}{d} \rfloor} 的上届应该是 max(\lfloor \frac{a}{d} \rfloor, \lfloor \frac{b}{d}\rfloor)，但是那个上面的位置太小了，写在这里

后面那一对东西可以直接整出分块算掉，所以最终的式子是：

\sum_{k=1}^{\lfloor \frac{a}{d} \rfloor} \mu(k) \times \lfloor {\frac{a}{dk} \rfloor} \times \lfloor \frac{b}{dk} \rfloor

Code

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

#define int long long

const int SIZE = 5e4 + 5;

int T, a, b, d, tot;
int prime[SIZE], isPrime[SIZE], mu[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 = 5e4;
    isPrime[1] = mu[1] = 1;
    for (int i = 2; i <= siz; ++ i) {
        if (!isPrime[i]) {
            prime[++ tot] = i, mu[i] = -1;
        }
        for (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();
    getMu();
    for (int i = 1; i <= 5e4; ++ i) sum[i] = sum[i - 1] + mu[i];
    while (T --) {
        a = read(), b = read(), d = read();
        a /= d, b /= d;
        int range = std::min(a, b), ans = 0;
        for (int l = 1, r; l <= range; l = r + 1) {
            if (!(range / l)) r = range;
            else r = std::min(range, std::min(a / (a / l), b / (b / l)));
            ans += (sum[r] - sum[l - 1]) * (a / l) * (b / l);
        }
        printf("%lld\n", ans);
    }
    return 0;
}