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

题意

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

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

(T leq 2 imes 10^4,1 leq n,mleq 10^5).

题解

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

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

然后使用常见套路变换式子：（下面默认(nleq 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}^{lfloorfrac{n}{d} floor} sum_{j=1}^{lfloorfrac{m}{d} floor}[gcd(i,j)=1] ]

[=sum_{d=1}^n sigma(d)[sigma(d)leq a] sum_{i=1}^{lfloorfrac{n}{d} floor} sum_{j=1}^{lfloorfrac{m}{d} floor} sum_{d'|i,d'j} mu(d') ]

[=sum_{d=1}^n sigma(d)[sigma(d)leq a]sum_{d'=1}^{lfloorfrac{n}{d} floor}mu(d') lfloorfrac{n}{dd'} floor lfloorfrac{m}{dd'} floor ]

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

[=sum_{i=1}^n f(d)sum_{d'=1}^{lfloorfrac{n}{d} floor}mu(d') lfloorfrac{n}{dd'} floor lfloorfrac{m}{dd'} floor ]

枚举乘积(k=dd')：

[=sum_{k=1}^nlfloorfrac{n}{k} floor lfloorfrac{m}{k} floorsum_{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} floor lfloor frac{m}{k} floor 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+Tsqrt 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
", ans[i]);
    return 0;
}