题目大意
有一个长度为无限的字符集只有T,H的随机字符串S,现在给出n个长度为m的匹配串T1..nT1..n,问每个匹配串作为S中最先成功匹配的串的概率。
解题思路
暴力怎么做?考虑使用AC自动机,设f[i][j]表示从节点i出发,最先走到匹配串j的终止节点,的概率。
那么每个点的转移是一个方程,把他们联立起来高斯消元即可。
这样的复杂度最坏是O(n33)O(n33)的,不过AC自动机节点一般难以达到n2n2。
看起来是要优化了。在暴力中,我们最后其实只需要f[0][1..n]这n个变量的值。考虑能不能只搞这些变量,即P(A)表示串A最先匹配到的概率。
我们需要寻找变量之间的关系,似乎找不到。
引入一个新变量P(N),N是一个随机的字符串,而且T1..nT1..n不能匹配它,即不在终止点的状态。P(N)则表示到达这种状态的概率。
这有什么用呢?考虑两个T中字符串A=TTH,B=HTTA=TTH,B=HTT。如果在N后面加上TTH,无论N是什么,游戏一定结束。考虑N的每一种情况,我们可以得到
P(N+TTH)=P(N′+HT+TTH)+P(N′+HH+TTH)+P(N′+TT+TTH)+P(N′+TH+TTH)=P(A)+P(B+H)+P(B+TH)P(N+TTH)=P(N′+HT+TTH)+P(N′+HH+TTH)+P(N′+TT+TTH)+P(N′+TH+TTH)=P(A)+P(B+H)+P(B+TH)。
那么P(N)∗2−3=P(A)+P(B)∗2−1+P(B)∗2−2P(N)∗2−3=P(A)+P(B)∗2−1+P(B)∗2−2,也就是说,考虑在N后面加上任意一个TiTi,如果某个TjTj的某个后缀为TiTi前缀,那么P(Tj)P(Tj)在等式里会出现。系数用KMP算出来。
这样能得到N个方程,再加上∑iP(Ti)=1∑iP(Ti)=1,就有N+1个方程,又有N+1个系数,可以高斯消元。
注意填系数的时候都乘上2m2m,不然m大的时候double会认为系数为0。
代码
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<map>
#include<set>
using namespace std;
typedef long long ll;
typedef double db;
#define fo(i,j,k) for(i=j;i<=k;i++)
#define fd(i,j,k) for(i=j;i>=k;i--)
#define cmax(a,b) (a=(a>b)?a:b)
#define cmin(a,b) (a=(a<b)?a:b)
const int N=300+5,M=5e6+5,mo=1e9+7;
const db eps=1e-7;
int nxt[N],len,n,m,i,j,K;
char s[N][N];
db a[N][N],xs;
int kmp(char *s1,char *s2)
{
int i,j;
j=0;
fo(i,2,m)
{
while (j&&s1[i]!=s1[j+1]) j=nxt[j];
if (s1[i]==s1[j+1]) j++;
nxt[i]=j;
}
j=0;
fo(i,1,m)
{
while (j&&s1[j+1]!=s2[i]) j=nxt[j];
if (s1[j+1]==s2[i]) j++;
}
return j;
}
void gauss(int n)
{
int i,j,k;
fo(i,0,n-1)
{
fo(j,i,n)
if (fabs(a[j][i])>eps) break;
fo(k,0,n+1) swap(a[j][k],a[i][k]);
fo(j,i+1,n)
{
xs=a[j][i]/a[i][i];
fo(k,i,n+1) a[j][k]-=a[i][k]*xs;
}
}
fd(i,n,0)
{
fd(j,i-1,0)
{
xs=a[j][i]/a[i][i];
a[j][i]-=a[i][i]*xs;
a[j][n+1]-=a[i][n+1]*xs;
}
}
fo(i,0,n)
a[i][n+1]/=a[i][i],a[i][i]/=a[i][i];
}
int main()
{
freopen("t4.in","r",stdin);
//freopen("t4.out","w",stdout);
scanf("%d %d",&n,&m);
fo(i,1,n)
scanf("%s",s[i]+1);
K=m;
fo(i,1,n)
{
fo(j,1,n)
{
len=kmp(s[i],s[j]);// both s[i]'s prefix and s[j]'s suffix
while (len)
{
a[i][j]-=pow(0.5,m-len-K);
len=nxt[len];
}
}
a[i][0]=pow(0.5,m-K);
}
fo(i,1,n) a[0][i]=1,a[0][n+1]=1;
gauss(n);
fo(i,1,n)
printf("%.6lf\n",a[i][n+1]);
}