首页 > 代码库 > 关于PE458(project euler 458 Permutations of Project)的思考

关于PE458(project euler 458 Permutations of Project)的思考

本文回顾了PE458的解题过程中遇到的问题,介绍了trie,AC自动机,自动机化简算法.

题意:给定project 这7个字母(等价于互不相同的7个字母),问长度为10^12的由这7个字母组成的串且串中任意连续7个字母中不包含所有的7个字母(即:任意连续7个字母不是project的排列)有多少个,答案对1e9取模.

尝试用AC自动机和矩阵二分解决问题:

首先想到的构造一个自动机,能够判断一个串是否满足要求.假定构造出来,我们可以构造转移矩阵tran[i][j],当在状态机j状态下输入某个字母并转移到i状态,tran[i][j]的值增加1.给定列向量x[0] = (1, 0, 0, ...)‘, x[i] = tran * x[i-1], 第一个分量的状态表示空串的状态.于是x[i][j]就表示长度为i,且状态是j的串的个数.考虑其中需要的状态,相加即可.

怎么构造这个自动机呢,当然用AC自动机了.也就是说,先构造一个trie,然后把7!个不合法的串插入到这个trie中.在构造的trie中,还有一些状态在接受输入时没有下一个状态,所以下一步是把这些状态填满.

考虑一个非法状态,无论输入是什么,下一个状态还是本身,因为只要某处出现project的排列了,无论后面的输入何如,这个串都是非法的.

考虑其它中间状态,这个状态对应一个字符串s,|s|>0 .我们只需要找出一个串t,|t| < |s|, t是s的一个前缀,同时t也是s的后缀,并且|t|最大化. 显然t是存在的,因为空串是满足要求的.t同时也是唯一的. 当s在输入字母l时,如果没有下一个状态,那么等价于t在输入字母l时的状态转移. 而关于t可以递归地构造.

int p = 0, q = 0;
sfx[1] = 1;
for (int i = 0; i < 7; ++i) if (trie[1][i])
{
	que[q++] = trie[1][i];
	sfx[trie[1][i]] = 1;
}
else
{
	trie[1][i] = 1;
}
while (p < q)
{
	int curr = que[p++];
	assert(danger[sfx[curr]] == 0);
	danger[curr] = danger[curr] || danger[sfx[curr]];
	if (danger[curr])
	{
		for (int i = 0; i < 7; ++i)
		trie[curr][i] = curr;
	}
	else
	{
		for(int i = 0; i < 7; ++i) if (trie[curr][i])
		{
			que[q++] = trie[curr][i];
			sfx[trie[curr][i]] = trie[sfx[curr]][i];
		}
		else
		{
			trie[curr][i] = trie[sfx[curr]][i];
		}
	}
}

其中状态1表示空串. 如果状态i对应的串是s那么状态sfx[i]的状态就是t.danger[i]表示该状态是非法的. trie[i][j] == 0的时候表示在输入j的时候trie没有边.可以看到这就是一个简单的bfs,再深入看看其实ac自动机就是kmp算法在多串上的扩展.

然后构造tran矩阵,矩阵二分可以计算出来x[1e12],把合法分量相加即可.

但是杯具的是,状态有13700个,n^3的矩阵乘法在一般机器上不可行.(在我机器上计算量为6.5204 * 10^12 的简单计算需要5个小时. 姑且认为矩阵乘法的每一个操作和我这个参考操作的计量时间一样,单个矩阵乘法的计算量是我的参数计算的40%不妨设为2小时.再考虑1e12对应的二进制大概是40位,总共需要80个小时的计算.如果再考虑两种基本操作的差异,缓存神马的,时间更不好估计了)

直接构造的自动机由于状态数太多,不可行,所以得想办法改.


尝试减少自动机的状态:

我总共采用了三种减少状态的方法.

