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\n",re);       
	return 0;  
}

　　

转载于:https://www.cnblogs.com/guangheli/p/11170862.html