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][n]；然后，个人理解，我们在这个时候从每个点都走到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;
}