bzoj [SDOI2014]数表 莫比乌斯反演 BIT

Posted 2021-11-27 dsrdsr

bzoj [SDOI2014]数表 莫比乌斯反演 BIT

链接

bzoj
luogu
loj

思路

\[ \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}a*[f[gcd(i,j)]<=a] \]

\[ f[]可以O(n)预处理出来 \]

\[ \sum\limits_{k=1}^{n}f[k]*\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{m}[gcd(i,j)==k] \]

\[ \sum\limits_{k=1}^{n}f[k]*\sum\limits_{i=1}^{\frac{n}{k}}{\frac{n}{ki}}{\frac{m}{ki}}\mu(i) \]

\[ \sum\limits_{i=1}^{n}f[i]*\sum\limits_{k|d}{\frac{n}{d}\frac{m}{d}}\mu(\frac{d}{i}) \]

d替换k*i
\[ \sum\limits_{d=1}^{n} \frac{n}{d} \frac{m}{d} \sum\limits_{k|d} \mu(\frac{d}{k})f(k) \]
\[ \sum\limits_{d=1}^{n} \frac{n}{d} \frac{m}{d} g(d) \]

\[ g(d)=\sum\limits_{k|d} \mu(\frac{d}{k})f(k) \]

询问按照a排序，每次加入f(k)时候影响的只是能k|d的g(d)
每次修改就是\(O(\sqrt{n}logn)\)
查询也是一样，
总的就是\(O(T\sqrt{n}logn)\)
注意，ll+取模的话，loj会超时，用int的自然溢出就快了三倍(300ms)，是int，不是unsigned int。

其他

线性筛约数和

\[ x=p_{1}^{w_1}p_{2}^{w_2}…p_{k}^{w_k} \]

那么

\[ SD(x)=约数和=(1+p_1^1+p_1^2+…+p_1^{w_1})(1+p_2^1+p_2^2+…+p_2^{w_2})(1+p_k^1+p_k^2+…+p_k^{w_k}) \]

0x00 是个素数

显然\(SD(pri)=pri+1\)

0x01 两两互质

是个积性函数
因为x,y两两互质，所以他们质因子互不相交，所以显然脑补公式

\[ SD(x*y)=SD(x)*SD(y)(gcd(x,y)==1) \]

0x02 两两不互质（i%pri[j]!=0）

再开个数组tmp，记录最小质因子因子的贡献\((1+p_1^1+p_1^2+…+p_1^{w_1})\)
因为pri[j]是他的最小质因子（因为这是线性筛）
我们之前求出的i的

\[ SD(i)=(1+p_1^1+p_1^2+…+p_1^{w_1})(1+p_2^1+p_2^2+…+p_2^{w_2})(1+p_k^1+p_k^2+…+p_k^{w_k}) \]

现在的\(i*pri[j]\)的SD显然就是

\[ SD(i*pri[j])=(1+p_1^1+p_1^2+…+p_1^{w_1}+p_1^{w_1+1})(1+p_2^1+p_2^2+…+p_2^{w_2})(1+p_k^1+p_k^2+…+p_k^{w_k}) \]

改变的只有最小因子的贡献，tmp的作用就来了
tmp[i]我们已经求出来了，那么
\[tmp[i*pri[j]]=tmp[i]*pri[j]+1（这是个等差数列）\]

\[ SD[i*pri[j]]=SD[i]/tmp[i]*tmp[i*pri[j]] \]

代码

#include <bits/stdc++.h>
using namespace std;
const int N = 1e5 + 6;
int read() {
    int x = 0, f = 1; char s = getchar();
    for (; s > '9' || s < '0'; s = getchar()) if (s == '-') f = -1;
    for (; s >= '0' && s <= '9'; s = getchar()) x = x * 10 + s - '0';
    return x * f;
}
int pri[N], tot, vis[N], mu[N];
int f[N], tmp[N];
struct node {
    int n, m, a, id;
    int ans;
    bool operator < (const node &b) const {
        return a < b.a;
    }
} Q[N];
bool cmp(node a, node b) {
    return a.id < b.id;
}
pair<int,int> F[N];
void Euler(int limit) {
    f[1] = tmp[1] = mu[1] = 1;
    for (int i = 2; i <= limit; ++i) {
        if (!vis[i]) {
            mu[i] = -1;
            pri[++tot] = i;
            f[i] = i + 1;
            tmp[i] = i + 1;
        }
        for (int j = 1; j <= tot && i * pri[j] <= limit; ++j) {
            vis[i * pri[j]] = 1;
            if (i % pri[j] == 0) {
                tmp[i * pri[j]] = tmp[i] * pri[j] + 1;
                f[i * pri[j]] =  f[i] / tmp[i] * tmp[i * pri[j]];
                mu[i * pri[j]] = 0;
                break;
            }
            mu[i * pri[j]] = -mu[i];
            f[i * pri[j]] = f[i] * f[pri[j]];
            tmp[i * pri[j]] = pri[j] + 1;
        }
    }
    for (int i = 1; i <= limit; ++i) {
        F[i].first = f[i], F[i].second = i;
    }
    sort(F + 1, F + 1 + limit);
}
namespace BIT {
    int sum[N];
    int lowbit(int x) {return x & (-x);}
    void add(int x, int ad) {
        for (int i = x; i <= 100000; i += lowbit(i)) sum[i] = (sum[i] + ad);
    }
    int query(int x) {
        int ans = 0;
        for (int i = x; i >= 1; i -= lowbit(i)) ans = (ans + sum[i]);
        return ans;
    }
}
int main() {
    Euler(100000);
    int T = read();
    for (int i = 1; i <= T; ++i) {
        Q[i].n = read(),Q[i].m = read(),Q[i].a = read(), Q[i].id = i;
    }
    sort(Q + 1, Q + 1 + T);
    int now = 0;
    for (int i = 1; i <= T; ++i) {
        while (now + 1 <= 100000 && Q[i].a >= F[now + 1].first) {
            now++;
            for (int j = F[now].second; j <= 100000; j += F[now].second) {
                BIT::add(j, mu[j / F[now].second] * F[now].first);
            }
        }
        int ans = 0;
        if (Q[i].n > Q[i].m) swap(Q[i].n, Q[i].m);
        for (int l = 1, r; l <= Q[i].n; l = r + 1) {
            r = min(Q[i].n / (Q[i].n / l), Q[i].m / (Q[i].m / l));
            ans += 1LL * (Q[i].n / l) * (Q[i].m / l) * (BIT::query(r) - BIT::query(l - 1));
        }
        Q[i].ans = ans;
    }
    sort(Q + 1, Q + 1 + T, cmp);
    for (int i = 1; i <= T; ++i) printf("%u\n", Q[i].ans < 0 ? Q[i].ans + 2147483648 : Q[i].ans);
    return 0;
}