【题目描述】
给定两个DNA串A, B,如果说A的某个连续子串和B长度相同,且对应位置上的字符不同的数量小于3,那我们说A的这个子串和B近似匹配。
求A中有多少个连续子串和B近似匹配。两个A的子串当且仅当起始位置和/或结束为止不同时视作不同子串。
【输入】
T组数据 每组数据两行,表示A,B
【输出】
对于每组数据输出一行为答案。
【思路】
这道题是经典的字符串匹配问题,方法很多。如后缀数组,二分+hash(O(k∗n∗lognk*n*lognk∗n∗logn)(k为允许的失配次数)),以及快速傅里叶变换(O(n∗logn∗∑n*logn*\sumn∗logn∗∑)∑\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=0∑m−1(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=0m−1a[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>=m−3。
代码:
#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);
}
}