反演与狄利克雷卷积
修订日期：2021.12.3. 更改了很多本质上的错误。
建议首先阅读组合数学相关 Part 3. 容斥原理。
1. 本质与第一反演公式
1.1. 定义
反演是 \(f\) 表示 \(g\) \(\to\) \(g\) 表示 \(f\)。
已知 \(g(n) = \sum_{\\i = 1} ^ n c_{n, i} f(i)\)，根据系数 \(c_{i, j}\) 可以反推出用 \(g\) 表示 \(f\) 的方法：\(f(n) = \sum_{\\i = 1} ^ n d_{n, i} g(i)\)。这样，就算 \(f\) 无法直接求，也可以通过所有 \(g(i)\) 推出 \(f(n)\)。上述二式互为反演公式。
1.2. 第一反演公式
直接给出定理：如果存在 \(n\) 次多项式 \(p, q\) 分别满足以下公式：\([x ^ n]p = \sum_{\\i = 1} ^ n c_{n, i} \times [x ^ i]q\) 以及 \([x ^ n]q = \sum_{\\i = 1} ^ n d_{n, i} \times [x ^ i]p\)，且 \(c_{i, i}\) 和 \(d_{i, i} \neq 0\)，那么对于任意序列 \(u, v\)，都有：
如果将 \(c, d\) 写成矩阵 \(\textbf{C}, \textbf{D}\)，那么它是下三角矩阵。将 \(p, q\) 写成列向量 \(\textbf{p}, \textbf{q}\)，那么条件中的两个公式可以表示为 \(\textbf{p} = \textbf{Cq}\)，\(\textbf{q} = \textbf{Dp}\)，即 \(\textbf{p} = \textbf{CDp} \Rightarrow \textbf{CD} = \textbf{I}\)。由于 \(\textbf{C}\) 与 \(\textbf{D}\) 互逆，所以它们的行列式不为 \(0\)，而因为它们是下三角矩阵，因此一定有对角线元素不为 \(0\)，这解释了为什么 \(c_{i, i}\) 和 \(d_{i, i} \neq 0\)。
只要我们找到任何一组符合条件的 \(p, q\) 与 \(c, d\)，就可以将 \(p, q\) 推广至任意 \(n\) 次多项式，因为 \(\textbf{C}\) 与 \(\textbf{D}\) 恒互逆，不随 \(p, q\) 的改变而改变。这说明为了证明某个反演公式成立，我们只需考虑一种方便的特殊非平凡情况即可。
2. 二项式反演
2.1. 引入
在把玩二项式定理的时候，尝试把 \(x\) 变成 \(1+(x - 1)\) 可能会有新的收获：
嗯？和反演公式的形式长得很像！那么试试能不能用 \(x\) 表示 \(x - 1\)：
很好，令 \(p_i = x ^ i\)，\(q_i = (x - 1) ^ i\)，\(c_{i, j} = \dbinom i j\)，\(d_{i, j} = \dbinom i j ( - 1) ^ {i - j}\)，那么我们就成功找到了一组反演公式，称为形式一。
代数化证明：
众所周知，右边的 \(\sum\) 当且仅当 \(n - j = 0\) 即 \(j = n\) 时为 \(1\)，否则为 \(0\)。因此 \(f_n = \sum_{\\j = 0} ^ n f_n \dbinom n j [n = j] = f_n\)，证毕。
如果将矩阵转置，那么还可以得到另一种更常用表示，形式二：
2.2. 组合意义
二项式反演的本质是容斥原理。\(f_i\) 表示钦定选 \(i\) 个的方案数，\(g_i\) 表示恰好选 \(i\) 个的方案数。对于任意的 \(i\geq n\) ，\(g_i\) 在 \(f_n\) 中被计算了 \(\dbinom{i}{n}\) 次，故 \(f_n=\sum_{\\i=n} ^ m\dbinom i n g_i\) ，其中 \(m\) 是数目上界，恰好对应形式二。这说明可以使用 \(\sum_{\\i = n} ^ m (-1) ^ {i - n} \dbinom i n f_i\) 求出 \(g_n\)。
注意：在定义中，\(f_n\) 表示先钦定 \(n\) 个，再统计钦定情况如此的方案数，其中会包含重复的方案，因为一个方案可以有多种钦定情况。具体地，对于恰好选择 \(i\) 个的方案，钦定情况数为 \(\dbinom i n\) ，故 \(g_i\) 在 \(f_n\) 中被计算了 \(\dbinom n i\) 次。切勿将 \(f_n\) 理解为普通的 \(g_i\) 的后缀和。
很多初学者最大的误区在于将 \(f_n\) 理解为至少选 \(n\) 个的方案数，这是完全错误的（如果真是这样那么 \(g_n\) 不就等于 \(f_n-f_{n+1}\) 么，还要二项式反演干啥 = . =）。
2.3. 应用
OI 界二项式反演的应用常结合动态规划：DP 求出钦定选择 \(i\) 个元素时的总方案，然后就可以通过二项式反演得到恰好选择 \(i\) 个元素时的方案数。综上所述，我们总结出二项式定理的核心：通过二项式系数的容斥进行 “钦定” 和 “恰好” 的转换。
2.3.1. 错排问题
\(n\) 个人排列，求每个人都恰好站错的方案数。即求 \(\sum_{\\p} [\forall i, \ p_i\neq i]\)。
设 \(g_n\) 表示 \(n\) 个人的错排方案数，而 \(f_n\) 表示所有排列方案数即 \(n!\)。那么有 \(f_n = \sum_{\\i = 0} ^ n \dbinom n i g_i\)，\(i\) 的意义是枚举站错的人的个数。反演得到 \(g_n = \sum_{\\i = 0} ^ n (-1) ^ {n-i} \dbinom n i f_i = \sum_{i=0} ^ n (-1) ^ {n - i} \dfrac {n!} {(n-i)!}\)（组合数分母上的 \(i!\) 与 \(f_i\) 抵消了）。
2.3.2. 染色问题 I（一维）
\(n\) 个格子，\(k\) 种颜色，相邻格子颜色不同，每个颜色至少出现一次，求方案数。
设 \(g_k\) 表示恰好出现 \(k\) 个颜色的方案数，\(f_k\) 表示 \(k\) 个颜色的总方案数，即 \(k \times (k - 1) ^ {n - 1}\)。有
反演得到 \(g_k = \sum_{\\i = 0}^k(-1) ^ {k - i} \dbinom k i f_i\)。可以 \(\mathcal{O}(k\log n)\) 计算。
2.3.3. 染色问题 II（二维）
\(n\times m\) 个格子，\(1\) 种颜色，每行和每列至少有一个格子被涂色，求方案数。
设 \(g_{i, j}\) 表示恰好有 \(i\) 行 \(j\) 列被涂色的方案数，\(f_{i, j}\) 表示 \(i\) 行 \(j\) 列被涂色的方案数即 \(2 ^ {ij}\)，有
对 \(i\) 和 \(j\) 分别二项式反演，有
2.3.4. 染色问题 III（三维）
\(n\times m\) 个格子，\(k\) 种颜色，每行和每列至少一个格子被涂色，每个颜色至少出现一次，格子可不被涂色，求方案数。
请读者结合染色问题 I. 与 II. 自行思考。设 \(g_{n, m, k}\) 表示恰有 \(n\) 行 \(m\) 列至少有一个格子被涂上颜色，且恰有 \(k\) 种颜色的方案数，\(f_{n, m, k}\) 表示总方案数 \((k + 1) ^ {nm}\)。
直接计算 \(n ^ 3\)，可以通过二项式定理化简减小复杂度。
2.4. 参考文章
- GXZlegend —— 二项式反演及其应用（推荐）。
- liu_jiangwen ——111111 组合数学——二项式反演。
- Icyfrog —— 二项式反演入门。
- WAautomaton —— 反演定理及其应用。
2.5. 例题
*I. P6478 [NOI Online #2 提高组] 游戏
很不错的题目，可以用来入门二项式反演。设 \(f_x\) 为钦定 \(x\) 对祖先关系的情况总数，\(g_x\) 为答案，那么根据二项式反演有 \(g_x = \sum_{\\i = x} ^ n (-1) ^ {i - x} \dbinom i xf_i\)。\(f_x\) 可以树形背包求，别忘了最后乘上 \((m-x)!\)，因为剩下 \(m-x\) 对点对顺序任意。时间复杂度 \(\mathcal{O}(n^2)\)。
*II. P4859 已经没有什么好害怕的了
和例题 I. 极度类似。DP 方法几乎一样，最后也都要乘以阶乘作为系数。实际上，上一题就是把本题搬到了树上。
const int mod = 1e9 + 9;
const int N = 2e3 + 5;
int n, k, ans, a[N], b[N], num[N], f[N], fc[N], ifc[N];
void add(int &x, int y) {x += y, x >= mod && (x -= mod);}
int ksm(int a, int b) {
int s = 1;
while(b) {
if(b & 1) s = 1ll * s * a % mod;
a = 1ll * a * a % mod, b >>= 1;
} return s;
}
int main(){
cin >> n >> k;
if((n & 1) != (k & 1)) return puts("0"), 0;
k = n + k >> 1, fc[0] = 1;
for(int i = 1; i <= n; i++) fc[i] = 1ll * fc[i - 1] * i % mod;
ifc[n] = ksm(fc[n], mod - 2);
for(int i = n - 1; ~i; i--) ifc[i] = 1ll * ifc[i + 1] * (i + 1) % mod;
for(int i = 1; i <= n; i++) cin >> a[i];
for(int i = 1; i <= n; i++) cin >> b[i];
sort(a + 1, a + n + 1), sort(b + 1, b + n + 1);
for(int i = 1; i <= n; i++) {
num[i] = num[i - 1];
while(num[i] < n && a[i] > b[num[i] + 1]) num[i]++;
} f[0] = 1;
for(int i = 1; i <= n; i++)
for(int j = min(num[i], i); j; j--)
add(f[j], 1ll * f[j - 1] * (num[i] - j + 1) % mod);
for(int i = k; i <= n; i++) {
int bin = 1ll * fc[i] * ifc[k] % mod * ifc[i - k] % mod;
int coef = 1ll * bin * f[i] % mod * fc[n - i] % mod;
add(ans, (i - k) & 1 ? mod - coef : coef);
} cout << ans << endl;
return 0;
}
*III. CF1228E Another Filling the Grid
设 \(g_{i, j}\) 表示恰有 \(i\) 行 \(j\) 列的最小值为 \(1\)。若钦定最小值为 \(1\) 无法计算，但钦定最小值不为 \(1\) 就很好算了。因此我们设 \(f_{i, j}\) 表示钦定 \(i\) 行 \(j\) 列使得其它行列的最小值非 \(1\) 的方案数，即 \((k - 1) ^ {(n - i)(n - j)} \times k ^ {n ^ 2 - (n - i)(n - j)}\)，那么有
二项式反演得到
可以 \(\mathcal{O}(n^2\log k)\) 计算，时间复杂度足够优秀。当然，精细实现可以做到平方甚至线性对数。
IV. [BZOJ2839]集合计数
考虑钦定 \(i\) 个元素在交集里面，那么方案数为 \(\dbinom n i 2^{2^{n-i}}\)。二项式反演掉，得到答案为 \(\sum_{\\i=k}^n(-1)^{k-i}\dbinom i k\dbinom n i 2^{2^{n-i}}\)。
V. [BZOJ4665]小 w 的喜糖
钦定有 \(k\) 组不同的方案数不好算，考虑钦定 \(k\) 组相同的方案数为 \(f_k\)，恰好 \(k\) 组相同的方案数为 \(g_k\)，有
二项式反演即可，\(f_k\) 可以通过简单背包乘以二项式系数计算（\(f_i=\sum_{\\j=0}^c f_{i-j}\dbinom c j c^{\underline{j}}\)）。注意最后要除掉 \(\prod (cnt_i)!\)，因为种类相同的糖本质相同。时间复杂度平方。
*VI. [BZOJ2162]男生女生
本题可以看做割裂的两部分，一部分是求使某部点尽量大的二分图最大团，见 网络流与二分图 Part 2.5. 二分图的最大团。另一部分则需要用到二项式反演：
考虑求出男生个数 \(B\) 和女生个数 \(G\) 后如何计算答案：钦定 \(i\) 个男生和 \(j\) 个女生不能被选，方案数为 \(f_{i,j}= \dbinom{(B-i)(G-j)} k\)。设恰好 \(i\) 男 \(j\) 女没有被选到的方案数为 \(g_{i,j}\)，那么 \(f_{i,j}=\sum_{\\i'=i}^B\sum_{j'=j}^G\dbinom {i'} i\dbinom {j'}jg_{i,j}\)，二项式反演掉就做完了。由于模数不是质数所以预处理组合数的复杂度是瓶颈，为 \(\mathcal{O}(n^4)\)。
3. 单位根反演
3.1. 公式
常规证法是用等比数列求和公式，我不太喜欢，因为这样不能使初学单位根反演的同学感受到它的自然与优美。
一个比较直观的我自己的理解是 \(\omega_n^x\) 在 \(n\mid x\) 时等于 \(1\)，\(n\) 个求和后再除以 \(n\) 就是 \(1\)；而在 \(n\nmid x\) 时不妨设 \(d=\gcd(x,n)\)。把单位根放到平面直角坐标系里感受一下，\(\sum_{i=0}^{n-1}(\omega_n^x)^i\) 就是把 \(\omega_n^x\) 旋转 \(n\) 次，每次旋转相同的 \(\dfrac{2\pi x}{n}\) 角度。那么最终形成了相同的 \(d\) 组 “向量”，每组内部有 \(\dfrac{n}{d}\) 个在单位圆上均匀分布的 “向量”，根据高中知识可知每一组向量的和都是 \(0\)。
它还有另外一种形式，在做题时经常会遇到：
3.2. 应用
我们为什么要用单位根反演？钻研这个问题，不仅可以帮助我们深入理解单位根反演，在做题时我们也能更容易看出它的蛛丝马迹。
3.2.1. 将 mod 转化为求和
首先看一道经典的单位根反演入门题：给出 \(n,s,a_0,a_1,a_2,a_3\leq 10^{18}\)，求 \(\sum_{\\i=0}^n\dbinom{n}{i}s^ia_{i\bmod 4}\) 对 \(998244353\) 取模。
看到组合数和幂结合，自然想到二项式定理，可是 \(a_{i\bmod 4}\) 的部分断绝了我们的出路 …… 吗？并没有。如果我们将 \(a_{i\bmod 4}\) 写作 \(\sum_{j=0}^3a_j[i\equiv j\pmod 4]\)，再使用单位根反演得到 \(\sum_{j=0}^3a_j\times \dfrac{1}{4}\sum_{k=0}^3\omega_4^{(i-j)k}\)，揉进原来的柿子：
交换求和符号并整理，得
逆用二项式定理，得
如果你学过 NTT 就知道在模 \(998244353\) 意义下，原根和单位根可以互换，因此这里的四次单位根等价于 \(g^{\frac{p-1}4}\)，快速幂解决。时间复杂度 \(\mathcal{O}(\log n)\)。
代码见例题 I.
3.3. 例题
*I. LOJ#6485. LJJ 学二项式定理
代码。
*II. P5591 小猪佩奇学数学
推柿子：根据 \(\left\lfloor\dfrac{i}k\right\rfloor=\dfrac{i-i\bmod k}k\) 得
略去 \(\dfrac 1 k\)，括号拆开分别考虑。
-
前半部分：\(\sum_{\\i=0}^n\dbinom n i p^i\times i\)。
根据 \(\dbinom{n}{m}=\dfrac n m\dbinom {n-1}{m-1}\)，得 \(np\sum_{\\i=0}^n\dbinom {n-1}{i-1} p^ {i-1}\)。
即 \(np(p+1)^{n-1}\)。
-
后半部分：略去负号，得 \(\sum_{\\i=0}^{n}\dbinom{n}{i}p^i\sum_{\\j=0}^{k-1}j[i\equiv j\pmod k]\)。
套入单位根反演，得 \(\sum_{\\i=0}^n\dbinom n ip^i\sum_{\\j=0}^{k-1}\dfrac j k\sum_{\\x=0}^{k-1}\omega_k^{(i-j)x}\)。
略去 \(\dfrac 1 k\)，交换求和符号，得 \(\sum_{\\j=0}^{k-1}j\sum_{\\x=0}^{k-1}\omega_k^{-jx}\sum_{\\i=0}^n\dbinom ni(p\omega_k^x)^i\)。
套入二项式定理，得 \(\sum_{\\j=0}^{k-1}j\sum_{\\x=0}^{k-1}\omega_k^{-jx}(p\omega_k^x+1)^n\)。
然后发现不会了。交换求和顺序，得 \(\sum_{\\x=0}^{k-1}(p\omega_k^x+1)^n\sum_{\\j=0}^{k-1}j(\omega_k^x)^{-j}\)。略去前面的柿子，将 \(\omega_k^{-x}\) 记做 \(c\)，仅关注 \(\sum_{\\j=0}^{k-1}jc^{j}\)。由于 \(k\) 是常数，因此将其记做 \(f(c)\)。看到这里，你发现了什么？对！是我们小学二年级就学过的错位相减法！
-
如果 \(c\neq 1\)，求 \(cf(c)-f(c)\)：比对每一项的系数，有 \(f(c)-cf(c)=\sum_{\\j=1}^{k-1}c^j-(k-1)c^k\)。
先特判掉 \(k=1\) 的情况，然后用等比数列求和公式得到 \(\sum_{\\j=1}^{k-1}c^j=-1+\sum_{\\j=0}^{k-1}c^j=-1+\dfrac{c^k-1}{c-1}\)，注意到 \(c^k\) 恒为 \(1\)，所以 \(f(c)=\dfrac{-1-(k-1)}{1-c}=\dfrac{k}{c-1}\)。
-
如果 \(c=1\)，那么原式等于 \(\sum_{\\j=0}^{k-1}j=\dfrac{k(k-1)}2\)。
因此直接计算即可！时间复杂度 \(\mathcal{O}(k\log n)\)。
-
ll n, p, k, w, ans, sub;
int main() {
cin >> n >> p >> k, w = ksm(3, (mod - 1) / k);
ans = n * p % mod * ksm(p + 1, n - 1) % mod;
for(int x = 0; x < k; x++) {
ll om = ksm(w, x), pw = ksm(p * om + 1, n), c = inv(om);
if(c == 1) sub = (sub + pw * (k * (k - 1) / 2 % mod)) % mod;
else sub = (sub + pw * k % mod * inv(mod + c - 1)) % mod;
} cout << (ans + mod - sub * inv(k) % mod) * inv(k) % mod;
return 0;
}
4. 莫比乌斯反演
极其重要的一类反演，对于推式子的帮助很大。
4.1. 前置知识
4.1.1. 数论函数与积性函数
数论函数，就是定义域是正整数的函数。
积性函数，就是对于任意 \(x,y\in \N^{+}\)，若 \(\gcd(x,y)=1\)，则 \(f(xy)=f(x)f(y)\) 的函数 \(f\)。若完全积性则不需要 \(\gcd(x,y)=1\) 的条件。
积性函数举例（下文中会出现）：
- 单位函数 \(\epsilon(n)=[n=1]\)。
- 常数函数 \(1(n)=1\)。
- 恒等函数 \(\mathrm{id}_k(n)=n^k\)，\(\mathrm{id}_1(n)\) 通常记作 \(\mathrm {id}(n)\)。
- 除数函数 \(\sigma_k(n)=\sum_{\\d\mid n}d^k\)。\(\sigma_0(n)\) 就是约数个数，记作 \(\tau(n)\) 或 \(\mathrm{d}(n)\)。\(\sigma_1(n)\) 就是约数和，记作 \(\sigma(n)\)。
- 大名鼎鼎的欧拉函数 \(\varphi(n)=\sum_{\\i=1}^n[\gcd(i,n)=1]\)。关于欧拉函数可以看 基础数论学习笔记。
- 本章节的核心莫比乌斯函数 \(\mu(n)=\begin{cases}0&\exist\ d>1,d^2\mid n \\(-1)^{\omega(n)}&\mathrm{otherwise}\end{cases}\)，其中 \(\omega(n)\) 表示 \(n\) 本质不同质因子个数。它的积性是显然的。
4.1.2. 狄利克雷卷积
由于学习莫比乌斯反演之前需要学会狄利克雷卷积，所以就放在这了。狄利克雷卷积竟沦落到给莫比乌斯反演作铺垫（大雾）。
定义：设 \(f,g\) 是数论函数，定义它们的狄利克雷卷积为 \(h(x)=\sum_{\\a\times b=x}f(a)g(b)\)，写作 \(h=f*g\)。显然的，狄利克雷卷积具有交换律，结合律和乘法分配律。它的另一种表达形式为 \(h(x)=\sum_{\\d\mid x}f(d)g\left(\dfrac x d\right)\)，非常有用。
FMT 和 FWT 叫位运算卷积，FFT 和 NTT 叫加法卷积，不难发现狄利克雷卷积实际上是乘法卷积。它有如下性质：两个积性函数的狄利克雷卷积也是积性函数，积性函数的逆元也是积性函数。证明略。
4.1.3. 莫比乌斯函数
因为莫比乌斯函数是积性函数，所以它可以由线性筛求得。莫比乌斯函数有以下几条性质：
首先 \(n=1\) 时显然，否则我们令 \(n\) 的本质不同质因数分别为 \(p_1,p_2,\cdots,p_k\)，从中选取 \(i\) 个质数的方案数为 \(\dbinom ki\)。由选奇数个方案数之和等于选偶数个方案数之和可知和为 \(0\)：\(k>0\) 时 \(\sum_{i=0}^k1^{k-i}\times\dbinom k i\times (-1)^{i}=(1+(-1))^k=0\)（二项式定定理）。因而证明了 \(\mu*1=\epsilon\)。
该式子在推式子时比较有用。可以先证明 \(\varphi*1=\mathrm{id}\)，两边同时卷积 \(\mu\) 有 \(\varphi*(1*\mu)=\varphi*\epsilon=\varphi=\mathrm{id}*\mu\)。得证。
变式：
令上式两端同除 \(n\) 即得。
极其重要的一个柿子，是莫比乌斯反演的核心。将 \(\gcd(i,j)\) 带入式 \(4.1\) 即得。
线性筛莫比乌斯函数：根据 \(\mu(n)\times \mu(p)=-\mu(n)=\mu(np)\ (\gcd(n,p)=1)\) 可得。
mu[1] = 1;
for(int i = 2; i < N; i++) {
if(!vis[i]) pr[++cnt] = i, mu[i] = -1;
for(int j = 1; j <= cnt && i * pr[j] < N; j++) {
vis[i * pr[j]] = 1;
if(i % pr[j] == 0) break;
mu[i * pr[j]] = -mu[i];
}
}
4.1.4. 整除分块
哪里有莫反，哪里就有整除分块。详见 基础数论学习笔记 Part 5 整除分块。
4.1.5. 线性筛积性函数
见本文 Part 6.1。
4.2. 公式
若 \(f,g\) 是两个数论函数，那么有
证明：已知 \(g*1=f\)，则 \(g*1*\mu=g=f*\mu\)。
证明（来自 OI - Wiki）：
4.3. 直观理解
只会套皮是不行的，要想熟练运用，不仅要多刷题，也要理解莫反的本质：容斥原理。其中莫比乌斯函数就是容斥系数。也许你会问：容斥系数不应该是 \((-1)^k\) 吗？其实在某种程度上没错，毕竟对比 \((-1)^k\) 与莫比乌斯函数的定义，它们的形式非常相近：只有当 \(x\) 是完全平方数的倍数时容斥系数才为 \(0\)（这其实也暗示了这样的 \(x\) 不同于别的数）。但我们是对于全体自然数做容斥，会有一些奇怪的事情发生：一个自然数被另一个自然数完全包含（这也说明它一定是完全平方数的倍数，接下来会解释），而在容斥原理中，一个集合被另一个集合包含这样的事情是不会发生的。
不妨设接下来有 \(f_n=\sum_{d\mid n}g_d\)。
4.3.1. \(n=2\)
当 \(n=2\) 时，为了求出 \(g_2\)，我们只需要用 \(f_2-f_1=(g_2+g_1)-g_1\)。这里已经可以看见容斥的影子了。
4.3.2. \(n=6\) 及 \(\prod p_i\)
当 \(n=6\) 时，为了求出 \(g_6\)，我们首先需要 \(f_6\)，但是这样会多加 \(g_1+g_2+g_3\)，所以减去 \(f_2+f_3\)，但是这样又多减了 \(g_1\)，所以我们还应加上 \(f_1\)。
如果 \(n\) 是若干不同质数 \(p_1,p_2,\cdots,p_k\) 的乘积，画出能够表示 \(n\) 的所有约数的韦恩图（\(k>3\) 时画不出来，不过可以理解一下）。\(k\) 个集合相交形成的所有子集就是 \(n\) 的所有约数。比如说，如果有一块区域落在集合 \(p_1\) 和集合 \(p_3\) 中，那么它所表示的约数就是 \(p_1p_3\)。显然，如果 \(x\mid y\)，那么区域 \(y\) 包含于区域 \(x\)。
结合上图，不妨将 \(f_d\) 看做整块区域 \(d\)（例如上图中区域 \(30\) 也属于区域 \(6\)），而 \(g_d\) 看做标上 \(d\) 的小区域，然后将所有区域表示的数从 \(i\) 替换成 \(\dfrac n i\)，那么 \(f_n\) 就表示最外层的大区域，即所有 \(g_d\) 的和，而 \(f_1\) 就表示最里层的小区域，即 \(g_1\)，刚好符合 \(f_n=\sum_{\\d\mid n} g_d\) 的条件。因此，由于 \(f\) 已知，为了得到仅属于区域 \(n\) 的部分所表示的数 \(g_n\)，我们只需要进行容斥：加上 \(f_n\)，减去 \(f_{p_1},\cdots, f_{p_k}\)，加上 \(f_{p_1p_2}, f_{p_1,p_3},\cdots,f_{p_{k-1}p_k}\) ，以此类推。不难发现每个数前的容斥系数就是 \(-1\) 的 \(\dfrac n i\) 所含质因子个数次幂，即 \(\mu\left(\dfrac n i\right)\)（\(n\) 由互不相同的质因子相乘得到，那么 \(\dfrac n i\) 亦然）。这不就是 \(g_n=\sum_{d\mid n}f_d\mu\left(\dfrac n d\right)\) 吗！
4.3.3. \(n=4\) 及任意自然数
有了上面的经验，你应该知道表示约数 \(4\) 的区域被表示约数 \(2\) 的区域 “完全包含”。画出韦恩图就是一大一小两个圆，其中一个圆完全落在另一个里面。因此，如果计算了 \(2\) 的贡献，那么完全没有必要再算 \(4\) 了，因为 \(4\) 被 \(2\) “完全包含”，所以它们的 “多加/多减状态” 应该是相同的。也就是说算完 \(2\) 的贡献，\(4\) 的贡献也恰好算完，既不多加，也不多减，因而 \(\mu(4)=0\)。不难发现，这种情况仅在一个数是非 \(1\) 的完全平方数的倍数时出现，从而得以将莫比乌斯函数推广至任意自然数。
综上，莫比乌斯函数可以看作对 \(\N^{+}\) 做容斥时的容斥系数。结合 \(1 *\mu=\epsilon\) 这一性质，\(\mu(n)\) 的定义变得十分自然：\(n\) 含有完全平方因子 \(p\) 时可以直接忽略不计，因为 \(n\) 的贡献在 \(\dfrac n p\) 处已经计算过。否则 \(n\) 含有的质因子个数就是其在容斥时所表示的集合大小，因为当 \(n=1\) 时需定义为 \(1\) 故 \(\mu(n)=(-1)^{\omega(n)}\)。
4.4. 应用
一般用来化简带有 \(\gcd\) 的柿子，几乎所有莫反题都需要用到整除分块。具体应用可以看 OI-Wiki，抄在博客里其实意义不大。
4.5. 例题
下面所有分式都是下取整。有部分整除分块的题目。
I. P2522 [HAOI2011]Problem b
二维差分将下界化为 \(1\)，然后推式子：
把 \(k\) 搞掉：
莫反：
枚举约数 \(d\)，记 \(c=\min\left(\dfrac n k,\dfrac m k\right)\)：
由于 \(1\sim x\) 中 \(y\) 的倍数有 \(\dfrac x y\) 个，故原式化为：
整除分块即可，时间复杂度 \(\mathcal{O}(n+T\sqrt n)\)。long long 改成 int，11s 变成 5s，离谱。
II. SP5971 LCM Sum
来推式子。
\(\gcd\) 在分母的位置，很不爽，尝试枚举 \(\gcd\)：
单独算 \(F(d)=\sum_{\\i=1}^di[\gcd(i,d)=1]\)：
至此，我们可以 \(n\ln n\) 枚举 \(d\) 及其倍数 \(kd\) 算出所有 \(F(n)\)，再根据 \(F(i)\) 枚举 \(d\) 及其倍数算出所有数的答案。时间复杂度 \(\mathcal{O}(T+n\ln n)\)。实际上可以继续化简：
这里用到了 \(\mu*1=\epsilon\) 以及 \(\mu*\mathrm{id}=\varphi\)。由于 \(F=\dfrac{\mathrm{id}\times(\epsilon+\varphi)}2\)，两部分都是积性函数，那么 \(1*F\) 的两部分仍是积性函数，可以线性筛求出。时间复杂度 \(\mathcal{O}(T+n)\)。
const int N = 1e6 + 5;
ll T, pr[N], G[N], low[N], cnt;
bool vis[N];
void sieve() {
G[1] = 1;
for(int i = 2; i < N; i++) {
if(!vis[i]) pr[++cnt] = i, G[i] = 1 + 1ll * i * (i - 1), low[i] = i;
for(int j = 1; j <= cnt && i * pr[j] < N; j++) {
vis[i * pr[j]] = 1;
if(i % pr[j] == 0) {
low[i * pr[j]] = low[i] * pr[j];
if(i == low[i]) G[i * pr[j]] = G[i] + low[i] * (pr[j] - 1) * low[i] * pr[j];
else G[i * pr[j]] = G[i / low[i]] * G[low[i] * pr[j]]; break;
} G[i * pr[j]] = G[i] * G[pr[j]], low[i * pr[j]] = pr[j];
}
}
}
int main(){
cin >> T, sieve();
while(T--) {
int n = read();
print((G[n] + 1) * n / 2), pc('\n');
} return flush(), 0;
}
*III. P4318 完全平方数
莫比乌斯函数作容斥系数的经典例子。答案满足可二分性，于是考虑如何求 \(x\) 以内不含有平方因子的数的个数。仍然是容斥：我们只需加上所有由偶数个质数相乘的数的平方的倍数，减去由奇数个质数相乘的数的平方的倍数。发现就是 \(\sum_{\\i^2\leq n}\mu(i)\dfrac{n}{i^2}\)，直接暴力求。时间复杂度 \(\mathcal{O}(T\log k\sqrt k)\)。
const int N = 1e5 + 5;
int T, n, cnt, pr[N], vis[N], mu[N];
int main() {
cin >> T, mu[1] = 1;
for(int i = 2; i < N;i++) {
if(!vis[i]) pr[++cnt] = i, mu[i] = -1;
for(int j = 1; j <= cnt && i * pr[j] < N; j++) {
vis[i * pr[j]] = 1;
if(i % pr[j] == 0) break;
mu[i * pr[j]] = -mu[i];
}
} while(T--){
cin >> n; ll l = n, r = n << 1;
while(l < r) {
ll m = l + r >> 1, res = 0;
for(int i = 1; i * i <= m; i++) res += mu[i] * (x / i / i);
res < n ? l = m + 1 : r = m;
} cout << l << endl;
} return 0;
}
IV. P2257 YY 的 GCD
考虑枚举 \(\gcd(i,j)=p\in \mathbb{P}\cap [2,\min(n,m)]\)：
做到这里发现推不动了，因为我们不得不枚举 \(p\in \mathbb P\)。一个 trick 是将 \(pd\) 写成 \(T\)，并将 \(pd\) 写成狄利克雷卷积的形式：
为什么要这样做：分母上的 \(pd\) 与两个变量有关，换种写法仅与一个变量有关，于是可以提前，但并没有改变求 \(\mu\) 的本质，只是调整了计算顺序使得可通过乘法分配律提出向下取整的柿子。或者说，对于使 \(pd\) 相同的 \(p,d\) 取值，\(\sum \mu(d)\) 的形式可以预处理避免重复计算。这个大概要写过很多莫反题目才能来感觉。
*另一种推导方式：考虑对 \([\gcd(i,j)=p]\) 进行容斥并用莫比乌斯函数作容斥系数：\(\sum_{i=1}^n\sum_{j=1}^m\sum_{\\d=kp}\mu(k)[d\mid \gcd(i,j)]\)，相当于 将 \(\mu\) 乘上 \(\epsilon_p(i)\) 即 \([i=p]\ (p\in \mathbb P)\)，然后求和。稍作化简得到 \(\sum_{\\d=1}^{\min(n,m)}\dfrac n d \dfrac m d\sum_{\\p\mid d\land p\in\mathbb P}\mu\left(\dfrac{n}{p}\right)\)，发现与上述柿子等价。
\(f(T)=\sum_{\\p\mid T\land p\in \mathbb P} \mu\left(\dfrac T p\right)\) 可以通过类似埃氏筛的筛法 \(n\log\log n\) 求出，前缀和之后整除分块即可。总时间复杂度 \(\mathcal{O}(T\sqrt n+n\log \log n)\)。\(f(T)\) 还可以线性筛，具体推导过程略去，复杂度 \(\mathcal{O}(T\sqrt n+n)\)。
const int N = 1e7 + 5;
int T, n, m, cnt, vis[N], pr[N >> 3], f[N], mu[N];
int main(){
for(int i = 2; i < N; i++) {
if(!vis[i]) pr[++cnt] = i, f[i] = 1, mu[i] = -1;
for(int j = 1; j <= cnt && i * pr[j] < N; j++) {
vis[i * pr[j]] = 1;
if(i % pr[j] == 0) {f[i * pr[j]] = mu[i]; break;}
f[i * pr[j]] = -f[i] + mu[i], mu[i * pr[j]] = -mu[i];
}
} for(int i = 2; i < N; i++) f[i] += f[i - 1]; cin >> T;
while(T--) {
nll ans = 0 cin >> n >> m;
for(int l = 1, r; l <= min(n, m); l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans += 1ll * (n / l) * (m / l) * (f[r] - f[l - 1]);
} cout << ans << "\n";
}
return 0;
}
V. P3455 [POI2007]ZAP-Queries
双倍经验 P2522。
VI. P2568 GCD
双倍经验 P2257.
VII. P1829 [国家集训队]Crash的数字表格 / JZPTAB
已经可以 \(n\ln n\) 做了，但是还不够。到这里我们可以看到和 IV. 差不多的套路，令 \(T\gets kd\)：
注意到前面一坨可以整除分块，后面相当于求 \(f(T)\) 其中 \(f=(\mu\times \mathrm{id}^2)*\mathrm{id}\)（化简一下也可以看成 \(T\times f(T)\) 其中 \(f = \mu*\mathrm{id}\)），是积性函数，可以线筛求出。时间复杂度 \(\mathcal{O}(n)\)。实际上可以多组询问。
const int N = 1e7 + 5;
const int mod = 20101009;
int n, m, c, pr[N >> 3], cnt, low[N];
bool vis[N]; ll f[N], ans;
ll S(ll x) {return x * (x + 1) / 2 % mod;}
int main(){
cin >> n >> m, c = min(n, m), f[1] = 1;
for(ll i = 2; i <= c; i++) {
if(!vis[i]) low[i] = pr[++cnt] = i, f[i] = (i - i * i % mod + mod) % mod;
for(int j = 1; j <= cnt && i * pr[j] <= c; j++) {
int k = i * pr[j]; vis[k] = 1;
if(i % pr[j] == 0) {
low[k] = low[i] * pr[j];
if(low[i] == i) f[k] = (k - pr[j] * pr[j] % mod * i % mod + mod) % mod;
else f[k] = f[i / low[i]] * f[low[k]] % mod; break;
} low[k] = pr[j], f[k] = f[i] * f[pr[j]] % mod;
}
}
for(int i = 2; i <= c; i++) f[i] = (f[i] + f[i - 1]) % mod;
for(int l = 1, r; l <= c; l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans = (ans + S(n / l) * S(m / l) % mod * (f[r] - f[l - 1] + mod)) % mod;
} cout << ans << endl;
return 0;
}
VIII. AT5200 [AGC038C] LCMs
原题的两个枚举上界都写到 \(n\) 方便求答案，那么需要将结果减掉 \(\sum a_i\) 后除以 \(2\)。记 \(b_i\) 表示 \(i\) 在 \(\{a_j\}\) 中的出现次数，有：
剩下来推式子的过程和上一题差不多，但是因为有 \(b_i\) 的存在所以无法写成一个简单式子的形式。注意推到 \(\ \sum_{\\d=1}^{V}\dfrac{1}{d}\sum_{\\k=1}^{\frac V d}\mu(k)S^2(kd)\) 其中 \(S(i)\) 表示 \(\sum ki\times b_{ki}\ (k\in \Z)\) 的时候可以 \(V\ln V\) 做，时间复杂度 \(\mathcal{O}(n+V\ln V)\)。运用 Part 5 所讲解的狄利克雷前 / 后缀和并改变枚举顺序线性筛 \(\dfrac{1}{\mathrm{id}}* \mu\) 可以 \(\mathcal{O}(n+V\log\log V)\)，即 \(\sum_{\\T=1}^V S^2(T)\sum_{d\mid T}\dfrac{1}{d}\mu\left(\dfrac{T}{d}\right)\)。为了卡常可以写成 \(\sum_{\\T=1}^VS^2(T)\dfrac{1}{T}\sum_{d\mid T}\dfrac{T}{d}\mu\left(\dfrac{T}{d}\right)\)，这样就只需要最后求答案时取模。
const int N = 1e6 + 5;
const int mod = 998244353;
int n, f[N], cnt, mx, pr[N >> 3];
ll b[N], inv[N], ans; bool vis[N];
int main(){
cin >> n, f[1] = inv[1] = 1;
for(int i = 1, x; i <= n; i++) cmax(mx, x = read()), ans -= x, b[x] += x;
for(int i = 2; i <= mx; i++) {
inv[i] = mod - mod / i * inv[mod % i] % mod;
if(!vis[i]) pr[++cnt] = i, f[i] = 1 - i;
for(int j = 1; j <= cnt && i * pr[j] <= mx; j++) {
int k = i * pr[j]; vis[k] = 1;
if(i % pr[j] == 0) {f[k] = f[i]; break;}
f[k] = f[i] * f[pr[j]];
}
} for(int i = 1; i <= cnt; i++) for(int j = mx / pr[i]; j; j--) b[j] += b[j * pr[i]];
for(int i = 1; i <= mx; i++) b[i] %= mod, ans += b[i] * b[i] % mod * f[i] % mod * inv[i] % mod;
cout << (ans % mod + mod) % mod * inv[2] % mod << endl;
return 0;
}
IX. P3911 最小公倍数之和
双倍经验。
X. P6156 简单题
其中 \(F(n)\) 表示 \(\sum_{\\i=1}^n\sum_{j=1}^n(i+j)^k\)。直接线性筛 \((\mu^2\times \mathrm{id})*\mu\) 以及 \(p^k\) 即可，时间复杂度 \(\mathcal{O}(n+\dfrac{n}{\ln n}\times \log k)\approx\mathcal{O}(n)\)。
ll k;
int low[N], kp[N], F[N >> 1], G[N];
int n, cnt, ans, pr[N / 15];
bitset <N> vis;
int main(){
cin >> n >> k, k %= mod - 1, kp[1] = 1, G[1] = 1;
for(int i = 2; i <= n << 1; i++) {
if(!vis[i]) kp[i] = ksm(i, k), pr[++cnt] = low[i] = i, G[i] = i - 1;
for(int j = 1; j <= cnt && i * pr[j] <= n << 1; j++) {
int c = i * pr[j]; vis[c] = 1, kp[c] = 1ll * kp[i] * kp[pr[j]] % mod;
if(i % pr[j] == 0) {
low[c] = low[i] * pr[j];
if(i == low[i]) G[c] = i == pr[j] ? mod - pr[j] : 0;
else G[c] = 1ll * G[i / low[i]] * G[low[c]] % mod; break;
} low[c] = pr[j], G[c] = 1ll * G[i] * G[pr[j]] % mod;
}
} for(int i = 2; i <= n; i++) G[i] = (G[i - 1] + 1ll * G[i] * kp[i]) % mod;
for(int i = 1; i <= n << 1; i++) kp[i] = (kp[i - 1] + kp[i]) % mod;
for(int i = 1, j = 2; i <= n; i++, j += 2) F[i] = (F[i - 1] + 2ll * (kp[j] - kp[i] + mod) - (kp[j] - kp[j - 1]) + mod) % mod;
for(int l = 1, r; l <= n; l = r + 1) r = n / (n / l), ans = (ans + 1ll * F[n / l] * (G[r] - G[l - 1] + mod)) % mod;
cout << ans << endl;
return 0;
}
XI. P6222 「P6156 简单题」加强版
双倍经验。
5. Dirichlet 卷积
6. 高级筛法
6.1. 线性筛积性函数
前置知识：对线性筛的工作原理有一定认识。
众所周知，欧拉筛在最小质因子处筛去所有合数，但需要讨论：1. 最小质因子次数 \(=1\)。2. 最小质因子次数 \(>1\)。了解两种情况进入的不同分支有助于理解线筛积性函数。
for(int i = 2; i < N; i++) {
if(!vis[i]) pr[++cnt] = i;
for(int j = 1; j <= cnt && i * pr[j] < N; j++) {
vis[i * pr[j]] = 1;
if(i % pr[j] == 0) break; // i * pr[j] 最小质因子 pr[j] 次数 > 1
// i * pr[j] 最小质因子 pr[j] 次数 = 1
}
}
并不是所有积性函数 \(f\) 都可以线筛求出，它需要满足：根据 \(f(p^k)\) 可快速求出 \(f(p^{k+1})\)。当 \(i\times p\) 最小质因子 \(p\) 的次数为 \(1\) 时，方便求得 \(f(ip)=f(i)f(p)\)。但不为 \(1\) 时需要另辟蹊径：若 \(i=p^k\) 则根据 \(f(i)\) 推得 \(f(ip)\) 即 \(f(p^{k+1})\)。 否则设 \(i\) 含有 \(c\) 个质因子 \(p\)，根据积性函数的性质与已经求得的 \(f(p^{c+1})\)（因为 \(i>p^c\)，在 \(i=p^c\) 时已经求过 \(f(p^{c+1})\)）容易得到 \(f(ip)=f\left(\dfrac i{p^c}\right)f(p^{c+1})\)。对此还需记录 \(low_i\) 表示 \(i\) 最小质因子 \(p\) 的 \(p^c\)。
for(int i = 2; i < N; i++) {
if(!vis[i]) pr[++cnt] = i, f[i] = ..., low[i] = i;
for(int j = 1; j <= cnt && i * pr[j] < N; j++) {
vis[i * pr[j]] = 1;
if(i % pr[j] == 0) { // 此时 i 与 pr 不互质, 需要另辟蹊径
low[i * pr[j]] = low[i] * pr[j];
if(i == low[i]) f[i * pr[j]] = ...; // 根据 f(p ^ k) 算 f(p ^ (k + 1))
else f[i * pr[j]] = f[i / low[i]] * f[low[i * pr[j]]]; break;
} f[i * pr[j]] = f[i] * f[pr[j]], low[i * pr[j]] = pr[j]
// i * pr 最小质因子次数 = 1, 可以通过 f[i] * f[pr] 算 f[i * pr] (积性函数的性质)
}
}
