51nod 1565 模糊搜索 fft

题意

有两个基因串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;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值