DNA Sequence 【POJ - 2778】【AC自动机+矩阵快速幂】

题目链接

---

  哇哇哇！！！过啦，美滋滋，捋了半天题意，不解为何最后答案为0，然后就一直一直的debug，后来发现竟然是在写字典树的时候写彪了……呜呜呜，好气啊QAQ……

  题目的思路倒也却是很有奥妙，也特别的有意思（尤其是想到之后），我们看到有长达"2e9"的步数得去走，就想到了用矩阵快速幂，[i][j]指的是从i到j所需要的路径可能的数量。然后，我们找到N步之后的所有从[0]状态到达的全体点[i]的总和就是我么需要的ans。其中，这里用到的就是矩阵快速幂的优化了。

  那么，用AC自动机求的是什么？我们用AC自动机处理的是状态，所有的状态利用AC自动机来写，把所有冲突的状态在矩阵中赋值为0，其余的能走就直接"++"找齐全体方案数。

  还有啊，这道题就不用释放内存了，因为我们要取的是所有点的方案数，也就是说从一个点出发可以达到的全体{A、C、G、T}点，那么就不会有停止的时候，就不要释放内存操作了，会RE的（亲测，还加了vis记录的，也是会RE），那么，或许会有人闻到为什么要开那么多的无用的状态，其实，这都是有用的，所有表面不能到达的状态实际上在矩阵中得表示为可以达到的状态，这也就是为什么在矩阵快速幂结束后，我们的ans的值得从[0]状态开始向上加了，就是为了得到那些不受限状态的可能路径了。

---

#include <iostream>
#include <cstdio>
#include <cmath>
#include <string>
#include <cstring>
#include <algorithm>
#include <limits>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <map>
#define lowbit(x) ( x&(-x) )
#define pi 3.141592653589793
#define e 2.718281828459045
#define INF 0x3f3f3f3f
#define HalF (l + r)>>1
#define lsn rt<<1
#define rsn rt<<1|1
#define Lson lsn, l, mid
#define Rson rsn, mid+1, r
#define QL Lson, ql, qr
#define QR Rson, ql, qr
#define myself rt, l, r
#define MP(a, b) make_pair(a, b)
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
const ll mod = 100000;
map<char, int> dir;
inline void pre_did() { dir['A'] = 0;   dir['C'] = 1;   dir['G'] = 2;   dir['T'] = 3; }
int M, cnt, N;
struct node
{
    node *next[4];
    node *fail;
    bool flag;
    int id;
    node()
    {
        memset(next, NULL, sizeof(next));
        fail = NULL;
        flag = false;
        id = 0;
    }
};
node *root;
char virus[15];
void update(char *s)
{
    node *temp = root;
    int len = (int)strlen(s);
    for(int i=0; i<len; i++)
    {
        int x = dir[s[i]];
        if(temp->next[x] == NULL)
        {
            temp->next[x] = new node();
            temp->next[x]->id = cnt++;
        }
        temp = temp->next[x];
    }
    temp->flag = true;
}
inline void build_fail()
{
    queue<node *> Q;
    Q.push(root);
    node *temp, *p;
    while(!Q.empty())
    {
        temp = Q.front();   Q.pop();
        for(int i=0; i<4; i++)
        {
            if(temp->next[i])
            {
                if(temp == root) temp->next[i]->fail = root;
                else
                {
                    p = temp->fail;
                    while(p)
                    {
                        if(p->next[i]) { temp->next[i]->fail = p->next[i];  break; }
                        p = p->fail;
                    }
                    if(!p) temp->next[i]->fail = root;
                    if(temp->next[i]->fail->flag) temp->next[i]->flag = true;
                }
                Q.push(temp->next[i]);
            }
            else temp->next[i] = (temp == root? root:(temp->fail->next[i]));
        }
    }
}
struct Multi
{
    ll a[105][105];
    int d;
    Multi(int ff=0)
    {
        memset(a, 0, sizeof(a));
        d = ff;
    }
    friend Multi operator * (Multi e1, Multi e2)
    {
        Multi qq(e1.d);
        for(int i=0; i<e1.d; i++)
        {
            for(int j=0; j<e1.d; j++)
            {
                for(int k=0; k<e1.d; k++)
                {
                    qq.a[i][j] = ( e1.a[i][k] * e2.a[k][j] + qq.a[i][j] ) % mod;
                }
            }
        }
        return qq;
    }
}mp;
Multi fast_mi(Multi e1, int x)
{
    Multi tt = Multi(e1.d);
    for(int i=0; i<cnt; i++) tt.a[i][i] = 1;
    while(x)
    {
        if(x & 1) tt = tt * e1;
        e1 = e1 * e1;
        x >>= 1;
    }
    return tt;
}
bool vis[105];
void Creat_Multi()
{
    node *temp;
    memset(vis, false, sizeof(vis));
    queue<node *> Q;
    Q.push(root);
    vis[root->id] = true;
    while(!Q.empty())
    {
        temp = Q.front();   Q.pop();
        for(int i=0; i<4; i++)
        {
            if(!temp->flag && !temp->next[i]->flag)
            {
                mp.a[temp->id][temp->next[i]->id]++;
            }
            if(!vis[temp->next[i]->id])
            {
                vis[temp->next[i]->id] = true;
                Q.push(temp->next[i]);
            }
        }
    }
}
inline void init()
{
    mp = Multi();
    cnt = 0;
    root = (node *)malloc(sizeof(node));
    for(int i=0; i<4; i++) root->next[i] = NULL;
    root->fail = NULL;
    root->flag = false;
    root->id = cnt++;
}
inline void FREE(node *temp)
{
    for(int i=0; i<4; i++) if(temp->next[i] && !vis[temp->next[i]->id])
    {
        vis[temp->next[i]->id] = true;
        FREE(temp->next[i]);
    }
    free(temp);
}
int main()
{
    pre_did();
    while(scanf("%d%d", &M, &N)!=EOF)
    {
        init();
        for(int i=0; i<M; i++)
        {
            scanf("%s", virus);
            update(virus);
        }
        build_fail();
        mp.d = cnt;
        Creat_Multi();
        ll ans = 0;
        mp = fast_mi(mp, N);
        for(int i=0; i<cnt; i++)
        {
            ans += mp.a[0][i];
            if(ans >= mod) ans -= mod;
        }
        printf("%lld\n", ans);
//        memset(vis, false, sizeof(vis));        //释放内存！别释放了…… 会RE的……
//        vis[root->id] = true;
//        FREE(root);
    }
    return 0;
}

