POJ2778 DNA Sequence AC自动机+矩阵快速幂

来源:互联网 发布:win7 删除网络驱动器 编辑:程序博客网 时间:2024/05/17 22:27

题目描述:给出m(m <= 10)个长度不超过10的仅由A T C G组成的模式串,问长度为n(1 <= n <= 2e9)的仅由ATCG组成的文本串中,不出现模式串的有多少个?


       这道题目深层次地利用了AC自动机,但是并没有用到AC自动机中的find(query)函数,而是用到了BFS序建立失败指针走向,以及为字典树中的空节点加上侧向边的BFS函数(或者有的叫getFail函数),所以,在继续读之前,确保已经非常清楚AC自动机为字典树加上侧向边是怎么回事。

      看下面的例子,模式串的集合为{ATC , T},那么建完侧向边之后的字典树应该是这个样子的。每个节点都应该有ATCG四个指针,这四个指针都应该指向别的节点,为了便于区分,用同样颜色的边表示同一个节点发出的侧向边,每条侧向边旁边都注明了这条侧向边是代表的是哪一个字母的指针。

      再次强调一遍,这些花花绿绿的边都是侧向边,不是每个节点fail函数的指向,fail函数只是间接用来辅助构建侧向边

       现在图上有0、1、2、3、4五个节点,每个节点都表示字符串当前的匹配状态,比如如果当前位于2节点,那么现在字符串的最后两位就是AT,想象我们正在构造一个字符串,每次我们向字符串末尾加入ATCG中的一个字母,比如如当前我们位于2状态,如果我们在字符串结尾加入字母A,那么我们就转移到了1状态。

      我们可以看到,对于任意一个状态,加入ATCG中的一个字母之后,都会转移到一个新的状态,我们用一个矩阵M来表示这种可达关系,在矩阵中,如果M[i][j]表示从节点i到节点j有多少种转移方式,比如M[1][0]等于2,因为从状态1可以加入C或G两种方式转到状态0。

      离散数学中学过,假设A等于M矩阵的n次方,那么A[i][j]表示通过n次转移之后,从状态i到状态j有几种转移方法。n次转移,具体的反映出来,就是从一个末尾的状态是i的字符串开始,加入n个字母之后,转移成为末尾状态为j字符串的方法。

     可是,这道题的条件是字符串中不能出现模式串,那么该怎么理解这个条件呢?其实很简单,不出现模式串,就是不出现模式串代表的状态,以上图为例,状态3、4分别代表了末尾是ATC、T的状态,而ATC、T都是模式串,那么我们要做的很简单,就是把矩阵M[i][j]中i为1、3或j为1、3的位置全部置为0,表示这种转移不允许发生。那么,只有状态3、4是不允许出现的吗?不是,其实状态2也不允许出现,因为状态2出现就意味着状态4出现了,所以要把2所在的行列置为0,寻找诸如2状态的具体方法,相信做过一些基础AC自动机题目之后都是会的。

    如此一来,就得到了M矩阵,令A等于M矩阵的n次方,根据上面的分析,A[0][1]就表示从0状态开始,加入n个字母之后末状态是1的情况数,也就是长度为n且不出现模式串且结尾状态为n的字符串的个数,由此一来,最终答案ans = ΣA[0][i]  (0 <= i < 字典树中节点总数),具体实现时,因为n比较大,所以算矩阵的n次方用矩阵快速幂。

#pragma warning(disable:4786)#pragma comment(linker, "/STACK:102400000,102400000")#include<iostream>#include<cstdio>#include<cstring>#include<algorithm>#include<stack>#include<queue>#include<map>#include<set>#include<vector>#include<cmath>#include<string>#include<sstream>#include<bitset>#define LL long long#define FOR(i,f_start,f_end) for(int i=f_start;i<=f_end;++i)#define mem(a,x) memset(a,x,sizeof(a))#define lson l,m,x<<1#define rson m+1,r,x<<1|1using namespace std;const int INF = 0x3f3f3f3f;const int mod = 1e5;const double PI = acos(-1.0);const double eps = 1e-6;const int MAXNODE = 105;const int SIGMA_SIZE = 4;const int MAXS = 12;const int NUM = 105;int MAXN;struct Matrix{LL a[NUM][NUM];void init(){mem(a, 0);for (int i = 0; i < MAXN; i++){a[i][i] = 1;}}}mat;Matrix mul(Matrix a, Matrix b){Matrix ans;for (int i = 0; i < MAXN; i++){for (int j = 0; j< MAXN; j++){ans.a[i][j] = 0;for (int k = 0; k < MAXN; k++){ans.a[i][j] += a.a[i][k] * b.a[k][j];ans.a[i][j] %= mod;}}}return ans;}Matrix qpow(Matrix a, int n){Matrix ans;ans.init();while (n){if (n & 1)       ans = mul(ans, a);a = mul(a, a);n >>= 1;}return ans;}struct AhoCorasickAutomata{int ch[MAXNODE][SIGMA_SIZE];int f[MAXNODE];int last[MAXNODE];int val[MAXNODE];int ok[MAXNODE];int sz;void init(){mem(ch[0], 0);sz = 1;}int idx(char c){if (c == 'A')        return 0;else if (c == 'T')       return 1;else if (c == 'C')       return 2;else                  return 3;}void insert(char * s){int u = 0, n = strlen(s);for (int i = 0; i<n; i++){int c = idx(s[i]);if (!ch[u][c]){mem(ch[sz], 0);val[sz] = 0;ch[u][c] = sz++;}u = ch[u][c];}++val[u];}void print(int j){if (j){printf("%d: %d\n", j, val[j]);print(last[j]);}}void getFail(){queue<int>Q;f[0] = 0;for (int i = 0; i < SIGMA_SIZE; i++){int u = ch[0][i];if (u){f[u] = 0;Q.push(u);last[u] = 0;}}while (!Q.empty()){int r = Q.front();          Q.pop();for (int c = 0; c< SIGMA_SIZE; c++){int u = ch[r][c];if (!u){ch[r][c] = ch[f[r]][c];         continue;}Q.push(u);int v = f[r];f[u] = ch[v][c];last[u] = val[f[u]] ? f[u] : last[f[u]];}}}void inimat(){MAXN = sz;mem(ok, 1);for (int i = 0; i < sz; i++){if (val[i] || last[i])      ok[i] = 0;}for (int i = 0; i < sz; i++){if (!ok[i])      continue;for (int j = 0; j < SIGMA_SIZE; j++){int u = ch[i][j];if (!ok[u])          continue;++mat.a[i][u];}}}}ac;char p[15];int main(){int m, n;while (scanf("%d %d", &m, &n) != EOF){ac.init();for (int i = 1; i <= m; i++){scanf("%s", p);ac.insert(p);}ac.getFail();ac.inimat();Matrix ans;ans = qpow(mat, n);LL res = 0;for (int i = 0; i<MAXN; i++){res = (res + ans.a[0][i]) % mod;}printf("%lld\n", res);}return 0;}


   

0 0