BZOJ 3994 [SDOI2015]约数个数和 莫比乌斯反演

题目链接

题意

求

$$\sum_{i=1}^{N}\sum_{j=1}^{M}d(ij)$$ 其中 $d(x)$ 为 $x$ 的约数个数

结论

$$d(ij)=\sum_{ii|i}\sum_{jj|j}[gcd(ii,jj)=1]$$
证明参见 PoPoQQQ
还有一版 iwtwiioi

推导

$$\begin{align}原式 &=\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{ii|i}\sum_{jj|j}[gcd(ii,jj)=1]\\ &=\sum_{ii=1}^{N}\sum_{jj=1}^{M}[gcd(ii,jj)=1]\sum_{i=1}^{\lfloor\frac{N}{ii}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{jj}\rfloor}\\ &=\sum_{i=1}^{N}\sum_{j=1}^{M}\lfloor\frac{N}{i}\rfloor\lfloor\frac{M}{j}\rfloor[gcd(i,j)=1]\\ &=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\lfloor\frac{N}{id}\rfloor\sum_{i=1}^{\lfloor\frac{M}{d}\rfloor}\lfloor\frac{M}{jd}\rfloor \end{align}$$
接下来的只要处理出 $$f(k)=\sum_{i=1}^{k}\lfloor\frac{k}{i}\rfloor$$ 即可。
这里
1. 可以直接 分块 $O(N\sqrt N)$预处理.
2. 注意到 $$\sum_{i=1}^{k}d(i)=\sum_{i=1}^{k}\sum_{j=1}^{i}[j|i]=\sum_{j=1}^{k}\sum_{i=1}^{\lfloor\frac{k}{j}\rfloor}=\sum_{j=1}^{k}\lfloor\frac{k}{j}\rfloor$$
即 $f(k)=\sum_{i=1}^{k}d(i)$，这样就又可以 线性筛 $O(N)$解决了.

Code

注释部分为法一分块求 $f(k)$

#include <bits/stdc++.h>
#define maxn 50000
#define maxm maxn + 10
typedef long long LL;
using namespace std;
bool check[maxm];
int prime[maxm], mu[maxm], cnt[maxn];
LL pre[maxm], sum[maxm], d[maxn];
//void init() {
//    int tot = 0; mu[1] = 1;
//    for (int i = 2; i <= maxn; ++i) {
//        if (!check[i]) {
//            prime[tot++] = i;
//            mu[i] = -1;
//        }
//        for (int j = 0; j < tot; ++j) {
//            if (i * prime[j] > maxn) break;
//            check[i * prime[j]] = true;
//            if (i % prime[j] == 0) {
//                mu[i * prime[j]] = 0;
//                break;
//            }
//            mu[i * prime[j]] = -mu[i];
//        }
//    }
//    for (int i = 1; i <= maxn; ++i) {
//        int lim = maxn / i;
//        for (int j = 0; j <= lim; ++j) {
//            int w = i * j;
//            sum[w] += j;
//            if (w + i <= maxn) sum[w + i] -= j;
//        }
//    }
//    for (int i = 0; i <= 20; ++i) printf("%d ", sum[i]); printf("\n");
//
//    for (int i = 1; i <= maxn; ++i) {
//        pre[i] = pre[i - 1] + mu[i];
//        sum[i] += sum[i - 1];
//    }
//}
void init() {
    int tot = 0; mu[1] = d[1] = 1;
    for (int i = 2; i <= maxn; ++i) {
        if (!check[i]) {
            prime[tot++] = i;
            mu[i] = -1, d[i] = 2, cnt[i] = 1;
        }
        for (int j = 0; j < tot; ++j) {
            if (i * prime[j] > maxn) break;
            check[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                cnt[i * prime[j]] = cnt[i] + 1;
                d[i * prime[j]] = d[i] / (cnt[i] + 1) * (cnt[i * prime[j]] + 1);
                break;
            }
            mu[i * prime[j]] = -mu[i];
            cnt[i * prime[j]] = 1;
            d[i * prime[j]] = d[i] << 1;
        }
    }

    for (int i = 1; i <= maxn; ++i) {
        pre[i] = pre[i - 1] + mu[i];
        sum[i] = sum[i - 1] + d[i];
    }
}
void work() {
    int n, m;
    scanf("%d%d", &n, &m);
    int lim = min(n, m), le, ri;
    LL ans = 0;
    for (int i = 1; i <= lim; i = ri + 1) {
        le = i, ri = min(n / (n / i), m / (m / i));
        ans += (pre[ri] - pre[le - 1]) * sum[n / i] * sum[m / i];
    }
    printf("%lld\n", ans);
}
int main() {
    init();
    int T;
    scanf("%d",&T);
    while (T--) work();
    return 0;
}