1. 从一道题看字符串匹配的进阶玩法第一次看到DNA Sequence这个题名很多人会以为是个生物信息学的模拟题结果点进去发现是一道关于字符串模式匹配的算法题。题目大意是这样的给定若干致病DNA片段模式串问长度为n的DNA序列中不包含任何致病片段的序列有多少种。DNA字符集只有A、T、C、G四种n可以大到20亿。这个题最直观的想法是动态规划设dp[i][j]表示长度为i、当前匹配到AC自动机上节点j时的合法方案数转移就是枚举下一个字符沿着AC自动机的fail指针走。但n到了20亿这个量级O(n * 节点数 * 4)的复杂度直接爆炸。这时候就需要祭出矩阵快速幂这个武器——把DP的转移写成矩阵形式然后用快速幂把n次转移压缩到O(log n)次矩阵乘法。我当年第一次做这道题的时候卡了整整两天。不是因为不会AC自动机也不是因为不会矩阵快速幂而是卡在怎么把AC自动机上的状态转移正确地映射成矩阵这个环节上。后来才发现这里面有几个非常容易踩的坑比如fail指针的传递、危险节点的标记、矩阵下标的对应关系等等。这篇文章就把整个思路从头到尾捋一遍把那些文档里不会写的细节全部摊开讲清楚。这道题适合已经掌握AC自动机基础、想进一步学习矩阵快速幂优化DP的读者。如果你连AC自动机的build和query都还没写过建议先去补一下基础否则后面的内容会有点吃力。但如果你已经会AC自动机只是不知道怎么和矩阵快速幂结合那这篇文章应该能帮你把思路彻底打通。2. 为什么朴素的AC自动机DP会超时2.1 先明确DP的状态定义在AC自动机上做DP状态定义非常关键。我们设dp[i][j]表示当前构造的DNA序列长度为i并且在AC自动机上走到了节点j且整个过程中没有经过任何危险节点的方案数。所谓危险节点就是那些代表某个致病DNA片段结尾的节点以及通过fail指针能间接到达这些结尾节点的节点。为什么后者也是危险的因为AC自动机的fail指针指向的是当前字符串的最长后缀如果这个后缀是某个致病片段那当前字符串本身就包含了致病片段自然不合法。转移方程很直接dp[i1][next[j][c]] dp[i][j]其中next[j][c]表示从节点j出发走字符c之后到达的节点如果走不了就沿着fail指针一直跳直到能走或者回到根。前提是next[j][c]不是危险节点。2.2 复杂度分析为什么n2e9直接劝退假设AC自动机的节点数是m最多不超过所有模式串长度之和字符集大小是4。那么每一层DP需要枚举m个节点每个节点枚举4个字符总复杂度是O(n * m * 4)。题目中n最大可以到2×10^9m最多大概100左右取决于模式串的数量和长度。算一下2×10^9 × 100 × 4 8×10^11次操作。这个量级就算给你一台超级计算机也得跑很久更别说普通的评测机了。所以朴素的DP思路虽然正确但在n很大的时候完全不可行。这时候就需要观察DP转移的结构——每一层的转移方式其实是完全相同的都是通过那个固定的next数组来转移。这种转移方式固定、只是重复执行n次的结构天然适合用矩阵快速幂来加速。2.3 矩阵快速幂的切入角度把dp[i]看成一个长度为m的行向量或者列向量那么从dp[i]到dp[i1]的转移可以写成一个m×m的矩阵Mdp[i1] dp[i] * M其中M[j][k]表示从节点j转移到节点k的方案数。因为每个节点走4个字符可能有多个字符到达同一个节点所以M[j][k]的值可能是0、1、2、3、4。那么dp[n] dp[0] * M^n。初始状态dp[0]是长度为0的序列此时在根节点节点0所以dp[0][0] 1其余为0。M^n用矩阵快速幂可以在O(m^3 * log n)的时间内算出来。m100时m^310^6log n≈31总共大约3×10^7次操作完全可以在1秒内跑完。注意这里的矩阵乘法是普通的矩阵乘法不是模2意义下的。因为方案数可能很大通常题目会要求对某个数取模比如100000所以矩阵乘法中要记得取模。3. 构建AC自动机时那些容易翻车的细节3.1 危险节点的标记不能只标记结尾很多人写AC自动机的时候习惯在insert的时候把每个模式串的结尾节点标记为危险。但这样做是不够的。考虑这种情况模式串是AT和TC现在有一个节点代表ATC它的fail指针指向TC的结尾节点。虽然ATC这个节点本身不是任何模式串的结尾但它包含了TC这个致病片段所以它也应该是危险的。正确的做法是在build fail指针的时候如果某个节点的fail指针指向的节点是危险的那么这个节点也要标记为危险。代码大概长这样void build() { queueint q; for (int c 0; c 4; c) { if (next[0][c]) { fail[next[0][c]] 0; q.push(next[0][c]); } else { next[0][c] 0; } } while (!q.empty()) { int u q.front(); q.pop(); // 如果fail指向的节点是危险的当前节点也危险 if (danger[fail[u]]) danger[u] true; for (int c 0; c 4; c) { if (next[u][c]) { fail[next[u][c]] next[fail[u]][c]; q.push(next[u][c]); } else { next[u][c] next[fail[u]][c]; } } } }注意danger[u] | danger[fail[u]]这一行必须在处理u的子节点之前执行否则子节点可能会漏掉这个危险标记的传递。3.2 根节点的自环处理在构建next数组的时候如果根节点节点0没有某个字符的出边通常会把next[0][c]设为0也就是指向自己。这其实是一个自环。这个自环在DP中的含义是如果当前在根节点走了一个字符之后还是回到根节点因为没有匹配到任何模式串的前缀这是合法的。但要注意如果根节点本身被标记为危险比如某个模式串是空串虽然一般不会这么出题那这个自环就不能走了。不过正常情况下根节点不会是危险的所以这个自环可以放心使用。3.3 矩阵下标的对应关系这是最容易出错的地方。AC自动机的节点编号是从0到m-1矩阵的下标也是从0到m-1看起来是一一对应的。但问题是哪些节点应该出现在矩阵里答案是只有非危险节点才应该出现在矩阵中。因为一旦走到危险节点这个方案就非法了不应该继续转移。所以我们在构建矩阵的时候要跳过所有危险节点。但这里有个细节如果我们在矩阵中跳过了危险节点那么矩阵的维度就不再是m×m而是(m - 危险节点数)×(m - 危险节点数)。这时候就需要重新建立节点编号到矩阵下标的映射。我当时的做法是先遍历一遍所有节点把非危险节点重新编号存到一个数组id[]中id[u]表示节点u对应的矩阵下标如果是危险节点id[u] -1。然后构建矩阵的时候只考虑id[u] ! -1的节点。int id[maxn], cnt 0; for (int i 0; i m; i) { if (!danger[i]) id[i] cnt; else id[i] -1; }这样矩阵的维度就是cnt×cnt比m×m小一些也能稍微优化一下常数。3.4 取模的时机矩阵快速幂中矩阵乘法的每一层都要取模。但取模的时机也有讲究。如果是在累加的时候取模可能会因为加法溢出导致结果错误。通常的做法是用long long来累加加完之后再取模。或者用一些技巧比如每次加法之后都判断是否超过模数超过就减去模数。我一般习惯用long long累加然后最后统一取模。这样代码简洁也不容易出错。但要注意如果模数很大比如10^97long long累加可能会溢出这时候就需要用一些更复杂的技巧比如快速乘或者__int128。不过这道题的模数通常是100000比较小用long long完全够用。4. 矩阵快速幂在AC自动机上的具体实现4.1 矩阵的结构设计矩阵的大小是cnt×cnt其中cnt是非危险节点的数量。矩阵M[i][j]表示从非危险节点i转移到非危险节点j的方案数。构建矩阵的代码如下struct Matrix { long long a[105][105]; int n; Matrix(int n) : n(n) { memset(a, 0, sizeof(a)); } Matrix operator*(const Matrix other) const { Matrix res(n); for (int i 0; i n; i) { for (int k 0; k n; k) { if (a[i][k] 0) continue; for (int j 0; j n; j) { res.a[i][j] (res.a[i][j] a[i][k] * other.a[k][j]) % MOD; } } } return res; } };注意这里用了if (a[i][k] 0) continue;这个优化因为矩阵中很多元素是0跳过它们可以节省不少时间。构建矩阵的代码Matrix buildMatrix() { Matrix M(cnt); for (int u 0; u m; u) { if (danger[u]) continue; for (int c 0; c 4; c) { int v next[u][c]; if (danger[v]) continue; M.a[id[u]][id[v]]; } } return M; }这里M.a[id[u]][id[v]]是因为可能有多个字符从u走到v每多一个字符就多一种方案。4.2 快速幂的实现细节矩阵快速幂的写法和普通的整数快速幂类似只是把乘法换成了矩阵乘法Matrix power(Matrix base, long long exp) { Matrix res(base.n); for (int i 0; i base.n; i) res.a[i][i] 1; while (exp) { if (exp 1) res res * base; base base * base; exp 1; } return res; }初始向量是dp[0]其中dp[0][id[0]] 1假设根节点不是危险的其余为0。最终答案是dp[n]中所有元素的和因为长度为n的合法序列可以停在任何一个非危险节点上。4.3 字符集映射的处理DNA的字符集是{A, T, C, G}我们需要把它们映射到0到3。这个映射方式可以随意只要保证一致性就行。我一般用int getChar(char c) { if (c A) return 0; if (c T) return 1; if (c C) return 2; return 3; // G }注意不要搞混顺序否则会导致next数组构建错误最终结果也会错。4.4 完整代码框架把上面的部分拼起来完整的代码框架大概是这样的#include bits/stdc.h using namespace std; const int MOD 100000; const int maxn 105; int next[maxn][4], fail[maxn], danger[maxn]; int id[maxn], m 0; // ... AC自动机的insert和build函数 ... // ... Matrix结构体和快速幂 ... int main() { int n, k; cin n k; for (int i 0; i k; i) { string s; cin s; insert(s); } build(); // 重新编号 int cnt 0; for (int i 0; i m; i) { if (!danger[i]) id[i] cnt; else id[i] -1; } Matrix M buildMatrix(); Matrix Mn power(M, n); long long ans 0; for (int i 0; i cnt; i) { ans (ans Mn.a[id[0]][i]) % MOD; } cout ans endl; return 0; }这个框架看起来简单但实际写的时候有很多细节需要注意比如数组大小、取模、边界条件等等。5. 实测中踩过的坑与排查过程5.1 答案总是比预期大危险节点漏标记我第一次提交的时候答案总是比预期大。比如样例输入是2个模式串n3预期答案是36我跑出来是40。多了4个。排查过程是这样的我先手动枚举了所有长度为3的DNA序列一共4^364种然后逐个检查哪些包含致病片段。发现我多算的那4个序列都是包含了某个致病片段的后缀的序列。也就是说这些序列本身不直接以致病片段结尾但它们的某个后缀是致病片段。这就对应到前面说的危险节点传递问题。我在build的时候没有把fail指针指向的危险节点传递下去导致一些间接危险的节点没有被标记。加上danger[u] | danger[fail[u]]之后答案就对了。经验AC自动机的危险标记一定要在build的过程中传递不能只在insert的时候标记结尾。5.2 矩阵快速幂结果全为0初始向量搞错了还有一次我矩阵快速幂跑出来的结果全是0。检查了半天发现是初始向量搞错了。我一开始把初始向量设成了dp[0][0] 1但这里的0是AC自动机的节点编号不是矩阵的下标。如果节点0恰好是危险的虽然一般不会或者节点0在重新编号后对应的矩阵下标不是0那就会出错。正确的做法是用id[0]来获取根节点对应的矩阵下标Matrix dp(1, cnt); // 1×cnt的矩阵 dp.a[0][id[0]] 1;或者直接在最后取答案的时候用Mn.a[id[0]][i]这样就不需要单独维护初始向量了。5.3 取模导致的负数问题有一次我把模数改成了1000000007结果发现答案有时候是负数。排查后发现是矩阵乘法中累加的时候溢出了。long long的最大值是9×10^18而两个模数接近10^9的数相乘是10^18再加上之前的累加结果很容易就溢出了。解决办法是每次加法之后立即取模或者用__int128来累加。不过这道题的原题模数是100000比较小不会溢出。但如果遇到模数很大的题目就需要注意这个问题。5.4 节点数估算错误导致数组越界AC自动机的节点数最多是所有模式串长度之和。如果模式串有10个每个长度10那最多100个节点。但如果模式串有100个每个长度100那就是10000个节点。这时候矩阵就是10000×10000矩阵乘法就是10^12次操作肯定会超时。所以这道题通常模式串的数量和长度都不会太大节点数一般在100以内。如果遇到节点数很大的情况可能需要考虑其他的优化方法比如稀疏矩阵乘法或者BM线性递推。注意做题之前一定要估算一下节点数的上限确保矩阵乘法的复杂度在可接受范围内。6. 从这道题延伸出的通用套路6.1 AC自动机矩阵快速幂的适用场景这个套路并不是只能用在DNA序列这道题上。它的核心适用场景是在一个AC自动机上做DPDP的转移方式固定但转移次数非常大比如10^9级别。常见的变体包括求长度为n且不包含某些模式串的字符串数量求长度为n且包含至少一个模式串的字符串数量用总数减去不包含的求长度为n且恰好包含k个模式串的字符串数量需要扩展状态这些问题的共同点是状态转移可以用矩阵表示且转移次数很大。只要满足这两个条件就可以考虑用矩阵快速幂来优化。6.2 矩阵维度的优化思路矩阵乘法的复杂度是O(m^3)所以m的大小直接决定了算法的可行性。如果m太大可以考虑以下优化稀疏矩阵乘法如果矩阵中大部分元素是0可以只存储非零元素乘法的时候只计算非零部分。这道题的矩阵通常比较稠密因为每个节点走4个字符大部分节点都有出边所以稀疏优化效果有限。状态压缩如果某些节点的行为完全相同可以把它们合并成一个状态从而减小矩阵维度。BM线性递推如果矩阵的维度还是太大可以考虑用Berlekamp-Massey算法求出答案的线性递推式然后用线性递推的快速算法比如Kitamasa算法来求解。这个方法的复杂度是O(m^2 * log n)比矩阵快速幂的O(m^3 * log n)要快。不过对于这道题来说m通常不超过100矩阵快速幂完全够用不需要这些高级优化。6.3 类似的题目推荐如果你已经掌握了这道题可以尝试以下类似的题目来巩固求长度为n且不包含某些模式串的字符串数量字符集更大求长度为n且包含至少一个模式串的字符串数量求长度为n且恰好包含k个不同模式串的字符串数量需要扩展状态矩阵维度会变大在AC自动机上做概率DP然后用矩阵快速幂求期望这些题目的核心思路都是一样的只是在状态定义和转移方程上有些变化。多练几道就能把这个套路彻底吃透。6.4 一些实用的调试技巧最后分享几个我在调试这类题目时常用的技巧先写一个暴力的DP版本用小数据比如n10跑出答案然后用矩阵快速幂的版本跑同样的数据对比结果。如果一致再逐步增大n观察是否有溢出或越界。把矩阵打印出来检查每一行每一列的含义是否正确。特别是危险节点是否被正确排除自环是否存在。如果答案不对先检查AC自动机的build是否正确。可以手动构造一些简单的模式串比如AA和A看看危险节点是否被正确标记。注意取模的时机和方式避免溢出和负数。这些技巧看起来简单但在实际调试中非常有用。我当年就是靠这些方法一步步把这道题从WA调到AC的。
