抱歉，您的浏览器无法访问本站

本页面需要浏览器支持（启用）JavaScript


了解详情 >

七つの海

今日、海を見た。もう怖くない

博客

随笔

朋友

归档

关于

简单数学：从约数个数和到莫比乌斯反演

七海

算法数学

发布于：2021年8月5日

次浏览

差不多算是基础的程度的知识的数论能力，就算已经不参加弹幕神乐的巫女想必也是应该要掌握的吧。

首先来点题目链接：P3327【SDOI2015】约数个数和 - 洛谷 | 2185. 「SDOI2015」约数个数和 - LibreOJ

\(\text d(x)\) 是 \(x\) 的约数个数，给定 \(N, M \le 50000\)，求： \[ \sum_{i=1}^N\sum_{j=1}^M \text{d}(ij) \]

莫比乌斯反演

这个题要用到它，但是我还不会。所以先随便学学吧！

前置名词

数论函数

定义域为正整数的函数称为数论函数。

积性函数

若 \(f(x)\) 是数论函数，且对于 \(\forall a,b \and (a,b)=1\) 有 \(f(ab)=f(a)f(b)\)，那么这样的数论函数称为积性函数。

常见的积性函数有欧拉函数（\(\varphi\)）、莫比乌斯函数（\(\mu\)）和除数函数（\(\sigma_k\)）

完全积性函数

若积性函数不需要 \((a,b)=1\) 也能有 \(f(ab)=f(a)f(b)\)，那样的 \(f\) 就是完全积性函数。

常见的完全积性函数：\(f(x)=x\)、\(1(x)=1\)

Dirichlet 卷积

两个数论函数 \(f,g\) 的 Dirichlet 卷积的定义如下： \[ (f*g)(n)=\sum_{d\mid n}f(d)g(\frac nd) \] 卷积运算符 \(*\) 满足运算的交换律结合律以及对于 \(+\) 的分配律；此外，对于满足 \(h(1)\ne0\) 的数论函数 \(h\)，卷积运算符还满足了等式的性质，即：\(f = g \Leftrightarrow f*h=g*h\)。

Dirichlet 卷积具有单位元，其单位元 \(e\) 定义如下： \[ \varepsilon:e(n)=\left\{\begin{matrix} 1\; (n=1)\\ 0\; (n\neq 1) \end{matrix}\right.=[n=1] \]

这个单位元也被称为单位函数，它是一个完全积性函数。

结论 1：若 \(f, g\) 是积性函数，那么 \(f*g\) 也是积性函数。

结论 2：积性函数的逆元（对于单位元）也是积性函数。

莫比乌斯函数

莫比乌斯函数是积性函数。

定义

记作 \(\mu\)；对于一个整数 \(n\)，令其标准分解形式为 \(\prod p_i^{a_i},\ i\in[1, k]\)，则莫比乌斯函数可以如下定义： \[ \mu(n)=\begin{cases} 0 &,\ \exists i: a_i > 1 \\ (-1)^k &,\ \text{else} \end{cases} \]

简单地说：如果 \(n\) 有平方因子，那么 \(\mu(n)=0\)；否则是 \((-1)^k\)，\(k\) 是 \(n\) 互不相同的质因子个数。

性质

莫比乌斯函数最重要的性质是： \[ \sum_{d\mid n}\mu(d) = \begin{cases} 1 &,\ n=1 \\ 0 &,\ n\ne1 \end{cases} = e(n) \] 用 Dirichlet 卷积的形式表示，就是 \(\mu * 1 = e\).

这个性质的证明

令 \(n\) 有标准分解形式为 \(\prod p_i^{a_i},\ i\in[1, k]\)，定义 \(n' = \prod p_i\)；那么显然有： \[ \sum_{d\mid n}\mu(d) = 0 + \sum_{d\mid n'}\mu(d) \] 而 \(n'\) 中是没有重复的因子的，所以 \(d\) 只需要在这 \(k\) 个本质不同的因子中随意选取组合即可： \[ \sum_{d\mid n'}\mu(d) = \sum_{i=0}^k \mathbf C_k^i\cdot(-1)^i \] 为每一项增加一个 \(1^{k-i}\)，由二项式定理可得： \[ \sum_{i=0}^k \mathbf C_k^i\cdot(-1)^i\cdot1^{k-i} = (-1 + 1)^k = 0^k \] 显然，当且仅当 \(k=0\) 时上式取值 \(1\)；其他情况下均为 \(0\)。而 \(k=0\) 时，\(n\) 没有任何素因子，故 \(n=1\).

