POJ 2778 DNA Sequence (AC自己主动机 + dp)

DNA Sequence

题意：DNA的序列由ACTG四个字母组成，如今给定m个不可行的序列。问随机构成的长度为n的序列中。有多少种序列是可行的（仅仅要包括一个不可行序列便不可行）。个数非常大。对100000取模。

思路：推荐一个博客，讲的很清楚。

这样的题目。n非常大，首先想到的就是用矩阵来优化。那么怎样构造转移方程呢：首先建立一棵Trie，然后依照AC自己主动机的方式构造fail指针，然后会发现。当一个状态分别加入ACTG之后，会得到还有一个状态。

（详细解释见代码）

代码：

/*
ID: wuqi9395@126.com
PROG:
LANG: C++
*/
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<cmath>
#include<cstdio>
#include<vector>
#include<string>
#include<fstream>
#include<cstring>
#include<ctype.h>
#include<iostream>
#include<algorithm>
#define INF (1<<30)
#define PI acos(-1.0)
#define mem(a, b) memset(a, b, sizeof(a))
#define rep(i, n) for (int i = 0; i < n; i++)
#define debug puts("===============")
typedef long long ll;
using namespace std;
const int maxn = 110;
const int maxm = 110;
ll mod = 100000;
struct Matrix {
    int n, m;
    ll a[maxn][maxm];
    void clear() {
        n = m = 0;
        memset(a, 0, sizeof(a));
    }
    Matrix operator * (const Matrix &b) const { //实现矩阵乘法
        Matrix tmp;
        tmp.clear();
        tmp.n = n;
        tmp.m = b.m;
        for (int i = 0; i < n; i++)
            for (int j = 0; j < m; j++) {
                if (!a[i][j]) continue;
                for (int k = 0; k < b.m; k++)
                    tmp.a[i][k] += a[i][j] * b.a[j][k], tmp.a[i][k] %= mod;
            }

        return tmp;
    }
}A, res;

const int maxnode = 11 * 11;
const int charset = 4;
struct ACAutomaton {
    int ch[maxnode][charset];
    int fail[maxnode];
    int Q[maxnode];
    int val[maxnode];
    int sz;
    int id(char ch) {
        if (ch == 'A') return 0;
        else if (ch == 'C') return 1;
        else if (ch == 'T') return 2;
        return 3;
    }
    void init() {
        fail[0] = 0;
        //for (int i = 0; i < charset; i++) ID[i] = i;
    }
    void reset() {
        sz = 1;
        memset(ch[0], 0, sizeof(ch[0]));
    }
    void Insert(char* s, int key) {
        int u = 0;
        for (; *s; s++) {
            int c = id(*s);
            if (!ch[u][c]) {
                memset(ch[sz], 0, sizeof(ch[sz]));
                val[sz] = 0;
                ch[u][c] = sz++;
            }
            u = ch[u][c];
        }
        val[u] = key;
    }
    void Construct () {
        int *s = Q, *e = Q;
        for (int i = 0; i < charset; i++) {
            if (ch[0][i]) {
                *e++ = ch[0][i];
                fail[ch[0][i]] = 0;
            }
        }
        while(s != e) {
            int u = *s++;
            if (val[fail[u]]) val[u] = 1;
            for (int i = 0; i < charset; i++) {
                int &v = ch[u][i];
                if (v) {
                    *e++ = v;
                    fail[v] = ch[fail[u]][i];
                } else {
                    v = ch[fail[u]][i];
                }
            }
        }
    }
    /*
        dp[i][j]表示长度为i。后缀为j的状态 最多就仅仅有10*10个后缀
        所以可以通过dp[n][j] = a0 * dp[n-1][0] + ... + ak * dp[n - 1][k]得到状态转移的矩阵
    */
    void work() {
        for (int i = 0; i < sz; i++) {
            for (int j = 0; j < charset; j++) {
                //对于i状态，通过加入ACTG可以得到新的状态（且之前已经构造过AC自己主动机，ch[i][j]便表示新状态）
                if (!val[i] && !val[ch[i][j]]) { //两个状态都必须是可行的，转化才有意义
                    A.a[i][ch[i][j]]++;
                }
            }
        }
    }
} AC;

Matrix Matrix_pow(Matrix A, ll k, ll mod) {
    res.clear();
    res.n = res.m = AC.sz;
    for (int i = 0; i < AC.sz; i++) res.a[i][i] = 1;
    while(k) {
        if (k & 1) res = res * A;
        A = A * A;
        k >>= 1;
    }
    return res;
}
int main () {
    int m, n;
    A.clear();
    AC.init();
    AC.reset();
    char str[15];
    scanf("%d%d", &m, &n);
    for (int i = 0; i < m; i++) {
        scanf("%s", str);
        AC.Insert(str, 1);
    }
    A.n = A.m = AC.sz;
    AC.Construct();
    //之前的都是AC自己主动机构造部分
    AC.work(); //得到状态转移的矩阵
    res = Matrix_pow(A, n, mod);
    int ans = 0;
    rep(i, AC.sz) ans += res.a[0][i];
    printf("%d
", ans % mod);
    return 0;
}