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

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

http://poj.org/problem?id=2778

题意:

        所有串只由A,C,T,G字符构成。给你m个模板串，然后给你一个n，问你长为n的由ACTG字符构成的所有字符串中有多少个是不包含任意一个模板串的？其中n很大。

分析:本题类似于UVA11468:

http://blog.csdn.net/u013480600/article/details/23294375

        首先用模板建立AC自动机(match数组的),然后从0号节点开始在AC自动机里面走,不可以走到match=1的节点上,问你当走n步的时候有多少种行走方式?(仔细想想是不是这个问题)

        令f[i][j]=x表示当前在i节点，行走了j步(但是不经过单词节点)共有x种行走方式。

        那么有递推公式：

        f[i][n] = a0*f[0][n-1]+a1*f[1][n-1]+a2*f[2][n-1]+…+am*f[m-1][n-1]

        也就是说当走到第n-1步的时候在0,1,2…m-1号点的方法总数可以为走n步到i号点做出贡献，这个贡献度就看从它们上一个节点有多少合法后继字母能走到i号点上了。

        建立AC自动机,对所有的非后缀单词节点求它的合法后继,如果能从i走到j,那么矩阵mat[j][i]++。

        下面是个简单的实例:

        举个例子：{AG, CG} ，首先构造 AC 自动机：

        那么转换方程为:

        上面的方程没有写出f3(n)和f4(n)来,不过程序中的矩阵是包括所有节点的.而且f3(n)=f4(n)=0,因为没有节点能到达3和4节点(3和4是非法节点,我们只算节点到达合法节点).

转换矩阵为:

当n=3时,走3步,所以要求上面矩阵的3次幂:

大致思想就在上面了,下面是代码.

AC代码:

#include<cstdio>
#include<cstring>
#include<queue>
using namespace std;
const int MOD = 100000;
const int maxnode=100+10;
const int sigma_size=4;

struct AC_Automata
{
    int ch[maxnode][sigma_size];
    int match[maxnode];//表示当前单词节点的后缀是否存在单词节点
    int f[maxnode];
    int sz;
    int mp[128];
    void init()
    {
        mp['A']=0,mp['C']=1,mp['T']=2,mp['G']=3;//初始化
        sz=1;
        match[0]=f[0]=0;
        memset(ch[0],0,sizeof(ch[0]));
    }
    void insert(char *s)
    {
        int n=strlen(s),u=0;
        for(int i=0;i<n;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;
    }
    void getFail()
    {
        queue<int> q;
        for(int i=0;i<sigma_size;i++)
        {
            int u=ch[0][i];
            if(u)
            {
                f[u]=0;
                q.push(u);
            }
        }
        while(!q.empty())
        {
            int r=q.front();q.pop();
            for(int i=0;i<sigma_size;i++)
            {
                int u=ch[r][i];
                if(!u)
                {
                    ch[r][i]=ch[f[r]][i];
                    continue;
                }
                q.push(u);
                int v=f[r];
                while(v && ch[v][i]==0) v=f[v];
                f[u]=ch[v][i];
                match[u] |= match[f[u]];
            }
        }
    }
};
AC_Automata ac;
long long z[maxnode][maxnode];
long long mat[maxnode][maxnode];
long long ans[maxnode][maxnode];
void mutiply(long long x[maxnode][maxnode],long long y[maxnode][maxnode])
{
    for(int i=0;i<ac.sz;i++)
    {
        for(int j=0;j<ac.sz;j++)
        {
            z[i][j]=0;
            for(int k=0;k<ac.sz;k++)
                z[i][j] +=x[i][k]*y[k][j];
            z[i][j] %=MOD;
        }
    }
    for(int i=0;i<ac.sz;i++)
        for(int j=0;j<ac.sz;j++)
            y[i][j]=z[i][j];
}
int main()
{
    int m,n;
    while(scanf("%d%d",&m,&n)==2)
    {
        ac.init();
        memset(mat,0,sizeof(mat));
        memset(ans,0,sizeof(ans));
        for(int i=0;i<m;i++)
        {
            char str[20];
            scanf("%s",str);
            ac.insert(str);
        }
        ac.getFail();
        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[ac.ch[i][j]][i]++;
        for(int i=0;i<ac.sz;i++) ans[i][i]=1;//ans就是E单位矩阵
        while(n)//这里没乘列向量(1,0,0...,0),因为最后求res的时候只算第一列的
        {
            if(n&1) mutiply(mat,ans);
            mutiply(mat,mat);
            n>>=1;
        }
        long long res=0;
        for(int i=0;i<ac.sz;i++) res += ans[i][0];//只算第一列的
        printf("%I64d\n",res%MOD);
    }
    return 0;
}