您的位置:首页 > 其它

Poj 2778 [AC自动机,矩阵乘法]

2017-04-03 23:47 246 查看
It’s well known that DNA Sequence is a sequence only contains A, C, T and G, and it’s very useful to analyze a segment of DNA Sequence,For example, if a animal’s DNA sequence contains segment ATC then it may mean that the animal may have a genetic disease. Until now scientists have found several those segments, the problem is how many kinds of DNA sequences of a species don’t contain those segments.

Suppose that DNA sequences of a species is a sequence that consist of A, C, T and G,and the length of sequences is a given integer n.

题意:有m种DNA序列是有疾病的,问有多少种长度为n的DNA序列不包含任何一种有疾病的DNA序列。(仅含A,T,C,G四个字符)

dp[i][j]=sum(dp[i-1][k])

枚举第i位填了’A’,’C’,’G’,’T’中的哪一个,自动机从k走到了哪里。

把i看作行,j看做列,那么dp[i][]和dp[i-1][]之间的关系显然和dp[i+1][]和dp[i-1][]的关系相同。

因此,可以用构造转移矩阵方法优化。。。矩阵就是dp【i][j]表示从I到j结点走一步(加一个字母)有几种方法,注意不能走到任何模式串的尾巴也不能它来转移,所以相当于尾结点的行列都不要,然后n次方,就是加了n个字母(题目要求了长度),初始化时候,只把0 0设置为1(没有字母),最后所有(0,i)的值加起来,就是0号结点转移到I的。

注意m=0特判。。。

#include <cstdio>
#include <iostream>
#include <cstring>
#include <queue>
using namespace std;
int n,m,node;
const int maxn = 15;
const long long MOD = 100000;
char ss[maxn][maxn];
int val[maxn*10],last[maxn*10],ch[maxn*10][maxn],fail[maxn*10],kmax;

int getpos(char a){
if(a == 'A')    return 0;
if(a == 'T')    return 1;
if(a == 'C')    return 2;
if(a == 'G')    return 3;
}

void Insert(int v){
int len = strlen(ss[v]);
int tmp = 0;
for(int i = 0; i < len ; i++){
int pos = getpos(ss[v][i]);
if(ch[tmp][pos] == 0){
memset(ch[node],0,sizeof(ch[node]));
val[node] = 0;
ch[tmp][pos] = node++;
}
tmp = ch[tmp][pos];
}
val[tmp] = v;
}

void getfail(){
queue<int> q;
int u = 0;
for(int i = 0; i < 4; i++){
if(ch[u][i]){
int tmp = ch[u][i];
fail[tmp] = 0;
q.push(tmp);
last[tmp] = 0;
}
}

while(!q.empty()){
int tmp = q.front();
q.pop();
for(int i = 0; i < 4; i++){
if(ch[tmp][i]){
int v = fail[tmp];
while(v && ch[v][i]==0) v = fail[v];
fail[ch[tmp][i]] = ch[v][i];
q.push(ch[tmp][i]);
// if(ch[tmp][i]==2)   printf("asd   ch[v][i] = %d \n",ch[v][i]);
last[ch[tmp][i]] = val[ch[v][i]]?ch[v][i] : last[ch[v][i]];
}
}
}
}

void init(){
node = 1;
kmax = 0;
memset(ch[0],0,sizeof(ch[0]));
for(int i = 1; i <= m ; i++){
scanf("%s",ss[i]);
Insert(i);
}

}
struct Matrix{
long long m[maxn*10][maxn*10];
}dp,ini;

void print(Matrix a){
for(int i = 0 ; i < 15 ; i++){
for(int j = 0 ; j < 15 ; j++)
printf("%lld ",a.m[i][j]);
printf("\n");
}
cout <<"\n";
}

void make_matrix(){
// cout <<"node = " << node <<endl;
memset(dp.m,0,sizeof(dp.m));
for(int i = 0 ; i < node; i++){
if(val[i]||val[last[i]])  continue;
for(int j = 0 ; j < 4; j++){
//printf("i = %d  j = %d \n",i,j);
//printf("last[%d]  = %d    val = %d\n",ch[i][j],last[ch[i][j]],val[ch[i][j]]);
if(last[ch[i][j]]||val[ch[i][j]])  continue;
int v = i;

if(ch[v][j]!=0){
dp.m[i][ch[v][j]] += 1;
kmax = max(kmax,max(i,ch[v][j]));
}
else{
while(v && ch[v][j] == 0 )    v = fail[v];
if(last[ch[v][j]]||val[ch[v][j]])  continue;
dp.m[i][ch[v][j]] += 1;
kmax = max(kmax,max(i,ch[v][j]));
}
}
}
// print(dp);
}

Matrix Mul(Matrix a,Matrix b){
Matrix p;
for(int i = 0 ;i <= kmax ; i++)
for(int j = 0 ; j <= kmax; j++)
p.m[i][j] = 0;
for(int i =  0 ;i <= kmax ;i++){
for(int j = 0 ; j <= kmax ; j++){
for(int k = 0 ; k <= kmax ; k++)
p.m[i][j] = (p.m[i][j]+a.m[i][k]*b.m[k][j])%MOD;
}
}
return p;
}

Matrix pow(Matrix a, int n){
Matrix p;
for(int i = 0 ;i <= kmax ; i++)
for(int j = 0 ; j <= kmax; j++)
p.m[i][j] = 0;
for(int i = 0 ; i <= kmax ; i++)  p.m[i][i] = 1;
while(n){
if(n & 1)   p = Mul(p,a);
n >>= 1;
a = Mul(a,a);
}
// print(p);
return p;
}

void sov(){
for(int i = 0 ;i <= kmax ; i++)
for(int j = 0 ; j <= kmax; j++)
ini.m[i][j] = 0;
ini.m[0][0] = 1;
//cout <<"n = "<<n<<endl;
ini =  Mul(ini,pow(dp,n));
int ans = 0;
for(int i = 0 ; i < node; i++)
ans = (ans + ini.m[0][i])%MOD;
printf("%d\n",ans);
}
/*
3 10
AT
GATC
TGAC
*/
int qpow(int n){
int a = 1,p = 4;
while(n){
if(n & 1)   a = (a*p)%MOD;
n >>= 1;
p = (p*p)%MOD;
}
return a;
}

int main(){
while(~scanf("%d%d",&m,&n)){
if(m == 0){
// cout <<"Asd"<<endl;
printf("%d\n",qpow(n));
continue;
}
init();
getfail();
make_matrix();
sov();
}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息