[BZOJ]4503 两个串：我的第一次FFT尝试

序言

接触FFT半年了，水平一直停留在只会写个高精度乘法那个层次，所以就试着做了一道FFT的题。

要是我没感觉错的话，这是一道权限题，然而我并没有权限号，我只是写了个暴力对拍了一下…

题目

[BZOJ]4503 两个串

Time Limit: 20 Sec；Memory Limit: 256 MB；Submit: 398 Solved: 190；

Description

兔子们在玩两个串的游戏。给定两个字符串S和T，兔子们想知道T在S中出现了几次，
分别在哪些位置出现。注意T中可能有“?”字符，这个字符可以匹配任何字符。

Input

两行两个字符串，分别代表S和T

Output

第一行一个正整数k，表示T在S中出现了几次
接下来k行正整数，分别代表T每次在S中出现的开始位置。按照从小到大的顺序输出，S下标从0开始。

Sample Input

bbabaababaaaaabaaaaaaaabaaabbbabaaabbabaabbbbabbbbbbabbaabbbababababbbbbbaaabaaabbbbbaabbbaabbbbabab

a?aba?abba

Sample Output

0

HINT

S 长度不超过 10^5， T 长度不会超过 S。 S 中只包含小写字母， T中只包含小写字母和“?”

Source

数学 FFT 字符串 脑洞题

思路

首先，因为有通配符，所以不能KMP。而又因为数据范围1e5，所以时间复杂度被锁定在O(nlogn)左右，暴力匹配显然也是不可能的，考虑通过FFT强行优化。

（下文中，l1表示S串长度，l2表示T串长度。S表示文本串，T表示模式串。）

然后，在S中的某一个位置x，T得到成功匹配的充要条件是什么。先考虑没有通配符的情况：
(c1是S中的某个字符，c2是T中的某个字符。)

c1=c2⇔(c1−c2)2=0

Sx...x+l2−1=T⇔∑i=0l2−1(Si+x−Ti)2=0

如果有通配符，只要是通配符，那么就一定可以匹配，不妨先把通配符都改成0。

c1=c2⇔c2(c1−c2)2=0

Sx...x+l2−1=T⇔∑i=0l2−1Ti(Si+x−Ti)2=0

并没发现什么问题，平方展开试一试。

Sx...x+l2−1=T⇔∑i=0l2−1Ti(S2i+x+T2i−2Si+xTi)=0

Sx...x+l2−1=T⇔∑i=0l2−1TiS2i+x+T3i−2Si+xT2i=0

请大家强行脑补，回忆一下卷积的定义是什么？

C=A∗B⇔Ck=∑i+j=kAi×Bj

在计算卷记得某一个位置k时，被使用到的i和j的和始终是一个定值。但本题中，两个相乘的数的下标始终是，i和i+x，他们的差是一个定值，但和不是。怎么办呢？为什么不把T字符串翻转一下！

令T′i=Tl2−i−1

这样，T’就是T翻转之后的结果。我们可以试着用T’去替代原式中的T,并且令其为f(x)。

f(x)=∑i=0l2−1T′l2−i−1S2i+x+T′3l2−i−1−2Si+xT′2l2−i−1=0

f(x)=(∑i=0l2−1T′l2−i−1S2i+x)+(∑i=0l2−1T′3l2−i−1)−2(∑i=0l2−1Si+xT′2l2−i−1)=0

这次是不是非常和谐了，当x一定时，被相乘的两个数的下标分别为l2-i-1和i+x，它们的和始终等于l2+x-1，是一个定值，满足卷积的性质。

把上式看成三个部分：

中间的那部分：

∑i=0l2−1T′3l2−i−1

是一个定值，可以直接求出。

左边的那部分：

∑i=0l2−1T′l2−i−1S2i+x

构造一个数组S’令它为S的平方，则有：

S′i=S2i

令S’和T的卷积为A，则有：

Ak=∑i+j=kS′jT′i

根据定义有：

Al2+x−1=∑i=0l2+x−1S′l2+x−1−iT′i=(∑i=0l2−1S′l2+x−1−iT′i)+(∑i=l2l2+x−1S′l2+x−1−iT′i)

因为T’的长度只有l2，我们默认T’在l2及以后的位置的值都是零（也就是假定都是通配符，可以和任何字符匹配）。那么这个式子的后半部分一定等于零，即：

Al2+x−1=∑i=0l2+x−1S′l2+x−1−iT′i=(∑i=0l2−1S′l2+x−1−iT′i)+(∑i=l2l2+x−1S′l2+x−1−iT′i)=∑i=0l2−1S′l2+x−1−iT′i

尝试用l2-i-1替换i：i=0时，l2-i-1=l2-1：i=l2-1时，l2-i-1=0。

Al2+x−1=∑i=0l2−1S′l2+x−1−iT′i=∑l2−i−1=0l2−i−1=l2−1S′l2+x−1−(l2−i−1)T′(l2−i−1)=∑l2−i−1=0l2−i−1=l2−1S′l2+x−1−(l2−i−1)T′l2−i−1

=∑i=l2−1i=0T′l2−i−1S′i+x=∑i=0l2−1T′l2−i−1S2i+x

惊人的事情出现了：T’和S’的卷积A中的第l2-x+1位的值，恰好与f(x)中的左边部分相等。

