POJ2778 DNA Sequence AC自动机+DP+矩阵快速幂

Problem Address：http://poj.org/problem?id=2778

【前言】

这道题算是出自Trie主题的。

算是一道比较综合的题目。

【思路】

（1）首先用所给的字符串构造出一个AC自动机。

          这一步很简单的，也不多说了。

          注意自动机的fail以及标记，如果子串为非法，则该串需标记为非法。

（2）用AC自动机构造一个DP矩阵。

          该矩阵表示到达某个状态的合法转移次数。

（3）求出矩阵的n次幂，第一行的合法结点之和即为答案。

          简单地矩阵快速幂。相当于普通的二分快速幂。

          

【代码】

#include <iostream>
#include <cstring>
using namespace std;

char map[1005][1005];
char word[1005];

const int kind = 4;
const int mod = 100000;
const int size = 110;

__int64 g[size][size];

struct node
{
	node *fail;
	node *next[kind];
	int end;
	int index;
}trie[size];

int total;
node *root;

node *q[size];
int head, tail;

inline node* new_node()
{
	node *p = &trie[total];
	p->index = total;
	total++;
	p->fail = NULL;
	p->end = 0;
	memset(p->next, 0, sizeof(p->next));
	return p;
}

inline int code(char c)
{
	switch (c)
	{
	case 'A': return 0;
	case 'T': return 1;
	case 'C': return 2;
	case 'G': return 3;
	}
	return -1;
}

void insert(node *root, char *s)
{
	node *p = root;
	int i=0, t;
	while(s[i]!='\0')
	{
		t = code(s[i]);
		if (p->next[t]==NULL)
			p->next[t] = new_node();
		p = p->next[t];
		i++;
	}
	p->end++;
}

void build_ac_automation(node *root)
{
	int i;
	node *temp;
	root->fail = NULL;
	head = 1;
	tail = 0;
	q[0] = root;
	while(head!=tail)
	{
		temp = q[tail];
		tail++;
		for (i=0; i<kind; i++)
		{
			if (temp->next[i]!=NULL)
			{
				if (temp==root)
				{
					temp->next[i]->fail = root;
				}
				else
				{
					temp->next[i]->fail = temp->fail->next[i];
					if (temp->next[i]->fail->end!=0)
						temp->next[i]->end++;
				}
				q[head] = temp->next[i];
				head++;
			}
			else
			{
				if (temp==root)
					temp->next[i] = root;
				else
					temp->next[i] = temp->fail->next[i];
			}
		}
	}
}

void MatrixMultiply(__int64 b[][size], __int64 c[][size], int sz)
{
	int i,j,k;
	__int64 temp[size][size] = {0};
	for (i=0; i<sz; i++)
	{
		for (j=0; j<sz; j++)
		{
			for (k=0; k<sz; k++)
			{
				temp[i][j] += b[i][k]*c[k][j];
				if (temp[i][j]>=mod)
					temp[i][j] %= mod;
			}
		}
	}
	for (i=0; i<sz; i++)
	{
		for (j=0; j<sz; j++)
		{
			b[i][j] = temp[i][j];
		}
	}
}

void MatrixPow(__int64 t[][size], __int64 a[][size], int sz, int n)
{
	while(n>0)
	{
		if (n&1) MatrixMultiply(t, a, sz);
		MatrixMultiply(a, a, sz);
		n >>= 1;
	}
}

int main()
{
	int m,n;
	char word[15];
	int i,j;
	__int64 tot[size][size] = {0};
	__int64 ans;
	scanf("%d %d", &m, &n);
	total = 0;
	root = new_node();
	for (i=0; i<m; i++)
	{
		scanf("%s", word);
		insert(root, word);
	}
	build_ac_automation(root);//构造AC自动机
	memset(g, 0, sizeof(g));
	for (i=0; i<total; i++)//构造DP转移矩阵
	{
		if (trie[i].end==0)
		{
			for (j=0; j<kind; j++)
			{
				if (trie[i].next[j]->end==0)
				{
					g[i][trie[i].next[j]->index]++;
				}
			}
		}
	}
	memset(tot, 0, sizeof(tot));
	for (i=0; i<total; i++)//初始化为单位矩阵
		tot[i][i] = 1;
	MatrixPow(tot, g, total, n);//快速幂
	ans = 0;
	for (i=0; i<total; i++)
	{
		if (trie[i].end==0)
		{
			ans += tot[0][i];
			if (ans>=mod)
				ans %= mod;
		}
	}
	printf("%I64d\n", ans);
	return 0;
}

【P.S】

最近在做POJ的contest。

只可惜功力不足，结局惨不忍睹。

然后跟别人聊了一些话，又想了很多。

于是开始糊涂起来。

不过这个暑假还是要好好努力下去，尽管和别人的差距是如此之大。

总之，不顾一切。