51nod 1565模糊搜索(FFT)

题目大意就是字符串匹配，不过有一个门限k而已

之前有提到过fft做字符串匹配，这里和之前那种有些许不同

因为只有A,C,G,T四种字符，所以就考虑构造4个01序列

例如，模板串a关于'A'的01序列中，1代表这个位置可以匹配，而0则代表不能匹配。

这样构造出4个序列后，再对匹配串b做同样的处理

下面用a['A']代表a关于'A'的01序列，b同理

然后可以知道a['A'][i]&b['A'][i]如果是1则代表可以匹配，如果是0则代表不能匹配。

那么在位置i两个串能否匹配就可以写做

for(x = i ~ i+lenb) ans += a['A'][i]&b['A'][i]

如果ans恰好等于b中'A'的个数，那么就说明'A'匹配成功，接下来做'C','G','T'的情况

如果都可以匹配成功，那么这个位置就可以匹配了

如何用fft做这个呢？实际上也很简单

把b串颠倒一下就变成了a['A'][i]&b['A'][lenb-i]

就是一个卷积形式，于是就可以fft了

（测试感觉stl里的complex比较慢，非递归比递归快很多）

#include <iostream>

#include <cstdio>

#include <cstring>

#include <bitset>

#include <cmath>

#include <complex>

using namespace std;

const double pi = acos(-);

const int maxn = ;

int A[][maxn], B[][maxn], ANS[maxn];

struct cp {

    double a,b;

    cp() { a = b = ; }

    cp(double _a, double _b):a(_a), b(_b) {}

    cp operator+(const cp&y)const { return cp(a+y.a, b+y.b); }

    cp operator-(const cp&y)const { return cp(a-y.a, b-y.b); }

    cp operator*(const cp&y)const { return cp(a*y.a-b*y.b,a*y.b+b*y.a); }

    cp operator!()const { return cp(a,-b); }

};

int T[];

int la, lb, k;

char str[maxn];

class FFT {

    int n, L, R[maxn];

    cp a[maxn], b[maxn];

    void DFT(cp *a,int n,int f) {

        for(int i = ; i < n; i++) if(i < R[i]) swap(a[i], a[R[i]]);

        for(int i = ; i < n; i <<= ) {

            cp t, x, wn(cos(pi/i), sin(pi*f/i));

            for(int j = ; j < n; j += (i<<)) {

                cp w(, );

                for(int k = ; k < i; k++,w = w*wn) {

                    x = a[j+k];

                    t = w*a[j+k+i];

                    a[j+k] = x+t;

                    a[j+k+i] = x-t;

                }

            }

        }

    }

public:

    int c[maxn];

    FFT() { memset(R, , sizeof(R)); }

    void init(int *A, int *B, int n1, int n2) {

        memset(a, , sizeof(a));

        memset(b, , sizeof(b));

        for(int i = ; i <= n1; i++) a[i] = cp(A[i], );

        for(int i = ; i <= n2; i++) b[i] = cp(B[i], );

        n2 += n1;

        for(n = , L = ; n <= n2; n <<= ) L++;

        for(int i = ; i < n; i++) R[i] = (R[i>>]>>)|((i&)<<(L-));

    }

    void calc() {

        DFT(a, n, );

        DFT(b, n, );

        for(int i = ; i <= n; i++) a[i] = a[i]*b[i];

        DFT(a, n, -);

        for(int i = ; i <= n; i++) c[i] = (a[i].a/n + 0.5);

    }

} fft;

int main() {

    cin>>la>>lb>>k;

    T['A'] = ; T['C'] = ; T['G'] = ; T['T'] = ;

    scanf("%s", str);

    for(int i = ; i < la; i++) {

        int ty = T[str[i]];

        A[ty][i] = ;

        for(int j = i+; j <= min(i+k, la-); j++) {

            if(ty == T[str[j]]) break;

            A[ty][j] = ;

        }

        for(int j = i-; j >= max(, i-k); j--) {

            if(A[ty][j]) break;

            A[ty][j] = ;

        }

    }

    scanf("%s", str);

    for(int i = ; i < lb; i++) {

        int ty = T[str[i]];

        B[ty][lb-i-] = ;

    }

    for(int i = lb-; i <= la+lb-; i++) ANS[i] = ;

    for(int i = ; i < ; i++) {

        fft.init(A[i], B[i], la, lb);

        fft.calc();

        int t = ;

        for(int j = ; j < lb; j++) if(B[i][j]) t++;

        for(int j = lb-; j <= la+lb-; j++)

            ANS[j] &= (fft.c[j] == t);

    }

    int ans = ;

    for(int i = lb-; i <= la+lb-; i++) ans += ANS[i];

    cout<<ans<<endl;

}

巴特西

51nod 1565模糊搜索(FFT)

最新文章

热门文章