除此之外，我们还有一个结论可以帮助我们将莫比乌斯函数和欧拉函数联系起来： \[ \sum_{d\mid n}\frac{\mu(d)}d = \frac{\varphi(n)}n \] 这个待会再想办法证（

线性筛求解

根据定义，稍微修改一下线性筛，我们可以写出下面的代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
template <int n> auto &mobius_sieve() {
    vector<int> prime;
    static array<int, n + 1> mu;
    bitset<n + 1> vis;
    mu[1] = 1;
    for (int i = 2; i <= n; ++ i) {
        if (!vis[i]) prime.push_back(i), mu[i] = -1;
        for (auto & pp : prime) {
            auto ipp = (llong)i * pp;
            if (ipp > n) break;
            vis[ipp] = true;
            if (i % pp == 0) { 
                mu[ipp] = 0; 
                break; 
            } else mu[ipp] = -mu[i];
        }
    }
    return mu;
}

事实上，线性筛基本可以求所有的积性函数。

莫比乌斯反演

对于一些函数 \(f(n)\)，若它本身难求但其倍数/约数和好求，那么就可以用莫比乌斯反演来简化其运算。

公式

设有数论函数 \(f,g\)，那么有如下公式： \[ \begin{align} f(n)=\sum_{d\mid n}g(d) &\iff g(n) = \sum_{d\mid n}\mu(d)f(\frac nd) \\ f(n)=\sum_{n\mid d}g(d) &\iff g(n) = \sum_{n\mid d}\mu(\frac dn)f(d) \end{align} \] 形式一是标准形式。

证明

我们先利用卷积知识简易证明上述的形式一：

首先，将上式左边看作 \(f = g*1\)；然后，有 Dirichlet 卷积的运算性质，我们在等式两侧卷积 \(\mu\)，有： \[ f * \mu = g * 1 * \mu \] 又因为莫比乌斯函数的性质有 \(e = 1 * \mu\)，所以上式的右侧消去单位元，得：\(g = f * \mu\)，即上式右侧。

或者通过数论变换的方式证明形式一：

数论变换就是反着推导的过程；对于上述关系右侧等式的右侧，代入关系左侧的条件可得： \[ \sum_{d\mid n}\mu(d)f(\frac nd) = \sum_{d\mid n}\mu(d)\sum_{k\mid \frac nd}g(k) = \sum_{k\mid n}g(k)\sum_{d\mid \frac nk}\mu(d) \] 因为 \(dk\) 是 \(n\) 的因数，\(d\) 和 \(k\) 是对于 \(dk\) 的进一步划分，这里交换求和顺序可以得到上式最右侧的形态；观察这个式子的右侧，利用莫比乌斯函数的主要性质，可以将其转化为单位函数的形式： \[ \sum_{k\mid n}g(k)\sum_{d\mid \frac nk}\mu(d) = \sum_{k\mid n}g(k)[\frac nk = 1] = g(n) \] 也就是左侧的求和只有在 \(k = n\) 时取值 \(g(n)\)，结果也是 \(g(n)\)，和形式一右侧等式的左侧一致。

我们也可以如法炮制的证明形式二：

将可被 \(n\) 整除的 \(d\) 表示成 \(kn\)，将形式二的前提代入关系右侧的等式的右侧可得： \[ \begin{align} g(n) &= \sum_{n\mid d}\mu(\frac dn)f(d) = \sum_{k=1}^{+\infty}\mu(k)f(kn)\\ &= \sum_{k=1}^{+\infty}\mu(k)\sum_{kn\mid q}g(q) = \sum_{n\mid q}g(q)\sum_{k\mid \frac qn}\mu(k) \end{align} \] 然后再进行如法炮制的交换求和顺序：\(kn\) 是无穷的 \(q\) 的因数，再进一步划分 \(kn\)；再次观察最右侧的式子，并转化为单位函数： \[ \sum_{n\mid q}g(q)\sum_{k\mid \frac qn}\mu(k) = \sum_{n\mid q}g(q)[\frac qn=1] = g(n) \] 显然，整个求和式子只有在 \(q = n\) 时才能取到值，且此时的值是 \(g(n)\)，和形式二右侧等式的左侧一致。

