COGS2216 你猜是不是KMP
第一道自己写的FFT......
不知为啥这题在网上找不到题解......真是麻烦，害得我推了半天......
还是写个简要题解吧......
首先把S和T拆成序列，a~z分别对应成1~26，?是0，设拆成的两个序列分别为A和B。
如果从i开始可以S匹配上T，那么就有
\begin{align}\sum_{j=0}^{m-1}[|A_{i+j}-B_j|=0 \space or \space B_j=0]=m\end{align}
换一个等价写法：
\begin{align}\sum_{j=0}^{m-1}(A_{i+j}-B_j)^2 B_j=0\end{align}
求出所有i对应的值就可以得知每一位是否匹配了。
展开之后发现并没有什么用，这个式子只能$O(n^2)$计算，还不如朴素匹配呢，这玩意儿有个卵用啊(╯‵□′)╯︵┻━┻
尝试把T反转，原式变为
\begin{align}\sum_{j=0}^{m-1}(A_i-B_{m-j-1})^2 B_{m-j-1}\end{align}
展开得
\begin{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}^3\end{align}
有没有发现这个式子看着有点眼熟......
前两项下标之和都是i+m-1，因此可以看成一个卷积形式，而令最后一项再卷上一个各项全为1的序列(显然不影响结果是吧......)，也是卷积，而下标之和是i+m-1，因此定义我们得到的结果为
\begin{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}^3\end{align}
令$D_i=A_i^2,E_i=B_i^2,F_i=B_i^3,I_i=1$，再写成卷积形式就是
\begin{align}C=D*B-2A*E+F*I\end{align}
分别求出三个卷积之后加一下就行了，FFT加速卷积即可，复杂度$O(nlogn)$。
然而常数大如狗，跑的比bitset慢到不知哪儿去了
求出三个卷积的和之后扫一遍判断哪些i对应的$C_{i+m-1}$是0，是的话说明两串在i位置匹配，否则不匹配。
1 #include<cstdio> 2 #include<cstring> 3 #include<cmath> 4 #include<algorithm> 5 using namespace std; 6 const int maxn=270010; 7 const double pi=acos(-1.0),eps=1e-4; 8 struct Complex{ 9 double a,b; 10 Complex(double a=0.0,double b=0.0):a(a),b(b){} 11 Complex operator+(const Complex &x)const{return Complex(a+x.a,b+x.b);} 12 Complex operator-(const Complex &x)const{return Complex(a-x.a,b-x.b);} 13 Complex operator*(const Complex &x)const{return Complex(a*x.a-b*x.b,a*x.b+b*x.a);} 14 Complex &operator*=(const Complex &x){return *this=*this*x;} 15 }A[maxn],B[maxn],C[maxn],D[maxn],E[maxn],F[maxn],I[maxn];//D是A每项取平方，E是b每项取平方，F是b每项取立方，I是伪单位元 16 void FFT(Complex*,int,int); 17 char S[maxn],T[maxn]; 18 int n,m,N=1,ans=0; 19 int main(){ 20 freopen("guess.in","r",stdin); 21 freopen("guess.out","w",stdout); 22 scanf("%s%s",S,T); 23 n=strlen(S); 24 m=strlen(T); 25 while(N<n+m)N<<=1; 26 for(int i=0;i<N;i++)I[i].a=1.0; 27 for(int i=0;i<n;i++){ 28 A[i].a=S[i]-'a'+1; 29 D[i].a=A[i].a*A[i].a; 30 } 31 reverse(T,T+m); 32 for(int i=0;i<m;i++){ 33 B[i].a=(T[i]=='?')?0:T[i]-'a'+1; 34 E[i].a=B[i].a*B[i].a; 35 F[i].a=B[i].a*E[i].a; 36 } 37 FFT(A,N,1); 38 FFT(B,N,1); 39 FFT(D,N,1); 40 FFT(E,N,1); 41 FFT(F,N,1); 42 FFT(I,N,1);//woc，这么多FFT慢不慢啊 43 for(int i=0;i<N;i++){ 44 D[i]*=B[i]; 45 A[i]*=E[i]; 46 F[i]*=I[i];//printf("I[%d]=(%lf,%lf)\n",i,I[i].a,I[i].b); 47 } 48 FFT(A,N,-1); 49 FFT(D,N,-1); 50 FFT(F,N,-1); 51 for(int i=0;i<n;i++){ 52 //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)\n",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); 53 C[i].a=D[i].a-A[i].a*2+F[i].a; 54 } 55 for(int i=0;i+m<=n;i++)if(C[i+m-1].a<eps)ans++; 56 printf("%d\n",ans); 57 for(int i=0;i+m<=n;i++)if(C[i+m-1].a<eps)printf("%d\n",i); 58 return 0; 59 } 60 void FFT(Complex *A,int n,int tp){ 61 for(int i=1,j=0,k;i<n-1;i++){ 62 k=N; 63 do{ 64 k>>=1; 65 j^=k; 66 }while(j<k); 67 if(i<j)swap(A[i],A[j]); 68 } 69 for(int k=2;k<=n;k<<=1){ 70 Complex wn(cos(-tp*2*pi/k),sin(-tp*2*pi/k)); 71 for(int i=0;i<n;i+=k){ 72 Complex w(1.0,0.0); 73 for(int j=0;j<(k>>1);j++,w*=wn){ 74 Complex a=A[i+j],b=w*A[i+j+(k>>1)]; 75 A[i+j]=a+b; 76 A[i+j+(k>>1)]=a-b; 77 } 78 } 79 } 80 if(tp<0)for(int i=0;i<n;i++)A[i].a/=n; 81 } 82 /* 83 把T反转，把S和T变成多项式a和b，令 84 c[i+m-1]=sum{(a[i+j]-b[m-j-1])^2*b[m-j-1]} 85 显然只有c是0的位置两串才会匹配，展开得 86 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} 87 前两项是卷积形式，第三项乘上一个伪单位元也是卷积形式，FFT加速即可 88 */
ps：其实对$I$跑的那个DFT完全没必要，手动逐项赋成1就行了，然后发现$F$卷完了$I$根本没有什么变化，所以直接对$F$做一下DFT再IDFT回去即可，把$I$写进去只是为了理解式子方便......
UPD：脑残了……那个ps是错的，忽视掉吧……
