1. 从一道题看字符串匹配的进阶玩法DNA Sequence POJ - 2778 这道题在算法竞赛圈子里算是字符串专题的一座小山头。第一次见到它的人多半会愣一下给一堆致病DNA片段问长度为n的DNA串里有多少个不包含任何致病片段的串。n可以大到20亿致病片段数量不超过10个每个片段长度不超过10。这个数据规模一摆出来暴力枚举4的n次方直接出局动态规划按位推也会因为n太大而超时。真正上手做过之后你会发现它把AC自动机和矩阵快速幂这两个看起来八竿子打不着的工具拧在了一起形成了一套非常漂亮的解题范式。这篇文章适合已经了解KMP或者Trie树基础、想往AC自动机方向进阶的读者也适合刷题刷到瓶颈、想看看字符串和矩阵怎么联手的同学。我会从为什么这么设计讲起把状态转移矩阵的构造过程一步步拆开再补上实际编码时容易踩的坑。整篇内容基于我自己的做题和教学经验整理代码用C写但思路对任何语言都通用。核心关键词就三个AC自动机、矩阵快速幂、状态转移。这三个词贯穿全文理解了它们之间的关系这道题就通了。2. 问题本质与方案选型为什么是AC自动机加矩阵快速幂2.1 暴力思路为什么走不通先明确题意。给定m个致病DNA片段每个片段由A、T、C、G四种字符组成长度不超过10。问长度为n的DNA串中有多少个不包含任何致病片段作为子串。n最大到20亿m不超过10。最直接的想法是枚举所有长度为n的串逐个检查是否包含致病片段。4的20亿次方这个数字大到没有意义直接放弃。退一步用动态规划设dp[i][j]表示长度为i、当前匹配到AC自动机上节点j时的合法串数量。转移就是dp[i1][next] dp[i][j]前提是next节点不是致病终止节点。这个DP的时间复杂度是O(n × 节点数 × 4)节点数最多是m×10100左右n是20亿算下来是8×10的11次方级别照样超时。问题的瓶颈在于n太大而状态数很少。这种“状态少、步数多”的结构天然适合用矩阵快速幂来加速线性递推。2.2 AC自动机在这里扮演什么角色AC自动机本质上是一棵Trie树加上失配指针。它解决的问题是在一个文本串中同时匹配多个模式串。放到这道题里我们把所有致病片段插入Trie树构建失配指针然后标记出所有“危险节点”——也就是那些本身是某个致病片段结尾的节点以及通过失配指针能跳到危险节点的节点。为什么需要AC自动机而不是普通的Trie树因为DNA串在逐字符生成的过程中匹配位置不是简单地沿着Trie树往下走。当当前字符在Trie树中没有对应子节点时需要沿着失配指针回跳找到最长的可行后缀继续匹配。这正是AC自动机的核心价值。没有失配指针我们就无法正确判断“当前生成的串是否已经包含了某个致病片段”。注意标记危险节点时不能只标记模式串结尾节点。如果一个节点的失配指针指向危险节点那么这个节点也是危险的。这一步必须在BFS构建失配指针的过程中同步完成否则会漏判。2.3 矩阵快速幂为什么能加速DP的转移方程可以写成状态向量乘以一个转移矩阵得到下一步的状态向量。设状态向量为V其中V[j]表示当前在节点j的合法串数量。转移矩阵M是一个size×size的矩阵M[j][k]表示从节点j能否一步转移到节点k通过某个字符能则值为1否则为0。但这里要注意从节点j通过不同字符可能转移到同一个节点k所以M[j][k]的值应该是转移方案数可能是0、1、2、3、4。那么长度为n的串对应的状态向量就是初始向量乘以M的n次方。矩阵快速幂把O(n)的乘法次数降到O(log n)每次矩阵乘法是O(size³)size最多100左右log n约31总算下来大约3×10的7次方完全可以在时限内跑完。这个方案的精妙之处在于它把“逐位递推”转化成了“矩阵幂运算”利用了矩阵乘法满足结合律这一性质。你可以类比成从A点到B点走n步如果每步的走法固定那总走法就是转移矩阵的n次方。这个类比虽然简化了细节但核心思想是一致的。3. 核心细节拆解从Trie树到转移矩阵的完整链路3.1 Trie树的构建与节点编号先定义节点结构。每个节点有四个子节点指针对应A、T、C、G一个失配指针fail一个布尔标记danger表示该节点是否危险。struct Node { int next[4]; int fail; bool danger; Node() { memset(next, -1, sizeof(next)); fail 0; danger false; } };插入模式串时从根节点出发逐字符往下走没有子节点就新建。最后一个字符对应的节点标记dangertrue。这里有一个细节字符到索引的映射。A→0T→1C→2G→3用一个map函数转换即可。节点总数最多是m×L1m≤10L≤10所以最多101个节点。这个规模非常小矩阵乘法的开销完全可以接受。3.2 失配指针的构建与危险标记传播失配指针的构建用BFS。根节点的所有存在的子节点入队它们的fail指向根。然后逐层处理对于当前节点u遍历它的四个子节点。如果子节点v存在那么v的fail等于u的fail对应的那个子节点如果存在的话否则继续沿着fail往上找。如果子节点不存在通常有两种处理方式一种是直接跳过另一种是补全成Trie图。这道题里我推荐补全成Trie图因为后续构造转移矩阵时会更方便。补全Trie图的做法是对于节点u的字符c如果next[c]不存在就令next[c] fail的next[c]。这样每个节点都有四个出边转移矩阵的构造就变成了直接遍历四个字符。危险标记的传播在BFS过程中如果节点v的fail指向的节点是danger的那么v也是danger的。这一步必须在BFS中完成因为fail指针是逐层构建的处理到v时fail指向的节点的danger标记已经确定。void build() { queueint q; for (int c 0; c 4; c) { if (nodes[0].next[c] ! -1) { nodes[nodes[0].next[c]].fail 0; q.push(nodes[0].next[c]); } else { nodes[0].next[c] 0; } } while (!q.empty()) { int u q.front(); q.pop(); if (nodes[nodes[u].fail].danger) nodes[u].danger true; for (int c 0; c 4; c) { int v nodes[u].next[c]; if (v ! -1) { nodes[v].fail nodes[nodes[u].fail].next[c]; q.push(v); } else { nodes[u].next[c] nodes[nodes[u].fail].next[c]; } } } }提示补全Trie图之后原来的子节点指针会被覆盖。如果你后续还需要区分“真实子节点”和“补全边”需要额外记录。但在这道题里不需要区分因为转移矩阵只关心“从当前节点走某个字符会到哪个节点”补全边正好满足这个需求。3.3 转移矩阵的构造逻辑转移矩阵的维度等于非危险节点的数量。为什么排除危险节点因为一旦进入危险节点就意味着当前生成的串已经包含了某个致病片段后续无论怎么走都不合法。所以危险节点不参与状态转移。构造方法对于每个非危险节点u遍历四个字符c找到next[c]对应的节点v。如果v不是危险节点那么矩阵M[u][v]加1。注意这里u和v都是重新编号后的索引不是原始节点编号。// 收集所有非危险节点 vectorint valid; for (int i 0; i nodeCount; i) { if (!nodes[i].danger) valid.push_back(i); } int sz valid.size(); // 建立原始编号到矩阵索引的映射 vectorint id(nodeCount, -1); for (int i 0; i sz; i) id[valid[i]] i; // 构造矩阵 Matrix M(sz); for (int i 0; i sz; i) { int u valid[i]; for (int c 0; c 4; c) { int v nodes[u].next[c]; if (!nodes[v].danger) { M.a[i][id[v]]; } } }这里有一个容易忽略的点M[i][j]的值可能大于1。因为从节点u走不同字符可能到达同一个节点v。比如节点u的A边和T边都指向v那么M[u][v]就是2。这个计数必须准确否则最终答案会偏小。3.4 矩阵快速幂的实现细节矩阵快速幂的核心是二分幂。设矩阵为M要求M的n次方。当n为偶数时M^n (M^(n/2))²当n为奇数时M^n M × M^(n-1)。递归或迭代实现都可以。struct Matrix { int sz; long long a[105][105]; Matrix(int s) : sz(s) { memset(a, 0, sizeof(a)); } Matrix operator*(const Matrix other) const { Matrix res(sz); for (int i 0; i sz; i) for (int k 0; k sz; k) { if (a[i][k] 0) continue; for (int j 0; j sz; j) { res.a[i][j] (res.a[i][j] a[i][k] * other.a[k][j]) % MOD; } } return res; } }; Matrix power(Matrix base, int exp) { Matrix res(base.sz); for (int i 0; i base.sz; i) res.a[i][i] 1; while (exp 0) { if (exp 1) res res * base; base base * base; exp 1; } return res; }取模的模数通常是100000因为POJ 2778的答案要求对100000取模。这个模数不是质数但矩阵快速幂不涉及除法所以无所谓。注意矩阵乘法中的累加要用long long因为两个小于100000的数相乘可能达到10的10次方超出int范围。取模操作要在每次累加后立即执行避免溢出。3.5 初始向量与最终答案初始状态是空串对应根节点。所以初始向量V是一个长度为sz的向量其中根节点对应的索引位置为1其余为0。最终答案是V乘以M的n次方后所有分量的和。为什么是求和因为长度为n的合法串可能停在任何一个非危险节点上每一种停法都对应一个合法串。所以把所有非危险节点的计数加起来就是答案。Matrix result power(M, n); long long ans 0; for (int i 0; i sz; i) { ans (ans result.a[id[0]][i]) % MOD; }这里result.a[id[0]][i]表示从根节点出发走n步到达节点i的方案数。把所有i加起来就是总方案数。4. 完整实操流程与编码实现4.1 整体代码框架把前面的模块拼起来整体流程是读入m和n读入m个模式串构建Trie树构建失配指针并标记危险节点收集非危险节点并构造转移矩阵矩阵快速幂求和输出。#include cstdio #include cstring #include queue #include vector using namespace std; const int MOD 100000; const int MAXN 105; int charId(char c) { if (c A) return 0; if (c T) return 1; if (c C) return 2; return 3; } struct Node { int next[4]; int fail; bool danger; Node() { memset(next, -1, sizeof(next)); fail 0; danger false; } } nodes[MAXN]; int nodeCount 1; void insert(const char* s) { int u 0; for (int i 0; s[i]; i) { int c charId(s[i]); if (nodes[u].next[c] -1) { nodes[u].next[c] nodeCount; } u nodes[u].next[c]; } nodes[u].danger true; } void build() { queueint q; for (int c 0; c 4; c) { if (nodes[0].next[c] ! -1) { nodes[nodes[0].next[c]].fail 0; q.push(nodes[0].next[c]); } else { nodes[0].next[c] 0; } } while (!q.empty()) { int u q.front(); q.pop(); if (nodes[nodes[u].fail].danger) nodes[u].danger true; for (int c 0; c 4; c) { int v nodes[u].next[c]; if (v ! -1) { nodes[v].fail nodes[nodes[u].fail].next[c]; q.push(v); } else { nodes[u].next[c] nodes[nodes[u].fail].next[c]; } } } } struct Matrix { int sz; long long a[MAXN][MAXN]; Matrix(int s 0) : sz(s) { memset(a, 0, sizeof(a)); } Matrix operator*(const Matrix o) const { Matrix res(sz); for (int i 0; i sz; i) for (int k 0; k sz; k) { if (a[i][k] 0) continue; for (int j 0; j sz; j) { res.a[i][j] (res.a[i][j] a[i][k] * o.a[k][j]) % MOD; } } return res; } }; Matrix power(Matrix base, int exp) { Matrix res(base.sz); for (int i 0; i base.sz; i) res.a[i][i] 1; while (exp 0) { if (exp 1) res res * base; base base * base; exp 1; } return res; } int main() { int m, n; while (scanf(%d%d, m, n) ! EOF) { nodeCount 1; for (int i 0; i MAXN; i) nodes[i] Node(); char buf[15]; for (int i 0; i m; i) { scanf(%s, buf); insert(buf); } build(); vectorint valid; for (int i 0; i nodeCount; i) if (!nodes[i].danger) valid.push_back(i); int sz valid.size(); vectorint id(nodeCount, -1); for (int i 0; i sz; i) id[valid[i]] i; Matrix M(sz); for (int i 0; i sz; i) { int u valid[i]; for (int c 0; c 4; c) { int v nodes[u].next[c]; if (!nodes[v].danger) { M.a[i][id[v]]; } } } Matrix R power(M, n); long long ans 0; for (int i 0; i sz; i) { ans (ans R.a[id[0]][i]) % MOD; } printf(%lld\n, ans); } return 0; }4.2 关键步骤的参数计算与选择节点数上限的计算m≤10每个模式串长度≤10所以Trie树最多有10×101101个节点。矩阵维度最多101矩阵乘法的复杂度是101³≈10的6次方log n≈31总运算量约3×10的7次方在POJ的时限内绰绰有余。模数的选择题目要求对100000取模。这个模数比较小矩阵元素在累加过程中容易超过int范围所以用long long存储。每次乘法后立即取模保证数值始终在可控范围内。n的范围n可以到20亿用int存储没问题因为int上限约21亿。但为了保险可以用long long读入n再转成int传给power函数。如果n超过int范围需要改用long long版本的快速幂但POJ 2778的n不会超过int。4.3 实操现场记录从WA到AC的调试过程我第一次写这道题的时候WA了三次。第一次是因为忘记传播危险标记只标记了模式串结尾节点导致一些通过失配指针能跳到危险节点的节点没有被排除。第二次是因为矩阵乘法没有用long long中间结果溢出。第三次是因为初始向量处理错了把根节点的索引搞混了。调试的时候我建议先写一个暴力版本用小数据对拍。比如n≤10的时候暴力枚举所有4的n次方个串逐个检查是否包含致病片段和矩阵版本的结果对比。这样能快速定位是AC自动机的问题还是矩阵的问题。还有一个实用的调试技巧把转移矩阵打印出来人工检查几行。比如根节点走A、T、C、G分别到哪个节点这些转移是否符合预期。如果发现某个转移指向了危险节点但矩阵里还有值说明危险标记有问题。提示POJ的输入可能有多组测试数据所以要用while(scanf(...)!EOF)循环处理。每次循环都要重置nodes数组和nodeCount否则上一组的数据会污染下一组。5. 常见问题与排查技巧实录5.1 危险标记传播遗漏这是最常见的错误。很多人只把模式串结尾节点标记为danger忘记了在BFS构建失配指针时传播危险标记。考虑一个例子模式串是AT和TC。Trie树中AT的结尾节点是危险的TC的结尾节点也是危险的。但是节点T本身不是任何模式串的结尾然而从T出发走C会到达TC的结尾所以T节点在匹配过程中如果遇到了C就会形成致病片段。更隐蔽的是如果某个节点的失配指针指向危险节点那么这个节点也是危险的。正确的做法是在BFS的每一层检查当前节点的fail指向的节点是否danger如果是把当前节点也标记为danger。这一步必须在处理当前节点的子节点之前完成因为子节点的fail可能指向当前节点。5.2 矩阵维度与节点编号的映射错误收集非危险节点后需要建立原始节点编号到矩阵索引的映射。这个映射容易搞混。比如原始节点编号是0、2、5、7矩阵索引是0、1、2、3。构造矩阵时M[i][j]中的i和j都是矩阵索引不是原始编号。如果直接用原始编号去访问矩阵会越界或者错位。我的习惯是单独开一个vector idid[原始编号]矩阵索引。构造矩阵和最终求和时都用这个映射转换。这样不容易出错。5.3 矩阵乘法溢出与取模时机矩阵乘法的三重循环中res.a[i][j] a[i][k] * o.a[k][j]。a[i][k]和o.a[k][j]都是小于100000的数乘积小于10的10次方用int存储会溢出。必须用long long。而且每次累加后要立即取模否则多次累加后long long也可能溢出。另一个细节是如果a[i][k]为0可以跳过内层循环节省时间。这个优化在矩阵稀疏的时候效果明显但在本题中矩阵比较稠密优化效果有限。5.4 常见问题速查表问题现象可能原因排查方法解决方案答案偏大危险节点未完全排除打印危险节点列表检查是否有遗漏在BFS中传播危险标记答案偏小矩阵转移计数错误打印转移矩阵检查每行之和确保M[i][j]累加而非赋值运行超时矩阵维度太大或快速幂实现低效打印矩阵维度检查循环跳过零元素用long long但及时取模结果不稳定多组数据未重置检查每次循环是否重置nodes和nodeCount在循环开头重新初始化编译错误数组越界或类型不匹配检查MAXN是否够大节点数上限101MAXN设为105足够5.5 独家避坑技巧第一个技巧在构造转移矩阵之前先把所有危险节点打印出来人工确认一遍。特别是那些通过失配指针传播得到的危险节点很容易漏掉。我通常会写一个debug函数输出每个节点的编号、fail指针、danger标记然后对着Trie树图检查。第二个技巧矩阵快速幂的指数n如果是0应该返回单位矩阵。虽然本题n≥1但养成处理边界条件的习惯没坏处。单位矩阵的对角线为1其余为0。第三个技巧如果答案要求对100000取模而最终结果是0不要惊讶。这是完全可能的因为模数很小。不要因为结果是0就怀疑代码有问题。第四个技巧POJ的评测机比较老cin/cout可能比scanf/printf慢。建议用scanf读入字符串和整数printf输出结果。字符串用char数组而不是string避免不必要的开销。6. 从这道题延伸出去的思考这道题最值得玩味的地方是它展示了“状态机矩阵加速”这个范式的通用性。AC自动机本质上是一个有限状态自动机状态就是Trie树上的节点转移就是字符驱动的边。一旦把问题抽象成状态机上的路径计数矩阵快速幂就自然而然地登场了。同样的思路可以迁移到很多问题上。比如给定一个图问从起点走n步到终点的方案数这就是邻接矩阵的n次方。再比如某些动态规划问题中状态转移只依赖于前一个状态且转移规则固定也可以用矩阵快速幂加速。识别这类问题的关键特征是状态数少、步数多、转移规则不随步数变化。我在实际做题时遇到n特别大而状态数在100以内的DP第一反应就是能不能写成矩阵。这个直觉帮我解决了不少看起来很难的题目。当然矩阵快速幂不是万能的如果转移规则随步数变化或者状态数太大就得另想办法。另一个值得注意的点是AC自动机的Trie图补全。很多教程在构建失配指针时只处理存在的子节点不补全不存在的边。这样做在匹配文本串时需要沿着fail指针回跳效率较低。补全成Trie图后每个节点都有完整的四个出边匹配时不需要回跳直接走就行。这个优化在本题中尤其重要因为转移矩阵的构造依赖于“从当前节点走某个字符到哪个节点”这个信息补全Trie图让这个信息变得直接可得。最后说一个编码习惯上的建议写这类涉及多个模块的题目时先把每个模块单独测试通过再拼在一起。比如先测试Trie树的插入和查询再测试失配指针的构建再测试矩阵乘法最后组装。这样调试的时候能快速定位问题出在哪个模块。我见过太多人把所有代码写在一起一运行就错然后花几个小时找bug其实如果分模块测试可能十分钟就搞定了。