FFT处理字符串匹配类问题
起因
hdu多校3遇到了一题带通配符*的字符串匹配题，赛中想用KMP写，但是发现是假算法，赛后用Bitset补了一发，优化常数后也T，题解是用FFT写的字符串匹配问题。
问题1:普通字符串匹配,可用KMP做的那种
参考：
(27条消息) FFT处理字符串匹配问题_CCCCTong的博客-CSDN博客
题解 P4173 【残缺的字符串】 - Ebola 的博客 - 洛谷博客 (luogu.com.cn)
数学原理
定义匹配函数$C(x,y)=A_x-B_y$，如果两个字符相同，则满足$C(x,y)=0$，要推广到字符串上：
再定义完全匹配函数$P = \sum^n_{i=1}C(i,i)=\sum^n_{i=1}(A_i-B_i)$，即每一位上的匹配函数求和。
故设A为文本串，B为模式串，取A串的子串A[k,k+m-1]（与B等长），与模式串B进行匹配，定义得到的$P(x) = \sum^{m-1}_{i=0}C(i,k+i)$=$\sum^{m-1}_{i=0} (A\_i-B\_{k+i})$
- 如果B的子串B[k~k+m-1]和A完全匹配，则他们的全匹配函数为0
但其实离上面的还差一步，考虑串A=ab和B=ba，在第一个位置有$C(1,1)=1-2=-1$;在第二个位置上$C(2,2)=2-1=1$;这样以来，虽然每个位上的匹配函数C不为0，但全匹配函数P为0，所以不行。
- 这里很容易就能想到使用绝对值来定义匹配函数，但是绝对值需要分类讨论，所以我们使用平方，即定义$C(i,i)=(A_i-B_i)^2$。
这样一来，$P(k) = \sum^{m-1}_{i=0} C(i,k+i)=\sum^{m-1}_{i=0} (A_i-B_{k+i})^2$
再拆括号：$\sum^{m-1}_{i=0} (A_i-B_{k+i})^2=\sum_{i=0}^{m-1}(A_i^2+B_{k+i}^2-2A_iB_{k+i})$
分开来求和就有$P(x)=\sum_{i=0}^{m-1}A_i^2+\sum_{i=0}^{m-1}B_{k+i}^2-2\sum_{i=0}^{m-1}A_iB_{k+i}$
因为是匹配问题，要求出所有k的P(k),第一项与k无关，即对$P(k)$为常数
第二项通过前缀和可以$O(n)$得到。
第三项则因为要枚举k和i，需要$O(n^2)$,可以怎么样优化？
我们可以观察$A_iB_{k+i}$，可以发现$(k+i)-i=k$。但是FFT在x的次数上是+，这里要-才能满足规律，咋办？
- reverse(A,A+m)反转B串，记为S串，则有$A_i=S_{m-1-i}$，则第三项$\sum_{i=0}^{m-1}A_iB_{k+i}$变为了$\sum_{i=0}^{m-1}S_{m-1-i}B_{k+i}$。两个下标之和刚好为$m-1+k$
这里k是左端点，为了符合字符串匹配类问题记右端点的习惯，我们用$x=k+m-1$即k开头子串的右端点来代换k：$k=x-m+1$
我们发现$x=k+m-1$正好是卷积结果中以k为左端点的结果的位置。
实现细节
注意右端点x代换左端点k有$k=x-m+1$
设第一项为$T=\sum_{i=0}^{m-1}A_i^2$
第二项为$\sum_{i=0}^{m-1}B_{(x-m+1)+i}^2$，我们预处理f[x]数组$f(x)=\sum_{j=0}^{j=x}B_j^2$（f为前缀和），即第二项为f[x]-f[x-m]
对于第三项，A和B在进行FFT之后得到G数组，有第三项$\sum_{i=0}^{m-1}A_iB_{k+i}=G_{m-1+k}=G_{x}$刚刚好就是x!
代码
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
int cansel_sync=(ios::sync_with_stdio(0),cin.tie(0),0);
#define ll long long
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
#define INF 0x3f3f3f3f
//FFT板子
const int MAXN = 214514;
const double Pi=acos(-1.0);
struct complex{
double x,y;
complex (double xx=0,double yy=0){x=xx,y=yy;}
}a[MAXN],b[MAXN];
complex operator + (complex a,complex b){ return complex(a.x+b.x , a.y+b.y);}
complex operator - (complex a,complex b){ return complex(a.x-b.x , a.y-b.y);}
complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);}//不懂的看复数的运算那部分
int l,r[MAXN];
int limit=1;
void FFT(complex *A,int type){
for(int i=0;i<limit;i++)
if(i<r[i]) swap(A[i],A[r[i]]);//求出要迭代的序列
for(int mid=1;mid<limit;mid<<=1){//待合并区间的中点
complex Wn( cos(Pi/mid) , type*sin(Pi/mid) ); //单位根
for(int R=mid<<1,j=0;j<limit;j+=R){//R是区间的右端点，j表示前已经到哪个位置了
complex w(1,0);//幂
for(int k=0;k<mid;k++,w=w*Wn){//枚举左半部分
complex x=A[j+k],y=w*A[j+mid+k];//蝴蝶效应
A[j+k]=x+y;
A[j+mid+k]=x-y;
}
}
}
}
//FFT板子
string s1,s2;
int n,m;
complex A[MAXN],B[MAXN];
double f[MAXN];
double g[MAXN];
const double eps = 0.0001;
void solve(){
limit = 1,l=0;
cin>>n>>m;
cin>>s1>>s2;
rep(i,0,n-1) B[i].x = s1[i]-'a'+1;
rep(i,0,m-1) A[i].x = s2[i]-'a'+1;
double T = 0;
//T = sigma A[i]^A[i] i=0~m-1
rep(i,0,m-1) T += A[i].x*A[i].x;
//f[x] = sigma B[i]^B[i] i=0~x
f[0] = B[0].x*B[0].x;
rep(i,1,n-1) f[i] = f[i-1]+B[i].x*B[i].x;
//g[x] = S[i]*B[j] i+j==x
reverse(A,A+m);//S = A.reverse
//FFT预处理
while(limit<=n+m-2) limit<<=1,l++;
rep(i,0,limit-1)
r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) );
FFT(A,1);FFT(B,1);
rep(i,0,limit) A[i]=A[i]*B[i];
FFT(A,-1);
rep(i,0,n-1) g[i] = (int)(A[i].x/limit+0.5);//四舍五入
//T + f(x) - f(x-m) - 2g(x);
double tmp;
rep(x,m-1,n-1){
tmp = T+f[x]-2*g[x];
if(x!=m-1) tmp -= f[x-m];
//cout<<tmp<<' ';
if(fabs(tmp)<eps) cout<<x-(m-1)+1<<endl;//输出匹配上的位置
}
cout<<endl;
}
int main(){
freopen("in.txt","r",stdin);
int z;
cin>>z;
while(z--) solve();
}
//洛谷P4173 残缺的字符串（带通配符的字符串匹配问题）
//这题的简化问题：没有通配符的情况如何用FFT进行O(nlogn)的字符串匹配
//https://www.luogu.com.cn/problem/P4173问题2:带通配符字符串匹配
题目:在普通字符串匹配的基础上，加上了可以匹配任何字符的通配符’*‘
P4173 残缺的字符串 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)
参考:
题解 P4173 【残缺的字符串】 - Ebola 的博客 - 洛谷博客 (luogu.com.cn)
数学原理
在普通的字符串匹配中，我们设字母a到z为1-26。
- 我们设定通配符为0（原因之后会讲
这样一来，先前的匹配函数$C(i,i)=(A_i-B_i)^2$在出现通配符时，不能满足$C(i,i)=0$
故构造$C(i,i)=(A_i-B_i)^2A_iB_i$，这样只要Ai或者Bi中有一个通配符(值为0)，匹配函数的值都为0！非常奇妙！
对应这个新的匹配函数$C$,有$P(k) = \sum^{m-1}_{i=0} C(i,k+i)=\sum^{m-1}_{i=0} (A_i-B_{k+i})^2A_iB_{k+i}$
考虑拆平方化简为:$\sum^{m-1}_{i=0} (A_i^3B_{k+i}+A_iB_{k+i}^3-2A_i^2B_{k+i}^2)$
其实和没有通配符的情况非常类似，依旧可以拆分为三个项：
$\sum^{m-1}_{i=0} (A_i^3B_{k+i})+\sum^{m-1}_{i=0}(A_iB_{k+i}^3)-2\sum^{m-1}_{i=0}(A_i^2B_{k+i}^2)$
同样把A给reverse$\sum^{m-1}_{i=0} (S_{m-1-i}^3B_{k+i})+\sum^{m-1}_{i=0}(S_{m-1-i}B_{k+i}^3)-2\sum^{m-1}_{i=0}(S_{m-1-i}^2B_{k+i}^2)$
在上一题中，前两项还可以用比较简单的方法去处理，第三项才需要FFT，这里三项都需要FFT来处理了😅😅😅
可以进行6次FFT，把三项都加到P上，最后对P进行一次逆变换。
踩的坑
- 在进行完一次A和B的卷积之后，记得要把A和B全部清空，包括n和m之后的数据，因为在多项式转化为点值的时候，n和m往后的位置也会被填上。
代码
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
using namespace std;
int cansel_sync=(ios::sync_with_stdio(0),cin.tie(0),0);
#define ll long long
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
#define INF 0x3f3f3f3f
//FFT板子
const int MAXN = 114514<<4;
const double Pi=acos(-1.0);
struct complex{
double x,y;
complex (double xx=0,double yy=0){x=xx,y=yy;}
};
complex operator + (complex a,complex b){ return complex(a.x+b.x , a.y+b.y);}
complex operator - (complex a,complex b){ return complex(a.x-b.x , a.y-b.y);}
complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);}//不懂的看复数的运算那部分
int l,r[MAXN];
int limit=1;
void FFT(complex *A,int type){
for(int i=0;i<limit;i++)
if(i<r[i]) swap(A[i],A[r[i]]);//求出要迭代的序列
for(int mid=1;mid<limit;mid<<=1){//待合并区间的中点
complex Wn( cos(Pi/mid) , type*sin(Pi/mid) ); //单位根
for(int R=mid<<1,j=0;j<limit;j+=R){//R是区间的右端点，j表示前已经到哪个位置了
complex w(1,0);//幂
for(int k=0;k<mid;k++,w=w*Wn){//枚举左半部分
complex x=A[j+k],y=w*A[j+mid+k];//蝴蝶效应
A[j+k]=x+y;
A[j+mid+k]=x-y;
}
}
}
}
//FFT板子
string s1,s2;
int n,m;
int Aa[MAXN],Bb[MAXN];
complex A[MAXN],B[MAXN],P[MAXN];
double f[MAXN];
double g[MAXN];
vector<int> vec;
const double eps = 0.0001;
inline int idx(char c){
if(c=='*') return 0;
else return c-'a'+1;
}
void solve(){
limit = 1,l=0;
cin>>m>>n;//注意读入顺序，这里是先读入模式串B
cin>>s2>>s1;
rep(i,0,n-1) Bb[i] = idx(s1[i]);
rep(i,0,m-1) Aa[i] = idx(s2[i]);
reverse(Aa,Aa+m);
//S = A.reverse,这里因为三项都要reverse,我直接把原数组给翻了
//FFT预处理
while(limit<=n+m-2) limit<<=1,l++;
rep(i,0,limit-1)
r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) );
//第一项g1[x] = sigma S[i]^3*B[j] i+j=x
rep(i,0,limit-1) A[i]=B[i]=complex();//n和m之后的数据也要清空,之前wa这里
rep(i,0,m-1) A[i] = complex(Aa[i]*Aa[i]*Aa[i],0);
rep(i,0,n-1) B[i] = complex(Bb[i],0);
FFT(A,1);FFT(B,1);
rep(i,0,limit) P[i]=A[i]*B[i];
//第二项g2[x] = sigma S[i]*B[j]^3 i+j=x
rep(i,0,limit-1) A[i]=B[i]=complex();
rep(i,0,m-1) A[i] = complex(Aa[i],0);
rep(i,0,n-1) B[i] = complex(Bb[i]*Bb[i]*Bb[i],0);
FFT(A,1);FFT(B,1);
rep(i,0,limit) P[i]=P[i]+A[i]*B[i];
//第三项g3[x] = sigma S[i]^2*B[j]^2 i+j=x
rep(i,0,limit-1) A[i]=B[i]=complex();
rep(i,0,m-1) A[i] = complex(Aa[i]*Aa[i],0);
rep(i,0,n-1) B[i] = complex(Bb[i]*Bb[i],0);
FFT(A,1);FFT(B,1);
rep(i,0,limit) P[i]=P[i]-complex(2,0)*A[i]*B[i];
FFT(P,-1);
vec.clear();
double tmp;
rep(x,m-1,n-1){
tmp = (int)(P[x].x/limit+0.5);
if(fabs(tmp)<eps) vec.push_back(x-(m-1)+1);//存左端点位置
}
cout<<vec.size()<<endl;
for(auto x:vec) cout<<x<<' ';
cout<<endl;
}
int main(){
//freopen("in.txt","r",stdin);
solve();
}
//洛谷P4173 残缺的字符串（带通配符的字符串匹配问题）
//https://www.luogu.com.cn/problem/P4173
//洛谷开O2才能过杭电多校2021（3） 1003.Forgiving Matching
Problem - 6975 (dingbacode.com)
概要
时间给8s，n和m的总和1e6，每组n和m为1e5。给长为n的字符串B和长为m的字符串A。串中包含通配符*，可以匹配任何符号。
定义一个原谅数k，即允许有k位失配，问对[0,m]的不同k，A有多少个连续子串与B匹配。
分析
这题其实比上面的题都简单。它对时间复杂度的要求很低。故可以分成0-9共10次FFT。
- 分析题目，可以转化为：求位置x的子串$B[x-m-2,x]$与A的不匹配的数量$tmp$，加到差分数组$cf[m-tmp]$上去即可。
我们维护10个数组$A[i][j]$，代表A串的j位置是否为数字i（*存在A[10]中）。
- 我们定义匹配函数$C(i,j)=A_i\cdot B_i$，故全匹配函数为$P(x)=\sum C(i,x+i)$,同样通过反转后卷积来维护，卷积后的数组$G[i][j]$代表$A$和$B[j-m-2,j]$有多少个相同的位置(都是数字i)。
- 通配符另外讨论，有匹配字符数 = 0-9匹配字符数 + A通配符数 + B子串通配符数 - AB位置重合通配符数量
A通配符数扫一遍可以得到，B子串通配符数通过前缀和得，AB位置重合的通配符数量，则和之前的其他字母一样，可以通过FFT维护后查询。
代码
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
int cansel_sync=(ios::sync_with_stdio(0),cin.tie(0),0);
#define ll long long
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
const int MAXN = 114514<<1;
const double Pi=acos(-1.0);//FFT板子
struct complex{
double x,y;
complex (double xx=0,double yy=0){x=xx,y=yy;}
};
complex operator + (complex a,complex b){ return complex(a.x+b.x , a.y+b.y);}
complex operator - (complex a,complex b){ return complex(a.x-b.x , a.y-b.y);}
complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);}//不懂的看复数的运算那部分
int l,r[MAXN<<2];
int limit=1;
void FFT(complex *A,int type){
for(int i=0;i<limit;i++)
if(i<r[i]) swap(A[i],A[r[i]]);//求出要迭代的序列
for(int mid=1;mid<limit;mid<<=1){//待合并区间的中点
complex Wn( cos(Pi/mid) , type*sin(Pi/mid) ); //单位根
for(int R=mid<<1,j=0;j<limit;j+=R){//R是区间的右端点，j表示前已经到哪个位置了
complex w(1,0);//幂
for(int k=0;k<mid;k++,w=w*Wn){//枚举左半部分
complex x=A[j+k],y=w*A[j+mid+k];//蝴蝶效应
A[j+k]=x+y;
A[j+mid+k]=x-y;
}
}
}
}//FFT板子
string s1,s2;
int n,m;
inline int idx(char c){if(c=='*')return 10;return c-'0';}
complex A[12][MAXN<<2],B[12][MAXN<<2];
int cnt[MAXN];//记录数字相同位置数量
int cnttp[MAXN];//通配符位置相同
int sum1[MAXN],sum2;
int cf[MAXN];//差分数组
void solve(){
l=0;limit=1;
cin>>n>>m;
cin>>s1>>s2;
rep(i,0,m) cf[i]=0;//差分数组
while(limit<=n+m-2) limit<<=1,l++;
rep(i,0,limit-1)
r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) );
rep(k,0,10){
rep(i,0,limit-1) A[k][i] = complex();
rep(i,0,limit-1) B[k][i] = complex();
}
rep(i,0,n-1){
B[idx(s1[i])][i]=complex(1);
cnt[i]=cnttp[i]=0;
if(i==0) sum1[i] = (s1[i]=='*');
else sum1[i]=sum1[i-1]+(s1[i]=='*');
}
sum2=0;
rep(i,0,m-1){
A[idx(s2[i])][i]=complex(1);
sum2 += (s2[i]=='*');
}
//cnt[x] = sigma A[i]*B[j] i+j=x
rep(i,0,10){//普通字符匹配的位置
reverse(A[i],A[i]+m);
FFT(A[i],1);FFT(B[i],1);
rep(j,0,limit-1) A[i][j] = A[i][j]*B[i][j];
FFT(A[i],-1);
if(i!=10) rep(j,m-1,n-1) cnt[j] += (int)(A[i][j].x/limit+0.5);
else rep(j,m-1,n-1) cnttp[j] += (int)(A[i][j].x/limit+0.5);
}
//匹配字符数 = cnt[j] + A通配 + B通配 - AB重合通配(FFT)
int tmp;
rep(i,m-1,n-1){
tmp = cnt[i] + sum2 + sum1[i] - cnttp[i];
if(i!=m-1) tmp -= sum1[i-m];//其实就是sum1[i]-sum1[i-m]，防止越界用
cf[m-tmp]++;//更新差分数组
}
int now = 0;
rep(i,0,m){
now += cf[i];
cout<<now<<endl;
}
}
int main(){
int z;
cin>>z;
while(z--) solve();
}
