UESTC 1709 DNA序列 AC自动机+dp+矩阵快速幂优化

DNA序列

Time Limit: 1000/300MS (Java/Others)     Memory Limit: 65535/65535KB (Java/Others)


Submit  Status
众所周知，DNA序列是仅含有A，C，T和G的序列，分析DNA序列的一段非常有用，例如，如果动物的DNA序列含有片段ATC，则可能意味着可能有 一种遗传性疾病。
到目前为止，科学家已经发现了几个片段，问题是一个物种的DNA序列有多少种不包含这些片段。
假设物种的DNA序列是由A，C，T和G组成的序列，序列长度为给定的整数n。

Input

第一行包含两个整数m(0≤m≤10)，n(1≤n≤2000000000)。 这里，m是遗传疾病段的数量，n是序列的长度。
接下来m行每行表示一个DNA遗传疾病段，并且这些段长度不为零且不大于10。

Output

一个整数，DNA序列数，mod 100000。

Sample input and output

Sample Input Sample Output
4 3                                   36
AT
AC
AG

AA

Source

2017 UESTC Training for Search Algorithm & String

UESTC 1709 DNA序列

My Solution

题意：给出m(0<=m<=10)个模式串(0<len<=10)，用AGTC构造长度为n的字符串，
要求每个串的子串都不出现给定的n个串中的任一个，求满足要求的字符串的个数。

AC自动机+dp+矩阵快速幂优化
因为构成的最终串是由一个字符一个字符添加到字符串尾部构成的，
那么如果一个串的后缀如果恰好是某个给定串的前缀时，这个串就可能最终成为非法串。
用k个给定串建立AC自动机，然后从根节点开始递推，
dpij表示递推到第j个字符当前在自动机上的i号节点时的方案数，如果下一个节点是k，
且不是危险节点，则把dpij加到dp[k][i+1]里，
跑一遍，然后答案就是所有非危险节点的方案数的和(其实危险节点上都是0)。
因为危险节点是给定串的终点或者其后缀节点是危险节点的点，
遍历到危险节点的点上的方案必定是包含了给定串的方案，故不能记录这些。
然后这里n比较大，所以要用矩阵快速幂优化，
邻接矩阵，矩阵i行j列的权值是结点i转移到结点j的方案数,所以可以转换成一个(m*len)*(m*len)的矩阵。
而进行k次转移，从结点i转移到结点j的方案数是这个矩阵的k次幂，这个结论离散数学的图论有。
时间复杂度 O((m*len)^3logn)
空间复杂度 O((m*len)^2)

#include <iostream>#include <cstdio>#include <cstring>#include <queue>using namespace std;typedef long long LL;const int CHAR_SIZE = 4;const int MAX_SIZE = 100 + 8;inline int mp(char ch){    if(ch == 'A') return 0;    if(ch == 'G') return 1;    if(ch == 'C') return 2;    if(ch == 'T') return 3;}struct AC_Machine{    int ch[MAX_SIZE][CHAR_SIZE], danger[MAX_SIZE], fail[MAX_SIZE];    int sz;    void init(){        sz = 1;        memset(ch[0], 0, sizeof ch[0]);        memset(danger, 0, sizeof danger);    }    void _insert(char *s){        int n = strlen(s);        int u = 0, c;        for(int i = 0; i < n; i++){            c = mp(s[i]);            if(!ch[u][c]){                memset(ch[sz], 0, sizeof ch[sz]);                danger[sz] = 0;                ch[u][c] = sz++;            }            u = ch[u][c];        }        danger[u] = 1;    }    void _build(){        queue<int> Q;        fail[0] = 0;        for(int c = 0, u; c < CHAR_SIZE; c++){            u = ch[0][c];            if(u){Q.push(u); fail[u] = 0;}        }        int r;        while(!Q.empty()){            r = Q.front();            Q.pop();            danger[r] |= danger[fail[r]];            for(int c = 0, u; c < CHAR_SIZE; c++){                u = ch[r][c];                if(!u){ch[r][c] = ch[fail[r]][c]; continue; }                fail[u] = ch[fail[r]][c];                Q.push(u);            }        }    }}ac;//n^3*log(m)const LL MOD = 100000;const LL M_SIZE = MAX_SIZE;LL mod(LL &x){    x -= x / MOD * MOD;}struct Matrix{    LL m[M_SIZE][M_SIZE];    LL N; //¾ØÕóµÄ½×Êý    void init(){        memset(m, 0, sizeof m);    }    void setOne(){        init();        for(int i = 0; i < M_SIZE; i++) m[i][i] = 1;    }    Matrix(){        init();    }    Matrix operator*(const Matrix &rhs) const{        Matrix ret;        ret.N = N;        int i, j, k;        for(k = 0; k <= N; k++){            for(i = 0; i <= N; i++){                for(j = 0; j <= N; j++){                    mod(ret.m[i][j] += m[i][k] * rhs.m[k][j]);                }            }        }        return ret;    }    void print(){        for(int i = 0; i <= N; i++){            for(int j = 0; j <= N; j++)                cout << m[i][j] << " ";            cout << endl;        }        cout << endl;    }};Matrix res, b;void quickPow(int index){    res.setOne();    //res.print();    //b.print();    while(index){        if(index&1) res = res * b;        index >>= 1;        b = b * b;    }}char s[MAX_SIZE];int main(){    #ifdef LOCAL    freopen("g.txt", "r", stdin);    //freopen("g.out", "w", stdout);    int T = 1;    while(T--){    #endif // LOCAL    //ios::sync_with_stdio(false); cin.tie(0);    int n, m;    scanf("%d%d", &m, &n);    ac.init();    while(m--){        scanf("%s", s);        ac._insert(s);    }    ac._build();    res.N = ac.sz; b.N = ac.sz;    int i, j, k;    //for(i = 1; i <= n; i++){        for(j = 0; j < ac.sz; j++){            if(ac.danger[j]) continue;            for(k = 0; k < CHAR_SIZE; k++){                if(ac.danger[ac.ch[j][k]]) continue;                b.m[j][ac.ch[j][k]]++;            }        }    //}    //cout << "?" << endl;    quickPow(n);    LL ans = 0;    for(i = 0; i < ac.sz; i++){        mod(ans += res.m[0][i]);    }    printf("%lld\n", ans);    #ifdef LOCAL    cout << endl;    }    #endif // LOCAL    return 0;}

  Thank you!

                                                                                                                                            ------from ProLights