不知道看到这里，是否对于”莫比乌斯函数是一个和容斥系数相关的函数“这句话有了什么新的理解。

应用

但是一般来说构造一个 \(f(n) = \sum_{d\mid n}g(d)\) 颇有难度，那个公式很多时候都意义不明。所以通常的做法都是想办法整出一个 \([\gcd(i, j)=1]\) 也就是 \(e(\gcd(i,j))\)，然后通过 \(\sum_{d\mid \gcd(i,j)}\mu(d)\) 来计算它；而这实际上就是莫比乌斯函数的性质——这样说也是十分意义不明，所以我们说一类相对比较常见的问题作为例子：

给定 \(n\le m\)，求解： \[ \sum_{i=1}^n\sum_{j=1}^m f(\gcd(i, j)) \]

我们考虑枚举 \(d = \gcd(i, j)\) 的取值，于是有了： \[ \sum_{i=1}^n\sum_{j=1}^m f(\gcd(i, j)) = \sum_{d=1}^n \sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor} f(d)e(\gcd(i,j)) \] 那么我们可以对于单位函数的部分代入莫比乌斯函数的性质——有： \[ \sum_{d=1}^{\min(n,m)} \sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor} f(d)[\gcd(i,j)=1] = \sum_{d=1}^n \sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor} f(d)\sum_{r\mid\gcd(i,j)}\mu(r) \] 对右侧式子进行反演——或者说调换求和顺序，去枚举 \(r\) 的值并将其挪到外层，有： \[ \sum_{d=1}^n \sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}f(d)\sum_{r\mid\gcd(i,j)}\mu(r) = \sum_{d=1}^n \sum_{r=1}^{\min(\lfloor\frac nd\rfloor,\lfloor\frac md\rfloor)}\mu(r)\sum_{i=1}^{\lfloor\frac n{dr}\rfloor}\sum_{j=1}^{\lfloor\frac m{dr}\rfloor}f(d) \] 可以注意到这个时候右侧式子中的 \(i,j\) 已经不会影响到所需要求和的东西了，相当于对 1 求和： \[ \sum_{d=1}^n \sum_{r=1}^{\lfloor\frac nd\rfloor}\mu(r)f(d)\sum_{i=1}^{\lfloor\frac n{dr}\rfloor}\sum_{j=1}^{\lfloor\frac m{dr}\rfloor}1 = \sum_{d=1}^n \sum_{r=1}^{\lfloor\frac nd\rfloor}\mu(r)f(d)\lfloor\frac n{dr}\rfloor\lfloor\frac m{dr}\rfloor \] 此时，我们假设一个 \(g(t)\)： \[ g(t) = \sum_{d\mid t}f(d)\mu(\frac td) = f * \mu \] 观察题设的公式和我们的推到结果： \[ \sum_{i=1}^n\sum_{j=1}^m f(\gcd(i, j)) = \sum_{d=1}^n \sum_{r=1}^{\lfloor\frac nd\rfloor}\mu(r)f(d)\lfloor\frac n{dr}\rfloor\lfloor\frac m{dr}\rfloor \] 我们把 \(dr\) 看作一个整体，它是 \(n\) 的因数，\(d,r\) 是对其进一步的划分；那么令 \(t=dr\)： \[ \sum_{i=1}^n\sum_{j=1}^m f(\gcd(i, j)) = \sum_{t=1}^ng(t)\lfloor\frac nt\rfloor\lfloor\frac mt\rfloor \]

至此，我们完成了题设公式的转换，即反演；但是这样做的意义，是什么呢？当然是推出来的这个式子相对比较好算了！一重求和不比二重求和容易？

钢筋（bushi

