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
ATCGCCCTA
CTTCA
Sample Output
2
题解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;
}
