CodeForces 528 D.Fuzzy Search（FFT）

最新推荐文章于 2021-10-13 21:54:39 发布

v5zsq

最新推荐文章于 2021-10-13 21:54:39 发布

阅读量560

点赞数

分类专栏： Code Forces FFT

本文链接：https://blog.csdn.net/v5zsq/article/details/77261420

版权

Code Forces 同时被 2 个专栏收录

606 篇文章 6 订阅

订阅专栏

FFT

25 篇文章 0 订阅

订阅专栏

Description
给出一个长度为 $n$ 的串 $S$ 和一个长度为 $m$ 的串 $T$ ，两个串均只由 $'A','G','C','T'$ 组成，定义 $T$ 和 $S$ 在 $i$ 位置匹配为对于任意 $1\leq j\leq m$ ，存在 $p(1\leq p\leq n)$ 使得 $|(i+j-1)-p|\leq k$ 且 $S[p]=T[j]$ ，问 $T$ 可以在多少位置和 $S$ 匹配
Input
第一行三个整数 $n,m,k$ 分别表示两个串的串长和匹配位置差上限，之后输入两个串 $S$ 和 $T$
$(1\leq m\leq n\leq 2\times 10^{5},0\leq k\leq 2\times 10^{5})$
Output
输出 $T$ 和 $S$ 可以匹配的位置数
Sample Input
10 4 1
AGCAATTCAT
ACAT
Sample Output
3
Solution

对四种字符单独去计算固定起点后的匹配长度，如果四个匹配长度加一起是m说明可以匹配
对某个字符，用一个 $01$ 序列 $a$ 表示 $S$ 串这个字符的存在情况， $a[i]=1$ 表示 $S$ 串的 $[i-k,i+k]$ 这个段中出现过该字符，用前缀和优化可以 $O(n)$ 得到 $a$ 序列，再用一个 $01$ 序列 $b$ 表示 $T$ 串的反串中该字符的存在情况， $b[i]=1$ 表示 $T[m-i]$ 是该字符，那么 $T$ 从 $S$ 的第 $i$ 个位置开始往前匹配的该种字符数量为 $c[i]=\sum\limits_{i=0}^{i}a[j]\times b[i-j]$ ，用 $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+5
const 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;
}