上面的部分主要说的是对于含有 \(\gcd\) 的式子的处理方法：弄出一个 \([\gcd(i,j)=1]\) 然后再利用它的等价式子 \(\sum_{d\mid \gcd(i,j)}\mu(d)\) 代换它，并求出了一个比较通用的“公式”；那么对于这种出现了 \(\gcd\) 的式子，我就是想要使用反演公式设函数套怎么办呢？套路在此：

首先，我们再写一遍莫比乌斯反演公式的某种形态： \[ f(n)=\sum_{n\mid d}g(d) \iff g(n) = \sum_{n\mid d}\mu(\frac dn)f(d) \] 一般来说，我们设 \(g(d)\) 为范围内满足 \(\gcd(i, j) = d\) 的数对个数，\(f(n)\) 为满足 \(n \mid \gcd(i, j)\) 的数对个数，那么它们就满足了： \[ \begin{align} g(d)&=\sum_{i=1}^N\sum_{j=1}^M[\gcd(i, j)=d]\\ f(n)&=\sum_{n\mid d}g(d)=\lfloor\frac Nn\rfloor\lfloor\frac Mn\rfloor \end{align} \] 关于 \(f(n)\) 的两种表达式：第一个表达式是利用我们定义的 \(g(d)\) 来定义；第二个表达式是根据我们的定义直接得到的——当 \(i,j\) 都有确定范围的时候，满足 \(\gcd(i, j) = kn\) 的数对数量当然是这个。

这样，我们就可以利用前面写的那个莫比乌斯反演公式，得到： \[ g(n) = \sum_{n\mid d}\mu(\lfloor\frac dn\rfloor)f(d) \] 当然，裸求 \(g(d)\) 是不好求的；但是 \(\mu\) 可以线性筛加维护前缀和，\(f(n)\) 可以数论分块；于是我们就利用反演公式将不太好求的 \(g(d)\) 简化的好求了一些。

例题

上面的简介仍然是比较抽象，所以还是看几个题：

【HAOI2011】Problem b

求 \(\sum_{i=a}^b \sum_{j=c}^d [\gcd(i,j)=k]\)

首先原式的求和有区间，所以很显然地将他转化为一个二维差分的形式： \[ \begin{align} &记\ A_{n,m}=\sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j)=k]\\ &则\ 原式 = A_{b,d} - A_{b,c-1} - A_{a-1,d} + A_{a-1,c-1} \end{align} \] 那么原来的问题转化成求 \(A_{n,m}\)：

参考上面运用小节推的那个比较具有代表性的”公式“——这里 \(f(x)=[x=k]\)；那么我们可以根据”公式“来套路的得到： \[ g(t) = \sum_{d\mid t}f(d)\mu(\frac td) \] 然后，将它代回那个公式里，可以得到： \[ \begin{align} &\sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j)=k]\\ =&\sum_{t=1}^{\min(n,m)}\sum_{d\mid t}[d=k]\mu(\frac td)\lfloor\frac nt\rfloor\lfloor\frac mt\rfloor\\ =&\sum_{t=1,k\mid t}^{\min(n,m)}\mu(\frac tk)\lfloor\frac nt\rfloor\lfloor\frac mt\rfloor\\ =&\sum_{d=1}^{\lfloor\frac {\min(n,m)}k\rfloor}\mu(d)\lfloor\frac n{dk}\rfloor\lfloor\frac m{dk}\rfloor \end{align} \] 现在这个公式已经可以通过 \(\mathcal O(n)\) 来计算了；但是注意到这里待求和的式子还可以使用数论分块计算，所以实际上上式的时间复杂度是 \(\mathcal O(\sqrt n)\) 的。如果不会数论分块可以看这个：数论分块入门 - 七海の参考書 (shiraha.cn)

然后，我们就能写出代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
constexpr int N = 5e4 + 50;
int mu_sum[N + 1], a, b, c, d, k;

int mu_seq(int l, int r) {
    return mu_sum[r] - mu_sum[l - 1];
}

template <int n> auto &mobius_sieve() {
    vector<int> prime;
    static array<int, n + 1> mu;
    bitset<n + 1> vis;
    mu[1] = 1;
    for (int i = 2; i <= n; ++ i) {
        if (!vis[i]) prime.push_back(i), mu[i] = -1;
        for (auto & pp : prime) {
            auto ipp = (llong)i * pp;
            if (ipp > n) break;
            vis[ipp] = true;
            if (i % pp == 0) {
                mu[ipp] = 0;
                break;
            } else mu[ipp] = -mu[i];
        }
    }
    return mu;
}

