[BZOJ4820]硬币游戏 - Blog

有一种暴力方式：对所有猜好的的序列建一个AC自动机，然后对这个有向图解方程求期望经过次数。
（然而 $O((nm)^3)$ 好像没什么用）

我们常试着少用一些中间状态去转移 $A$ 的期望。
设 $P(A)$ 为 $A$ 胜利的期望。
定义 $F$ 为不满足任何一种序列的情况，现在有两个同学进行了猜测，设序列分别为 $A,B$
为方便解释，这里给出一个例子。
设 $A=TTH,B=HTT$
把所有 $F$ 往后延伸 $TTH$ 这个串一定能结束，所有 $A$ 胜利的情况也可以由某个 $F$ 加上 $TTH$ 组成。
现在有一个问题：并不是所有 $F$ 加上 $TTH$ 都是 $A$ 胜利的情况。
$THTTH=T+HTT+H ->F+B+H$
$THTTH=TH+TTH->F+A$
这里 $B$ 先碰到底。
（ $F$ 的后缀加上 $A$ 的前缀构成了 $B$ ）
这里会发现，任何串的某一后缀如果与 $A$ 的某一前缀相同，都会有可能把 $A$ 的胜利从中途劫走。
那些后缀与 $A$ 的前缀相同的 $B$ ，他们胜利的情况可以用 $F+*$ 表示出来，但这里我们把 $*$ 定义成与 $A$ 相同的前缀部分。
这里的 $A$ 对所有 $F$ 都做了一次增加，所以对于 $B$ ，它所有可能的情况都能卡掉 $A$ 。
所以在计算 $A$ 的时候，需要减去那些出现了后缀的那些字符串的期望。
设 $B$ 的前缀与 $A$ 的后缀的公共部分长度为 $k$ 。
这里如果统计时出现 $B$ 的某种前缀卡掉了 $A$ ，那卡掉 $A$ 之后必须沿着 $A$ 添加的字符继续添加下去，否则不会被计入 $A$ 的统计。
所以实际上减去的是 $B$ 加上 $A$ 与 $B$ 不同的那一段，而每走一步期望都要减少到原来的一半。
这样要减去的期望就是 $B$ 的期望的 $\frac{1}{2^{m-k}}$
一个很有趣的事情是：要在统计里排除的情况只需要考虑某个字符串的后缀，前缀会由 $F$ 来解决。
这样的话一个 $A$ 可能减去很多次 $B$ 。
现在把AC自动机中非结束节点的期望到达次数之和定为 $Q$ (其实就是 $F$ 的期望啦)
这样对于每个 $A$ 都有方程:

P(A)=\frac{1}{2^m}Q- \sum_ B \frac{1}{2^ {m-k}}P(B)

（这里 $Q$ 要乘 $\frac{1}{2^m}$ 是因为我们钦定所有 $F$ 按照 $A$ 序列走 $m$ 步）

最后加上 $\sum_ i ^ n P(i)=1$ 就可以凑齐 $n+1$ 个式子高斯消元了。
这里实现的时候有一些问题要注意一下：

$\frac{1}{2^m}Q$ 最好把它当作一个变量看，要不然会出现爆精度导致系数全0得到 $n$ 未知数 $n+1$ 方程
GCC在处理超过字长的位移时直接把位移膜了字长，而不是变成0，所以对于超过字长的位移需要特判。

code:

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
const int N=302;
const long double eps=1e-12;
int n,m;
char ch[N][N];
long double mp[N][N];

inline int nxi(){
	int x=0;
	char c;
	while((c=getchar())>'9'||c<'0');
	while(x=x*10-48+c,(c=getchar())>='0'&&c<='9');
	return x;
}

inline void get_mp(int t){
	static int nx[N];
	for(int j=0,i=2;i<=m;++i){
		while(j&&ch[t][i]!=ch[t][j+1]) j=nx[j];
		if(ch[t][i]==ch[t][j+1]) ++j;
		nx[i]=j;
	}
	mp[t][t]=-1;
	for(int i=1;i<=n;++i){
		int j=1,k=0;
		for(;j<=m;++j){
			while(k&&ch[i][j]!=ch[t][k+1]) k=nx[k];
			if(ch[i][j]==ch[t][k+1]) ++k;
		}
		if(k==m) k=nx[k];
		for(;k;k=nx[k]){
			if(m-k>63) break;
			mp[t][i]-=1.0/(1ull<<(m-k));
		}
	}
}

inline void gauss(){
	for(int i=1;i<=n+1;++i){
		int j=i;
		while(j<=n+1&&fabs(mp[j][i])<eps) ++j;
		if(j>n+1) return;
		std::swap(mp[i],mp[j]);
		for(int k=1;k<=n+1;++k){
			if(k==i) continue;
			const long double p=mp[k][i]/mp[i][i];
			for(int l=0;l<=n+1;++l){
				mp[k][l]-=mp[i][l]*p;
			}
		}
	}
	for(int i=1;i<=n+1;++i){
		if(fabs(mp[i][0])>eps) mp[i][i]=mp[i][0]/mp[i][i];
		else mp[i][i]=0;
	}
}

inline void gauss_init(){
	for(int i=1;i<=n;++i){
		mp[n+1][i]=1;
		get_mp(i);
		mp[i][n+1]=1;
		//we don't need to know men.
		//1/(2^m)*men instead of men is also ok.
	}
	mp[n+1][0]=1;
}

int main(){
#ifndef ONLINE_JUDGE
//	freopen("a.in","r",stdin);
#endif
	n=nxi(),m=nxi();
	for(int i=1;i<=n;++i){
		scanf("%s",ch[i]+1);
	}
	gauss_init();
	gauss();
	for(int i=1;i<=n;++i){
		printf("%.10lf\n",(double)mp[i][i]);
	}
	return 0;
}