poj2778 DNA Sequence

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. 

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. 

Output

An integer, the number of DNA sequences, mod 100000.

Sample Input

4 3
AT
AC
AG
AA

Sample Output

36

Source

POJ Monthly--2006.03.26,dodo

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;
}