P3763 [TJOI2017]DNA (FFT)

原题题面

加里敦大学的生物研究所，发现了决定人喜不喜欢吃藕的基因序列 $S$,有这个序列的碱基序列就会表现出喜欢吃藕的性状，但是研究人员发现对碱基序列 $S$，任意修改其中不超过 3 个碱基，依然能够表现出吃藕的性状。现在研究人员想知道这个基因在 DNA 链 $S_0$上的位置。所以你需要统计在一个表现出吃藕性状的人的 DNA 序列 $S_0$上，有多少个连续子串可能是该基因，即有多少个$S_0$的连续子串修改小于等于三个字母能够变成 $S$。

输入格式

第一行有一个整数$T$，表示有几组数据。
每组数据第一行一个长度不超过 $10^5$的碱基序列$S_0$。
每组数据第二行一个长度不超过$10^5$的吃藕基因序列 $S$。

输出格式

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

输入样例

1
ATCGCCCTA
CTTCA

输出样例

2

题面分析

这本是一道后缀自动机的题目，但也有题解用了FFT，本篇博客采用了FFT的方法。
我们记第一个串为$s1$，第二个串为$s2$。
对于字符集{A,C,G,T}\{A,C,G,T\}{A,C,G,T}，我们对于每个字符依次处理，并记匹配到了该字符时为1，未匹配到该字符时为0，得到两个数组$A,B$。
那么对于$s1$串中的$[k,k+len(s2)]$和$s2$的$[0,len(s2)-1]$可以成功匹配的位数为
$$\sum _{i=0}^{len(s2)-1}A[k+i]B[i]$$
这个格式长得就很像卷积，我们设$C[i]=B[len(s2)-i-1]$,得到原式等价于
$$\sum _{i=0}^{len(s2)-1}A[k+i]C[len(s2)-i-1]$$
即
$$\sum _{i+j=len(s2)+k-1}A[i]C[j]$$
用FFT算出四个字符的结果后相加，再统计$[len(s2)-1,len(s1)]$的位置的结果就是符合题意的串个数即可。
假设我们记FFt的结果中某个符合条件的为$ans[i]$,那么$ans[i]+3\ge len(s2)$。

AC代码(900ms)

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define cp complex<double>
const int maxn=4e5;
cp a[maxn+10];
cp b[maxn+10];
char s1[maxn+10];
char s2[maxn+10];
ll ans[maxn+10];
char charset[]={'A', 'C', 'G', 'T'};
namespace FFT
{
    const double pi=acos(-1.0);
    ll rev[maxn+10];
    int getBit(int k)
    {
        int bit=0;
        while((1<<bit)<k)
            bit++;
        return bit;
    }
    void solve(cp *a, int n, int inv)
    //inv是1时是系数转点值,-1是点值转系数
    {
        int bit=getBit(n);
        for(int i=0; i<n; ++i)
        {
            rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
            if (i<rev[i]) swap(a[i], a[rev[i]]);
        }
        for(int mid=1; mid<n; mid*=2)
        {
            cp w(cos(pi*1.0/mid), inv*sin(pi*1.0/mid));//单位根
            for(int i=0; i<n; i+=mid*2)
            {
                cp omega(1, 0);
                for(int j=0; j<mid; j++, omega*=w)
                {
                    cp x=a[i+j];
                    cp y=omega*a[i+j+mid];
                    a[i+j]=x+y;
                    a[i+j+mid]=x-y;//蝴蝶变换
                }
            }
        }
        if (inv==-1)
        {
            for(int i=0; i<n; i++)
            {
                a[i].real(a[i].real()/n);
            }
        }
    }
}
void solve()
{
    int t;
    scanf("%d", &t);
    while(t--)
    {
        scanf("%s%s", s1, s2);
        int len1=strlen(s1);
        int len2=strlen(s2);
        int len=1<<FFT::getBit(len1+len2);
        for(int i=0; i<len; ++i)
            ans[i]=0;
        for(int z=0; z<4; ++z)
        {
            for(int i=0; i<len1; ++i)
            {
                a[i].real(s1[i]==charset[z] ? 1 : 0);
                a[i].imag(0);
            }
            for(int i=len1;i<len;++i)
            {
                a[i]=cp(0,0);
            }
            for(int i=0; i<len2; ++i)
            {
                b[len2-i-1]=(s2[i]==charset[z] ? 1 : 0);
                b[i].imag(0);
            }
            for(int i=len2;i<len;++i)
            {
                b[i]=cp(0,0);
            }
            FFT::solve(a, len, 1);
            FFT::solve(b, len, 1);
            for(int i=0; i<len; ++i)
            {
                a[i]*=b[i];
            }
            FFT::solve(a, len, -1);
            for(int i=0; i<len; ++i)
            {
                ans[i]+=(ll) (a[i].real()+0.5);
//                printf("%lld ", ans[i]);
            }
//            printf("\n");
        }
        ll cnt=0;
        for(int i=len2-1; i<len1; ++i)
        {
            if (ans[i]+3>=len2)
                cnt++;
        }
        printf("%lld\n", cnt);
    }
}
signed main()
{
//    ios_base::sync_with_stdio(false);
//    cin.tie(0);
//    cout.tie(0);
#ifdef ACM_LOCAL
    freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);
    long long test_index_for_debug=1;
    char acm_local_for_debug;
    while(cin>>acm_local_for_debug)
    {
        cin.putback(acm_local_for_debug);
        if (test_index_for_debug>100)
        {
            throw runtime_error("Check the stdin!!!");
        }
        auto start_clock_for_debug=clock();
        solve();
        auto end_clock_for_debug=clock();
        cout<<"\nTest "<<test_index_for_debug<<" successful"<<endl;
        cerr<<"Test "<<test_index_for_debug++<<" Run Time: "
            <<double(end_clock_for_debug-start_clock_for_debug)/CLOCKS_PER_SEC<<"s"<<endl;
        cout<<"--------------------------------------------------"<<endl;
    }
#else
    solve();
#endif
    return 0;
}

后记

FFT居然可以处理部分字符串问题，绝了…
DrGilbert 2020.10.5