第一道自己写的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$写进去只是为了理解式子方便......
原文地址：http://www.cnblogs.com/hzoier/p/6351822.html