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