int ybb(int n, int m) {
    if (n > m) swap(n, m);
    int ret = 0, nk = n / k;
    for (int l = 1, r; l <= nk; l = r + 1) {
        int nlk = n / l / k, mlk = m / l / k;
        r = min({nk, n / (n / l), m / (m / l)});
        ret += mu_seq(l, r) * nlk * mlk;
    }
    return ret;
}

signed main() {
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif
    const auto mu = mobius_sieve<N>();
    for (int i = 1; i <= N; ++ i)
        mu_sum[i] = mu_sum[i - 1] + mu[i];
    for (auto T = read(); T --;) {
        a = read(), b = read(), c = read();
        d = read(), k = read();
        auto ans = ybb(b, d) - ybb(b, c - 1) - ybb(a - 1, d) + ybb(a - 1, c - 1);
        cout << ans << endl;
    }
    return 0;
}

遇事不决开 long long 是吧？long long 不是你的电子宠物（

YY的GCD

求 \(\sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j)\ \text{is}\ \text{prime}]\)

令 \(n\le m\)，还是很套路地把上面的式子变成枚举 \(\gcd\) 的值： \[ 原式 = \sum_{k=1}^{\min(n,m)}\sum_{i=1}^n \sum_{j=1}^m [\gcd(i, j)=k] = \sum_{k=1}^{n}\sum_{i=1}^{\lfloor\frac nk\rfloor} \sum_{j=1}^{\lfloor\frac mk\rfloor} [\gcd(i, j)=1],\ k\ \text{is}\ \text{prime} \] 还是运用“应用”部分得到的公式；我们轻而易举地发现 \(f=e\)，那么套路地设 \(g\)： \[ g(t) = f * \mu = \sum_{d\mid t}e(d)\mu(\frac td) \] 然后还是代回原来的那个包含 \(f\) 的表达式中，可以得到： \[ \begin{align} &\sum_{k=1}^n\sum_{i=1}^{\lfloor\frac nk\rfloor} \sum_{j=1}^{\lfloor\frac mk\rfloor} e(\gcd(i, j)),\ \ k\ \text{is}\ \text{prime}.\\ =&\sum_{k=1}^n\sum_{t=1}^{\lfloor\frac nk\rfloor}\sum_{d\mid t}e(d)\mu(\frac td)\lfloor\frac n{tk}\rfloor\lfloor\frac m{tk}\rfloor\\ =&\sum_{k=1}^n\sum_{t=1}^{\lfloor\frac nk\rfloor}\mu(t)\lfloor\frac n{tk}\rfloor\lfloor\frac m{tk}\rfloor \end{align} \] 令 \(x = tk\)，那么： \[ 上式 = \sum_{x=1}^n\sum_{k\mid x}\mu(\frac xk)\lfloor\frac nx\rfloor\lfloor\frac mx\rfloor = \sum_{x=1}^n\lfloor\frac nx\rfloor\lfloor\frac mx\rfloor\sum_{k\mid x}\mu(\frac xk),\ \ k\ \text{is}\ \text{prime}. \] 那么式子就推到这里；上面的式子左边可以数论分块，右边那个东西虽然看起来玄乎但是总归还是可以预处理的：只需要对于所有的质数在范围内的倍数“对数筛”即可，复杂度不明不会素数定理。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
constexpr int N = 1e7 + 50;
int f[N + 1], f_sum[N + 1];

int f_sec(llong l, llong r) {
    return f_sum[r] - f_sum[l - 1];
}

template <int n> auto mobius_sieve() {
    static vector<int> prime;
    static array<int, n + 1> mu;
    bitset<n + 1> vis;
    mu[1] = 1;
    for (int i = 2; i <= n; ++ i) {
        if (!vis[i]) prime.push_back(i), mu[i] = -1;
        for (auto & pp : prime) {
            auto ipp = (llong)i * pp;
            if (ipp > n) break;
            vis[ipp] = true;
            if (i % pp == 0) {
                mu[ipp] = 0;
                break;
            } else mu[ipp] = -mu[i];
        }
    }
    return make_pair(ref(prime), ref(mu));
}

