CodeForces 528D Fuzzy Search 多项式 FFT

时间：2022-02-01 04:15:49

原文链接http://www.cnblogs.com/zhouzhendong/p/8782849.html

题目传送门 - CodeForces 528D

题意

　　给你两个串$A,B(|A|\geq|B|)$，以及一个$k$。

　　其中$A_i$与$B_j$匹配的条件是$A_{i-k\dots i+k}$中至少有一个与$B_j$相同。

　　问$B$能在$A$中匹配多少次。

　　字符集：$\{'A','C','G','T'\}$。

　　$|B|\leq|A|\leq 2\times 10^5,k\leq 2\times 10^5$

题解

　　由于我太弱了做这题又百度了真难受。

　　首先我们先学习一个用$FFT$快速匹配字符串的办法$\longrightarrow$链接。

　　首先注意到字符集很小。

　　我一开始的想法：

令

$$f[c]_i=[A串的第i位可以与字符c匹配]$$

$$g[c]_i=[B串的第i位可以与字符c匹配]$$

于是很容易写出卷积式子：

$$h_i=\sum_{j=0}^{i}\prod_{c\in\{'A','C','G','T'\}}(f[c]_j-g[c]_{i-j})^2$$

如果$h_i=0$则从$A$的第$i-|B|+1$个位置开始可以匹配。

但是马上发现这个式子的锅太重了。

一个是你展开之后……（你先把他展开吧QAQ）

第二个是他会因为他硕大的常数而最终GG。

　　然后就求助度娘了。

　　发现4个字符可以分开贡献。主要是因为$B$串只能死板匹配。

　　对于一个字符$c$，

　　如果$A$串的第$i$个位置可以放$c$，那么$f_i=1$，否则$f_i=0$。

　　如果$B$串的第$i$个位置可以放$c$，那么$g_i=1$，否则$g_i=0$。

　　然后卷积：

　　$$h_i=\sum_{j=0}^i f_jg_{i-j}$$

　　然后把$4$个字符的解累加起来。

　　$$res_i=\sum_{c\in\{'A','C','G','T'\}}h[c]_i$$。

　　显然，如果$res_i=|B|$，那么从$A$串的第$i-|B|+1$位开始可以匹配。

　　然后时间复杂度还是$O((|A|+|B|)\log(|A|+|B|))$。常数小了很多，也不用自己展开了。

代码

#include <bits/stdc++.h>

using namespace std;

const int N=1<<19;

double PI=acos(-1.0);

int n,L,R[N];

struct C{

	double r,i;

	C(){}

	C(double a,double b){r=a,i=b;}

	C operator + (C x){return C(r+x.r,i+x.i);}

	C operator - (C x){return C(r-x.r,i-x.i);}

	C operator * (C x){return C(r*x.r-i*x.i,r*x.i+i*x.r);}

}w[N],X[N],Y[N],Z[N];

double x[N],y[N],z[N];

char str1[N],str2[N];

int cti[300];

int k,a[N],b[N],A,B,res[N];

void FFT(C a[]){

	for (int i=0;i<n;i++)

		if (i>R[i])

			swap(a[i],a[R[i]]);

	for (int t=n>>1,d=1;d<n;d<<=1,t>>=1)

		for (int i=0;i<n;i+=(d<<1))

			for (int j=0;j<d;j++){

				C tmp=w[t*j]*a[i+j+d];

				a[i+j+d]=a[i+j]-tmp;

				a[i+j]=a[i+j]+tmp;

			}

}

void FFT_times(double x[],double y[],double z[]){

	for (int i=0;i<n;i++)

		X[i]=C(x[i],0),Y[i]=C(y[i],0);

	FFT(X),FFT(Y);

	for (int i=0;i<n;i++)

		Z[i]=X[i]*Y[i],w[i].i*=-1.0;

	FFT(Z);

	for (int i=0;i<n;i++)

		z[i]=Z[i].r/n,w[i].i*=-1.0;

}

int main(){

	cti['A']=0,cti['C']=1,cti['G']=2,cti['T']=3;

	scanf("%d%d%d%s%s",&A,&B,&k,str1,str2);

	for (int i=0;i<A;i++)

		a[i]=cti[str1[i]];

	for (int i=0;i<B;i++)

		b[i]=cti[str2[i]];

	for (int i=0;i<B/2;i++)

		swap(b[i],b[B-i-1]);

	for (n=1,L=0;n<A+B;n<<=1,L++);

	for (int i=0;i<n;i++){

		R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));

		w[i]=C(cos(2*i*PI/n),sin(2*i*PI/n));

	}

	memset(res,0,sizeof res);

	for (int i=0;i<4;i++){

		for (int j=0;j<n;j++)

			x[j]=y[j]=0;

		int cnt=0;

		for (int j=0;j<k;j++)

			cnt+=a[j]==i;

		for (int j=0;j<A;j++){

			if (j+k<A)

				cnt+=a[j+k]==i;

			if (j-k-1>=0)

				cnt-=a[j-k-1]==i;

			x[j]=cnt>0?1:0;

		}

		for (int j=0;j<B;j++)

			y[j]=b[j]==i?1:0;

		FFT_times(x,y,z);

		for (int j=B-1;j<A;j++)

			res[j]+=(int)(z[j]+0.5);

	}

	int ans=0;

	for (int i=B-1;i<A;i++)

		if (res[i]==B)

			ans++;

	printf("%d",ans);

	return 0;

}

　　

相关文章