CodeForces 528 D.Fuzzy Search(FFT)

来源:互联网 发布:如何提高淘宝店信誉 编辑:程序博客网 时间:2024/06/07 05:56

Description
给出一个长度为n的串S和一个长度为m的串T,两个串均只由A,G,C,T组成,定义TSi位置匹配为对于任意1jm,存在p(1pn)使得|(i+j1)p|kS[p]=T[j],问T可以在多少位置和S匹配
Input
第一行三个整数n,m,k分别表示两个串的串长和匹配位置差上限,之后输入两个串ST
(1mn2×105,0k2×105)
Output
输出TS可以匹配的位置数
Sample Input
10 4 1
AGCAATTCAT
ACAT
Sample Output
3
Solution

对四种字符单独去计算固定起点后的匹配长度,如果四个匹配长度加一起是m说明可以匹配
对某个字符,用一个01序列a表示S串这个字符的存在情况,a[i]=1表示S串的[ik,i+k]这个段中出现过该字符,用前缀和优化可以O(n)得到a序列,再用一个01序列b表示T串的反串中该字符的存在情况,b[i]=1表示T[mi]是该字符,那么TS的第i个位置开始往前匹配的该种字符数量为c[i]=i=0ia[j]×b[ij],用FFT加速该卷积即可,为减少DFT次数,可以先把处理四种字符时的a,b序列DFT后结果乘积先加一起,最后做一次IDFT即可

Code

#include<cstdio>#include<iostream>#include<cstring>#include<algorithm>#include<cmath>#include<vector>#include<queue>#include<map>#include<set>#include<ctime>using namespace std;#define maxn 200001#define maxfft 524288+5const double pi=acos(-1.0);struct cp{    double a,b;    cp operator +(const cp &o)const {return (cp){a+o.a,b+o.b};}    cp operator -(const cp &o)const {return (cp){a-o.a,b-o.b};}    cp operator *(const cp &o)const {return (cp){a*o.a-b*o.b,b*o.a+a*o.b};}    cp operator *(const double &o)const {return (cp){a*o,b*o};}    cp operator !() const{return (cp){a,-b};}}w[maxfft];int pos[maxfft];void fft_init(int len){    int j=0;    while((1<<j)<len)j++;    j--;    for(int i=0;i<len;i++)        pos[i]=pos[i>>1]>>1|((i&1)<<j);}void fft(cp *x,int len,int sta){    for(int i=0;i<len;i++)        if(i<pos[i])swap(x[i],x[pos[i]]);    w[0]=(cp){1,0};    for(unsigned i=2;i<=len;i<<=1)    {        cp g=(cp){cos(2*pi/i),sin(2*pi/i)*sta};        for(int j=i>>1;j>=0;j-=2)w[j]=w[j>>1];        for(int j=1;j<i>>1;j+=2)w[j]=w[j-1]*g;        for(int j=0;j<len;j+=i)        {            cp *a=x+j,*b=a+(i>>1);            for(int l=0;l<i>>1;l++)            {                cp o=b[l]*w[l];                b[l]=a[l]-o;                a[l]=a[l]+o;            }        }    }    if(sta==-1)for(int i=0;i<len;i++)x[i].a/=len,x[i].b/=len;}cp x[maxfft],y[maxfft],z[maxfft];int n,m,k,sum[maxn];char S[maxn],T[maxn];void Solve(char s,int len){    for(int i=0;i<len;i++)x[i]=y[i]=(cp){0,0};    for(int i=0;i<n;i++)sum[i+1]=(S[i]==s);    for(int i=1;i<=n;i++)sum[i]+=sum[i-1];    for(int i=0;i<n;i++)    {        int L=max(0,i-k),R=min(n-1,i+k)+1;        if(sum[R]-sum[L])x[i]=(cp){1,0};    }    for(int i=0;i<m;i++)        if(T[i]==s)y[i]=(cp){1,0};    fft(x,len,1),fft(y,len,1);    for(int i=0;i<len;i++)z[i]=z[i]+x[i]*y[i];}int main(){    while(~scanf("%d%d%d",&n,&m,&k))    {        scanf("%s%s",S,T);        reverse(T,T+m);        int len=1;        while(len<n+m)len<<=1;        fft_init(len);        for(int i=0;i<len;i++)z[i]=(cp){0,0};        Solve('A',len),Solve('G',len),Solve('C',len),Solve('T',len);        fft(z,len,-1);        int ans=0;        for(int i=m-1;i<n;i++)            if((int)(z[i].a+0.5)>=m)ans++;        printf("%d\n",ans);     }    return 0;}
原创粉丝点击