Bzoj4503 两个串 FFT

最新推荐文章于 2021-03-22 10:46:30 发布

原创最新推荐文章于 2021-03-22 10:46:30 发布 · 480 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#FFT

FFT 专栏收录该内容

1 篇文章

订阅专栏

本文介绍了一种利用快速傅立叶变换（FFT）进行模糊字符串匹配的方法，解决了含有通配符'?'的问题，通过数学公式和代码实现详细阐述了该算法的工作原理。

问题描述

给定两个串s1，s2，问你s2在s1中出现的次数和位置，s2中有‘？’字符，可以匹配任何字符。字符串长度为 ${10}^{5}$ 。

题解

我们考虑没有‘？’时，两个长度相等的字符串相等，当且仅当 $\sum^{n-1}_{i=0}{(s1[i]-s2[i])}^{2}$ 。那么，s2有‘？’怎么办呢。只需要把s2赋为0，式子就变成了 $\sum^{n-1}_{i=0}{s2[i]\times(s1[i]-s2[i])}^{2}$ 。考虑将s2翻转，这个式子就变成了一个卷积的形式，使用FFT即可。

代码

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define N 800010
using namespace std;
struct cp{double x,y;};
double pi=acos(-1);
int n,m,M,tot,Ans[N];
long long sum;
char s1[N],s2[N];
cp A[N],B[N],C[N],D[N],E[N],cur[N];
cp operator *(cp x,cp y){return (cp){x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x};}
cp operator +(cp x,cp y){return (cp){x.x+y.x,x.y+y.y};}
cp operator -(cp x,cp y){return (cp){x.x-y.x,x.y-y.y};}
cp operator *(int x,cp y){return (cp){x*y.x,x*y.y};}
void FFT(cp *a,int n,int fl)
{
    for (int i=n>>1,j=1;j<n;j++)
    {
        if (i<j) swap(a[i],a[j]);
        int k=n>>1;
        for (;k&i;i^=k,k>>=1);i^=k;
    }
    for (int m=2;m<=n;m<<=1)
    {
        cp w=(cp){cos(2*pi*fl/m),sin(2*pi*fl/m)};
        cur[0]=(cp){1,0};
        for (int i=1;i<m;i++) cur[i]=cur[i-1]*w;
        for (int i=0;i<n;i+=m)
            for (int j=i;j<i+(m>>1);j++)
            {
                cp u=a[j],v=a[j+(m>>1)]*cur[j-i];
                a[j]=u+v;
                a[j+(m>>1)]=u-v;
            }
    }
}
int main()
{
    scanf("%s",s1+1);n=strlen(s1+1);
    scanf("%s",s2+1);m=strlen(s2+1);
    for (int i=1;i<=n;i++)
    {
        B[i].x=s1[i];
        C[i].x=s1[i]*s1[i];
    }
    for (int i=1;i<=m;i++) if (s2[i]!='?')
    {
        A[i].x=s2[i];
        D[i].x=s2[i]*s2[i];
    }
    reverse(A+1,A+m+1);
    reverse(D+1,D+m+1);
    for (int i=1;i<=m;i++) sum+=A[i].x*A[i].x*A[i].x;
    M=1;while (M<=(n+m)*2) M<<=1;
    FFT(A,M,1);FFT(B,M,1);FFT(C,M,1);FFT(D,M,1);
    for (int i=0;i<M;i++) E[i]=A[i]*C[i]-2*B[i]*D[i];
    FFT(E,M,-1);
    for (int i=0;i<M;i++) E[i].x=fabs(E[i].x/M+sum);
    for (int i=m+1;i<=n+1;i++) if (E[i].x<=0.5) Ans[++tot]=i-m;
    printf("%d\n",tot);
    for (int i=1;i<=tot;i++) printf("%d\n",Ans[i]-1);
}