51nod 1565 模糊搜索 fft

来源:互联网 发布:台达触摸屏编程视频 编辑:程序博客网 时间:2024/06/03 19:31

题意

有两个基因串S和T,他们只包含AGCT四种字符。现在你要找出T在S中出现了几次。
有一个门限值k≥0。T在S的第i(1≤i≤|S|-|T|+1)个位置中出现的条件如下:把T的开头和S的第i个字符对齐,然后T中的每一个字符能够在S中找到一样的,且位置偏差不超过k的,那么就认为T在S的第i个位置中出现。也就是说对于所有的 j (1≤j≤|T|),存在一个 p (1≤p≤|S|),使得|(i+j-1)-p|≤k 和[p]=T[j]都成立。
例如,根据这样的定义”ACAT”出现在”AGCAATTCAT”的第2,3和6的位置。
这里写图片描述
如果k=0,那么这个就是经典的字符串匹配问题。
现在给定门限和两个基因串S,T,求出T在S中出现的次数。
1≤|T|≤|S|≤200000, 0≤k≤200000

分析

对于字符串匹配次数问题有经典的fft做法。
这题的话多了一个门限的限制,那么我们可以分四次来做。假设现在做A,若第一个串的第i位可以匹配A则x[i]=1否则x[i]=0。若第二个串的第i位为A则y[i]=1否则y[i]=0。然后对两个数组做fft,若卷积后第i位上的数字等于第二个串中A出现的次数,则表示i这个位置的所有A都可以匹配掉,反之则不行。

代码

#include<iostream>#include<cstdio>#include<cstdlib>#include<cstring>#include<algorithm>#include<complex>#include<cmath>#define pi acos(-1)using namespace std;typedef complex<double> com;const int N=200005;int n,m,k,a[N],b[N],num[5],L,rev[N*3],p[26];bool vis[N];com x[N*3],y[N*3];char ch[N];void fft(com *a,int f){    for (int i=0;i<L;i++) if (i<rev[i]) swap(a[i],a[rev[i]]);    for (int i=1;i<L;i<<=1)    {        com wn(cos(pi/i),f*sin(pi/i));        for (int j=0;j<L;j+=i<<1)        {            com w(1,0);            for (int k=0;k<i;k++)            {                com u=a[j+k],v=w*a[j+k+i];                a[j+k]=u+v;a[j+k+i]=u-v;                w*=wn;            }        }    }    if (f==-1) for (int i=0;i<L;i++) a[i]/=L;}int main(){    scanf("%d%d%d",&n,&m,&k);    scanf("%s",ch);    p[0]=1;p[2]=2;p[6]=4;p[19]=3;    for (int i=0;i<n;i++) a[i]=p[ch[i]-'A'];    scanf("%s",ch);    for (int i=0;i<m;i++) b[m-i-1]=p[ch[i]-'A'],num[b[m-i-1]]++;    for (int i=m-1;i<n;i++) vis[i]=1;    int lg=0;    for (L=1;L<=n;L<<=1,lg++);    for (int i=0;i<L;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));    for (int i=1;i<=4;i++)    {        int mx=-k-1;        for (int j=0;j<L;j++) x[j]=0;        for (int j=0;j<n;j++)        {            mx=a[j]==i?j:mx;            if (j-mx<=k) x[j]=1;        }        mx=n+k;        for (int j=n-1;j>=0;j--)        {            mx=a[j]==i?j:mx;            if (mx-j<=k) x[j]=1;        }        for (int j=0;j<m;j++) y[j]=b[j]==i?1:0;        for (int j=m;j<L;j++) y[j]=0;        fft(x,1);fft(y,1);        for (int j=0;j<L;j++) x[j]*=y[j];        fft(x,-1);        for (int j=m-1;j<n;j++)            if ((int)(x[j].real()+0.1)<num[i]) vis[j]=0;    }    int ans=0;    for (int i=m-1;i<n;i++) ans+=vis[i];    printf("%d",ans);    return 0;}
原创粉丝点击