【快速傅里叶变换应用(可失配匹配问题)】【TJOI2017】DNA

本文探讨了在DNA序列中寻找近似匹配子串的方法,使用快速傅里叶变换优化匹配过程,适用于允许少量错配的情况。介绍了算法原理,包括字符串匹配、后缀数组和多项式乘法的应用。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

【题目描述】
给定两个DNA串A, B,如果说A的某个连续子串和B长度相同,且对应位置上的字符不同的数量小于3,那我们说A的这个子串和B近似匹配。
求A中有多少个连续子串和B近似匹配。两个A的子串当且仅当起始位置和/或结束为止不同时视作不同子串。

【输入】
T组数据 每组数据两行,表示A,B

【输出】
对于每组数据输出一行为答案。

【思路】

这道题是经典的字符串匹配问题,方法很多。如后缀数组,二分+hash(O(k∗n∗lognk*n*lognknlogn)(k为允许的失配次数)),以及快速傅里叶变换(O(n∗logn∗∑n*logn*\sumnlogn)∑\sum是字符集大小)等。
这里主要介绍快速傅里叶变换的做法。我们知道,对于一般的字符串匹配问题,我们定义c数组满足:
c[i]=∑j=0m−1(a[i+j]−b[j])2 c[i]=\sum_{j=0}^{m-1}(a[i+j]-b[j])^2 c[i]=j=0m1(a[i+j]b[j])2
然后翻转其中一个字符串以变为多项式乘法形式。
那么对于这道题,我们自然希望c数组能告诉我们失配或匹配次数以便统计答案。考虑枚举字符集,对于枚举的字符cc,我们定义a,b,c,d满足:
a[i]=[s1[i]==cc]a[i]=[s1[i]==cc]a[i]=[s1[i]==cc]
b[i]=[s2[i]==cc]b[i]=[s2[i]==cc]b[i]=[s2[i]==cc]
c[i]=∑j=0m−1a[i+j]∗b[j]c[i]=\sum_{j=0}^{m-1}a[i+j]*b[j]c[i]=j=0m1a[i+j]b[j]
d[i]=∑ccc[i]d[i]=\sum_{cc} c[i]d[i]=ccc[i]
那么我们发现,如果d数组中是枚举每个字符后c对应位置的和,那么一对字符对d数组对应位置有贡献,当且仅当这一对字符匹配,且贡献为1。所以我们只需要翻转其中一个字符串,求出d数组就行了。d数组中某一个位置对答案有贡献,当且仅当该位置的数x满足x>=m−3x>=m-3x>=m3

代码:

#include<bits/stdc++.h>
#define re register
using namespace std;
const int N=2.7e5+5;
int n,m;
inline int red()
{
    int data=0;int w=1; char ch=0;
    ch=getchar();
    while(ch!='-' && (ch<'0' || ch>'9')) ch=getchar();
    if(ch=='-') w=-1,ch=getchar();
    while(ch>='0' && ch<='9') data=(data<<3)+(data<<1)+ch-'0',ch=getchar();
    return data*w;
}
int lim=0,l=0,rev[N],t=0;
inline void pre(int mx)
{
	lim=1,l=0;
	while(lim<mx)lim<<=1,++l;rev[0]=0;
	for(int re i=0;i<lim;i++)
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
}
struct cp{
	double x,y;
	cp(double _x=0,double _y=0){x=_x,y=_y;}
	friend cp operator +(cp a,cp b){return cp(a.x+b.x,a.y+b.y);}
	friend cp operator -(cp a,cp b){return cp(a.x-b.x,a.y-b.y);}
	friend cp operator *(cp a,cp b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
};
double pi=3.141592653589793;
inline void fft(cp *f,int op)
{
	for(int re i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
	for(int re mid=1;mid<lim;mid<<=1)
	{
		cp tmp(cos(pi/mid),op*sin(pi/mid));
		for(int re j=0;j<lim;j+=(mid<<1))
		{
			cp w(1,0);
			for(int re k=0;k<mid;k++,w=w*tmp)
			{
				cp x=f[j+k],y=f[j+k+mid]*w;
				f[j+k]=x+y;f[j+k+mid]=x-y;
			}
		}
	}
	if(op==-1)for(int re i=0;i<lim;i++)f[i].x/=lim;
}
char s1[N],s2[N];cp a[N],b[N];
char cc[4]={'A','T','C','G'};
int d[N];char c;
int main()
{
	scanf("%d",&t);
	while(t--)
	{
		int ans=0;
		memset(d,0,sizeof(d));
		scanf("%s%s",s1,s2);
		n=strlen(s1),m=strlen(s2);
		reverse(s1,s1+n);pre(n+m-1);
		for(int re j=0;j<4;j++)
		{
			c=cc[j];memset(a,0,sizeof(a));memset(b,0,sizeof(b));
			for(int re i=0;i<n;i++)a[i].x=(s1[i]==c);
			for(int re i=0;i<m;i++)b[i].x=(s2[i]==c);
			fft(a,1);fft(b,1);
			for(int re i=0;i<lim;i++)a[i]=a[i]*b[i];fft(a,-1);
			for(int re i=m-1;i<n;i++)d[i]+=(int)(a[i].x+0.5);
		}
		for(int re i=m-1;i<n;i++)
			if(d[i]>=m-3)ans++;
		printf("%d\n",ans);
	}
}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值