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
",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
",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;
}