BZOJ2553: [BeiJing2011]禁忌 AC自动机 期望DP 矩阵

题目大意：给出n个禁忌串，定义任意一个字符串的禁忌伤害是不相交的禁忌子串的最大个数，求长度为len的随机串禁忌伤害期望。

所有的串只包含前alphabet个拉丁字母。

N ≤ 5，len ≤109，1 ≤ alphabet ≤ 26。

对所有禁忌串建AC自动机，问题转化成：每次在AC自动机上走一步，走到危险节点（即有危险标记的节点或者能沿fail走到危险标记的节点）就自动回到根，求危险节点的期望经过次数。

考虑期望的定义：sigma（权值 * 概率）。权值就是1，因此将概率累加即可。可以先不考虑危险节点，则其为一个普通的图，用矩阵乘法即可求出走n步后从i点走到j点的概率。考虑危险节点，每一次转移，都会有一些节点以一定概率转移到危险节点，走到这些点的概率*转移概率就是这次转移对答案的贡献，可以随着矩阵的递推而计算每一次的贡献，因此新建一个代表所有危险点的节点，能走到危险点的点同时向根和这个点连边，同时这个点自身连一条1边，代表着累加上一次的结果。自乘len次后f[0][danger]就表示了从0号节点走len步，danger节点的期望经过次数。

据说有精度问题，要long double才能过。

#include<cstdio>
#include<queue>
using std::queue;
typedef long double cnt_t;
int n,len,alpha;
struct node
{
node *s[26],*fail;
bool da;
node():s(),fail(),da(){}
inline void* operator new(size_t);
}a[100];
inline void* node::operator new(size_t)
{
static size_t t=0xffffffffu;
return a+ ++t;
}
struct trie
{
node *rt;
trie():rt(new node){}
inline void insert(const char *s)
{
node *now=rt;
while(*s)
{
node*& x=now->s[*s-'a'];
if(!x) x=new node;
now=x;
++s;
}
now->da=1;
}
};
size_t maxn;
struct martix
{
cnt_t p[80][80];
cnt_t* operator [](size_t x){return p[x];}
martix():p(){}
martix(size_t):p()
{
for(size_t i=0;i<=maxn;++i)
p[i][i]=1.0;
}
martix operator *(martix &b) const
{
martix res;
for(size_t i=0;i<=maxn;++i)
for(size_t j=0;j<=maxn;++j)
for(size_t k=0;k<=maxn;++k)
res[i][j]+=p[i][k]*b[k][j];
return res;
}
}ans;
struct ac_auto:public trie
{
ac_auto():trie(){}
inline size_t operator() ()
{
static queue<node*> q;
size_t res=1;
rt->fail=rt;
for(int i=0;i<alpha;++i)
{
node*& x=rt->s[i];
if(!x) x=rt;
else x->fail=rt,q.push(x);
}
while(!q.empty())
{
node *now=q.front();q.pop();
now->da|=now->fail->da;
for(int i=0;i<alpha;++i)
{
node*& x=now->s[i];
if(!x) {x=now->fail->s[i];continue;}
x->fail=now->fail->s[i];
q.push(x);
}
++res;
}
return res;
}
}ac;
template<typename T,typename Int>
T power(T x,Int a)
{
T ans(1);
while(a)
{
if(a&1) ans=ans*x;
x=x*x;
a>>=1;
}
return ans;
}
char s[20];
int main()
{
scanf("%d%d%d",&n,&len,&alpha);
while(n--) scanf("%s",s),ac.insert(s);
maxn=ac();
ans[maxn][maxn]=1.0;
cnt_t base=1.0/alpha;
for(size_t i=0;i<maxn;++i)
for(int k=0;k<alpha;++k)
{
size_t j=a[i].s[k]-a;
if(a[j].da) ans[i][0]+=base,ans[i][maxn]+=base;
else ans[i][j]+=base;
}
ans=power(ans,len);
double res=ans[0][maxn];
printf("%.10lf\n",res);
return 0;
}