[BZOJ4259][FFT]残缺的字符串

复习一下FFT
http://blog.csdn.net/u011542204/article/details/50708834
感觉这种字符串匹配的方法很棒啊
然而我的代码被卡常卡内存……

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cmath>
#define eps 1e-7
#define N 1100010

using namespace std;

const double pi=acos(-1);

struct E{
  double real,imag;
  E(double r=0,double i=0):real(r),imag(i){}
  friend E operator +(const E &a,const E &b){ return E(a.real+b.real,a.imag+b.imag); }
  friend E operator -(const E &a,const E &b){ return E(a.real-b.real,a.imag-b.imag); }
  friend E operator *(const E &a,const E &b){
    return E(a.real*b.real-a.imag*b.imag,a.real*b.imag+a.imag*b.real);
  }
  friend E operator /(const E &a,const double &b){ return E(a.real/b,a.imag/b); }
};

int n,m,L,M,p;
int rev[N],Ans[N];
char A[N],B[N];
E a[N],b[N],a1[N],b1[N],f[N];

inline void FFT(E *a,int r){
  for(int i=0;i<n;i++) if(rev[i]>i) swap(a[i],a[rev[i]]);
  for(int i=1;i<n;i<<=1){
    E wn(cos(pi/i),r*sin(pi/i));
    for(int j=0;j<n;j+=(i<<1)){
      E w(1,0);
      for(int k=0;k<i;k++,w=w*wn){
    E x=a[j+k],y=w*a[j+k+i];
    a[j+k]=x+y; a[j+k+i]=x-y;
      }
    }
  }
  if(r==-1) for(int i=0;i<n;i++) a[i]=a[i]/n;
}

int main(){
  scanf("%d %d",&n,&m); int in=n;
  scanf("%s%s",A,B);
  reverse(A,A+n);
  for(int i=0;i<n;i++) a[i]=(A[i]=='*')?0:(A[i]-'a'+1);
  for(int i=0;i<m;i++) b[i]=(B[i]=='*')?0:(B[i]-'a'+1);
  for(n=1,M=m*2;n<=M;n<<=1,L++);
  for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L-1);
  for(int i=0;i<m;i++) a1[i]=a[i]*a[i]*a[i],b1[i]=b[i];
  FFT(a1,1); FFT(b1,1);
  for(int i=0;i<n;i++) f[i]=a1[i]*b1[i];
  for(int i=0;i<n;i++) a1[i]=a[i]*a[i],b1[i]=b[i]*b[i];
  FFT(a1,1); FFT(b1,1);
  for(int i=0;i<n;i++) f[i]=f[i]-a1[i]*b1[i]-a1[i]*b1[i];
  for(int i=0;i<n;i++) a1[i]=a[i],b1[i]=b[i]*b[i]*b[i];
  FFT(a1,1); FFT(b1,1);
  for(int i=0;i<n;i++) f[i]=f[i]+a1[i]*b1[i];
  FFT(f,-1);
  for(int i=in-1;i<m;i++)
    if(int(0.1+f[i].real)==0) Ans[++p]=i;
  printf("%d\n",p);
  for(int i=1;i<=p;i++) printf("%d%c",Ans[i]-in+2,i==p?'\n':' ');
  return 0;
}