llong fuck(int n, int m) {
    if (n > m) swap(n, m);
    llong ret = 0;
    for (llong l = 1, r; l <= n; l = r + 1) {
        r = min({n / (n / l), m / (m / l), (llong) n});
        ret += (llong) f_sec(l, r) * (n / l) * (m / l);
    }
    return ret;
}

signed main() {
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif
    const auto [prime, mu] = mobius_sieve<N>();
    for (auto pp : prime)
        for (int i = 1;; ++ i) {
            auto ipp = (llong)i * pp;
            if (ipp > N) break;
            else f[ipp] += mu[i];
        }
    for (int i = 1; i <= N; ++ i)
        f_sum[i] = f_sum[i - 1] + f[i];
    for (auto T = read(); T --;) {
        auto n = read(), m = read();
        cout << fuck(n, m) << endl;
    }
    return 0;
}

所谓条条大路通罗马（罗马！），直接用莫比乌斯反演公式也能推出一样的式子，看各自喜好了（

思路

那么，你已经学会了莫比乌斯反演了，快上！（指做开篇的那个题

关于 d(x)

首先关于这个都不太熟的 \(\text d(x)\)，我们有一个结论： \[ \text d(NM) = \sum_{i\mid N}\sum_{j\mid M}[\gcd(i, j)=1] \] 这个式子为什么是正确的？首先考虑对于一个整数 \(X\) 的标准分解，即 \(X = \prod p_i^{a_i}\)，约数的个数为 \(\text d(X) = \prod (a_i + 1)\)；这非常的好理解，它的约数必定由它的质因子构成，每个质因子 \(p_i\) 有 \(a_i + 1\) 种不同的选法。合数因子本身就是对于这个数字所有的质因子进行这样的选择组合而成；对于 \(NM\)，我们也可以从这个角度入手考虑：

对于一个素数 \(p\)，有 \(N = n\times p^x\) 和 \(M=m\times p^y\)；那么显然 \(p\) 在 \(NM\) 中出现的次数是 \(x + y\) 次，关于这个质因子有 \(x + y + 1\) 种选法。那么怎么枚举选法呢？显然，若枚举 \(i\mid N\) 和 \(j\mid M\)，得到的 \(ij \mid NM\)。那么问题就在于两次不同的枚举得到的 \(ij\) 的乘积可能实际上是一样的；为了避免重复，我们定义下面的取法：

  • 令对于 \(NM\) 和素因子 \(p\)，我们要选取其中的 \(k\) 次，那么这 \(k\) 个 \(p\) 就要由 \(N\) 和 \(M\) 来提供。
  • 当 \(k \le x\) 时，我们要求 \(p\) 完全由 \(N\) 提供；此时可以选出 \(k = x + 1\) 种不同的有序数对 \((p^k, 1)\)
  • 当 \(k > x\) 时，我们要求超过 \(x\) 的部分由 \(M\) 提供，但不再在 \(N\) 中选择；此时有 \(y\) 种有序数对 \((1,p^{k-x})\)

因此对于每一个因子，在上述规则的限制下，在 \(NM\) 中只会选出 \(x + y + 1\) 种不同的选法，符合我们的要求；而实现这样的限制条件很显然可以通过 \(\gcd(i, j)=1\) 来达成目的，因为我们从来没有在 \(i\) 和 \(j\) 中同时选择 \(p\).

变换

综上所述，我们要求的式子可以写成下面的形式： \[ \sum_{n=1}^N\sum_{m=1}^M\text d(nm) = \sum_{n=1}^N\sum_{m=1}^M\sum_{i\mid n}\sum_{j\mid m}[\gcd(i, j)=1] \] 改变四层求和的枚举顺序，先枚举 \(i\) 和 \(j\)，那么可以得到下面的形式： \[ \sum_{n=1}^N\sum_{m=1}^M\text d(nm) = \sum_{i=1}^N\sum_{j=1}^M\lfloor\frac Ni\rfloor\lfloor \frac Mj\rfloor[\gcd(i,j)=1] \] 那么就转化了题设的公式，可以基于这个公式进行莫比乌斯反演了。

