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
[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\n", ans % mod);
return 0;
}