DNA Sequence POJ - 2778（AC自动机+矩阵快速幂）

It’s well known that DNA Sequence is a sequence only contains A, C, T and G, and it’s very useful to analyze a segment of DNA Sequence，For example, if a animal’s DNA sequence contains segment ATC then it may mean that the animal may have a genetic disease. Until now scientists have found several those segments, the problem is how many kinds of DNA sequences of a species don’t contain those segments.

Suppose that DNA sequences of a species is a sequence that consist of A, C, T and G，and the length of sequences is a given integer n.
Input
First line contains two integer m (0 <= m <= 10), n (1 <= n <=2000000000). Here, m is the number of genetic disease segment, and n is the length of sequences.

Next m lines each line contain a DNA genetic disease segment, and length of these segments is not larger than 10.
Output
An integer, the number of DNA sequences, mod 100000.
Sample Input
4 3
AT
AC
AG
AA
Sample Output
36

题意：
只由 A G C T组成的长度为n序列，有一些组合不能出现，求合法序列的数目。

思路：
这是一道很综合的题目！

1. 建图

在使用给定的字符串跑完AC自动机以后，实际上就可以把Trie树当成一个有向图。
题意可以转换为：从root节点（0号）开始，走n步且不会走到有序列结尾的节点的方案数。

举个栗子，当$n=0$的时候，那0的子树的fail指针全部指向自己，每次都有4种可能，所以答案是 $4^n$

当序列为 AG，CG，所得图为

其中2和4是结尾点，不需要考虑。

2. 状态定义

1. 定义 $F[i][j]$代表在第$i$个点，已经走了$j$步的方案数。
2. 易得 $F[0][0]=1,F[0][1]=F[0][3]=0$
3. 转移方程为 $F[i][j]=ai0\ast F[0][j-1]+ai1\ast F[1][j-1]+ai2\ast F[2][j-1]...$
4. 其中$aij$代表图中从$i$到$j$的边数。如上面那个例子的矩阵为：
5. 2 1 0 1 0
   1 1 0 1 0
   0 0 0 0 0
   1 1 0 1 0
   0 0 0 0 0

3. 矩阵递推

有没有发现，状态转移方程特别像矩阵乘法？

没看出来？
。
。
。
。
。
。
。
。
现在再看呢？

（灵魂画手.jpg…

这个矩阵方程，和上面的递推方程等价。
但表示成矩阵，则加速了这个递推的过程。

4. 总结

所以到上一步就结束了。

1. 我们先用AC自动机建好图，去除结尾点（包括Fail指针能到结尾点的，因为Fail指针对应有结尾点，意味着当前对应后缀为结尾，那当前节点肯定不能用了），
2. 然后算出点之间边数的矩阵，
3. 推出一个状态方程，这个转移方程可以转换为矩阵乘法，
4. 最后再用矩阵快速幂加速递推这个过程。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <queue>

using namespace std;

const int maxn = 1e6 + 7;
const int mod = 1e5;
typedef long long ll;

struct Trie {
    int son[30],val,fail;
}t[maxn];

struct Matrix {
    ll m[101][101];
}M;

ll m;
int n,cnt,id[505];
char s[maxn];

void init() {
    for(int i = 0;i <= cnt;i++) {
        memset(t[i].son,0,sizeof(t[i].son));
        t[i].val = t[i].fail = 0;
    }
    cnt = 0;
    id['A'] = 0;
    id['G'] = 1;
    id['C'] = 2;
    id['T'] = 3;
}

void insert(char *s) {
    int u = 0,len = strlen(s);
    for(int i = 0;i < len;i++) {
        int v = id[s[i]];
        if(!t[u].son[v]) t[u].son[v] = ++cnt;
        u = t[u].son[v];
    }
    t[u].val++;
}

void getFail() {
    queue<int>q;
    q.push(0);t[0].fail = 0;
    while(!q.empty()) {
        int u = q.front();q.pop();
        if(t[t[u].fail].val) {
            t[u].val = 1;
        }
        for(int i = 0;i < 4;i++) {
            int v = t[u].son[i];
            int Fail = t[u].fail;
            if(!v) {
                t[u].son[i] = t[Fail].son[i];continue; //直接让不存在的儿子指向失配位置
            }
            if(Fail == 0 && t[Fail].son[i] == v) t[v].fail = 0;
            else t[v].fail = t[Fail].son[i];
            q.push(v);
        }
    }
}

void build() {
    for(int i = 0;i <= cnt;i++) {
        for(int j = 0;j < 4;j++) {
            if(!t[i].val && !t[t[i].son[j]].val) {
                M.m[i][t[i].son[j]]++;
            }
        }
    }
}

Matrix mul(Matrix A,Matrix B) {
    Matrix C;
    memset(C.m,0,sizeof(C.m));
    for(int i = 0;i <= cnt;i++) {
        for(int j = 0;j <= cnt;j++) {
            C.m[i][j] = 0;
            for(int k = 0;k <= cnt;k++) {
                C.m[i][j] = (C.m[i][j] + A.m[i][k] * B.m[k][j]) % mod;
            }
        }
    }
    return C;
}

Matrix Pow(Matrix A,ll n) {
    Matrix res;
    memset(res.m,0,sizeof(res.m));
    for(int i = 0;i <= cnt;i++) res.m[i][i] = 1;
    while(n) {
        if(n & 1) {
            res = mul(res,A);
        }
        A = mul(A,A);
        n >>= 1;
    }
    return res;
}

int main() {
    init();
    scanf("%d %lld",&n,&m);
    for(int i = 1;i <= n;i++) {
        scanf("%s",s);
        insert(s);
    }
    getFail();
    build();
    Matrix res = Pow(M,m);
    ll ans = 0;
    for(int i = 0;i <= cnt;i++) {
        ans = (ans + res.m[0][i]) % mod;
    }
    printf("%lld\n",ans);
    return 0;
}