bzoj4503: 两个串

我又变成sb了

首先通过%可以得到一个这样的柿子：

设s1是模版串，s2是匹配串，j表示模版串以j为起始的子串

fj(0<j<=s1len-s2len) =sigema(i=0~s2len-1) (s2[i]-s1[j+i])^2*s2[i]

当fj=0时j是一个解。

既然要fft，这个加号很不爽，我们把s2翻转：

fj(0<j<=s1len-s2len) =sigema(i=0~s2len-1) (s2[s2len-i]-s1[j+i])^2*s2[s2len-i]

改为枚举s2len-i

fj(0<j<=s1len-s2len) =sigema(i=0~s2len-1) (s2[i]-s1[j-i+s2len])^2*s2[i]

这个时候发现s2len没法消啊，改一下定义，j表示以j为结尾的子串

fj(s2len-1<=j<=s1len-1) =sigema(i=0~s2len-1) (s2[i]-s1[j-i])^2*s2[i]

把它拆开跑三次fft就可以了。

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
using namespace std;
const double pi=acos(-1.0);

struct complex
{
    double r,i;
    complex(){}
    complex(double R,double I){r=R,i=I;}
    friend complex operator +(complex x,complex y){return complex(x.r+y.r,x.i+y.i);}
    friend complex operator -(complex x,complex y){return complex(x.r-y.r,x.i-y.i);}
    friend complex operator *(complex x,complex y){return complex(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);}
}A[410000],B[410000],C[410000]; int Re[410000];
void fft(complex *a,int n,int op)
{
    for(int i=0;i<n;i++)
        if(i<Re[i])swap(a[i],a[Re[i]]);
    
    for(int i=1;i<n;i<<=1)
    {
        complex wn(cos(pi/i),sin(op*pi/i));
        for(int j=0;j<n;j+=(i<<1))
        {
            complex w(1,0);
            for(int k=0;k<i;k++,w=w*wn)
            {
                complex t1=a[j+k],t2=a[j+k+i]*w;
                a[j+k]=t1+t2;a[j+k+i]=t1-t2;
            }
        }
    }
}

char ss[410000];
int s1len,s1[410000],s2len,s2[410000],f[410000];
int aslen,as[410000];
int main()
{
    scanf("%s",ss);s1len=strlen(ss);
    for(int i=0;i<s1len;i++)
        if('a'<=ss[i]&&ss[i]<='z')s1[i]=ss[i]-'a'+1;
        else s1[i]=0;
    scanf("%s",ss);s2len=strlen(ss);
    for(int i=0;i<s2len;i++)
        if('a'<=ss[i]&&ss[i]<='z')s2[s2len-1-i]=ss[i]-'a'+1;
        else s2[s2len-1-i]=0;
    //~~~~~~~~~~~~~~~~~~~~read~~~~~~~~~~~~~~~~~~~~~~~~~~~    
    
    int n,m=s1len+s2len-2,L=0;
    for(n=1;n<=m;n*=2)L++;
    for(int i=1;i<=n;i++)Re[i]=(Re[i>>1]>>1)|((i&1)<<(L-1));
    memset(f,0,sizeof(f));
    //~~~~~~~~~~~~~~~~~~~~~init~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    for(int i=0;i<=n;i++)
        A[i].r=s2[i]*s2[i]*s2[i],B[i].r=1;
    fft(A,n,1),fft(B,n,1);
    for(int i=0;i<n;i++)C[i]=A[i]*B[i];
    fft(C,n,-1);
    for(int i=s2len-1;i<s1len;i++)f[i]+=int(C[i].r/n+0.5);
    
    memset(A,0,sizeof(A));
    memset(B,0,sizeof(B));
    for(int i=0;i<=n;i++)
        A[i].r=2*s2[i]*s2[i],B[i].r=s1[i];
    fft(A,n,1),fft(B,n,1);
    for(int i=0;i<n;i++)C[i]=A[i]*B[i];
    fft(C,n,-1);
    for(int i=s2len-1;i<s1len;i++)f[i]-=int(C[i].r/n+0.5);
    
    memset(A,0,sizeof(A));
    memset(B,0,sizeof(B));
    for(int i=0;i<=n;i++)
        A[i].r=s2[i],B[i].r=s1[i]*s1[i];
    fft(A,n,1),fft(B,n,1);
    for(int i=0;i<n;i++)C[i]=A[i]*B[i];
    fft(C,n,-1);
    for(int i=s2len-1;i<s1len;i++)f[i]+=int(C[i].r/n+0.5);
    
    aslen=0;
    for(int i=s2len-1;i<s1len;i++)
        if(f[i]==0)as[++aslen]=i-s2len+1;
    printf("%d
",aslen);
    for(int i=1;i<=aslen;i++)printf("%d
",as[i]);
    return 0;
}