算法学习笔记(24): 狄利克雷卷积和莫比乌斯反演

合集 - 算法学习笔记(32)

1.算法学习笔记(1): 欧几里得算法及其扩展01-122.算法学习笔记(2): 欧拉定理与逆元01-133.算法学习笔记(3): 倍增与ST算法01-124.算法学习笔记(4): 并查集及其优化01-125.算法学习笔记(5): 最近公共祖先(LCA)01-126.算法学习笔记(6): 树链剖分01-127.算法学习笔记(7): 二分图01-138.算法学习笔记(8): 网络流01-139.算法学习笔记(8.0): 网络流前置知识01-1310.算法学习笔记(8.1): 网络最大流算法 EK, Dinic, ISAP01-1411.算法学习笔记(8.2): 上下界网络流01-1412.算法学习笔记(9): 中国剩余定理(CRT)以及其扩展(EXCRT)01-1513.算法学习笔记(10): BSGS算法及其扩展算法01-1614.算法学习笔记(11): 原根01-1615.算法学习笔记(12): 线性基02-0316.算法学习笔记(13): Manacher算法02-0717.算法学习笔记(14): 字符串哈希02-0618.算法学习笔记(15): Trie（字典树）02-0819.算法学习笔记(16): 组合数学基础02-0820.算法学习笔记(17): 快速傅里叶变换（FFT）02-1021.算法学习笔记(18): 平衡树（一）03-1022.算法学习笔记(19): 树上启发式合并（DSU on tree）03-1823.算法学习笔记(20): AC自动机03-2624.算法学习笔记(21): 平衡树（二）05-1325.算法学习笔记(22): 逆序对与原序列05-3126.算法学习笔记(23): 马尔可夫链中的期望问题06-01

27.算法学习笔记(24): 狄利克雷卷积和莫比乌斯反演07-27

28.算法学习笔记(25): 矩阵树定理06-1629.算法学习笔记(26): 计算几何07-2230.算法学习笔记(27): 后缀排序07-2631.算法学习笔记(28): 筛法07-2632.算法学习笔记(∞)：杂项08-17

收起

狄利克雷卷积和莫比乌斯反演

看了《组合数学》，再听了学长讲的……感觉三官被颠覆……

目录

• 狄利克雷卷积和莫比乌斯反演
  • 狄利克雷卷积
    • 特殊的函数
    • 函数之间的关系
      • 除数函数和幂函数
      • 欧拉函数和恒等函数
      • 莫比乌斯函数和欧拉函数
    • 狄利克雷卷积的逆元
  • 莫比乌斯函数与莫比乌斯反演
    • 求法
    • 数论分块（整除分块）
  • 莫比乌斯反演的经典结构
    • 结构1
    • 结构2
    • 结构3
    • 结构4
    • 结构总结
  • 莫比乌斯再认识
    • 二项式反演
    • 扩展到偏序集

---

狄利克雷卷积

如此定义：

(f∗g)(n)=∑xy=nf(x)g(y)

或者可以写为

(f∗g)(n)=∑d|nf(d)g(nd)

---

特殊的函数

• 单位根 ε：满足 f∗ε=ε∗f=f。

