POJ 2778 DNA Sequence (AC自动机+矩阵幂DP)

 看完了ac自动机的基本应用题目，进军自动机dp的道路很艰辛，第一个是个自动机+概率dp（UVA 11468），不算难，但是比起基本自动机来，有了新的方向：就是在bfs找失败指针时，遇到now的不存在的i子节点不是直接continue，而是将该子节点只想fail[now]的i子节点（也不存在即指向了0），也就是所谓的吧所有边都连上。

    然后，这一类题的意思，大概是给了n个模板，问你任意一个长为len的串不能包含这些模板的问题。做法是建完trie树后，build_AC_Automaton()里①连所有边 并 ②维护match[]，然后dp解决相应问题。

    这个题值得讲的地方在矩阵快速幂上，饶齐博客上给配了图，讲了很多但是很混乱，一个图不全，而是节点标号不对（虽说是为了好写好理解但是没有说明，开始我是搞混了的），另外代码mat[][]是反着的，别扭了，以至于在求ans的时候说不清楚，所以代码我选择了笑着走完自己的路的，他博客里也有一些解释，两人结合着看题还是很容易理解的。

饶齐：点击打开链接       笑着：点击打开链接         矩阵快速幂入门：点击打开链接

    定义意思都在饶齐博客上，mat[i][j]反过来，代表i能走到j。建完树后，一遍遍历sz个节点*其4个儿子，另mat[i][ch[i][j]]++。这样也就是求得了f[i][1]=∑mat[1...sz][i]（从任何一点x开始走，走1步到i点，的方案），然后经过（n-1）次，求得f[i]
；然后，个人理解，我们在这个时候从每个点都走到0点，再走一次，也就是n次幂了。这样相当于：从每个点开始走到根（0）点结束；反过来就是从根（0）开始走到每个点结束的结果了。也就可以解释最后对ans【0】【i】求和（笑着版），也就是饶齐的矩阵图上那个n次幂后突然多了一个×[1
0 0]矩阵的原因了。

#include<cstdio>
#include<iostream>
#include<queue>
#include<cstring>
using namespace std;
typedef long long LL;
const int MAXNode=110;
const int SIZE=4;
const int MOD=100000;
int vis[1010];
struct AC_Automaton
{
int ch[MAXNode][SIZE];
int match[MAXNode]; // 代表单词节点
int fail[MAXNode]; // fail数组
int sz;//节点个数
int mp[128];
void init()
{
sz=1;
memset(ch[0], 0, sizeof(ch[0]));
match[0]=fail[0]=0;
mp['A']=0; mp['C']=1;
mp['T']=2; mp['G']=3;
}
void Insert(char *s)
{
int len=strlen(s);
int u=0;
for(int i=0; i<len; i++)
{
int id=mp[s[i]];
if(ch[u][id]==0)
{
ch[u][id]=sz;
memset(ch[sz], 0, sizeof(ch[sz]));
match[sz++]=0;
}
u=ch[u][id];
}
match[u]=1;
}
//build_AC_Automaton函数负责构造fail和last数组
void build_AC_Automaton()
{
queue<int>q;
fail[0]=0;
for(int i=0; i<SIZE; i++)//初始化root的儿子
{
int u=ch[0][i];
if(u)
{
fail[u]=0;//
q.push(u);
}
}
while(!q.empty())// 按BFS顺序计算fail
{
int now=q.front();
q.pop();
for(int i=0; i<SIZE; i++)
{
int u=ch[now][i];//儿子u
if(u==0)//该儿子不存在，不考虑
{
ch[now][i]=ch[fail[now]][i];//连全部的边，Trie图
continue;
}
q.push(u);
int v=fail[now];
while(v && ch[v][i]==0)//一直找失败指针
v=fail[v];
fail[u]=ch[v][i];//求得 fail【u】
match[u]=match[u]|match[fail[u]];
}
}
}
}ac;
struct Matrix
{
LL a[MAXNode][MAXNode];
int n;
};
Matrix mat;//从i走到j
Matrix ans;
Matrix mutiply(Matrix x, Matrix y)
{
Matrix z;
z.n=x.n;
for(int i=0; i<x.n; i++)
{
for(int j=0; j<y.n; j++)
{
z.a[i][j]=0;
for(int k=0; k<z.n; k++)
z.a[i][j]+=x.a[i][k]*y.a[k][j];
z.a[i][j]%=MOD;
}
}
return z;
}
char s[20];
int main()
{
int N, M;
while(scanf("%d%d", &N, &M)==2)
{
ac.init();
memset(mat.a, 0, sizeof(mat.a));
memset(ans.a, 0, sizeof(ans.a));
for(int i=1; i<=N; i++)
{
scanf("%s", s);
ac.Insert(s);
}
ac.build_AC_Automaton();
for(int i=0; i<ac.sz; i++)//求走一遍
{
if(ac.match[i]==0)//非后缀单词
{
for(int j=0; j<4; j++)
{
if(ac.match[ac.ch[i][j]]==0)//非后缀单词
mat.a[ac.ch[i][j]][i]++;
}
}
}
mat.n=ans.n=ac.sz;
for(int i=0; i<ac.sz; i++)
ans.a[i][i]=1;//E单位矩阵
while(M)
{
if(M & 1)
ans=mutiply(mat, ans);
mat=mutiply(mat, mat);
M>>=1;
}
LL sum=0;
for(int i=0; i<ac.sz; i++)
sum+=ans.a[i][0];//只算第一列的
printf("%I64d\n", sum%MOD);
}
return 0;
}