[BZOJ4892][TJOI2017]DNA(后缀数组)
题目描述
加里敦大学的生物研究所,发现了决定人喜不喜欢吃藕的基因序列S,有这个序列的碱基序列就会表现出喜欢吃藕的性状,但是研究人员发现对碱基序列S,任意修改其中不超过3个碱基,依然能够表现出吃藕的性状。现在研究人员想知道这个基因在DNA链S0上的位置。所以你需要统计在一个表现出吃藕性状的人的DNA序列S0上,有多少个连续子串可能是该基因,即有多少个S0的连续子串修改小于等于三个字母能够变成S。
输入输出格式
输入格式:
第一行有一个数T,表示有几组数据 每组数据第一行一个长度不超过10^5的碱基序列S0
每组数据第二行一个长度不超过10^5的吃藕基因序列S
输出格式:
共T行,第i行表示第i组数据中,在S0中有多少个与S等长的连续子串可能是表现吃藕性状的碱基序列
输入输出样例
说明
对于20%的数据,S0,S的长度不超过10^4
对于20%的数据,S0,S的长度不超过10^5,0<T<=10
两个串连起来中间插个特殊字符然后后缀数组求四次LCP即可。
启示:永远不要低估SA模板的默写难度。
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #define rep(i,l,r) for (int i=l; i<=r; i++) 5 typedef long long ll; 6 using namespace std; 7 8 const int N=200100; 9 int n,m,T,ans,x[N],y[N],sa[N],log[N],c[N],rk[N],h[N],st[N][20]; 10 char s[N],s1[N],S[N]; 11 12 int Cmp(int a,int b,int l){ return y[a]==y[b] && y[a+l]==y[b+l]; } 13 14 void build(int m){ 15 memset(y,0,sizeof(y)); 16 rep(i,0,m) c[i]=0; 17 rep(i,1,n) c[x[i]]++; 18 rep(i,1,m) c[i]+=c[i-1]; 19 for (int i=n; i; i--) sa[c[x[i]]--]=i; 20 for (int k=1,p=0; p<n; k<<=1,m=p){ 21 p=0; 22 rep(i,n-k+1,n) y[++p]=i; 23 rep(i,1,n) if (sa[i]>k) y[++p]=sa[i]-k; 24 rep(i,0,m) c[i]=0; 25 rep(i,1,n) c[x[y[i]]]++; 26 rep(i,1,m) c[i]+=c[i-1]; 27 for (int i=n; i; i--) sa[c[x[y[i]]]--]=y[i]; 28 rep(i,1,n) y[i]=x[i]; p=1; x[sa[1]]=1; 29 rep(i,2,n) x[sa[i]]=Cmp(sa[i-1],sa[i],k) ? p : ++p; 30 } 31 } 32 33 void get(){ 34 int k=0; 35 rep(i,1,n) rk[sa[i]]=i; 36 rep(i,1,n){ 37 for (int j=sa[rk[i]-1]; i+k<=n && j+k<=n && S[i+k]==S[j+k]; k++); 38 h[rk[i]]=k; if (k) k--; 39 } 40 } 41 42 void rmq(){ 43 rep(i,1,n) st[i][0]=h[i]; 44 rep(i,1,log[n]) 45 rep(j,1,n-(1<<i)+1) st[j][i]=min(st[j][i-1],st[j+(1<<(i-1))][i-1]); 46 } 47 48 int ask(int l,int r){ 49 if (l>r) swap(l,r); 50 l++; int t=log[r-l+1]; 51 return min(st[l][t],st[r-(1<<t)+1][t]); 52 } 53 54 int main(){ 55 log[1]=0; rep(i,2,N-100) log[i]=log[i>>1]+1; 56 for (scanf("%d",&T); T--; ){ 57 scanf("%s",s+1); scanf("%s",s1+1); 58 int n0=strlen(s+1),n1=strlen(s1+1); 59 rep(i,1,n0) S[i]=s[i]; S[strlen(s+1)+1]='$'; 60 rep(i,1,n1) S[strlen(s+1)+i+1]=s1[i]; 61 n=n0+n1+1; 62 rep(i,1,n) x[i]=(int)S[i]; 63 build(300); get(); rmq(); ans=0; 64 rep(i,1,n0-n1+1){ 65 int a=i,b=n0+2,f=0; 66 rep(j,1,4){ 67 int t=ask(rk[a],rk[b]); 68 if (b+t>n) { f=1; break; } 69 if (i+t>n0) break; 70 a+=t+1; b+=t+1; 71 } 72 if (f) ans++; 73 } 74 printf("%d\n",ans); 75 } 76 return 0; 77 }