• ●codeforces 528D Fuzzy Search


    题链:

    http://codeforces.com/problemset/problem/528/D

    题解:

    FFT

    先解释一下题意:

    给出两个字符串(只含'A','T','C','G'四种字符),一个为文本串T(长度为n),一个为模式串S(长度为m)。

    要用模式串去匹配文本串。

    同时给出一个正整数k,表示允许的匹配误差范围为k,即:

    如果对于T[i]和S[j],只要在T[i-k~i+k]范围中存在一个字符与S[j]相同,那么T[i]和S[j]就匹配。

    求出T中有多少个位置i满足从该位置开始的长度为m的子串T[i~i+m-1]可以和S串匹配。


    由于字符集很小,我们可以对每种字符处理。

    假设现在只考虑'A'字符,

    我们标记出T串中有哪些位置可以和A字符匹配,得到数组f,(1表示该位置匹配,0表示无法匹配)

    同时也用0,1标记出S串中的A字符,得到数组g。

    然后不难发现,如果让T串的第i位和S串的第j位匹配,那么匹配是否成功就可以用$f[i]*g[j]$表示

    所以如果要让S串和T的第$l$位开始匹配,我们可以得到匹配的贡献$D_A(l)$:

    $$D_A(l)=sum_{k=0}^{m-1}f(l+k)g(k)$$

    这个式子可以写成卷积的形式,即只需要把S串翻转一下,就可以得到:

    $$D_A(l)=D_A'(l+m-1)=sum_{k=0}^{m-1}f(l+k)g(m-1-k)$$

    然后就可以用FFT求出所有的$D_A$。

    同理可以的到$D_T,D_C,D_G$。

    那么对于T的l位置开始的长度为m的子串是否于S串匹配,就只需要判断$D_A(l)+D_T(l)+D_C(l)+D_G(l)$是否等于m即可。

    代码:

    #include<bits/stdc++.h>
    #define MAXN 524289
    #define INF 0x3f3f3f3f
    using namespace std;
    const double Pi=acos(-1);
    struct Complex{
    	double real,image;
    	Complex(double _real=0,double _image=0):real(_real),image(_image){}
    	Complex operator - () const{return Complex(-real,-image);}
    	friend Complex operator + (const Complex &A,const Complex &B){return Complex(A.real+B.real,A.image+B.image);}
    	friend Complex operator - (const Complex &A,const Complex &B){return A+(-B);}
    	friend Complex operator * (const Complex &A,const Complex &B){return Complex(A.real*B.real-A.image*B.image,A.image*B.real+A.real*B.image);}
    }null(0,0);
    int cnt[MAXN],T[MAXN],S[MAXN],order[MAXN];
    Complex A[MAXN],B[MAXN];
    int idx(char ch){
    	switch(ch){
    		case 'A':return 1;	case 'T':return 2;
    		case 'C':return 3;	case 'G':return 4;
    	}
    	return 0;
    }
    void getstring(int *s,int len){
    	static char ch;
    	for(int i=0;i<len;i++)
    		scanf(" %c",&ch),s[i]=idx(ch);
    }
    void mark(int *s,int len,int id,int lim,Complex *Y,int n){
    	static int last; last=-INF;
    	for(int i=0;i<n;i++) Y[i]=null;
    	for(int i=0;i<len;i++){
    		if(s[i]==id) last=i;
    		if(i-last<=lim) Y[i].real=1;
    	} last=INF;
    	for(int i=len-1;i>-1;i--){
    		if(s[i]==id) last=i;
    		if(last-i<=lim) Y[i].real=1;
    	}
    }
    void FFT(Complex *Y,int n,int sign){
    	for(int i=0;i<n;i++) if(i<order[i]) swap(Y[i],Y[order[i]]);
    	for(int d=2;d<=n;d<<=1){
    		Complex dw(cos(2*Pi/d),sin(sign*2*Pi/d)),w,tmp;
    		for(int i=0;w=Complex(1,0),i<n;i+=d)
    			for(int k=i;k<i+d/2;w=w*dw,k++)
    				tmp=w*Y[k+d/2],Y[k+d/2]=Y[k]-tmp,Y[k]=Y[k]+tmp;
    	}
    }
    int main(){
    	int n,m,k,N,len,ans=0;
    	scanf("%d%d%d",&n,&m,&k);
    	getstring(T,n);
    	getstring(S,m);
    	reverse(S,S+m);
    	for(N=1,len=0;N<n+m-1;N<<=1) len++;
    	for(int i=1;i<N;i++) order[i]=(order[i>>1]>>1)|((i&1)<<(len-1));
    	for(int id=1;id<=4;id++){
    		mark(T,n,id,k,A,N); 
    		mark(S,m,id,0,B,N);
    		FFT(A,N,1); FFT(B,N,1);
    		for(int i=0;i<N;i++) A[i]=A[i]*B[i];
    		FFT(A,N,-1);
    		for(int l=0;l<n;l++) cnt[l]+=(int)((A[l+m-1].real+0.5)/N);
    	}
    	for(int l=0;l<=n;l++) if(cnt[l]==m) ans++;
    	printf("%d
    ",ans);
    	return 0;
    }
    

      

  • 相关阅读:
    【对拍√】
    hdu5791 TWO
    luogu P1220 关路灯
    【NOI2001】食物链
    【HAOI2016】食物链
    luogu P1006 传纸条
    可持久化平衡树
    可持久化并查集
    线段树合并(【POI2011】ROT-Tree Rotations)
    可持久化数组
  • 原文地址:https://www.cnblogs.com/zj75211/p/8341224.html
Copyright © 2020-2023  润新知