COGS2216 你猜是不是KMP

第一道自己写的FFT......

不知为啥这题在网上找不到题解......真是麻烦，害得我推了半天......

还是写个简要题解吧......

首先把S和T拆成序列，a~z分别对应成1~26，?是0，设拆成的两个序列分别为A和B。

如果从i开始可以S匹配上T，那么就有

egin{align}sum_{j=0}^{m-1}[|A_{i+j}-B_j|=0 space or space B_j=0]=mend{align}

换一个等价写法：

egin{align}sum_{j=0}^{m-1}(A_{i+j}-B_j)^2 B_j=0end{align}

求出所有i对应的值就可以得知每一位是否匹配了。

展开之后发现并没有什么用，这个式子只能$O(n^2)$计算，还不如朴素匹配呢，这玩意儿有个卵用啊(╯‵□′)╯︵┻━┻

尝试把T反转，原式变为

egin{align}sum_{j=0}^{m-1}(A_i-B_{m-j-1})^2 B_{m-j-1}end{align}

展开得

egin{align}sum_{j=0}^{m-1}A_{i+j}^2 B_{m-j-1}-2A_{i+j}B_{m-j-1}^2+B_{m-j-1}^3end{align}

有没有发现这个式子看着有点眼熟......

前两项下标之和都是i+m-1，因此可以看成一个卷积形式，而令最后一项再卷上一个各项全为1的序列(显然不影响结果是吧......)，也是卷积，而下标之和是i+m-1，因此定义我们得到的结果为

egin{align}C_{i+m-1}=sum_{j=0}^{m-1}A_{i+j}^2 B_{m-j-1}-2A_{i+j}B_{m-j-1}^2+1B_{m-j-1}^3end{align}

令$D_i=A_i^2,E_i=B_i^2,F_i=B_i^3,I_i=1$，再写成卷积形式就是

egin{align}C=D*B-2A*E+F*Iend{align}

分别求出三个卷积之后加一下就行了，FFT加速卷积即可，复杂度$O(nlogn)$。

然而常数大如狗，跑的比bitset慢到不知哪儿去了

求出三个卷积的和之后扫一遍判断哪些i对应的$C_{i+m-1}$是0，是的话说明两串在i位置匹配，否则不匹配。

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const int maxn=270010;
const double pi=acos(-1.0),eps=1e-4;
struct Complex{
    double a,b;
    Complex(double a=0.0,double b=0.0):a(a),b(b){}
    Complex operator+(const Complex &x)const{return Complex(a+x.a,b+x.b);}
    Complex operator-(const Complex &x)const{return Complex(a-x.a,b-x.b);}
    Complex operator*(const Complex &x)const{return Complex(a*x.a-b*x.b,a*x.b+b*x.a);}
    Complex &operator*=(const Complex &x){return *this=*this*x;}
}A[maxn],B[maxn],C[maxn],D[maxn],E[maxn],F[maxn],I[maxn];//D是A每项取平方，E是b每项取平方，F是b每项取立方，I是伪单位元
void FFT(Complex*,int,int);
char S[maxn],T[maxn];
int n,m,N=1,ans=0;
int main(){
    freopen("guess.in","r",stdin);
    freopen("guess.out","w",stdout);
    scanf("%s%s",S,T);
    n=strlen(S);
    m=strlen(T);
    while(N<n+m)N<<=1;
    for(int i=0;i<N;i++)I[i].a=1.0;
    for(int i=0;i<n;i++){
        A[i].a=S[i]-'a'+1;
        D[i].a=A[i].a*A[i].a;
    }
    reverse(T,T+m);
    for(int i=0;i<m;i++){
        B[i].a=(T[i]=='?')?0:T[i]-'a'+1;
        E[i].a=B[i].a*B[i].a;
        F[i].a=B[i].a*E[i].a;
    }
    FFT(A,N,1);
    FFT(B,N,1);
    FFT(D,N,1);
    FFT(E,N,1);
    FFT(F,N,1);
    FFT(I,N,1);//woc，这么多FFT慢不慢啊
    for(int i=0;i<N;i++){
        D[i]*=B[i];
        A[i]*=E[i];
        F[i]*=I[i];//printf("I[%d]=(%lf,%lf)
",i,I[i].a,I[i].b);
    }
    FFT(A,N,-1);
    FFT(D,N,-1);
    FFT(F,N,-1);
    for(int i=0;i<n;i++){
        //printf("i=%d i+m-1=%d  A[i+m-1]=(%lf,%lf)  D[i+m-1]=(%lf,%lf)  F[i+m-1]=(%lf,%lf)
",i,i+m-1,A[i+m-1].a,A[i+m-1].b,D[i+m-1].a,D[i+m-1].b,F[i+m-1].a,F[i+m-1].b);
        C[i].a=D[i].a-A[i].a*2+F[i].a;    
    }
    for(int i=0;i+m<=n;i++)if(C[i+m-1].a<eps)ans++;
    printf("%d
",ans);
    for(int i=0;i+m<=n;i++)if(C[i+m-1].a<eps)printf("%d
",i);
    return 0;
}
void FFT(Complex *A,int n,int tp){
    for(int i=1,j=0,k;i<n-1;i++){
        k=N;
        do{
            k>>=1;
            j^=k;
        }while(j<k);
        if(i<j)swap(A[i],A[j]);
    }
    for(int k=2;k<=n;k<<=1){
        Complex wn(cos(-tp*2*pi/k),sin(-tp*2*pi/k));
        for(int i=0;i<n;i+=k){
            Complex w(1.0,0.0);
            for(int j=0;j<(k>>1);j++,w*=wn){
                Complex a=A[i+j],b=w*A[i+j+(k>>1)];
                A[i+j]=a+b;
                A[i+j+(k>>1)]=a-b;
            }
        }
    }
    if(tp<0)for(int i=0;i<n;i++)A[i].a/=n;
}
/*
把T反转，把S和T变成多项式a和b，令
c[i+m-1]=sum{(a[i+j]-b[m-j-1])^2*b[m-j-1]}
显然只有c是0的位置两串才会匹配，展开得
c[i+m-1]=sum{a[i+j]^2*b[m-j-1]-2*a[i+j]*b[m-j-1]^2+b[m-j-1]^3}
前两项是卷积形式，第三项乘上一个伪单位元也是卷积形式，FFT加速即可
*/

ps：其实对$I$跑的那个DFT完全没必要，手动逐项赋成1就行了，然后发现$F$卷完了$I$根本没有什么变化，所以直接对$F$做一下DFT再IDFT回去即可，把$I$写进去只是为了理解式子方便......

UPD：脑残了……那个ps是错的，忽视掉吧……