zoukankan      html  css  js  c++  java
  • CodeForces 528D Fuzzy Search 多项式 FFT

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

    题目传送门 - CodeForces 528D

    题意

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

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

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

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

      $|B|leq|A|leq 2 imes 10^5,kleq 2 imes 10^5$

    题解

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

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

      首先注意到字符集很小。

      我一开始的想法:

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

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

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

    $$h_i=sum_{j=0}^{i}prod_{cin{'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_{cin{'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;
    }
    

      

  • 相关阅读:
    SQL后台分页三种方案和分析
    SQL分页查询语句
    SQL利用临时表实现动态列、动态添加列
    查询sybase DB中占用空间最多的前20张表
    敏捷软件开发之TDD(一)
    敏捷软件开发之开篇
    Sql Server 2012启动存储过程
    改变VS2013的菜单栏字母为小写
    Sql Server获得每个表的行数
    Sql Server trace flags
  • 原文地址:https://www.cnblogs.com/zhouzhendong/p/CF528D.html
Copyright © 2011-2022 走看看