---

静态版本

#include <iostream>
#include <cstdio>
#include <cmath>
#include <string>
#include <cstring>
#include <algorithm>
#include <limits>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <map>
#include <bitset>
//#include <unordered_map>
//#include <unordered_set>
#define lowbit(x) ( x&(-x) )
#define pi 3.141592653589793
#define e 2.718281828459045
#define INF 0x3f3f3f3f
#define HalF (l + r)>>1
#define lsn rt<<1
#define rsn rt<<1|1
#define Lson lsn, l, mid
#define Rson rsn, mid+1, r
#define QL Lson, ql, qr
#define QR Rson, ql, qr
#define myself rt, l, r
using namespace std;
typedef unsigned long long ull;
typedef unsigned int uit;
typedef long long ll;
const ll mod = 1e5;
int N, M, Q, tot;
struct matrice
{
    ll a[105][105];
    matrice() { memset(a, 0, sizeof(a)); }
    friend matrice operator * (matrice x, matrice y)       //重载函数，矩阵的乘积得到的新矩阵
    {
        matrice ans;
        for(int i=0; i<=tot; i++)
        {
            for(int j=0; j<=tot; j++)
            {
                ll tmp = 0;
                for(int k=0; k<=tot; k++)
                {
                    tmp=(tmp + x.a[i][k] * y.a[k][j]) % mod;       //最后得到的需要取mod
                }
                ans.a[i][j] = tmp;
            }
        }
        return ans;
    }
}Bas, Want;
matrice fast_mi(matrice x, ll ti)
{
    matrice ans;
    for(int i=0; i<=tot; i++) ans.a[i][i] = 1;  //构造单位阵
    while(ti)
    {
        if(ti & 1) ans = ans * x;
        x = x * x;
        ti >>= 1;
    }
    return ans;
}
map<char, int> dir;
char s[15];
struct node
{
    int nex[4], fail; bool flag;
    node() { memset(nex, 0, sizeof(nex)); fail = 0; flag = false; }
} a[105];
inline void Insert()
{
    int len = (int)strlen(s), u = 0;
    for(int i=0, id; i<len; i++)
    {
        id = dir[s[i]];
        if(!a[u].nex[id])
        {
            a[++tot] = node();
            a[u].nex[id] = tot;
        }
        u = a[u].nex[id];
    }
    a[u].flag = true;
}
inline void build_fail()
{
    queue<int> Q; Q.push(0);
    int tmp, p, son;
    while(!Q.empty())
    {
        tmp = Q.front(); Q.pop();
        for(int i=0; i<4; i++)
        {
            son = a[tmp].nex[i];
            if(son)
            {
                if(!tmp) a[son].fail = 0;
                else
                {
                    p = a[tmp].fail;
                    while(p && !a[p].nex[i]) p = a[p].fail;
                    a[son].fail = a[p].nex[i];
                }
                Q.push(son);
                if(a[a[son].fail].flag || a[tmp].flag) a[son].flag = true;
            }
            else a[tmp].nex[i] = a[a[tmp].fail].nex[i]; //这一步必须要有
        }
    }
}
int main()
{
    dir['A'] = 0; dir['C'] = 1; dir['G'] = 2; dir['T'] = 3;
    scanf("%d%d", &M, &N);
    for(int i=1; i<=M; i++)
    {
        scanf("%s", s);
        Insert();
    }
    build_fail();
    Bas = matrice();
    for(int i=0; i<=tot; i++)
    {
        for(int j=0; j<4; j++)
        {
            if(!a[i].flag && !a[a[i].nex[j]].flag)  //有可能本身是非法的但是却加一个新值变成合法项
            {
                Bas.a[i][a[i].nex[j]] ++;
            }
        }
    }
    Want = fast_mi(Bas, N);
    int ans = 0;
    for(int i=0; i<=tot; i++)
    {
        ans += Want.a[0][i];
        if(ans >= mod) ans -= mod;
    }
    printf("%d\n", ans);
    return 0;
}
/*
2 3
AT
GT
ans:48
*/