洛谷P4173 残缺的字符串【FFT】【带通配符单模式串匹配】

时空限制 1000ms / 128MB

题目描述

很久很久以前，在你刚刚学习字符串匹配的时候，有两个仅包含小写字母的字符串A和B，其中A串长度为m，B串长度为n。可当你现在再次碰到这两个串时，这两个串已经老化了，每个串都有不同程度的残缺。

你想对这两个串重新进行匹配，其中A为模板串，那么现在问题来了，请回答，对于B的每一个位置i，从这个位置开始连续m个字符形成的子串是否可能与A串完全匹配?

输入格式：

第一行包含两个正整数m,n，分别表示AA串和BB串的长度。
第二行为一个长度为m的字符串A。
第三行为一个长度为n的字符串B。

两个串均仅由小写字母和$\ast$号组成，其中*号表示相应位置已经残缺。

输出格式：

第一行包含一个整数k，表示BB串中可以完全匹配A串的位置个数。
若k > 0，则第二行输出k个正整数，从小到大依次输出每个可以匹配的开头位置（下标从1开始）

说明

100％的数据满足$1\le m\le n\le 300000$

---

题目分析

(以下记$S(x)$为字符串$S$第$x$位的ASCII码，且字符串都从第0位开始)

无通配符匹配单模式串匹配

假设先不考虑通配符$\ast$
定义$A$的第$x$位与$B$的第$y$位匹配当且仅当$A(x)-B(y)=0$
那么可以进一步定义针对一整个字符串的匹配函数$P(x)={\sum }_{i=0}^{m-1}[A(i)-B(x-m+1+i){]}^2$
当且仅当$P(x)=0$时，$B$中以第$x$位结束的连续$m$位，与$A$完全匹配

这里为什么还需要一个平方呢，因为可能有负数出现
例如$A="ab",B="ba"$，$P(x)=1+(-1)=0$
所以需要乘方全部转化为整正数保证正确性

现在问题转化为如何快速求解所有$P(x)$
尝试将$A$翻转，记为$A^{\prime }$，则有$A(x)=A^{\prime }(m-x-1)$
于是更改$P(x)={\sum }_{i=0}^{m-1}[A^{\prime }(m-i-1)-B(x-m+1+i){]}^2$
将其展开再稍作变换
$$P(x)=\sum _{i=0}^{m-1}{A}^{\prime }(m-i-1{)}^2+\sum _{i=0}^{m-1}B(x-m+1+i{)}^2-2\ast \sum _{i=0}^{m-1}{A}^{\prime }(m-i-1)B(x-m+1+i)$$
$$P(x)=\sum _{i=0}^{m-1}{A}^{\prime }(i{)}^2+\sum _{i=0}^{m-1}B(x-m+1+i{)}^2-2\ast \sum _{i+j=x}{A}^{\prime }(i)B(j)$$
我们发现第一项为定值，第二项前缀和即可，第三项是卷积形式
所以可以直接FFT优化

带通配符匹配单模式串匹配

受上述启发，我们可以很快修改出一个新的匹配函数
若$A(x){=}^{\prime }{\ast }^{\prime }$，令$A(x)=0$
$$P(x)=\sum _{i=0}^{m-1}[A^{\prime }(m-i-1)-B(x-m+1+i){]}^2{A}^{\prime }(m-i-1)B(x-m+1+i)$$
再次大力展开
$$P(x)=\sum _{i=0}^{m-1}{A}^{\prime }(m-i-1{)}^{3}B(x-m+1+i)+\sum _{i=0}^{m-1}B(x-m+1+i{)}^3{A}^{\prime }(m-i-1)-2\ast \sum _{i=0}^{m-1}{A}^{\prime }(m-i-1{)}^{2}B(x-m+1+i{)}^2$$
$$P(x)=\sum _{i+j=x}{A}^{\prime }(i{)}^{3}B(j)+\sum _{i+j=x}B(i{)}^3{A}^{\prime }(j)-2\ast \sum _{i+j=x}{A}^{\prime }(i{)}^{2}B(j{)}^2$$

我们惊讶的的发现这不就是三个卷积吗
再次FFT大力优化即可

---

#include<iostream>
#include<cmath>
#include<algorithm>
#include<map>
#include<cstring>
#include<cstdio>
using namespace std;
typedef long long lt;
typedef double dd;
#define eps 1e-7

int read()
{
    int f=1,x=0;
    char ss=getchar();
    while(ss<'0'||ss>'9'){if(ss=='-')f=-1;ss=getchar();}
    while(ss>='0'&&ss<='9'){x=x*10+ss-'0';ss=getchar();}
    return f*x;
}

const dd Pi=acos(-1.0);
const int maxn=1200010;
int n,m;
struct complex{
    dd x,y;
    complex(dd _x=0,dd _y=0){ x=_x; y=_y;}
}A[maxn],B[maxn],C[maxn];
int lim=1,L,R[maxn];
char pat[maxn],txt[maxn];
int a[maxn],b[maxn];
int ans[maxn],cnt;

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);}

void FFT(complex* a,int opt)
{
    for(int i=0;i<lim;++i)
    if(i<R[i]) swap(a[i],a[R[i]]);
    
    for(int i=1;i<lim;i<<=1)
    {
        complex wn(cos(Pi/i),opt*sin(Pi/i));
        for(int j=0;j<lim;j+=(i<<1))
        {
            complex w(1,0);
            for(int k=0;k<i;++k)
            {
                complex nx=a[j+k],ny=w*a[i+j+k];
                a[j+k]=nx+ny;
                a[i+j+k]=nx-ny;
                w=w*wn;
            }
        }
    }
    
    if(opt==-1)
    for(int i=0;i<=lim;++i) a[i].x/=(dd)lim;
}

void solve1()
{
    for(int i=0;i<lim;++i) 
    A[i]=complex(a[i]*a[i]*a[i],0),B[i]=complex(b[i],0);
    FFT(A,1); FFT(B,1);
    for(int i=0;i<lim;++i) C[i]=C[i]+A[i]*B[i];
}

void solve2()
{
    for(int i=0;i<lim;++i) 
    A[i]=complex(a[i]*a[i],0),B[i]=complex(b[i]*b[i],0);
    FFT(A,1); FFT(B,1);
    for(int i=0;i<lim;++i) C[i]=C[i]-(A[i]*B[i]*complex(2,0));
}

void solve3()
{
    for(int i=0;i<lim;++i) 
    A[i]=complex(a[i],0),B[i]=complex(b[i]*b[i]*b[i],0);
    FFT(A,1);FFT(B,1);
    for(int i=0;i<lim;++i) C[i]=C[i]+A[i]*B[i];
}

int main()
{
    n=read();m=read();
    scanf("%s%s",pat,txt);
    reverse(pat,pat+n);
    
    while(lim<=n+m-2) lim<<=1,L++;
    for(int i=0;i<lim;++i)
    R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    
    for(int i=0;i<n;++i) a[i]=(pat[i]!='*')?(pat[i]-'a'+1):0;
    for(int i=0;i<m;++i) b[i]=(txt[i]!='*')?(txt[i]-'a'+1):0;
    
    solve1(); solve2(); solve3();
    
    FFT(C,-1);
    for(int i=n-1;i<m;++i)
    if(fabs(floor(C[i].x+0.5))<eps)
    ans[++cnt]=i-n+2;
    
    printf("%d\n",cnt);
    for(int i=1;i<=cnt;++i) printf("%d ",ans[i]);
    return 0;
}