同理，你也可以令T”表示T’的平方，即：

T′′i=T′2i

令T”与S的卷积为B，用同样的方法也可以证得B中的第l2-i+1位的值，恰好与f(x)中的右边部分相等。

这样的话，你只需要先预处理出中间的部分记为W，然后求一下S’和T’的卷积A，再求一下S和T”的卷积B就可以得出所谓f(x)。

f(x)=Al2+x−1−2×Bl2+x−1+W

用FFT求两次卷积，时间复杂度为O(nlogn)。其他操作的时间复杂度都为O(n)，总时间复杂度为O(nlogn)。

代码

我的代码不能保证正确性，但是通过了我自己的对拍。

#include<cstdio>#include<cstdlib>#include<algorithm>#include<cmath>#include<queue>#include<complex>using namespace std;typedef complex<double> cd;const int maxl=262145;const double PI=acos(-1.0);int rev[maxl];void fft(cd* a,int n,int dft){    for(int i=0;i<n;i++){        if(i<rev[i])            swap(a[i],a[rev[i]]);    }    for(int step=1;step<n;step<<=1){        cd wn=exp(cd(0,dft*PI/step));        for(int j=0;j<n;j+=step<<1){            cd wnk(1,0);            for(int k=j;k<j+step;k++){                cd x=a[k];                cd y=wnk*a[k+step];                a[k]=x+y;                a[k+step]=x-y;                wnk*=wn;            }        }    }    if(dft==-1){        for(int i=0;i<n;i++)            a[i]/=n;    }}void get_rev(int bit){    for(int i=0;i<(1<<bit);i++){        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));    }}char S[maxl],T[maxl];//emptyinline double sqr(double x){return x*x;}inline double cub(double x){return x*x*x;}double cubeB=0;cd a[maxl],b[maxl];cd sqrA[maxl],sqrB[maxl];cd bSqrA[maxl],aSqrB[maxl];int tran_str(char* s){    for(int i=0;;i++){        if(s[i]==0)            return i;//return len        s[i]=(s[i]=='?')?0:(s[i]-'a'+1);    }}int l1,l2;bool check(int x){    return int((cubeB+bSqrA[l2+x-1]-2.0*aSqrB[l2+x-1]).real()+0.5)==0;}int main(){    scanf("%s%s",S,T);    l1=tran_str(S);l2=tran_str(T);    reverse(T,T+l2);    int bit=1,s=2;    for(bit=1;(1<<bit)<l1+l2-1;bit++)s<<=1;    get_rev(bit);    for(int i=0;i<l1;i++){        a[i]=S[i];sqrA[i]=sqr(S[i]);    }    for(int i=0;i<l2;i++){        b[i]=T[i];sqrB[i]=sqr(T[i]);        cubeB+=cub(T[i]);    }    fft(sqrA,s,1);fft(b,s,1);    for(int i=0;i<s;i++)bSqrA[i]=sqrA[i]*b[i];    fft(bSqrA,s,-1);//b*(a^2)    fft(a,s,1);fft(sqrB,s,1);    for(int i=0;i<s;i++)aSqrB[i]=a[i]*sqrB[i];    fft(aSqrB,s,-1);//a*(b^2)    int cnt=0;queue<int>ans;    for(int x=0;x<=l1-l2;x++){        if(check(x)){            cnt++;            ans.push(x);        }    }    printf("%d\n",cnt);    while(!ans.empty()){        printf("%d\n",ans.front());        ans.pop();    }    return 0;}

我的数据生成器的代码：

#include<cstdio>#include<cstdlib>#include<ctime>#include<algorithm>using namespace std;const int maxl=10001;char S[maxl],T[maxl];int rand(int L,int R){    return rand()%(R-L+1)+L;}int main(){    srand(time(NULL));    int l1=rand(100,4000),l2=rand(2,10);    for(int i=0;i<l1;i++)S[i]=rand('a','b');    //因为数据是随机的，如果字母的范围太大可能最后成功匹配的次数很少    //我认为只保留两种字母并不影响程序正确性的测试    int st=rand(0,l1-l2-1);    for(int i=0;i<l2;i++){        T[i]=S[st+i];        if(!rand(0,4))T[i]='?';    }    printf("%s\n%s\n",S,T);    return 0;}

我的暴力的代码：

#include<cstdio>#include<cstdlib>#include<queue>#include<algorithm>using namespace std;const int maxl=10001;char S[maxl],T[maxl];bool check(int x){    for(int i=0;T[i];i++){        if(!(T[i]=='?' || S[x+i]==T[i])){            return false;        }    }    return true;}int main(){    scanf("%s%s",S,T);    int cnt=0;queue<int>ans;    for(int i=0;S[i];i++){        if(check(i)){            cnt++;            ans.push(i);        }    }    printf("%d\n",cnt);    while(!ans.empty()){        printf("%d\n",ans.front());        ans.pop();    }    return 0;}

尾声

上文中的这种 通过反转一个向量 从而 把对应相乘操作转化成卷积的方法是一种使用FFT时比较常见的套路。其主要思想是把对应相乘位置的两个下标“差一定”的关系转化成两个下标“和一定”的关系，从而利用卷积处理，然后再利用FFT优化卷积。

上文中，如有谬误，敬请指正。