BZOJ4892: [Tjoi2017]dna

Description

加里敦大学的生物研究所，发现了决定人喜不喜欢吃藕的基因序列S,

有这个序列的碱基序列就会表现出喜欢吃藕的性状。

但是研究人员发现对碱基序列S,任意修改其中不超过3个碱基，依然能够表现出吃藕的性状。

现在研究人员想知道这个基因在DNA链S0上的位置。

所以你需要统计在一个表现出吃藕性状的人的DNA序列S0上，有多少个连续子串可能是该基因。

即有多少个S0的连续子串修改小于等于三个字母能够变成S。

Input

第一行有一个数T,表示有几组数据

每组数据第一行一个长度不超过10^5的碱基序列S0

每组数据第二行一个长度不超过10^5的吃藕基因序列S

Output

共T行，第i行表示第i组数据中，在S0中有多少个与S等长的连续子串可能是表现吃藕性状的碱基序列

Sample Input

1
ATCGCCCTA
CTTCA

Sample Output

题解Here!

这题解法比较多。

后缀数组：

先求出后缀数组，然后对$height$数组建立$RMQ$。

于是$LCP$就可以$O(1)$求了。

我们枚举原串的每一个字符作为起点，如果相同就用$lcp$求，不相同就给记录不相同的计数器$++$。

保证不相同的计数器不超过$3$。

如果在$3$以内我们就统计答案。

注：这个好像在BZOJ上被卡常了，我也不知道为什么。。。

据说$SAM$跑得飞快，啥时候去学一学。。。

附代码：

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<cmath>
#define MAXN 200010
using namespace std;
int n,m;
int top,sa[MAXN],rk[MAXN],height[MAXN],tax[MAXN],tp[MAXN],f[MAXN][20];
char str[MAXN];
inline int read(){int date=0,w=1;char c=0;while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}return date*w;
}
void radixsort(){for(int i=0;i<=top;i++)tax[i]=0;for(int i=1;i<=n;i++)tax[rk[i]]++;for(int i=1;i<=top;i++)tax[i]+=tax[i-1];for(int i=n;i>=1;i--)sa[tax[rk[tp[i]]]--]=tp[i];
}
void suffixsort(){top=128;for(int i=1;i<=n;i++){rk[i]=str[i];tp[i]=i;}radixsort();for(int w=1,p=0;p<n;top=p,w<<=1){p=0;for(int i=1;i<=w;i++)tp[++p]=n-w+i;for(int i=1;i<=n;i++)if(sa[i]>w)tp[++p]=sa[i]-w;radixsort();swap(tp,rk);rk[sa[1]]=p=1;for(int i=2;i<=n;i++)rk[sa[i]]=(tp[sa[i-1]]==tp[sa[i]]&&tp[sa[i-1]+w]==tp[sa[i]+w])?p:++p;}
}
void getheight(){for(int i=1,j,k=0;i<=n;i++){if(k)k--;j=sa[rk[i]-1];while(str[i+k]==str[j+k])k++;height[rk[i]]=k;}
}
void step(){for(int i=1;i<=n;i++)f[i][0]=height[i];for(int i=1;(1<<i)<=n;i++)for(int j=1;j+(1<<i)-1<=n;j++)f[j][i]=min(f[j][i-1],f[j+(1<<(i-1))][i-1]);
}
int query(int l,int r){int x=rk[l],y=rk[r],k;if(x>y)swap(x,y);x++;k=(int)log2(y-x+1);return min(f[x][k],f[y-(1<<k)+1][k]);
}
void work(){int ans=0;for(int i=1;i<=m*2-n;i++){int s=0;for(int j=1;j<=n-m&&s<=3;){if(str[i+j-1]!=str[m+j]){s++;j++;if(s>3)break;}else j+=query(i+j-1,m+j);}if(s<=3)ans++;}printf("%d\n",ans);
}
void init(){scanf("%s",str+1);m=strlen(str+1);str[++m]='#';scanf("%s",str+m+1);n=strlen(str+1);suffixsort();getheight();step();
}
int main(){int t=read();while(t--){init();work();}return 0;
}