首先是经典的自动机化简的算法.
用belong[i]表示i状态被化简后的状态.初始情况下,非法状态的值是2,合法状态状态的值是1.
给定一个无限循环:
用map<vector<int>, vector<int> > mem;来记录当前的转移情况.
对于状态i,key的第一个分量是belong[i],其它分量是belong[trie[i][0], belong[trie[i][1]], 
..., belong[i][6]. value对应的vector包含i. 这样就通过key把所有的状态区分出来了.
然后再遍历一次这个map,对于一个key,为value对应的所有状态分配相同的belong值.
如果在分配前和分配后,状态数不发生变化,则算法终止.

for (int i = 1; i < top; ++i) if (danger[i] == 0)
{
	belong[i] = 1;
}
else
{
	belong[i] = 2;
}
int size = 2;
for (;;)
{
	map<vector<int>, vector<int> > mem;
	for (int i = 1; i < top; ++i)
	{
		vi key;
		key.pb(belong[i]);
		for (int j = 0; j < 7; ++j)
		key.pb(belong[trie[i][j]]);
		mem[key].pb(i);
	}
	
	int id = 1;
	for (auto& it: mem)
	{
		for (auto& s: it.second)
		belong[s] = id;
		++id;
	}
	if (mem.size() == size)
		break;
	size = mem.size();
}

算法其实很快.在coursera上的dfa课程中也给了个算法.算法的复杂度是输入状态数的平方,好处是给了一个有确定复杂度的算法,坏处是平方级的太慢,不如这样迭代.

在化简后,状态数为8661,仍然不可行.

其次是在构造自动机的时候可以试图减少一些状态.
考虑到所有的trie的叶子结点都是不合法的,可以在构造trie的时候把这些点视为一个,然后再使用自动机化简的算法.化简后状态数还是8661,反而验证了自动机存在一个最小状态数,也一定程度上说明我写的化简算法大概是正确的.

最后,不考虑trie,试图直接构造自动机,这样可以充分利用题目条件,有可能达到减少状态的目的.

考虑某个状态,对应于长度为l的串s,在面对输入i的时候,如果i不在s中出现,则将i追加到l,然后转移到一个新的状态.如果i在s中出现,设出现位置是pos,那么把s中从pos+1的串取出来,将i追加到这个取出来的串后面,并转移到这个状态.

struct Pt
{
	int s;	// 表示自动机中的状态
	vector<int> v;	// 表示这个状态对应的串
};
queue<Pt> last;
map<vector<int>, int> mem; // key为串,value为对应的状态
// last中放入一个初始状态,Pt.s = 1, Pt.v 为空
// 标记状态1为访问过
// 将空串加入到mem中
for (int len = 0; len <= 5; ++len)
{
	queue<Pt> next;
	对于last中每一个状态now:
		枚举输入字母i
			如果i不在now.v中出现,则转移到的串B, 如果B出现过则找到对应的状态,否则就分配新的状态
				并将状态加入到next队列中. 标志状态转移.
			如果i在now.v中出现,也能构造一个串B, B一定出现过, 标记状态转移.
	last.swap(next);
}
// last中的长度为6的状态,特殊处理长度为6的状态到长度为7的状态的转移

现在看一下状态数:
长度为0的:1
长度为1的:1!*C(7,1)=P(7,1)=7
长度为2的:2!*C(7,2)=P(7,2)=42
长度为3的:3!*C(7,3)=P(7,3)=210
长度为4的:4!*C(7,4)=P(7,4)=840
长度为5的:5!*C(7,5)=P(7,5)=2520
长度为6的:6!*C(7,6)=P(7,6)=5040
长度为7的:1
1+7+42+210+840+2520+5040+1=8661
如果不特殊处理长度为7的串则:
1+7+42+210+840+2520+5040+5040=13700

在化简的三个思路中得到了相同的结果,也再次验证了这个自动机的状态数的下界是8661.说明了对于这个自动机,状态数的化简只对非法结点产生了效果.同是第三次直接给出自动机的算法也更进一步说明了trie对应的AC自动机的内涵.

[其实在实现这个算法前我错误地估计了状态数,把后面的组合数弄掉了,长度为6的才720,再加上其它的也不超过1000,1000^3一次的矩阵乘法还是可行的]


正解是什么:
在上述碰壁后,我恍然大悟,茅塞顿开,删除了很多代码,保留了一些,补了不到十行代码,解决了这个问题.由于PE的游戏规则这里就不细说了.

关于PE458(project euler 458 Permutations of Project)的思考