DNA Sequence
Time Limit: 1000MS | Memory Limit: 65536K | |
Total Submissions: 15538 | Accepted: 5984 |
Description
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.
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.
Input
First line contains two integer m (0 <= m <= 10), n (1 <= n <=2000000000). Here, m is the number of genetic disease segment, and n is the length of sequences.
Next m lines each line contain a DNA genetic disease segment, and length of these segments is not larger than 10.
Next m lines each line contain a DNA genetic disease segment, and length of these segments is not larger than 10.
Output
An integer, the number of DNA sequences, mod 100000.
Sample Input
4 3 AT AC AG AA
Sample Output
36
Source
For Myself
调了整整一天,还是说明AC自动机不熟,且码力不强,不过经过这道题,对AC自动机的理解又加深了很多。
一开始走了很多弯路,最后还是回到sth的正解上来。
一开始这么搞的,把所有的情况建成一颗AC自动机,有4^4个节点,超时。
Analysis
第一步先把所有的串读入并建成AC自动机,那么考虑一种爆搜:我每次枚举序列中的一个字母,就相当于在AC自动机上move了这个字母。
bad[i]是布尔数组,表示这个点是否是合法的,bad[]可以在求fail的同时求出。
令f[i][j]表示长度为i的串匹配到了j这个节点的方案数。
那么DP就出来了,对于每个节点node,枚举每一种字母,并且move一下,move到的点记为pos,则i的方案数就能够加给pos,即f[i][pos]+=f[i-1][node],即i能够更新pos。
初始状态f[1][1]=1,其余0。
这样我们就得到了一张更新表,map[i][j]表示i能够更新j几次,不难发现将f构造成一个1*节点数的矩阵,将它与map矩阵相乘,正好对应DP中的一次转移,于是,矩阵快速幂。
还有一种理解,map[i][j]表示从i到j走一步有几种走法,那么根据一个结论,map进行n次幂运算之后,的map[i][j]就表示从i到j有几种走法,那么显然∑map[1][i]就是答案,这和上面那个一样。
这里有一个常数优化,如果不加的话就会TLE。不合法的点我们应该去掉,也就是说所有和不合法的点相关的更新都是无效的,这在矩阵中是一些行和列,我们可以把这些行或列裁掉。
Code
//AC自动机 矩阵乘法加速DP
#include <cstdio>
#include <algorithm>
#include <queue>
#define ll long long
#define maxs 110
#define mod 100000
using namespace std;
ll n, m, trie[maxs][10], tot=1, bad[maxs], fail[maxs];
queue<ll> q;
struct matrix
{
ll m[70][70], r, c;
ll* operator[](ll x){return m[x];}
matrix operator*(matrix x)
{
ll i, j, k;
matrix t;
t.r=r,t.c=x.c;
for(i=1;i<=t.r;i++)
for(j=1;j<=t.c;j++)
{
t[i][j]=0;
for(k=1;k<=c;k++)t[i][j]=(t[i][j]+m[i][k]*x[k][j]%mod)%mod;
}
return t;
}
}one, map, f;
void insert(char *s)
{
ll pos;
for(pos=1;*s;s++)
pos=trie[pos][*s]?trie[pos][*s]:trie[pos][*s]=++tot;
bad[pos]=true;
}
void acabuild()
{
ll u, f, i;
q.push(1);
while(!q.empty())
{
u=q.front(),q.pop();
for(i=1;i<=4;i++)
if(trie[u][i])
{
for(f=fail[u];f and !trie[f][i];f=fail[f]);
if(f)fail[trie[u][i]]=trie[f][i];
else fail[trie[u][i]]=1;
bad[trie[u][i]]=bad[trie[u][i]] or bad[u] or bad[fail[trie[u][i]]];
q.push(trie[u][i]);
}
}
}
void init()
{
ll i, j, f, pos;
char s[100], *p;
scanf("%lld%lld",&m,&n);
for(i=1;i<=m;i++)
{
scanf("%s",s);
for(p=s;*p;p++)
{
if(*p=='A')*p=1;
if(*p=='T')*p=2;
if(*p=='G')*p=3;
if(*p=='C')*p=4;
}
insert(s);
}
acabuild();
for(i=1;i<=tot;i++)
for(j=1;j<=4;j++)
{
for(pos=i;pos and !trie[pos][j];pos=fail[pos]);
pos=pos?trie[pos][j]:1;
map[i][pos]++;
}
for(i=1;i<=tot;i++)
if(bad[i])for(j=1;j<=tot;j++)map[i][j]=map[j][i]=0;
}
void cut()
{
int i, j, k;
for(i=1;i<=tot;i++)map[i][tot+1]=map[tot+1][i]=1;
for(i=1;i<=tot;i++)
{
for(j=1;j<=tot;j++)if(map[i][j] or map[j][i])break;
if(j>tot)
{
tot--;
for(j=1;j<=tot;j++)
for(k=i;k<=tot;k++)map[j][k]=map[j][k+1];
for(j=1;j<=tot;j++)
for(k=i;k<=tot;k++)map[k][j]=map[k+1][j];
i--;
}
}
}
matrix pow(matrix a, ll b)
{
matrix tmp=a, ans=one;
for(;b;b>>=1,tmp=tmp*tmp)
if(b&1)ans=ans*tmp;
return ans;
}
int main()
{
ll i, ans=0, j;
init();
cut();
one.r=one.c=map.r=map.c=f.c=tot,f.r=1;
for(i=1;i<=tot;i++)one[i][i]=1;
f[1][1]=1;
map=pow(map,n);
f=f*map;
for(i=1;i<=tot;i++)ans=(ans+f[1][i])%mod;
printf("%lld\n",ans);
return 0;
}