后缀自动机：

（坑，未填。。。）

FFT：

$FFT$？

对！$FFT$！

你不禁想问，这题是字符串匹配，跟$FFT$八竿子打不着啊。。。

但是假如转换一下匹配的结果呢？

把两个串中的一个（为了方便，选择短的串）翻转，卷积的结果就是在文本串每一位置匹配的结果。

$A,C,T,G$分开算，以$A$为例：

让一个串中有$A$的位置为$1$，无$A$的位置为$0$，另一个串相反。

这样两串字符不同处会为结果的值贡献$1$ 。

对于$<=3$的结果个数和为答案。

玄妙的解法啊。。。

当然，$FFT$常数巨大，BZOJ日常卡常。

附代码：

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<cmath>
#define MAXN (1<<19)
using namespace std;
const char f[4]={'A','C','G','T'};
int l1,l2;
int sum[MAXN];
char str_one[MAXN],str_two[MAXN];
inline int read(){int date=0,w=1;char c=0;while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}return date*w;
}
namespace FFT{int n,l,rev[MAXN];struct node{double r,i;node operator +(const node &p){return (node){r+p.r,i+p.i};}node operator -(const node &p){return (node){r-p.r,i-p.i};}node operator *(const node &p){return (node){r*p.r-i*p.i,r*p.i+i*p.r};}}a[MAXN],b[MAXN],d[MAXN];void DFT(node *a,int f){for(int i=0;i<l;i++)d[i]=a[rev[i]];for(int i=0;i<l;i++)a[i]=d[i];for(int i=2;i<=l;i<<=1){node wi=(node){cos(2.00*M_PI/i),f*sin(2.00*M_PI/i)};for(int k=0;k<l;k+=i){node w=(node){1.00,0.00},x,y;for(int j=0;j<i/2;j++){x=a[k+j];y=w*a[k+j+i/2];a[k+j]=x+y;a[k+j+i/2]=x-y;w=w*wi;}}}if(f==-1)for(int i=0;i<l;i++)a[i].r/=l;}void init(int x){/*for(n=l=1;l<x;l<<=1)n++;l<<=1;for(int i=0;i<l;i++)for(int j=0,t=i;j<n;j++,t>>=1){rev[i]<<=1;rev[i]|=t&1;}*/l=1;while (l<l1+l2) l<<=1;for (int i=0;i<l;++i)rev[i]=((i&1)?rev[i^1]+(l>>1):rev[i>>1]>>1);}void FFT(){DFT(a,1);DFT(b,1);for(int i=0;i<l;i++)a[i]=a[i]*b[i];DFT(a,-1);}void DSC(int l1,int l2){memset(rev,0,sizeof(rev));init(l1+l2);for(int num=0;num<4;num++){for(int i=0;i<=l;i++)a[i].r=a[i].i=b[i].r=b[i].i=0.00;for(int i=0;i<l1;i++)a[i].r=(str_one[i+1]!=f[num]);for(int i=0;i<l2;i++)b[i].r=(str_two[i+1]==f[num]);FFT();for(int i=l2-1;i<l1;i++)sum[i]+=(int)(a[i].r+0.5);}}
}
void work(){int ans=0;for(int i=l2-1;i<l1;i++)ans+=(sum[i]<=3);printf("%d\n",ans);
}
void init(){memset(sum,0,sizeof(sum));scanf("%s%s",str_one+1,str_two+1);l1=strlen(str_one+1);l2=strlen(str_two+1);for(int i=1;i<=(l2>>1);i++)swap(str_two[i],str_two[l2-i+1]);FFT::DSC(l1,l2);
}
int main(){int t=read();while(t--){init();work();}return 0;
}

转载于:https://www.cnblogs.com/Yangrui-Blog/p/9446494.html