zoukankan      html  css  js  c++  java
  • CF528D Fuzzy Search 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中出现的次数。

    Input 单组测试数据。 第一行有三个整数 |S|,|T|,k (1≤|T|≤|S|≤200000, 0≤k≤200000),表示S的长度,T的长度,门限值。 第二行给出基因串S。 第三行给出基因串T。 两个串都只包含大写字母'A', 'T', 'G' 和'C'。 Output 输出一个整数,表示T在S中出现的次数。

    题解:如果题中给定的 $k$ 等于 $0$ 的话就是一道 $SAM/KMP/AC$自动机模板题了
     
    有了 $k$ 的限制就让问题的难度提高了一个层次
     
    发现字符集大小只有 $4$,我们枚举每种字符依次考虑
     
    对于字符 $h$
     
    令 $a_{i}$ 表示 $S$ 中第 $i$ 位是否能匹配上 $h$
     
    令 $b_{i}$ 表示 $T$ 中第 $i$ 位是否是 $h$ 
     
    令 $C_{j}=sum_{i=0}^{|T|-1}b_{i}a_{j+i}$,若 $C_{j}=T$ 中 $h$ 数量则 $h$ 是满足要求的
     
    然而上式 $C_{j}=sum_{i=0}^{|T|-1}b_{i}a_{j+i}$ 的下标和是不固定的,翻转 $b$
     
    得 $b_{i}Rightarrow b_{|T|-1-i}$ 则 $C_{j}=sum_{i=0}^{|T|-1}b_{|T|-1-i}a_{j+i}$
     
    下标和为 $|T|-1+j$ ,是固定的
     
    将上述过程作用于 $4$ 个字符,若 $sum C_{j}=|T|$ 则 $j$ 位置是合法的
    #include<bits/stdc++.h>
    #define setIO(s) freopen(s".in","r",stdin) 
    #define ll long long 
    #define maxn 800002  
    using namespace std;
    const double pi=acos(-1.0);  
    inline int get(char c)
    {
    	if(c=='A') return 1; 
    	if(c=='T') return 2; 
    	if(c=='G') return 3; 
    	if(c=='C') return 4; 
    } 
    struct cpx 
    {
    	double x,y; 
    	cpx(double a=0,double b=0) {x=a,y=b; } 
    	cpx operator+(const cpx b) { return cpx(x+b.x,y+b.y); } 
    	cpx operator-(const cpx b) { return cpx(x-b.x,y-b.y); } 
    	cpx operator*(const cpx b) { return cpx(x*b.x-y*b.y,x*b.y+y*b.x); }
    }A[maxn],B[maxn];      
    void FFT(cpx *a,int n,int flag) 
    {
    	for(int i=0,k=0;i<n;++i) 
    	{
    		if(i>k) swap(a[i],a[k]); 
    		for(int j=(n>>1);(k^=j)<j;j>>=1); 
    	}
        for(int mid=1;mid<n;mid<<=1) 
        { 
        	cpx wn(cos(pi/mid), flag*sin(pi/mid)),x,y; 
        	for(int i=0;i<n;i+=(mid<<1)) 
        	{
        		cpx w(1,0); 
        		for(int j=0;j<mid;++j) 
        		{
        			x=a[i+j], y=w*a[i+j+mid];  
        			a[i+j]=x+y,a[i+j+mid]=x-y;           
        			w=w*wn; 
        		}
        	}
        }
        if(flag==-1) for(int i=0;i<n;++i) a[i].x/=(double)n;     
    } 
    int len_s,len_t,k; 
    int a[10][maxn],b[10][maxn],answer[10][maxn];    
    int S[maxn],T[maxn];   
    char srr[maxn],trr[maxn];     
    inline void Initialize(int h)
    {
    	int cnt=0; 
    	for(int i=0;i<len_s;++i) 
    	{   
    		cnt+=(S[i]==h);   
    		if(i-k-1>=0) cnt-=(S[i-k-1]==h); 
    		if(cnt) a[h][i]=1;  
    	} 
    	cnt=0; 
    	for(int i=len_s-1;i>=0;--i) 
    	{
    		cnt+=(S[i]==h); 
    		if(i+k+1<len_s) cnt-=(S[i+k+1]==h); 
    		if(cnt) a[h][i]=1;  
    	}
    }
    inline void solve(int h,int len)
    { 
    	for(int i=0;i<len;++i) A[i].x=A[i].y=B[i].x=B[i].y=0;   
    	for(int i=0;i<len;++i) A[i].x=a[h][i];    
    	for(int i=0;i<len;++i) B[i].x=b[h][len_t-1-i];
    	FFT(A,len,1),FFT(B,len,1); 
    	for(int i=0;i<len;++i) A[i]=A[i]*B[i]; 
    	FFT(A,len,-1);  
    	for(int i=0;i<len_s;++i) answer[h][i]=(ll)(A[len_t-1+i].x+0.5);      
    }
    int main()
    {
    	//  setIO("input");   
    	scanf("%d%d%d",&len_s,&len_t,&k);     
    	scanf("%s%s",srr,trr);   
    	for(int i=0;i<len_s;++i) S[i]=get(srr[i]); 
    	for(int i=0;i<len_t;++i) T[i]=get(trr[i]);
    	for(int i=0;i<len_t;++i) b[T[i]][i]=1;        
    	for(int i=1;i<=4;++i) Initialize(i);            
    	int len; 
    	for(len=1;len<=(len_s+len_t);len<<=1);    
    	for(int i=1;i<=4;++i) solve(i,len);   
    	int re=0; 
    	for(int i=0;i<len_s;++i) if(answer[1][i]+answer[2][i]+answer[3][i]+answer[4][i]>=len_t) ++re; 
    	printf("%d
    ",re);       
    	return 0;  
    }
    

      

  • 相关阅读:
    2019总结及2020计划
    蓝牙BLE连接与操作
    Android蓝牙操作
    Cannot use JSX unless the '--jsx' flag is provided.
    PyQt打包可执行文件
    Springboot项目报错【java.base/jdk.internal.loader.ClassLoaders$AppClassLoader cannot be cast to java.base/java.net.URLClassLoader】
    typescript枚举字符串型不能使用函数问题
    beautifulsoap常用取节点方法
    numpy常用矩阵操作
    MYSQL 碎片查询
  • 原文地址:https://www.cnblogs.com/guangheli/p/11170862.html
Copyright © 2011-2022 走看看