DNA序列 AC自动机+dp+矩阵快速幂优化
Source
2017 UESTC Training for Search Algorithm & String
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
#include
#include
#include
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 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("%lldn", ans);
#ifdef LOCAL
cout << endl;
}
#endif // LOCAL
return 0;
}
非特殊说明,本博所有文章均为博主原创,未经许可不得转载。
https://www.prolightsfxjh.com/
Thank you!
------from ProLightsfx
非特殊说明,本博所有文章均为博主原创,未经许可不得转载。
如经许可后转载,请注明出处:https://prolightsfxjh.com/article/uestc-1709/
共有 0 条评论