反演

根据上面的讨论，我们已经得到了一个包含了 \(e(\gcd(i, j))\) 形态的公式，现在对它进行操作：

根据 Dirichlet 卷积的单位元的性质，也就是 \(e = 1*\mu\)，得： \[ \sum_{i=1}^N\sum_{j=1}^M\lfloor\frac Ni\rfloor\lfloor \frac Mj\rfloor[\gcd(i,j)=1] = \sum_{i=1}^N\sum_{j=1}^M\lfloor\frac Ni\rfloor\lfloor \frac Mj\rfloor\sum_{d\mid\gcd(i, j)}\mu(d) \] 对于最右边这个子式，我们很自然地想到枚举 \(d\)；令 \(\min(N,M)=N\)，调换枚举顺序： \[ \sum_{i=1}^N\sum_{j=1}^M\lfloor\frac Ni\rfloor\lfloor \frac Mj\rfloor\sum_{d=1}^{\min(N,M)}[d\mid\gcd(i,j)]\mu(d) = \sum_{d=1}^{N}\mu(d)\sum_{i=1}^N\sum_{j=1}^M\lfloor\frac Ni\rfloor\lfloor \frac Mj\rfloor[d\mid\gcd(i,j)] \] 右边的双层枚举子式又是典型的求 \(\gcd\) 是倍数的类型，进行套路地转换： \[ \sum_{d=1}^{N}\mu(d)\sum_{i=1}^{\lfloor\frac Nd\rfloor}\sum_{j=1}^{\lfloor\frac Md\rfloor}\lfloor\frac N{id}\rfloor\lfloor \frac M{jd}\rfloor = \sum_{d=1}^{N}\mu(d)\sum_{i=1}^{\lfloor\frac Nd\rfloor}\lfloor\frac N{id}\rfloor\sum_{j=1}^{\lfloor\frac Md\rfloor}\lfloor\frac M{jd}\rfloor \] 综上所述，我们通过反演将原式转化为了： \[ \sum_{n=1}^N\sum_{m=1}^M\text d(nm) = \sum_{d=1}^{\min(N,M)}\mu(d)(\sum_{i=1}^{\lfloor\frac Nd\rfloor}\lfloor\frac N{id}\rfloor\cdot\sum_{j=1}^{\lfloor\frac Md\rfloor}\lfloor\frac M{jd}\rfloor) \] 那么这个化简后的式子要怎么去求呢？容易发现两个子求和是近乎一致的，可以预先处理；因此我们定义函数 \(h(n)=\sum_{i=1}^n\lfloor\frac ni\rfloor\)，于是有： \[ \sum_{n=1}^N\sum_{m=1}^M\text d(nm) = \sum_{d=1}^{\min(N,M)}\mu(d)h(\lfloor\frac Nd\rfloor)h(\lfloor\frac Md\rfloor) \] 而实际上，\(h(n)\) 就是约数个数函数 \(\text d(x)\) 的前缀和，这个也可以用线性筛求出来维护前缀和；当然如果不想使用线性筛来维护这个，也可以直接分块计算后加起来——复杂度是 \(\mathcal O(n\sqrt n)\)；

代码实现

综上所述，我们可以用分块和分块求解上面反演得到的式子：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
constexpr int N = 5e4 + 50;
int h[N + 1], mu_sum[N + 1];

int mu_sec(llong l, llong r) {
    return mu_sum[r] - mu_sum[l - 1];
}

template <int n> auto mobius_sieve() {
    static vector<int> prime;
    static array<int, n + 1> mu;
    bitset<n + 1> vis;
    mu[1] = 1;
    for (int i = 2; i <= n; ++ i) {
        if (!vis[i]) prime.push_back(i), mu[i] = -1;
        for (auto & pp : prime) {
            auto ipp = (llong)i * pp;
            if (ipp > n) break;
            vis[ipp] = true;
            if (i % pp == 0) {
                mu[ipp] = 0;
                break;
            } else mu[ipp] = -mu[i];
        }
    }
    return make_pair(ref(prime), ref(mu));
}

llong partition(llong n) {
    llong ret = 0;
    for (llong l = 1, r; l <= n; l = r + 1) {
        r = min(n / (n / l), n);
        ret += (r - l + 1) * (n / l);
    }
    return ret;
}

