题解 P3763 【[TJOI2017]DNA】

· · 题解

题意：

给出两个长度小于 10^5 的字符串 S 和 T，求 S 中有多少个连续子串在不改变超过 3 个字符后可以和 T 匹配。多组数据，字符集大小为 4。

设 S 的长度为 n， T 的长度为 m。

字符集大小为 4，考虑对每个字符分开考虑。

样例：

S:\ \texttt{ATCGCCCTA} T:\ \texttt{CTTCA}

不妨先对字符 \texttt{C} 进行处理。将 S 和 T 中为 \texttt{C} 的位置变为 1 ，不为 \texttt{C} 的位置变为 0 再进行匹配。

S:\ \texttt{001011100} T:\ \texttt{10010}

那么对于某个位置 k， S[k\ ...\ k+m-1] 和 T[0\ ...\ m-1] 的匹配个数可以写为：

\sum\limits_{i=0}^{m-1}S[i+k]T[i]

然后我们发现这个式子特别像 FFT，于是进行一番处理：

设 A[m-i-1] = T[i]

上式变为：

\sum\limits_{i=0}^{m-1}S[i+k]A[m-i-1] \sum\limits_{i+j=m-k-1}S[i]A[j]

然后跑 FFT 就可以快速求出上式。

最后位置 S[k\ ...\ k+m-1] 和 T[0\ ...\ m-1] 的匹配个数就是这四次 FFT 结果相加。

因为可以修改三个位置，最后判断是否匹配可以有 3 的误差。

#include<bits/stdc++.h>
#define ll long long
using namespace std;
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
const int maxn=100010*4;
struct Complex{
    double x,y;
};
Complex operator +(Complex a,Complex b){return (Complex){a.x+b.x,a.y+b.y};}
Complex operator -(Complex a,Complex b){return (Complex){a.x-b.x,a.y-b.y};}
Complex operator *(Complex a,Complex b){return (Complex){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
const double Pi=acos(-1.0);
Complex get_w(int k,int len,bool op){
    int T=1;if(!op)T=-1;
    return (Complex){cos(2.0*k*Pi/len),T*sin(2.0*k*Pi/len)};
}
Complex A[maxn],B[maxn];
int i,j,k,n,m;
int R[maxn];
char C[4]={'A','T','C','G'};
void FFT(int len,Complex *a,bool op){
    for(int i=0;i<len;i++)
        if(i<R[i])swap(a[i],a[R[i]]);
    for(int i=1;i<len;i=i*2){
        Complex w,Wn=get_w(1,i*2,op);
        for(int j=0;j<len;j+=i*2){
            w=(Complex){1,0};
            for(int k=0;k<i;k++,w=w*Wn){
                int K=k+j;
                Complex S1=a[K],S2=w*a[K+i];
                a[K]=S1+S2;
                a[K+i]=S1-S2;
            }
        }
    }
}
int Ans[maxn],ans;
int main(){
    int T;
    cin>>T;
    string None_;
    getline(cin,None_);
    while(T--){
        memset(Ans,0,sizeof(Ans));ans=0;
        string S,T;
        cin>>S>>T;
        n=S.size();
        m=T.size();
        for(i=0;i<4;i++){
            for(j=0;j<n;j++)
                A[j]=(Complex){(S[j]==C[i]),0};
            for(j=0;j<m;j++)
                B[m-j-1]=(Complex){(T[j]==C[i]),0};
            int len=1,L=0;
            while(len<=n+m)len*=2,++L;
            for(j=0;j<len;j++)
                R[j]=R[j/2]/2+(j%2)*(1<<L-1);
            FFT(len,A,true);FFT(len,B,true);
            for(j=0;j<len;j++)A[j]=A[j]*B[j];
            FFT(len,A,false);
            for(j=m-1;j<=m+n-1;j++)
                Ans[j-(m-1)]+=((int)(A[j].x/len+0.5));
            for(j=0;j<len;j++)A[j]=B[j]=(Complex){0,0};
        }
        for(i=0;i<n-m+1;i++)
            if(Ans[i]+3>=m)++ans;
        cout<<ans<<endl;
    }
    return 0;
}