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

题目链接

题意

求

∑i=1N∑j=1Md(ij)

其中 d(x) 为 x 的约数个数

结论

d(ij)=∑ii|i∑jj|j[gcd(ii,jj)=1]

证明参见PoPoQQQ
还有一版iwtwiioi

推导

原式=∑i=1N∑j=1M∑ii|i∑jj|j[gcd(ii,jj)=1]=∑ii=1N∑jj=1M[gcd(ii,jj)=1]∑i=1⌊Nii⌋∑j=1⌊Mjj⌋=∑i=1N∑j=1M⌊Ni⌋⌊Mj⌋[gcd(i,j)=1]=∑d=1min(n,m)μ(d)∑i=1⌊Nd⌋⌊Nid⌋∑i=1⌊Md⌋⌊Mjd⌋

接下来的只要处理出

f(k)=∑i=1k⌊ki⌋

即可。
这里
1. 可以直接分块O(NN−−√)预处理.
2. 注意到

∑i=1kd(i)=∑i=1k∑j=1i[j|i]=∑j=1k∑i=1⌊kj⌋=∑j=1k⌊kj⌋

即f(k)=∑ki=1d(i)，这样就又可以线性筛O(N)解决了.

Code

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

#include <bits/stdc++.h>#define maxn 50000#define maxm maxn + 10typedef 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;}