「BZOJ 3529」「SDOI 2014」数表「莫比乌斯反演」

题意

有一张 \(n\times m\) 的数表，其第\(i\)行第\(j\)列的数值为能同时整除\(i\)和\(j\)的所有自然数之和。

\(T\)组数据，询问对于给定的 \(n,m,a\) , 计算数表中\(\leq a\) 的数之和。

\(T \leq 2\times 10^4,1 \leq n,m\leq 10^5\).

题解

令\(\sigma(x)\)表示\(x\)的约数和，容易写出答案的式子：

\[\sum_{i=1}^n \sum_{j=1}^m\sigma(\gcd(i,j))[\sigma(\gcd(i,j))\leq a] \]

然后使用常见套路变换式子：（下面默认\(n\leq m\)）

\[=\sum_{d=1}^n \sigma(d)[\sigma(d)\leq a] \sum_{i=1}^n \sum_{j=1}^m[\gcd(i,j)=d] \]

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

\[=\sum_{d=1}^n \sigma(d)[\sigma(d)\leq a] \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor} \sum_{j=1}^{\lfloor\frac{m}{d}\rfloor} \sum_{d'|i,d'j} \mu(d') \]

\[=\sum_{d=1}^n \sigma(d)[\sigma(d)\leq a]\sum_{d'=1}^{\lfloor\frac{n}{d} \rfloor}\mu(d') \lfloor\frac{n}{dd'}\rfloor \lfloor\frac{m}{dd'}\rfloor \]

我们记前面那部分\(f(d)=\sigma(d)[\sigma(d)\leq a]\)：

\[=\sum_{i=1}^n f(d)\sum_{d'=1}^{\lfloor\frac{n}{d} \rfloor}\mu(d') \lfloor\frac{n}{dd'}\rfloor \lfloor\frac{m}{dd'}\rfloor \]

枚举乘积\(k=dd'\)：

\[=\sum_{k=1}^n\lfloor\frac{n}{k}\rfloor \lfloor\frac{m}{k}\rfloor\sum_{d|k}f(d)\mu(\frac{k}{d}) \]

可以发现后面变成了狄利克雷卷积的形式！

记$$g(k)=\sum_{d|k} f(d)\mu(\frac{k}{d})$$，即\(g=f*\mu\)：

\[=\sum_{i=1}^n \lfloor \frac{n}{k} \rfloor \lfloor \frac{m}{k}\rfloor g(k) \]

可以看出我们对于每个\(n,m\)，可以数论分块

但是\(g(k)\)会随着\(a\)的变化而变化，不能每次计算一遍\(g\)

但是注意到，\(g\)变化是因为\(f\)，\(f\)只有\(10^5\)种取值。

于是我们可以把询问离线后按\(a\)从小到大排序，每次加入一些满足\(\sigma(d)\leq a\)的\(d\)

然后考虑加入一个\(d\)有什么影响：使得所有\(d|k\)的\(g(k)\)加上\(f(d)\mu(\frac{k}{d})\)

就可以每次枚举\(d\)的倍数\(k\)，然后更新\(g(k)\)，最多更新\([1,n]\)这些数，根据调和级数可得枚举的总复杂度为\(O(n \log n)\)

别忘了数论分块是需要前缀和的，

所以现在我们需要一种支持快速单点加、求前缀和的数据结构，树状数组就很合适

然后就做完了，复杂度\(O(n \log^2 n+T\sqrt n \log n)\)

一个小细节：\(\sigma\)函数怎么线性筛？

考虑质因数分解后\(\sigma (x)=(1+p_1+p_1^2+..+p_1^{c_1})..(1+p_k+p_k^2+..+p_k^{c_k})\)

因此线性筛的时候记录一个\(tmp(x)\)表示\(1+p+p^2..+p^c\)，其中\(p\)为\(x\)的最小质因子，\(c\)的\(p\)的指数；然后就可以做了

一个减小常数的小技巧：分块的时候记录一个变量，减少重复的询问，详见代码

#include <algorithm>
#include <cstdio>
using namespace std;

typedef long long ll;

const int R = 1e5;
const int N = 1e5 + 10;
const int M = 2e4 + 10;
const ll mo = 1ll << 31;

struct qs {
    int n, m, a, id;
    bool operator < (const qs &b) const {
        return a < b.a;
    }
} q[M];
int t, tot, p[N], mu[N], sig[N], tmp[N], ans[M];
bool tag[N];

struct num {
    int x;
    bool operator < (const num &b) const {
        return sig[x] < sig[b.x];
    }
} a[N];

void sieve(int n) {
    sig[1] = mu[1] = 1;
    for(int i = 2; i <= n; i ++) {
        if(!tag[i]) {
            p[tot ++] = i; mu[i] = -1;
            sig[i] = tmp[i] = i + 1;
        }
        for(int j = 0; j < tot; j ++) {
            if(p[j] * i > n) break ;
            tag[i * p[j]] = 1;
            if(i % p[j] == 0) {
                mu[i * p[j]] = 0;
                tmp[i * p[j]] = tmp[i] * p[j] + 1;
                sig[i * p[j]] = sig[i] / tmp[i] * tmp[i * p[j]];
                break ;
            }
            mu[i * p[j]] = - mu[i];
            tmp[i * p[j]] = 1 + p[j];
            sig[i * p[j]] = sig[i] * (1 + p[j]);
        }
    }
    for(int i = 1; i <= n; i ++) a[i].x = i;
    sort(a + 1, a + n + 1);
}

int bit[N];

void add(int x, int y) {
    for(; x <= R; x += x & (-x))
        bit[x] = ((ll) bit[x] + y) % mo;
}

int qry(int x) {
    int ans = 0;
    for(; x >= 1; x &= x - 1)
        ans = ((ll) ans + bit[x]) % mo;
    return ans;
}

void ins(int x) {
    for(int k = 1; x * k <= R; k ++)
        add(x * k, (ll) mu[k] * sig[x] % mo);
}

int query(int n, int m) {
    int ans = 0, la = 0, nw;
    for(int i = 1, j; i <= n; i = j + 1, la = nw) {
        j = min(n / (n / i), m / (m / i)); nw = qry(j);
        ans = (ans + ((ll) nw - la) % mo * (n / i) % mo * (m / i) % mo) % mo;
    }
    return ((ll)ans + mo) % mo;
}

int main() {
    sieve(R); scanf("%d", &t);
    for(int i = 1; i <= t; i ++) {
        scanf("%d%d%d", &q[i].n, &q[i].m, &q[i].a);
        if(q[i].n > q[i].m) swap(q[i].n, q[i].m);
        q[i].id = i;
    }
    sort(q + 1, q + t + 1);
    int j = 1;
    for(int i = 1; i <= t; i ++) {
        for(; j <= R && sig[a[j].x] <= q[i].a; j ++) ins(a[j].x);
        ans[q[i].id] = query(q[i].n, q[i].m);
    }
    for(int i = 1; i <= t; i ++)
        printf("%d\n", ans[i]);
    return 0;
}