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

题目链接:http://poj.org/problem?id=2778

 

题目大意:

  给出m个字符串,求长度为n且不包含给出字符串的字符串个数

 

题目思路:

  推荐博客:传送门
  这道题一共分为两步,第一步是得到字典树上每个节点走一步到达各个节点的合法方案数(起点和终点不能为非法节点,非法节点的定义下面介绍),这个方案数通过矩阵存储,第i行第j列表示从i节点表示的状态走一步到达第j步表示的状态的方案数。第二步是这个矩阵的n次方,就得到了走n步的情况,然后节点1是根节点,所以1的那一行相加就是答案。
  首先说明第一步的做法。难点在于确定非法节点,非法节点的概念定义很容易想到,就是表示的字符串包含给定字符串的节点。首先很明显,给定字符串在字典树中的结束位置一定是非法节点,但是对于ACG和C的情况,ACG中的C也是非法节点,但是没有标记上,这个时候就需要利用ac自动机的fail指针。ACG中的C的fail指针指向C,C是非法节点,所以ACG中的C也是非法节点。所以这里的处理方法就是反向建立fail指针,然后跑个dfs,如果父亲是非法节点,那么儿子一定都是非法节点,反向建边的定义是说,如果x->y,那么x是y的后缀,x是非法节点,那么y也一定是非法节点,因为y包含x。由于y已经被列为非法节点,所以它后面的节点就无法到达,因为它后面的节点需要先得到y才能获得,而y无法获得,所以一石二鸟,这一部分要做的就是反向建立fail边,只要父亲是非法节点,所有儿子都设为非法节点。然后通过遍历节点,枚举他们通过四个字母能够到达的状态,如果是非法状态就跳过,否则就给对应的情况++,就得到了我们需要的矩阵。
  第二步就是矩阵快速幂的裸题了。基础矩阵是我们需要的矩阵,需要乘的矩阵也是我们需要的矩阵,因为是我们需要的矩阵的n次,基础矩阵已经有一次了,所以他的次方只用n-1次就行。

 

以下是代码:

#include<iostream>
#include<cstdio>
#include<cstring>
#include<queue>
#include<vector>

using namespace std;
#define ll long long
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define per(i,a,b) for(ll i=a;i>=b;i--)
#define ll long long
const int MAXN = 105;
const int MOD = 100000;
char s[105];
int trie[MAXN][5],tot,pos[MAXN],fail[MAXN],num[MAXN];
int a[MAXN][MAXN],f[MAXN][MAXN];
void mul(int f[105][105],int a[105][105]){
    int c[105][105];
    memset(c,0,sizeof(c));
    rep(i,1,tot){
        rep(j,1,tot){
            rep(k,1,tot){
                c[i][j]=(c[i][j]+(ll)f[i][k]*a[k][j])%MOD;
            }
        }
    }
    memcpy(f,c,sizeof(c));
}
void mulself(int a[105][105]){
    int c[105][105];
    memset(c,0,sizeof(c));
    rep(i,1,tot){
        rep(j,1,tot){
            rep(k,1,tot){
                c[i][j]=(c[i][j]+(ll)a[i][k]*a[k][j])%MOD;
            }
        }
    }
    memcpy(a,c,sizeof(c));
}
void Insert(int x){
    int len=strlen(s),p=1;
    rep(i,0,len-1){
        int ch;
        if(s[i]=='A')ch=0;
        if(s[i]=='C')ch=1;
        if(s[i]=='T')ch=2;
        if(s[i]=='G')ch=3;
        if(!trie[p][ch])trie[p][ch]=++tot;
        p=trie[p][ch];
    }
    num[p]++;
}
queue<int>q;
vector<int>v[MAXN];
void dfs(int u,int x){
    int len=v[u].size();
    if(x==1)num[u]=1;
    rep(i,0,len-1){
        int y=v[u][i];
        if(num[u])dfs(y,1);
        else dfs(y,x);
    }
}
int main()
{
    int m,n;
    while(~scanf("%d%d",&m,&n)){
        tot=1;
        memset(trie,0,sizeof(trie));
        memset(num,0,sizeof(num));
        rep(i,1,m){
            scanf("%s",s);
            Insert(i);
        }
        rep(i,0,tot)v[i].clear();
        rep(i,0,3)trie[0][i]=1;
        q.push(1);
        while(!q.empty()){
            int u=q.front();
            q.pop();
            rep(i,0,3){
                if(trie[u][i]){
                    fail[trie[u][i]]=trie[fail[u]][i];
                    q.push(trie[u][i]);
                }
                else{
                    trie[u][i]=trie[fail[u]][i];
                }
            }
        }

        rep(i,2,tot)v[fail[i]].push_back(i);
        dfs(1,0);
        memset(a,0,sizeof(a));
        rep(i,1,tot){
            if(num[i])continue;
            rep(j,0,3){
                int p=trie[i][j];
                if(num[p])continue;
                a[i][p]++;
            }
        }
        rep(i,1,tot){
            rep(j,1,tot)f[i][j]=a[i][j];
        }
        n--;
        for(;n;n>>=1){
            if(n&1)mul(f,a);
            mulself(a);
        }
        int ans=0;
        rep(i,1,tot){
            ans=(ans+f[1][i])%MOD;
        }
        printf("%d\n",ans);
    }
    return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值