[POJ 2778] DNA Sequence （AC自动机+trie图+矩阵快速幂）

链接

POJ 2778

---

题意

给出一些关键词，问长度为n的不含有这些关键词中任何一个的字符串有多少个。
关键词和字符串均由’A’、’C’、’T’、’G’四个字符构成。

---

题解

一道关于trie图的题目。

首先补充一个关于“矩阵和有向图”的知识：
现在有一个有向图含有n个点，可以用一个n*n的邻接矩阵代表该图节点间的关系。现在问你从i节点到j节点走k步，有多少种走法。
事实上我们算出该矩阵的k次幂矩阵gk，gk[i][j]就是答案，这个结论研究一下矩阵乘法的定义就ok了，或者百度“矩阵乘法”。

然后再强调一下利用AC自动机求多文本匹配时候的过程：
一个直观的想法是，走到trie树的某节点k的时候看nd[k]是否是某模式串的结尾，如果是则计数。但事实上是有漏解的，考虑两个模式串：”acbd”、“ijkabcdabcdabcdeeeee”，显然匹配第二个模式串的时候会漏解。
实际上trie树上的任何一个节点代表的都是某模式串的一个前缀，fail指针的含义是“所有模式串中，与当前前缀的后缀匹配的最长前缀”，相当于是kmp扩展到多串的情况。
所以我们需要对trie上的每个节点求出所有前缀，检查这些前缀中是否含有其它模式串（其实任何一个子串都可以表示成原串某前缀的一个后缀），有则计数。

这道题目也需要做类似的事，在建自动机过程中，需要检查当前前缀的fail，如果nd[fail[i]]为1，那nd[i]当然也要为1。凡nd为1的节点均不可到达，建立trie图后，矩阵快速幂即可求解。

---

代码

#include <cstdio>
#include <cstring>
#include <queue>
#include <iostream>
using namespace std;
typedef long long lint;
const int mod = 100000;
#define maxn (1010)
struct Matrix
{
    int n;
    int a[111][111];
    void clear() { n = 0; memset(a, 0, sizeof(a)); }
    Matrix(int k) { clear(); n = k; }
};
Matrix operator* (const Matrix& a, const Matrix& b)
{
    Matrix tmp = Matrix(a.n);
    for(int i = 0, n = a.n; i < n; i++)
    for(int j = 0; j < n; j++)
    for(int k = 0; k < n; k++)
    {
        tmp.a[i][j] += ((lint)a.a[i][k] * b.a[k][j]) % mod;
        tmp.a[i][j] %= mod;
    }
    return tmp;
}
Matrix mat_pow(Matrix a, int k)
{
    Matrix ans = Matrix(a.n);
    for(int i = 0; i < ans.n; i++)
        ans.a[i][i] = 1;
    while(k)
    {
        if(k & 1)
            ans = ans * a;
        a = a * a;
        k >>= 1;
    }
    return ans;
}
int ID[1<<8];
struct Aho_Corasick
{
    #define degree (4)
    int next[maxn][degree], nd[maxn], fail[maxn];
    int root, L;
    int newnode()
    {
        for(int i = 0; i < degree; i++)
            next[L][i] = -1;
        nd[L] = 0;
        return L++;
    }
    void init()
    {
        L = 0;
        root = newnode();
    }
    void insert(char s[])
    {
        int now = root;
        for(int i = 0, key, sz = strlen(s); i < sz; i++)
        {
            key = ID[s[i]];
            if(next[now][key] == -1)
                next[now][key] = newnode();
            now = next[now][key];
        }
        nd[now] = 1;
    }
    void build()
    {
        queue<int> que;
        fail[root] = root;
        for(int i = 0; i < degree; i++)
        if(next[root][i] == -1)
            next[root][i] = root;
        else
        {
            fail[next[root][i]] = root;
            que.push(next[root][i]);
        }
        while(!que.empty())
        {
            int now = que.front();
            que.pop();
            if(nd[fail[now]]) nd[now] = 1;
            for(int i = 0; i < degree; i++)
            if(next[now][i] == -1)
                next[now][i] = next[fail[now]][i];
            else
            {
                fail[next[now][i]] = next[fail[now]][i];
                que.push(next[now][i]);
            }
        }
    }
    Matrix getans(int step)
    {
        Matrix triegraph = Matrix(L);
        for(int i = 0; i < L; i++)
        for(int j = 0; j < 4; j++)
        if(!nd[next[i][j]]) triegraph.a[i][next[i][j]]++;
        return mat_pow(triegraph, step);
    }
} actree;

char keyword[111];
int main()
{
    ID['A'] = 0, ID['C'] = 1, ID['T'] = 2, ID['G'] = 3;

    int m, n;
    while(cin >> m >> n)
    {
        actree.init();
        while(m--)
        {
            scanf("%s", keyword);
            actree.insert(keyword);
        }
        actree.build();
        Matrix ans = actree.getans(n);
        int o = 0;
        for(int i = 0; i < ans.n; i++)
        {
            o += ans.a[0][i];
            o %= mod;
        }
        cout << o << endl;
    }
    return 0;
}