题解：[BJWC2011] 禁忌
题意分析
记 \(N,\text{alphabet},T_i\) 分别为 \(n,v,t_i\)。
贪心选择
发现 \(s\) 分割成为若干个禁忌串实在太难了，考虑如何高效地对于确定的 \(s\) 求出能分割出的 \(t_i\) 的数量。
观察样例，aabb 的禁忌伤害为 \(1\)，划分为 aa bb 或 a abb。经过手玩，我们可以发现一种贪心策略是从左往右扫 \(s\)，如果匹配到一个禁忌串就立即匹配。
因为禁忌串对于禁忌伤害的贡献都是 \(1\)，我们要最大化划分出的禁忌串的数量。而把一个子串留到后面再划分，不能带来更多的禁忌串。
AC 自动机
所以建立 \(t_i\) 的 AC 自动机，从而快速匹配 \(s\) 中的禁忌串数量。
DP
暴力枚举 \(s\) 是无法接受的，考虑 DP 求解。
记 \(\textit{dp}_{i,j}\) 表示扫到 \(s_i\)，处于 AC 自动机节点 \(j\) 时即将到达禁忌串的终止节点，产生禁忌伤害的概率。
记 \(\operatorname{child}_k(j)\) 表示 AC 自动机节点 \(j\) 的代表字符 \(k\) 的子节点，\(0\leq k<v\)。设 \(\textit{id}_j\) 表示 AC 自动机节点 \(j\) 是否为一个禁忌串的终止节点。
-
若 \(\textit{id}_{\operatorname{child}_k(j)}=1\)，则有 \(\textit{dp}_{i,0}\leftarrow\textit{dp}_{i,0}+\dfrac{\textit{dp}_{i-1,j}}v\)。因为 AC 自动机匹配到禁忌串后会回到根节点 \(0\)，\(s_i\) 有 \(v\) 种选择。
此时还要记录期望伤害，\(\textit{ans}\leftarrow\textit{ans}+\dfrac{\textit{dp}_{i-1,j}}v\)。
-
否则有 \(\textit{dp}_{i,\operatorname{child}_k(j)}\leftarrow\textit{dp}_{i,\operatorname{child}_k(j)}+\dfrac{\textit{dp}_{i-1,j}}v\)。
记 \(L\) 为 AC 自动机节点数，得到了 \(\mathcal O(\textit{len}\cdot L)\) 的做法。因为 \(\textit{len}\) 可以达到 \(10^9\)，无法接受。
朴素 DP 参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
#define double long double
constexpr const int N=5,T=15,LEN=1e9,V=26,X=75;
int n,len,v,ans;
string s;
string t[N+1];
struct AC{
struct node{
int m[26];
int fail;
bool id;
}t[X+1];
int size;
void insert(string s,int id){
int p=0;
for(char i:s){
if(!t[p].m[i-'a']){
t[p].m[i-'a']=++size;
}
p=t[p].m[i-'a'];
}
t[p].id=true;
}
int q[N+1],front,rear;
void build(){
for(int i=0;i<V;i++){
if(t[0].m[i]){
q[rear++]=t[0].m[i];
}
}
while(front<rear){
int x=q[front++];
for(int i=0;i<v;i++){
if(t[x].m[i]){
t[t[x].m[i]].fail=t[t[x].fail].m[i];
t[t[x].m[i]].id|=t[t[t[x].m[i]].fail].id;
q[rear++]=t[x].m[i];
}else t[x].m[i]=t[t[x].fail].m[i];
}
}
}
double DP(){
static double dp[X+1][X+1];
double ans=0;
dp[0][0]=1;
for(int i=1;i<=len;i++){
for(int j=0;j<=size;j++){
for(int k=0;k<v;k++){
if(t[t[j].m[k]].id){
dp[i][0]+=dp[i-1][j]/v;
ans+=dp[i-1][j]/v;
}else{
dp[i][t[j].m[k]]+=1.0*dp[i-1][j]/v;
}
}
}
}
return ans;
}
}AC;
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int n;
cin>>n>>len>>v;
for(int i=1;i<=n;i++){
string t;
cin>>t;
AC.insert(t,i);
}
AC.build();
cout<<fixed<<setprecision(8)<<AC.DP()<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
矩阵加速转移
容易发现求解 \(\textit{dp}_{i,j}\) 的时候，与 \(i\) 无关，每次转移都是一样的线性转移。
因此考虑矩阵快速幂加速转移，因为 \(L\leq75\)，可以接受。
矩阵乘法加速 DP 其实就是枚举所有 DP 式中的元素，再把 DP 式里的系数填入矩阵，没有就是 \(0\)。不要忘了 \(\textit{ans}\) 对自己的贡献（右下角的 \(1\)）。
AC 代码
粘了一个矩阵板子。另外根据实现不同，可能会存在卡精度的问题，或许需要 long double。
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
constexpr const int N=5,T=15,LEN=1e9,V=26,X=75;
int n,len,v,ans;
string s;
string t[N+1];
struct Matrix{
int n,m;
double a[X+2+1][X+2+1];
Matrix(int nn=-1,int mm=-1){
if(mm==-1){
mm=nn;
}
if(nn!=-1){
n=nn,m=mm;
}else{
n=0,m=0;
}
for(int i=0;i<X+2+1;i++){
for(int j=0;j<X+2+1;j++){
a[i][j]=0;
}
}
}
Matrix(initializer_list<initializer_list<double>> x){
n=x.size();m=0;
for(auto &i:x){
m=max(m,(int)i.size());
}
int pi=0;
for(auto &i:x){
int pj=0;
for(auto &j:i){
a[pi][pj]=j;
pj++;
}
pi++;
}
}
double* operator [](const int &x){
return a[x];
}
void unit(){
for(int i=0;i<n;i++){
for(int j=0;j<m;j++){
a[i][j]=0;
}
}
for(int i=0;i<n;i++){
a[i][i]=1;
}
}
void print(){
for(int i=0;i<n;i++){
for(int j=0;j<m;j++){
cout<<a[i][j]<<' ';
}
cout<<endl;
}
}
};
Matrix operator *(Matrix A,Matrix B){
Matrix C(A.n,B.m);
for(int i=0;i<A.n;i++){
for(int j=0;j<B.m;j++){
for(int k=0;k<A.m;k++){
C[i][j]+=A[i][k]*B[k][j];
}
}
}
return C;
}
Matrix& operator *=(Matrix &A,Matrix B){
return A=A*B;
}
Matrix operator +(Matrix A,Matrix B){
Matrix C(A.n,A.m);
for(int i=0;i<A.n;i++){
for(int j=0;j<A.m;j++){
C[i][j]=A[i][j]+B[i][j];
}
}
return C;
}
Matrix& operator +=(Matrix &A,Matrix B){
for(int i=0;i<A.n;i++){
for(int j=0;j<A.m;j++){
A[i][j]+=B[i][j];
}
}
return A;
}
Matrix qpow(Matrix base,int n){
Matrix ans(base.n);
ans.unit();
while(n){
if(n&1){
ans=base*ans;
}
base=base*base;
n>>=1;
}
return ans;
}
struct AC{
struct node{
int m[26];
int fail;
bool id;
}t[X+1];
int size;
void insert(string s,int id){
int p=0;
for(char i:s){
if(!t[p].m[i-'a']){
t[p].m[i-'a']=++size;
}
p=t[p].m[i-'a'];
}
t[p].id=true;
}
int q[X+1],front,rear;
void build(){
for(int i=0;i<v;i++){
if(t[0].m[i]){
q[rear++]=t[0].m[i];
}
}
while(front<rear){
int x=q[front++];
for(int i=0;i<v;i++){
if(t[x].m[i]){
t[t[x].m[i]].fail=t[t[x].fail].m[i];
t[t[x].m[i]].id|=t[t[t[x].m[i]].fail].id;
q[rear++]=t[x].m[i];
}else t[x].m[i]=t[t[x].fail].m[i];
}
}
}
double DP(){
Matrix dp(size+2,1),f(size+2,size+2);
dp[0][0]=1;
for(int j=0;j<=size;j++){
for(int k=0;k<v;k++){
if(t[t[j].m[k]].id){
// dp[i][0]+=dp[i-1][j]/v;
f[0][j]+=1.0/v;
// ans+=dp[i-1][j]/v;
f[size+1][j]+=1.0/v;
}else{
// dp[i][t[j].m[k]]+=1.0*dp[i-1][j]/v;
f[t[j].m[k]][j]+=1.0/v;
}
}
}
f[size+1][size+1]=1;
dp=qpow(f,len)*dp;
return dp[size+1][0];
}
}AC;
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int n;
cin>>n>>len>>v;
for(int i=1;i<=n;i++){
string t;
cin>>t;
AC.insert(t,i);
}
AC.build();
cout<<fixed<<setprecision(8)<<AC.DP()<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
浙公网安备 33010602011771号