首页 > 代码库 > POJ 2778 DNA Sequence —— (AC自动机+矩阵快速幂)

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

  距离上次做AC自动机有很久了=。=,以前这题的思路死活看不懂,现在还是觉得很好理解的。

  思路参见:http://blog.csdn.net/morgan_xww/article/details/7834801#。

  我用cnt=1表示这个节点是危险的,然后再匹配fail指针的时候,如果一个节点的前缀是危险的,那么这个节点也是危险的,这么维护即可。

  顺便一提,我以前的AC自动机模板是没有build过程中失配时的nxt指针的(以前是在match的过程中体现),但是失败时候需要的nxt指针又是很好用的,因此以后的模板中在build中新增这个内容(其实上次的AC自动机DP中就已经有了)。

  另外两点可能不是很重要的是:1.我的矩阵模板统一是从1开始的,而这里有0节点;2.在结构体内似乎不能直接初始化字符串= =。

  代码如下(我的代码跑的有点慢。。):

  1 #include <stdio.h>
  2 #include <algorithm>
  3 #include <string.h>
  4 #include <vector>
  5 #include <queue>
  6 using namespace std;
  7 const int MAX_Tot = 100 + 50;
  8 const int mod = 100000;
  9 
 10 int m,n;
 11 
 12 void add(int &a,int b)
 13 {
 14     a += b;
 15     if(a < 0) a += mod;
 16     a %= mod;
 17 }
 18 
 19 struct matrix
 20 {
 21     int e[MAX_Tot][MAX_Tot],n,m;
 22     matrix() {}
 23     matrix(int _n,int _m): n(_n),m(_m) {memset(e,0,sizeof(e));}
 24     matrix operator * (const matrix &temp)const
 25     {
 26         matrix ret = matrix(n,temp.m);
 27         for(int i=1;i<=ret.n;i++)
 28         {
 29             for(int j=1;j<=ret.m;j++)
 30             {
 31                 for(int k=1;k<=m;k++)
 32                 {
 33                     add(ret.e[i][j],1LL*e[i][k]*temp.e[k][j]%mod);
 34                 }
 35             }
 36         }
 37         return ret;
 38     }
 39     matrix operator + (const matrix &temp)const
 40     {
 41         matrix ret = matrix(n,m);
 42         for(int i=1;i<=n;i++)
 43         {
 44             for(int j=1;j<=m;j++)
 45             {
 46                 add(ret.e[i][j],(e[i][j]+temp.e[i][j])%mod);
 47             }
 48         }
 49         return ret;
 50     }
 51     void getE()
 52     {
 53         for(int i=1;i<=n;i++)
 54         {
 55             for(int j=1;j<=m;j++)
 56             {
 57                 e[i][j] = i==j?1:0;
 58             }
 59         }
 60     }
 61 };
 62 
 63 matrix qpow(matrix temp,int x)
 64 {
 65     int sz = temp.n;
 66     matrix base = matrix(sz,sz);
 67     base.getE();
 68     while(x)
 69     {
 70         if(x & 1) base = base * temp;
 71         x >>= 1;
 72         temp = temp * temp;
 73     }
 74     return base;
 75 }
 76 
 77 char way[4] = {A,T,C,G};
 78 struct Aho
 79 {
 80     struct state
 81     {
 82         int nxt[4];
 83         int fail,cnt;
 84     }stateTable[MAX_Tot];
 85 
 86     int find(char c) {for(int i=0;i<4;i++) if(c == way[i]) return i;}
 87 
 88     int size;
 89     
 90     queue<int> que;
 91 
 92     void init()
 93     {
 94         while(que.size()) que.pop();
 95         for(int i=0;i<MAX_Tot;i++)
 96         {
 97             memset(stateTable[i].nxt,0,sizeof(stateTable[i].nxt));
 98             stateTable[i].fail = stateTable[i].cnt = 0;
 99         }
100         size = 1;
101     }
102 
103     void insert(char *s)
104     {
105         int n = strlen(s);
106         int now = 0;
107         for(int i=0;i<n;i++)
108         {
109             char c = s[i];
110             int to = find(c);
111             if(!stateTable[now].nxt[to])
112                 stateTable[now].nxt[to] = size++;
113             now = stateTable[now].nxt[to];
114         }
115         stateTable[now].cnt = 1;
116     }
117 
118     void build()
119     {
120         stateTable[0].fail = -1;
121         que.push(0);
122 
123         while(que.size())
124         {
125             int u = que.front();que.pop();
126             for(int i=0;i<4;i++)
127             {
128                 if(stateTable[u].nxt[i])
129                 {
130                     if(u == 0) stateTable[stateTable[u].nxt[i]].fail = 0;
131                     else
132                     {
133                         int v = stateTable[u].fail;
134                         while(v != -1)
135                         {
136                             if(stateTable[v].nxt[i])
137                             {
138                                 stateTable[stateTable[u].nxt[i]].fail = stateTable[v].nxt[i];
139                                 // 在匹配fail指针的时候顺便更新cnt
140                                 if(stateTable[stateTable[stateTable[u].nxt[i]].fail].cnt == 1)
141                                     stateTable[stateTable[u].nxt[i]].cnt = 1;
142                                 break;
143                             }
144                             v = stateTable[v].fail;
145                         }
146                         if(v == -1) stateTable[stateTable[u].nxt[i]].fail = 0;
147                     }
148                     que.push(stateTable[u].nxt[i]);
149                 }
150                 /*****建立自动机nxt指针*****/
151                 else
152                 {
153                     if(u == 0) stateTable[u].nxt[i] = 0;
154                     else
155                     {
156                         int p = stateTable[u].fail;
157                         while(p != -1 && stateTable[p].nxt[i] == 0) p = stateTable[p].fail;
158                         if(p == -1) stateTable[u].nxt[i] = 0;
159                         else stateTable[u].nxt[i] = stateTable[p].nxt[i];
160                     }
161                 }
162                 /*****建立自动机nxt指针*****/
163             }
164         }
165     }
166     
167     matrix build_matrix()
168     {
169         matrix ans = matrix(size,size);
170         for(int i=0;i<size;i++)
171         {
172             for(int j=0;j<4;j++)
173             {
174                 if(!stateTable[i].cnt && !stateTable[stateTable[i].nxt[j]].cnt)
175                     ans.e[i+1][stateTable[i].nxt[j]+1]++;
176             }
177         }
178         return ans;
179     }
180 }aho;
181 
182 int main()
183 {
184     while(scanf("%d%d",&m,&n) == 2)
185     {
186         aho.init();
187         char s[15];
188         for(int i=1;i<=m;i++)
189         {
190             scanf("%s",s);
191             aho.insert(s);
192         }
193         aho.build();
194         matrix p = aho.build_matrix();
195         p = qpow(p,n);
196         int ans = 0;
197         for(int i=1;i<=aho.size;i++) add(ans, p.e[1][i]);
198         printf("%d\n",ans);
199     }
200     return 0;
201 }

 

 

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