ε(n)={1,if n = 10,otherwise

• 幂函数 Idk(n)=nk。特殊的，Id1(n)=n 为恒等函数，Id0(n)=1 为常函数，简记为 I。
• 除数函数 σk(n)=∑d|ndk。特殊的，σ1(n) 为因数和函数，简记为 σ(n)，σ0(n) 为因数个数函数，简记为 τ(n)。
• 欧拉函数 φ(n)。质因数分解 n=pc11pc22...pckk，则 φ(n)=n∏ki=1pi−1pi。

这些函数都是积性函数，满足 gcd(i,j)=1⟹f(ij)=f(i)f(j)。

---

函数之间的关系

可能有不完整的地方……

其中可以通过 I 转化的，都可以通过 μ 转换回去。考虑 I∗μ=ε 的事实，后面会讲。

除数函数和幂函数

根据定义，我们有

(Idk∗I)(n)=∑d|nIdk(d)=∑d|ndk=σk(n)

即

Idk∗I=σk

欧拉函数和恒等函数

根据卷积：

(φ∗I)(n)=∑d|nφ(d)

在 n=pk 时（p 为质数），有：

∑d|nφ(d)=φ(1)+k∑i=1φ(pi)=1+k∑i=1pi−pi−1=pm=d

所以 (φ∗I)(pk)=pk

将 n 质因数分解为 ∏pk，根据积性函数的定义，可知：

(φ∗I)(n)=n=Id1(n)

莫比乌斯函数和欧拉函数

应该把后面看完了再回来看这个……

莫比乌斯函数定义为：

μ∗I=ε

根据： φ∗I=Id1，两边同时卷上 μ

有：

φ∗I∗μ=Id1∗μ⟺φ=Id1∗μ

---

狄利克雷卷积的逆元

对于一个函数 f，我们可以如下的定义一个函数 g。

首先设 g(1)=1f(1)。

然后令 g(x)=−1f(1)∑d|x,d>1g(d)f(xd)

于是 (f∗g)=ε

展开带入证明即可。

---

莫比乌斯函数与莫比乌斯反演

终于到这里了 QwQ

我们定义莫比乌斯函数是 I 的逆函数，也就是说 (μ∗I)=ε。

所以，在狄利克雷卷积中：

μ(n)={1if n=10if ∃x∃k,n=kx2(−1)mn=p1p2...pm

至于为什么强调狄利克雷卷积……后文会提及。

莫比乌斯函数常用于以下形式

g(n)=∑d|nf(d)⟺f(n)=∑d|nμ(d)g(nd)

或者可以写作：

f∗I=g⟺f=g∗μ

而这就是莫比乌斯反演的核心公式。

很简单的公式，根本无需记忆……

求法

和欧拉函数 φ 类似，可以通过线性筛的方法在 O(n) 内求出。

vector<int> prms;
int mob[N], notp[N];
void getMob(int n) {
    mob[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (!notp[i])
            mob[i] = -1, prms.push_back(i);

        for (int p : prms) {
            int ip = i * p;
            if (ip > n) break;
            notp[ip] = 1;
            if (i % p == 0) {
                mob[ip] = 0;
                break;
            } else mob[ip] = -mob[i];
        }
    }
}

---

数论分块（整除分块）

这部分虽然不属于莫比乌斯反演，但是在求很多相关问题的时候需要用到。

开篇膜拜 Pecco 大佬：# 算法学习笔记(36): 莫比乌斯反演

核心问题：求解 ∑ni=1⌊ni⌋,n≤1012。

不难得知， ⌊ni⌋ 不同的取值只有 O(√n) 个，并且是连续的。所以考虑对于每一种取值，求出有其总贡献。问题转化为，对于 i，需要求出满足 ⌊ni⌋=⌊nj⌋ 的最大的 j。

于是设 ⌊ni⌋=k

⌊nj⌋=k⟹k≤nj≤k+1⟹1k+1<jn≤1k⟹j≤nk⟹j≤⌊n⌊ni⌋⌋

也就是说，对于每一个取值 ⌊ni⌋，最大在 ⌊n⌊ni⌋⌋ 时满足。

于是可以这样写出代码：

for (int l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    ans += (r - l + 1) * (n / l);
}

练习题：[CQOI2007]余数求和 - 洛谷

---

补充：现在我们考虑更极端的情况，例如求：

f(n)=n∑i=1g(⌊ni⌋)

其中：

g(n)=n∑i=1⌊ni⌋

写出来 n=1e9 的时候可以轻松过掉。其复杂度为 O(n34)，证明见讨论。

膜拜 @jijidawang in 讨论。

所以见到类似的数论分块套数论分块，大胆的写下去吧！

---

其实数论分块不仅仅可以解决 ⌊ni⌋ 的问题，也可以套用在 √i 的问题上。写法十分类似。这里就不做过多讲解。

---

莫比乌斯反演的经典结构

莫比乌斯反演的经典套路其实就是 ε,φ,μ,I,Id 的灵活应用和转换，以及交换合式的小技巧。

其中有 ε(x)=[x=1]=∑d|xμ(d),Id(x)=x=∑d|xφ(x)。

结构1

[gcd(i,j)=1]=ε(gcd(i,j))=∑d|gcd(i,j)μ(d)

于是对于：

n∑i=1m∑j=1[gcd(i,j)=1]=n∑i=1m∑j=1∑d|gcd(i,j)μ(d)

在这里，有一个非常经典的处理方法：提取公因数

也就是枚举 gcd(i,j)

=min(n,m)∑d=1⌊ni⌋∑i=1⌊nj⌋∑j=1μ(d)=min(n,m)∑d=1⌊nd⌋⌊md⌋μ(d)

于是最终利用数论分块求即可。复杂度为 O(n+√min(n,m))。

但是代码需要注意，每一次取小步跳：

r = min(n / (n / l), m / (m / l));

结构2

原题：GCD SUM - 洛谷

gcd(i,j)=Id1(gcd(i,j))=(I∗φ)(gcd(i,j))=∑d|gcd(i,j)φ(d)

于是求 ∑ni=1∑mj=1gcd(i,j) 的方法与结构1类似即可。

结构3

τ(x)=∑k|x1τ(xy)=∑i|x∑j|y[gcd(i,j)=1]

于是求 ∑nx=1∑my=1τ(xy) 也就很简单了。

结构4

原题：YY的GCD - 洛谷

令 P 表示质数集合，求：

n∑i=1m∑j=1[gcd(i,j)∈P]

我们首先提取公因数：

=∑p∈P⌊np⌋∑i=1⌊mp⌋∑j=1[gcd(i,j)=1]

于是根据模型1：

=∑p∈P⌊min(n,m)p⌋∑x=1μ(x)⌊npx⌋⌊npx⌋

接下来是一个非常经典的套路：枚举 T=px

=min(n,m)∑T=1⌊nT⌋⌊mT⌋∑p|T,p∈Pμ(Tp)

于是问题转化为求 ∑p|T,p∈Pμ(Tp) 的前缀和，这样就可以单次询问 O(√n)。

这完全可以通过埃氏筛筛出，复杂度 O(nloglogn)，十分优秀。

不过也可以通过线性筛筛出（因为这个函数非积性函数，所以这里不展开）。

类似的题还有 [国家集训队]Crash的数字表格 / JZPTAB - 洛谷

问题：∑ni=1∑mj=1lcm(i,j)。

考虑转化为 gcd，有

lcm(i,j)=ijgcd(i,j)=gcd(i,j)×igcd(i,j)×jgcd(i,j)

带入原式中，枚举 gcd(i,j)，有：

=min(n,m)∑d=1dnd∑i=1md∑j=1ij[gcd(i,j)=1]=min(n,m)∑d=1dnd∑i=1md∑j=1ij∑t|gcd(i,j)μ(t)=min(n,m)∑d=1dmin(n,m)d∑t=1μ(t)t2ndt∑i=1imdt∑j=1j

后面部分可以通过 g(n)=n(n+1)2 简化。

=min(n,m)∑d=1dmin(n,m)d∑t=1μ(t)t2g(⌊ndt⌋)g(⌊mdt⌋)

于是只需要预处理 μ(t)t2 的前缀和即可。

但是显然数论分块套数论分块不够优秀，考虑继续优化。

所以枚举 T=dt。

=min(n,m)∑T=1g(⌊nT⌋)g(⌊mT⌋)T∑t|Tμ(t)t

也就是说只需要线性求出 ∑t|Tμ(t)t 即可（其满足积性）。

有 [SDOI2014]数表 - 洛谷

问题：∑ni=1∑mj=1σ(gcd(i,j))

不过这道题要难一些，因为涉及到了更多的操作。

主要就是将询问离线，按照限制的大小一一处理即可。

有 毒瘤之神的考验 - 洛谷

求 ∑ni=1∑mj=1φ(ij)。

需要知道转化式子：

φ(ij)=φ(i)φ(j)gcd(i,j)φ(gcd(i,j))

将 φ 展开即可证明。

于是可以得到玄妙的式子。

然后通过值域分治求解即可。

……说着简单

化简后有三个函数：

f(n)=∑d|ndμ(nd)φ(d)g(n,t)=n∑i=1φ(it)h(n,m,t)=f(t)∗g(n,t)∗g(m,t)

原式为：

min(n,m)∑t=1h(⌊nt⌋,⌊mt⌋,t)

显然，由于 h 的状态数很恨很多，所以考虑值域分治。

在 nt 比较大的时候暴力求解，到很小的时候再利用预处理的前缀和快速求解。

类似于 √n 分治吧。

于是你可以写出类下的代码：

lint solve(lint n, lint m) {
    if (n > m) swap(n, m);

    lint ans = 0;
    for (lint i = 1; i <= min(n, S); ++i)
        (ans += h(i, n / i, m / i)) %= mod;

    for (lint l = S + 1, r; l <= n; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans += t[n / l][m / l][r] - t[n / l][m / l][l - 1] + mod;
        ans %= mod;
    }

    return ans;
}

其中 S 表示分治的边界，t 表示对于 h 在第三个参数位置的前缀和。

预处理也是很简单：

inline lint h(int t, int n, int m) {
    return f[t] * g[t][n] % mod * g[t][m] % mod;
}

void init(lint n = 1e5) {
    for (lint i = 1; i <= n; ++i) {
        for (lint j = 1; i * j <= n; ++j)
            (f[i * j] += (mod + mob[j]) % mod * i % mod * iphi[i] % mod) %= mod;
    }

    for (lint t = 1; t <= n; ++t) {
        g[t].resize(n / t + 1);
        for (int i = 1; i * t <= n; ++i)
            g[t][i] = (g[t][i - 1] + phi[i * t]) % mod;    
    }

    for (int i = 1; i < B; ++i) for (int j = 1; j < B; ++j) {
        int len = n / max(i, j);
        t[i][j].resize(len + 2);
        t[i][j][0] = 0;
        for (int k = 1; k <= len; ++k)
            t[i][j][k] = (t[i][j][k - 1] + f[k] * g[k][i] % mod * g[k][j] % mod) % mod;
    }

    S = N / B; 
}

很烦写而已。

---

结构总结

在这类莫比乌斯反演中，经典的两个套路：

• 提取公因数

• 枚举 T=px

其实在 Pecco 的文章中，对于提取公因数这个方法有更加深刻的阐释。其不仅能应用在只有 i,j 两个变量的模型中，还可以扩展到更多的变量上。

再次膜拜大佬：# 算法学习笔记(36): 莫比乌斯反演

---

其实一般讲莫比乌斯反演到这里就没有了，但是我看了《离散数学》中的莫比乌斯反演一章。我发现两者根本不在同一个位阶上……这就是颠覆我认知的原因。

所以这里我还要把莫比乌斯反演扩展出来。

---

莫比乌斯再认识

我们考虑这么一个情况：

有集合 X 和偏序关系 (P(X),⊆)，设：

F:P(X)→RG:P(X)→R

其中：G(S)=∑T⊆SF(T)。

则可以求得：F(S)=∑T⊆S(−1)|S|−|T|G(T)

由 G 求的 F 的过程称为反解，其中，(−1)|S|−|T| 就是莫比乌斯函数在这种情况下的取值，也称为容斥系数。

顺便回顾一下基本容斥原理：

设 A1,A2,⋯,An 是有限集 S 的子集（代表 n 中属性？）定义 F(K) （K 为下标集合，⊆{1,2,⋯,n}）为集合 S 中 ∈Ai(i∉K) 的元素的个数。也就是对于 s∈S，计数 s 当且仅当：

∀i∈K,s∉Ai∀j∉K,s∈Aj

于是设 G(K)=∑L⊆KF(L)，

其计数 S 中属于 j 不在 K 中的所有 Aj 的元素，以及属于其他的一些集合的元素的个数。因而还有：

G(K)=|⋂i∉KAi|

根据上文，有

F(K)=∑L⊆K(−1)|K|−|L|G(L)

此时 F(Xn)(Xn={1,2,⋯,n}) 计数的是那些仅属于满足 i∉Xn 的集合 Ai 的元素，因此：

F(Xn)=|⋂i∈Xn¯Ai|

带入 (1) 中可以得到：

|¯A1∩¯A2∩⋯∩¯An|=∑L⊆Xn(−1)n−|L||⋂i∉LAi|

如过等价的利用 L 的补集，那么我们有：

|¯A1∩¯A2∩⋯∩¯An|=∑J⊆Xn(−1)J|⋂i∈JAi|

这就是基本的容斥原理。

---

二项式反演

为什么突然到这里了……

二项式反演可以说是上面内容的一种特殊形式。其满足：

|S|=|T|⟹F(S)=F(T),G(S)=G(T)

此时我们可以直接通过集合大小表示：F(S)=f(|S|),G(S)=g(|S|)

于是对于 G(K)=∑L⊆KF(L)，合并相同大小的子集，即可得到：

g(k) = \sum_{l = 0}^{k} {k \choose l} f(l)

根据反演，也就有：

f(k) = \sum_{l = 0}^k (-1)^{k - l} {k \choose l} g(l)

这也就是二项式反演在此的推导，这里的莫比乌斯函数 \mu，后文再说。

---

扩展到偏序集

在此，我们扩展到任意有限偏序集 (X, \le)。不过为了得到莫比乌斯函数，我们首先考虑二元变量。

设 \mathbb{F}(X) 是满足只要 x \not \le y 就有 f(x, y) = 0 的所有实数函数的集合。

f: X \times X \to \R

我们如此定义卷积 h = f * g：

h(x, y) = \begin{cases} \sum_{\{z: x \le z \le y\}} f(x, z)g(z, y) & (x \le y) \\ 0 & otherwise \end{cases}

不难证明卷积满足结合律，这部分留个读者思考。

于是，我们重新定义如下函数：

• 单位函数（克罗内克 delta 函数）：

\delta(x, y) = \begin{cases} 1 & \text{if } x = y \\ 0 & \text{otherwise } \end{cases}

• 常数函数（zeta function）：

\zeta(x, y) = \begin{cases} 1 & \text{if } x \le y \\ 0 & \text{otherwise} \end{cases}

至于莫比乌斯函数，与上文的定义类似，也就是 \zeta 的逆函数：

\mu * \zeta = \delta

于是通过卷积的定义，我们得到：

\sum_{\{z: x \le z \le y\}} \mu(x, z)\zeta(z, y) = \delta(x, y) \qquad (x \le y)

或等价的：

\sum_{\{z: x \le z \le y\}} \mu(x, z) = \delta(x, y) \qquad (x \le y) \tag{2.1}

而等式 (2.1) 意味着，对于所有的 x:

\mu(x, x) = 1

以及：

\mu(x, y) = -\sum_{\{z: x \le z \lt y\}} \mu(x, z) \qquad (x < y)

至于莫比乌斯反演，无非还是：

f * \zeta = g \iff f = g * \mu

---

于是我们重新思考二项式反演，其实就是偏序集 (P(X_n), \subseteq) 上的莫比乌斯反演。

设 A 和 B 是 X_n 的子集，且 A \subseteq B，有 \mu(A, B) = (-1)^{|B| - |A|}。

这可以通过归纳假设证明，这里不过多展开。

---

开篇讲的狄利克雷卷积上的莫比乌斯反演，其实就是偏序集 (\Z, |) 上的莫比乌斯反演。

这东西谁都知道，一元的莫比乌斯函数 \mu(x) 怎么求。不过我们的目标是计算该偏序集的 \mu(a, b)。

但是，由于如果 a | b 则 \mu(a, b) = \mu(1, \frac ba)。所以我们只需要算 \mu(1, x) 即可。

而 \mu(x) 其实就是 \mu(1, x) ……

---

考虑离散傅立叶变换。

越扯越远了……QwQ

不了解离散傅立叶变换的可以参考：算法学习笔记(17): 快速傅里叶变换（FFT）

我们不是有：

h(\omega^x) = \sum_{k = 0}^{n - 1} c_k \omega^{kx}

我们整理一下重新写出：

g(x) = \sum_{k = 0}^{n - 1} \omega^{kx} f(k)

根据离散傅立叶逆变换，则有：

f(x) = \frac 1n \sum_{k = 0}^{n - 1} \omega^{-kx} g(k)

其中，\frac 1n \omega^{-kx} 就是容斥系数，\mu(k, x)。

---

最后的最后，提一个题吧：[春季测试 2023] 幂次 - 洛谷。

其实也可以通过容斥（求 \mu）求解，但并非反演。

参见博客：幂次 - Jijidawang