51nod 1565 模糊搜索 FFT

这。。。好强啊(QwQ)


(FFT)?

(5)次

(QwQ)

题解:

对四种字符分开考虑:我们设(a[char][i])表示在第一个串(s)中,对于(char in {'A','C','G','T'})来说(i)位置是否能模糊匹配,换言之,若(s[i]==char),则(a[char][j]=1,jin [i-k,i+k])
而对于第二个串(t)不做特殊处理,直接(b[char][i]=[t[i]==char])
我们在不加优化时,计算答案是(O(n^2))

for(R i=0;i<=lens-lent;++i) for(R j=1;j<=lent;++j) if(a[char][i+j]&&b[char][j]) ++c[i];

最后需要统计(tmp=sum [t[i]==char]),若(c[i]==tmp)表示对(char)匹配成功。
所以我们要做(4)遍,对于每一个字符都做一遍。
考虑优化(O(n^2))的过程:我们发现把(t)(b)倒过来的话,上面的枚举相当于是一个卷积。
于是我们可以(FFT)

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#define R register int
using namespace std;
namespace Luitaryi {
template<class I> inline I g(I& x) { x=0;
  register I f(1); register char ch; while(!isdigit(ch=getchar())) f=ch=='-'?-1:f;
  do x=x*10+(ch^48); while(isdigit(ch=getchar())); return x*=f;
} const int N=524292; const double PI=acos(-1.0);
int n,m,anss,k,K=1,l,l1,l2,p[N],a[4][N],b[4][N],c[N];
char s[N];
bool ans[N];
struct complex { double x,y; complex() {}
  complex(double _x,double _y) {x=_x,y=_y;}
  complex operator + (const complex& that) {return complex(x+that.x,y+that.y);}
  complex operator - (const complex& that) {return complex(x-that.x,y-that.y);}
  complex operator * (const complex& that) {return complex(x*that.x-y*that.y,x*that.y+y*that.x);}
}f[N],h[N];
inline int cvt(const char& c) {
  if(c=='A') return 0; if(c=='C') return 1; 
  if(c=='G') return 2; if(c=='T') return 3;
}
inline void fft(complex* a,short op) {
  for(R i=0;i<K;++i) if(i<p[i]) swap(a[i],a[p[i]]);
  for(R l=1;l<K;l<<=1) { register complex w1(cos(PI/l),op*sin(PI/l));
    for(R len=l<<1,i=0;i<K;i+=len) { register complex wn(1,0);
      for(R j=0;j<l;++j,wn=wn*w1) { register complex x=a[i+j],y=a[i+j+l]*wn;
        a[i+j]=x+y,a[i+j+l]=x-y;
      }
    }
  } if(op==-1) for(R i=0;i<=n;++i) a[i].x=1.0*fabs(a[i].x)/K;
}
inline void main() { 
  g(l1),g(l2),g(k); scanf("%s",s); 
  for(R i=0;i<l1;++i) {
    R tmp=cvt(s[i]); a[tmp][i]=1;
    for(R j=i+1,lim=min(l1-1,i+k);j<=lim&&s[i]!=s[j];++j) a[tmp][j]=1;
    for(R j=i-1,lim=max(0,i-k);j>=lim&&!a[tmp][j];--j) a[tmp][j]=1;
  } scanf("%s",s); for(R i=0;i<l2;++i) b[cvt(s[i])][l2-i-1]=1; 
  n=l1+l2; while(K<=n) K<<=1,++l;
  for(R i=0;i<K;++i) p[i]=(p[i>>1]>>1)|((i&1)<<(l-1)); 
  for(R i=l2-1;i<=n-2;++i) ans[i]=1;
  for(R i=0;i<4;++i) {
    memset(f,0,sizeof(complex)*(K+2)),memset(h,0,sizeof(complex)*(K+2));
    for(R j=0;j<l1;++j) f[j]=complex(a[i][j],0.0);
    for(R j=0;j<l2;++j) h[j]=complex(b[i][j],0.0);
    fft(f,1),fft(h,1); for(R i=0;i<=K;++i) f[i]=f[i]*h[i];
    fft(f,-1); for(R i=0;i<=n-2;++i) c[i]=(int)(f[i].x+0.5);
    R tmp=0; for(R j=0;j<l2;++j) tmp+=b[i][j]; 
    for(R j=l2-1;j<=n-2;++j) ans[j]&=(tmp==c[j]);
  } for(R i=l2-1;i<=n-2;++i) anss+=ans[i]; printf("%d
",anss);
}
} signed main() {Luitaryi::main(); return 0;}

2019.08.12
88