ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

从容斥原理到NTT:解析组合计数问题的多项式解法

从容斥原理到NTT:解析组合计数问题的多项式解法 1. 项目概述从一道集训队题目看组合计数的核心思想最近在复盘一些经典的组合数学与多项式题目LOJ3395这道来自2020-2021年集训队的作业题“Yet Another Permutation Problem”又一次吸引了我的注意。题目名字听起来很普通——“又一个排列问题”但内核却一点也不简单它完美地将容斥原理、生成函数以及多项式技术这三块组合计数的基石串联在了一起。很多刚接触组合计数的朋友往往对容斥感觉“很玄”对生成函数觉得“形式复杂”对多项式操作更是望而却步。这道题恰恰是一个绝佳的案例它逼迫你不得不将这些工具融会贯通去解决一个关于排列的、带有约束条件的计数问题。简单来说这类问题通常描述为对于一个长度为 n 的排列我们需要计算满足某种特定约束比如某些位置不能出现某些数字或者某些数字之间必须满足特定的大小关系的排列数量。直接正面计数往往极其困难因为约束条件之间可能相互交织。这时容斥原理就成为了我们手中的利器它的核心思想是“正难则反”通过考虑违反约束条件的情况来间接求出满足所有条件的情况。而生成函数尤其是指数型生成函数则是将组合对象如排列进行代数化编码的优雅工具它能将复杂的组合关系转化为多项式的运算。最后多项式技术如卷积、求逆、快速傅里叶变换FFT/数论变换NTT则为我们高效处理这些代数运算提供了可能。所以这道题的价值远不止于得到一个答案。它更像一个引子引导我们深入理解如何将具体的组合约束抽象成数学模型如何用容斥进行问题转化以及如何用多项式工具进行高效计算。接下来我将结合自己的解题和教学经验详细拆解这道题的解决思路、核心细节以及实现中会遇到的各种“坑”。2. 问题核心与容斥框架的建立2.1 问题重述与约束分析首先我们需要把LOJ3395的题意明确下来原题描述可能较长这里提炼其组合核心。通常这类“Yet Another Permutation Problem”会定义一个由排列构成的集合并对其施加一系列限制。一个典型的限制模式是“禁止模式”或“位置-值”冲突。例如可能规定对于排列 p不允许出现 p[i] j 对于某些特定的 (i, j) 对。或者更复杂的规定某些数字必须出现在某些数字之前。为了进行一般性的讨论我们假设一个简化但不失一般性的模型给定一个 n * n 的 0/1 矩阵 M如果 M[i][j] 1则表示在排列 p 中允许 p[i] j如果 M[i][j] 0则表示禁止 p[i] j。我们需要计算有多少个排列 p 满足对于所有 i都有 M[i][p[i]] 1。这就是一个经典的“带禁位的排列”问题或称“棋盘问题”更专业的叫法是“积和式”计算问题。而积和式的精确计算是#P-Complete的对于大的 n 没有多项式时间算法。因此题目通常会给出矩阵 M 具有特殊的结构使得我们可以利用容斥和生成函数来求解。容斥原理的切入点非常直接。设全集 U 是所有 n! 个排列。对于每个禁止位置 (i, j)即 M[i][j]0我们定义一个坏性质 A_{i,j}排列 p 满足 p[i] j。我们的目标是计算不满足任何坏性质的排列数量即 |⋂ (A_{i,j}的补集)|。根据容斥原理 [ \text{答案} \sum_{S \subseteq \mathcal{B}} (-1)^{|S|} \cdot |\bigcap_{(i,j)\in S} A_{i,j}| ] 其中 (\mathcal{B}) 是所有坏性质的集合。(|\bigcap_{(i,j)\in S} A_{i,j}|) 表示同时满足 S 中所有坏性质的排列数。这相当于在排列中固定了一系列“位置-值”的对应关系。如果 S 中包含的两个坏性质要求同一个位置 i 取两个不同的值 j 和 j‘那么交集显然为空贡献为0。同样如果要求同一个值 j 出现在两个不同的位置 i 和 i’贡献也为0。因此只有那些 S 中所有 (i, j) 两两不同行且不同列时其贡献才非零。这样的一个集合 S 实际上对应了一个“部分匹配”或“部分置换”。2.2 从容斥到生成函数的转化假设我们选取了一个大小为 k 的集合 S {(i1, j1), (i2, j2), ..., (ik, jk)}并且它们两两不同行不同列。这意味着我们已经固定了 k 个位置的值。剩下的 n-k 个位置和 n-k 个值可以任意组成一个排列因此 (|\bigcap_{(i,j)\in S} A_{i,j}| (n-k)!)。于是容斥公式变为 [ \text{答案} \sum_{k0}^{n} (-1)^k \cdot (\text{大小为 k 且两两不同行不同列的禁止位置子集 S 的个数}) \cdot (n-k)! ] 这里的关键在于如何计算“大小为 k 且两两不同行不同列的禁止位置子集 S 的个数”。这等价于在禁止位置构成的二分图左侧 n 个位置右侧 n 个值边对应 M[i][j]0中选取 k 条互不相邻的边的方案数。这正好是二分图匹配的计数问题。生成函数在这里闪亮登场。我们为整个禁止位置矩阵定义一个权重生成函数。设 a_k 为在禁止位置中选取 k 个两两不同行不同列的点的方案数。那么我们可以构造一个关于 x 的普通生成函数 (A(x) \sum_{k0}^{n} a_k x^k)。根据容斥公式最终的答案可以写成 [ \text{答案} \sum_{k0}^{n} (-1)^k \cdot a_k \cdot (n-k)! \sum_{k0}^{n} a_k \cdot (-1)^k (n-k)! ] 这看起来像是 A(x) 的系数与另一个序列的卷积。更优雅的方式是使用指数型生成函数。回忆一下对于排列计数指数型生成函数EGF是更自然的工具。一个数列 {c_n} 的 EGF 是 (C(x) \sum_{n\ge0} c_n \frac{x^n}{n!})。我们定义 F(x) 为“带权”排列的 EGF其中每个排列的权重是 ((-1)^{\text{其包含的禁止位置匹配数}})这个思路有点绕。更标准的做法是利用“匹配多项式”或“棋盘多项式”。实际上有一个著名的结论设 r_k 为在禁止位置中放置 k 个互不攻击的车Rook的方案数这正好就是我们定义的 a_k那么满足禁位条件的排列数由以下公式给出 [ \text{答案} n! \sum_{k0}^{n} \frac{(-1)^k r_k}{(n-k)!} ] 这个形式已经非常接近我们上面的推导了。为了套用多项式技术我们考虑构造两个多项式 [ R(x) \sum_{k0}^{n} r_k x^k ] [ F(x) \sum_{k0}^{n} (-1)^k r_k \cdot (n-k)! ] 注意 F(x) 并不是一个标准的多项式形式因为系数里包含了阶乘。但我们可以将其视为 R(x) 的系数与另一个序列的卷积。观察发现 [ F \sum_{k0}^{n} r_k \cdot (-1)^k (n-k)! \sum_{k0}^{n} r_k \cdot g_{n-k} ] 其中 (g_m (-1)^m m!)。所以 F 实际上是序列 {r_k} 和序列 {g_m} 的卷积的第 n 项如果我们把序列看成多项式系数。但这里 n 是固定的我们最终要求的就是这个卷积结果。因此解题的核心步骤就清晰了根据禁止矩阵 M计算出所有 r_k即大小为 k 的互不攻击的车放置方案数。这对应计算二分图禁位构成的匹配数。利用容斥公式或卷积公式计算出最终答案。难点和算法的精髓全部落在了第一步如何高效计算 r_k 序列。3. 核心算法利用生成函数与多项式计算匹配数3.1 二分图匹配计数的多项式解法对于一个一般的二分图计算其匹配数是经典的#P-Complete问题。但题目中给出的禁位矩阵 M 通常具有特殊的结构否则题目将无法在多项式时间内求解。常见的特殊结构包括禁位构成一个“阶梯状”、“带状”、或者是由几个独立的“矩形”区域组成。这些结构使得其对应的二分图可能是若干个互不连通部分的并集或者具有树形、路径形等简单结构。对于结构简单的二分图其匹配数可以通过动态规划直接计算。但对于稍复杂的结构我们需要更强大的工具。一个关键观察是二分图的匹配数生成函数即匹配多项式可以通过其关联矩阵的积和式Pfaffian或行列式来联系但这通常适用于一般图。对于二分图有一个更直接的方法是利用行列式和容斥的另一种形式。考虑二分图 G(X,Y,E)其中 |X||Y|n。我们构造一个 n x n 的矩阵 B其中 [ B_{i,j} \begin{cases} x, \text{if } (i, j) \text{ 是禁止位置即 } M[i][j]0 \text{} \ 0, \text{otherwise} \end{cases} ] 这里 x 是一个形式变量。那么矩阵 B 的行列式 (\det(B)) 展开后每一项对应一个排列 σ 和其符号以及因子 (x^{\text{该排列中属于禁止位置的个数}})。但这并不是我们想要的因为行列式要求 σ 是一个完整的排列而我们想要的是部分匹配允许有些点不匹配。实际上计算部分匹配数的经典方法是利用图的邻接矩阵和塔特Tutte矩阵的思想但对于二分图可以简化。我们真正需要的是二分图的匹配生成函数(M_G(x) \sum_{k0}^{n} m_k x^k)其中 m_k 是大小为 k 的匹配数即我们的 r_k。计算 M_G(x) 有一个基于高斯消元/行列式的优美方法但它通常适用于一般图且计算复杂度较高约 O(n^3 * 2^n) 对于所有 k。在算法竞赛中更实用的方法是基于状态压缩的动态规划或利用生成函数的卷积。如果二分图可以分解成若干条链或者若干独立的部分那么整个图的匹配生成函数就是各部分匹配生成函数的卷积。因为不同部分之间的匹配是独立的。设图 G 可以分解为两个不相交的部分 G1 和 G2那么 G 中选一个 k-匹配的方案数等于在 G1 中选 i-匹配在 G2 中选 (k-i)-匹配的方案数之和对所有 i 求和。这正是多项式乘法的定义 [ M_G(x) M_{G1}(x) \cdot M_{G2}(x) ]因此如果题目中的禁位矩阵具有分块对角结构或者每一行/每一列的禁位是连续的区间形成一条链那么问题就简化为先计算每个简单子结构的匹配生成函数然后将它们卷积起来。3.2 子结构匹配生成函数的计算我们来看几个典型子结构的匹配生成函数如何手算或简单DP得出。情况一孤立点无边。对应的匹配生成函数是 (1 0*x)因为只能选 0-匹配方案数为1无法选1-匹配。情况二一条边。连接左点 i 和右点 j 的一条边。匹配生成函数是 (1 1*x)。解释不选这条边是 0-匹配1种方案选这条边是 1-匹配1种方案。情况三一条长度为 L 的链。这通常对应于一维的禁位区间。例如左点 i 禁止匹配右点区间 [a, b] 内的所有点并且这些禁止关系是连续的。更精确的模型是左点集 {1,2,...,L} 和右点集 {1,2,...,L}禁止边为 (i, i), (i, i1) 或其他特定模式。计算链的匹配生成函数可以用线性 DP。设 dp[i][0/1] 表示考虑前 i 个左点第 i 个左点是否被匹配的方案生成函数是一个多项式。转移时考虑第 i 个左点如果不匹配则从 dp[i-1][0] 和 dp[i-1][1] 转移如果匹配则需要枚举它能匹配的右点根据禁位规则并从 dp[i-1][0] 转移因为 i-1 必须未匹配否则右点冲突。最终将 dp[L][0] 和 dp[L][1] 相加即可得到整个链的匹配生成函数。这个 DP 的复杂度是 O(L^2)因为多项式乘法卷积在 DP 过程中进行。情况四一个完全二分图 K_{a,b} 去掉一些边。这对应一个矩形禁位区域。计算其匹配生成函数更为复杂通常需要更一般的二分图匹配计数算法。在LOJ3395这道题中经过对题目条件的深入分析通常需要仔细阅读原题描述禁位结构往往可以归结为若干个独立的链或者顺序结构。这是出题人为了考察生成函数与多项式卷积而设计的经典套路。例如条件可能是“对于每个 ip[i] 不能等于 i, i1, ..., id”。这就形成了一个“带状”禁位每一行的禁位是一个连续区间。这种结构可以通过巧妙的转化将其分解为相互独立的链式结构。3.3 多项式卷积与答案计算假设我们已经通过分析将原问题分解为 m 个独立的子结构并计算出了每个子结构的匹配生成函数 (P_1(x), P_2(x), ..., P_m(x))每个都是多项式其次数不超过该子结构的大小。那么整个禁位图的匹配生成函数 (R(x) \prod_{i1}^{m} P_i(x))。计算这个乘积就需要用到多项式乘法卷积。由于 n 可能很大10^5 量级我们必须使用快速傅里叶变换FFT或数论变换NTT来加速多项式乘法。这里通常采用分治策略将所有子多项式两两合并使用 NTT 进行卷积。得到 (R(x) \sum_{k0}^{n} r_k x^k) 后我们代入容斥公式计算答案 [ \text{ans} \sum_{k0}^{n} (-1)^k \cdot r_k \cdot (n-k)! ] 这个计算本身是 O(n) 的但前提是我们已经得到了所有 r_k。注意这里 n 可能很大直接计算 n! 和 (n-k)! 需要预处理阶乘和逆元。然而这个公式还可以进一步用多项式操作来优雅地表示。考虑构造两个多项式 [ A(x) \sum_{k0}^{n} r_k x^k \quad \text{(这就是 R(x))} ] [ B(x) \sum_{j0}^{n} (-1)^j j! \cdot x^j ] 那么观察容斥公式 [ \text{ans} \sum_{k0}^{n} r_k \cdot [x^{n-k}] B(x) [x^n] (A(x) \cdot B(x)) ] 这里 ([x^m]F(x)) 表示多项式 F(x) 中 (x^m) 项的系数。注意在 (A(x)B(x)) 的乘积中(x^n) 项的系数正是 (\sum_{k0}^{n} r_k \cdot b_{n-k})其中 (b_j (-1)^j j!)。这正是我们需要的。因此最终答案就是多项式 (A(x)) 和 (B(x)) 卷积后的第 n 次项系数。这给了我们一个统一的算法框架根据禁位结构分解并计算出每个独立子结构的匹配生成函数 (P_i(x))。使用分治 NTT 计算 (A(x) \prod P_i(x))。预处理阶乘构造多项式 (B(x) \sum_{j0}^{n} (-1)^j j! x^j)。计算 (C(x) A(x) * B(x))多项式乘法。输出 (C(x)) 中 (x^n) 项的系数。这个框架将组合意义的容斥完美地转化为了多项式的卷积运算体现了生成函数将组合问题“代数化”的强大力量。4. 实战实现细节与NTT应用4.1 子结构匹配生成函数的DP求法我们以一个最常见的子结构为例一条长度为 L 的链其中左点 i 只能匹配右点 i 和 i1或者类似的简单邻接关系。但实际上在禁位问题中往往是“禁止”某些边。我们考虑一个等价的正向模型左点 i 允许匹配右点集合 S_i。为了简化假设 S_i 是连续的比如 S_i {1,2,...,L} \setminus {i}即不能匹配自己。计算这个链的匹配生成函数。定义 dp[i][0/1] 为考虑前 i 个左点且第 i 个左点是否被匹配的情况下匹配方案的生成函数多项式。初始状态dp[0][0] 1一个空多项式常数项为1dp[0][1] 0。对于第 i 个左点i从1到L状态 dp[i][0]第 i 点不匹配那么它可以从前一个点的任意状态转移而来且不增加任何边。所以 dp[i][0] dp[i-1][0] dp[i-1][1]。状态 dp[i][1]第 i 点被匹配那么它必须从前一个点的“未匹配”状态转移而来即 dp[i-1][0]因为前一个点如果已匹配可能会占用第 i 点想匹配的右点。同时第 i 点有 |S_i| 种匹配右点的选择。但这里需要注意右点必须未被前 i-1 个点匹配过。在我们的链式假设和连续禁位下通常第 i 点可匹配的右点集合中可能包含已经被前 i-1 点匹配的右点这需要根据具体禁位规则仔细分析转移系数。更严谨的做法是采用“匹配 DP”的标准形式设 dp[i][mask] 表示考虑前 i 个左点当前右点的占用状态为 mask因为右点也只有 L 个可以用状态压缩。但这样复杂度是 O(L * 2^L)对于 L 稍大就不行。因此对于具有简单拓扑结构如链、区间的二分图其匹配生成函数往往有线性递推关系可以直接用线性常系数递推来求解或者通过构造一个转移矩阵其行列式就是匹配生成函数。这在一些更理论的题目中会出现。在实际编码中如果 L 不大比如 L 20我们可以用状态压缩 DP 暴力求出该子结构的匹配生成函数多项式 P(x)。具体地枚举所有左点的匹配情况检查合法性并统计匹配边数 k然后给 P(x) 的 x^k 项系数加1。这样得到的 P(x) 是精确的。4.2 分治NTT合并多项式假设我们通过上述方法得到了 m 个子多项式 (P_1, P_2, ..., P_m)。我们需要计算 (A(x) \prod_{i1}^{m} P_i(x))。直接顺序乘复杂度是 O(m * n log n)并且会因为多项式长度不断变长而效率低下。标准的方法是分治将多项式列表分成两半分别递归计算乘积然后再用 NTT 将两个结果乘起来。伪代码如下def divide_and_conquer(polys, l, r): if l r: return polys[l] # 单个多项式 mid (l r) // 2 left_poly divide_and_conquer(polys, l, mid) right_poly divide_and_conquer(polys, mid1, r) return multiply_by_ntt(left_poly, right_poly) # NTT 乘法其中multiply_by_ntt是封装好的 NTT 乘法函数。使用 NTT 需要选取合适的模数如 998244353其原根为 3并实现标准的 NTT 变换、点值乘法、逆变换流程。注意多项式乘法后次数会增长需要保留足够的长度通常是两个多项式长度之和减一。实操心得在分治乘的时候要注意多项式长度的管理。每次乘法前可以先将两个多项式裁剪到实际需要的长度即次数避免对很多零系数进行无谓的计算。此外对于非常大的 n如 10^5递归深度为 O(log m)整体复杂度约为 O(n log n log m)是可以接受的。4.3 构造B(x)与最终卷积多项式 (B(x) \sum_{j0}^{n} (-1)^j j! x^j) 的构造很简单直接预处理阶乘 fac[0..n]然后循环生成系数即可。注意系数可能是负数在模运算下要处理成非负整数。然后计算 (C(x) A(x) * B(x))。这里 A(x) 的次数最多为 n因为最多选 n 条边B(x) 的次数为 n所以 C(x) 的次数最多为 2n。我们只关心第 n 次项的系数。所以最终答案就是C[n]。4.4 模运算与实现注意事项这类题目通常要求对一个大质数如 998244353取模。预处理需要预处理阶乘 fac[i] 和阶乘的逆元 inv_fac[i]用于计算组合数或公式中的阶乘。同时需要预处理 NTT 所需的原根及其逆元。负数的处理容斥公式中有 ((-1)^k)。在模 p 下-1 等价于 p-1。所以 ((-1)^k) 在模 p 下等于1如果 k 为偶数或p-1如果 k 为奇数。在计算 B(x) 的系数(-1)^j * j! mod p时要注意先计算j! mod p再根据 j 的奇偶性决定是加还是减。NTT 的长度进行 NTT 卷积时需要将多项式长度扩展到大于等于len(A)len(B)-1的最小的 2 的幂。自己实现的 NTT 函数需要包含这个补零操作。常数优化如果子多项式数量 m 很多但每个多项式都很稀疏只有少数低次项系数非零可以考虑使用普通的多项式乘法O(n^2)来合并这些小多项式直到乘积达到一定规模后再切换为 NTT。这是一种实用的启发式优化。5. 常见问题与调试技巧实录5.1 答案总是0或负数模运算下这是最常遇到的问题。检查容斥符号确认公式中 ((-1)^k) 是否正确实现。特别是在计算 B(x) 系数和最终求和时奇偶性判断不能出错。检查匹配生成函数确认你计算的 r_k 是否是“在禁止位置中选 k 个互不冲突的位置”的方案数。最容易出错的地方在于原问题的约束可能是“允许位置”而你错误地计算了“允许位置”的匹配数。必须明确容斥的对象是“坏事件”违反约束。所以 r_k 应该是从禁止边中选 k 条互不相邻的边的方案数。验证小数据务必用手工或暴力程序对 n 很小的情况比如 n8进行验证。写一个暴力枚举所有排列的程序直接统计满足条件的排列数。然后与你多项式程序的结果对比。这是最可靠的调试方法。5.2 NTT 卷积后结果不对长度扩展确保 NTT 前多项式长度是 2 的幂且至少为lenAlenB-1。不足时要补零。模数与原根确认使用的模数如 998244353和原根如 3是否正确。逆变换时别忘了乘长度的逆元。数组越界多项式系数数组要开足够大通常是while(limit nm) limit 1;数组大小要大于 limit。清零进行多次 NTT 时要确保数组的剩余部分补零的位置是真正的零。可以在每次乘法后手动将多余位置清零。5.3 时间复杂度太高无法通过分析子结构确认你是否将原问题分解为了最简形式。也许有更高效的合并方法或者某些子结构的匹配生成函数有封闭表达式如斐波那契数列可以直接写入无需 DP 计算。分治合并的优化如果子多项式数量极多比如 n 个但每个多项式都是1 x的形式那么它们的乘积就是(1x)^n可以用二项式定理直接计算复杂度 O(n)远优于分治 NTT。要观察子结构是否相同。DP 求子多项式的优化如果每个子结构都需要 DP但结构相同只有大小不同可以考虑先预处理出不同大小 L 对应的生成函数避免重复计算。5.4 对原题条件转化的困惑这是解决此类问题的最大难点。LOJ3395 的原题条件可能需要一些巧妙的转化才能变成独立的链或简单结构。常见的转化技巧包括重新标号对排列的值域或定义域进行重新排序或映射使禁止关系呈现规律。补图思想考虑“允许位置”构成的图有时其补图禁止位置的结构更简单。模型转换将排列约束转化为棋盘上的车放置问题、转化为图的路径问题等。多画图将抽象的约束可视化是找到规律的关键。个人体会解决这类组合多项式问题就像在搭积木。容斥原理提供了设计图公式生成函数是标准化的积木块将组合对象多项式化而 NTT 就是高效的粘合剂快速合并积木。最重要的第一步永远是理解题意正确地将自然语言描述转化为组合模型。这一步错了后面所有精妙的算法都是徒劳。我自己的习惯是在编码前先用最小的例子n3,4在纸上完全演算一遍整个容斥和多项式过程确保思路闭环这能避免很多后期调试的折磨。
返回列表