llong fuck(int n, int m) {
    if (n > m) swap(n, m);
    llong ret = 0;
    for (llong l = 1, r; l <= n; l = r + 1) {
        r = min({n / (n / l), m / (m / l), (llong)n});
        ret += (llong) mu_sec(l, r) * h[n / l] * h[m / l];
    }
    return ret;
}

signed main() {
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif
    const auto [prime, mu] = mobius_sieve<N>();
    for (int i = 1; i <= N; ++ i)
        mu_sum[i] = mu_sum[i - 1] + mu[i];
    for (int i = 1; i <= N; ++ i)
        h[i] = (int) partition(i);
    for (auto T = read(); T --;) {
        auto n = read(), m = read();
        cout << fuck(n, m) << endl;
    }
    return 0;
}

如果用线性筛维护 \(\text d(x)\) 再求其前缀和 \(h(x)\)，则需要引入新的线性筛：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
constexpr int N = 5e4 + 50;
int h[N + 1], mu_sum[N + 1];

int mu_sec(llong l, llong r) {
    return mu_sum[r] - mu_sum[l - 1];
}

template <int n> auto mobius_sieve() {
    static vector<int> prime;
    static array<int, n + 1> mu;
    bitset<n + 1> vis;
    mu[1] = 1;
    for (int i = 2; i <= n; ++ i) {
        if (!vis[i]) prime.push_back(i), mu[i] = -1;
        for (auto & pp : prime) {
            auto ipp = (llong)i * pp;
            if (ipp > n) break;
            vis[ipp] = true;
            if (i % pp == 0) {
                mu[ipp] = 0;
                break;
            } else mu[ipp] = -mu[i];
        }
    }
    return make_pair(ref(prime), ref(mu));
}

template <int n> auto divisor_sieve() {
    static vector<int> prime;
    static array<int, n + 1> d, c;
    bitset<n + 1> vis;
    d[1] = 1, c[1] = 1;
    for (int i = 2; i <= n; ++ i) {
        if (!vis[i]) {
            prime.push_back(i);
            d[i] = 2, c[i] = 1;
        }
        for (auto & pp : prime) {
            auto ipp = (llong)i * pp;
            if (ipp > n) break;
            else vis[ipp] = true;
            if (i % pp == 0) {
                c[ipp] = c[i] + 1;
                d[ipp] = d[i] / c[ipp] * (c[ipp] + 1);
                break;
            } else c[ipp] = 1, d[ipp] = d[i] * 2;
        }
    }
    return make_pair(ref(prime), ref(d));
}

llong fuck(int n, int m) {
    if (n > m) swap(n, m);
    llong ret = 0;
    for (llong l = 1, r; l <= n; l = r + 1) {
        r = min({n / (n / l), m / (m / l), (llong)n});
        ret += (llong) mu_sec(l, r) * h[n / l] * h[m / l];
    }
    return ret;
}

signed main() {
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif
    const auto [prime, mu] = mobius_sieve<N>();
    for (int i = 1; i <= N; ++ i)
        mu_sum[i] = mu_sum[i - 1] + mu[i];
    const auto &d = divisor_sieve<N>().second;
    for (int i = 1; i <= N; ++ i)
        h[i] = h[i - 1] + d[i];
    for (auto T = read(); T --;) {
        auto n = read(), m = read();
        cout << fuck(n, m) << endl;
    }
    return 0;
}

当然，肯定是要把两个线性筛写在一起的；像上面那样写的人脑子多半是有点大病（

于是，这个问题就解决了！

后记

莫比乌斯反演还是比较有趣的；这里列举的也仅仅是最基础最基础的板子题，用来加深对这个公式的推导以及这种方法的理解——也就是说有趣的题还有很多……之后有时间做了再整理一篇吧（

参考资料

博客内容遵循 署名-非商业性使用-相同方式共享 4.0 国际 (CC BY-NC-SA 4.0) 协议

本文永久链接是：https://shiraha.cn/2021/sdoi-2015-the-summation-of-divisors/

更新于：2021年8月8日

算法

ACM

数论

莫比乌斯反演

评论