【xsy1154】 DNA配对 FFT

题目大意:给你一个字符串$s$和字符串$w$,字符集为${A,T,C,G}$,你要在字符串$s$中选出一个与$w$长度相同的子串,使得这两个串的差异度最小。

两个字符$c1$,$c2$的差异度为给定的$c[c1][c2]$。

字符串长度$≤2*10^5$。

$FFT$套路题。

我们将串$w$翻转。

设$p[i]$为$s$中子串$s[i-|w|+1.......i]$与$w$的差异度。

显然$p[i]=\sum_{j=0}^{i} c[s[j]][w[i-j]]$。(此处的$w$是翻转后的)

显然的卷积形式。

五次$FFT$即可。

 

 1 #include<bits/stdc++.h>
 2 #define M (1<<19)
 3 #define PI acos(-1)
 4 #define INF 19890604
 5 using namespace std;
 6 
 7 struct cp{
 8     double i,r; cp(){i=r=0;}
 9     cp(double rr,double ii){i=ii;r=rr;}
10     friend cp operator +(cp a,cp b){return cp(a.r+b.r,a.i+b.i);}
11     friend cp operator -(cp a,cp b){return cp(a.r-b.r,a.i-b.i);}
12     friend cp operator *(cp a,cp b){return cp(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
13     friend cp operator /(cp a,double b){return cp(a.r/b,a.i/b);}
14     int num(){return (int)(i+0.499)/2;}
15 }a[M],b[M],c[M],d[M],ans[M];
16 
17 void change(cp a[],int len){
18     for(int i=0,j=0;i<len-1;i++){
19         if(i<j) swap(a[i],a[j]);
20         int k=len>>1;
21         while(j>=k) j-=k,k>>=1;
22         j+=k;
23     }
24 }
25 void FFT(cp a[],int len,int on){
26     change(a,len);
27     for(int h=2;h<=len;h<<=1){
28         cp wn=cp(cos(2*on*PI/h),sin(2*on*PI/h));
29         for(int j=0;j<len;j+=h){
30             cp w=cp(1,0);
31             for(int k=j;k<(j+(h>>1));k++){
32                 cp u=a[k],t=w*a[k+(h>>1)];
33                 a[k]=u+t; a[k+(h>>1)]=u-t;
34                 w=w*wn;
35             }
36         }
37     }
38     if(on==-1) for(int i=0;i<len;i++) a[i]=a[i]/len;
39 }
40 
41 int D[4][4]={0};
42 int get(char c){
43     if(c=='A') return 0;
44     if(c=='T') return 1;
45     if(c=='C') return 2;
46     if(c=='G') return 3;
47 }
48 
49 char s[M]={0},w[M]={0};
50 int n,m;
51 int main(){
52     scanf("%s%s",s,w); n=strlen(s); m=strlen(w);
53     for(int i=0;i<16;i++) scanf("%d",D[0]+i);
54     for(int i=0;i<n;i++) s[i]=get(s[i]);
55     for(int i=0;i<m;i++) w[i]=get(w[i]);
56     reverse(w,w+m);
57     for(int i=0;i<n;i++){
58         if(s[i]==0) a[i].i=1;
59         if(s[i]==1) b[i].i=1;
60         if(s[i]==2) c[i].i=1;
61         if(s[i]==3) d[i].i=1;
62     }
63     for(int i=0;i<m;i++){
64         a[i].r=D[0][w[i]];
65         b[i].r=D[1][w[i]];
66         c[i].r=D[2][w[i]];
67         d[i].r=D[3][w[i]];
68     }
69     int len=1; while(len<n+m) len<<=1;
70     FFT(a,len,1); FFT(b,len,1); FFT(c,len,1); FFT(d,len,1);
71     for(int i=0;i<len;i++){
72         ans[i]=a[i]*a[i]+b[i]*b[i]+c[i]*c[i]+d[i]*d[i];
73     }
74     FFT(ans,len,-1);
75     int minn=INF;
76     for(int i=m-1;i<n;i++) 
77     minn=min(minn,ans[i].num());
78     cout<<minn<<endl;
79 }

 

posted @ 2018-10-17 20:49  AlphaInf  阅读(183)  评论(0编辑  收藏  举报