zoukankan      html  css  js  c++  java
  • 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(-1);
    const int maxn = 800050;
    int A[4][maxn], B[4][maxn], ANS[maxn];
    struct cp {
        double a,b;
        cp() { a = b = 0; }
        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[128];
    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 = 0; i < n; i++) if(i < R[i]) swap(a[i], a[R[i]]);
            for(int i = 1; i < n; i <<= 1) {
                cp t, x, wn(cos(pi/i), sin(pi*f/i));
                for(int j = 0; j < n; j += (i<<1)) {
                    cp w(1, 0);
                    for(int k = 0; 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, 0, sizeof(R)); }
        void init(int *A, int *B, int n1, int n2) {
            memset(a, 0, sizeof(a));
            memset(b, 0, sizeof(b));
            for(int i = 0; i <= n1; i++) a[i] = cp(A[i], 0);
            for(int i = 0; i <= n2; i++) b[i] = cp(B[i], 0);
            n2 += n1;
            for(n = 1, L = 0; n <= n2; n <<= 1) L++;
            for(int i = 0; i < n; i++) R[i] = (R[i>>1]>>1)|((i&1)<<(L-1));
        }
        void calc() {
            DFT(a, n, 1);
            DFT(b, n, 1);
            for(int i = 0; i <= n; i++) a[i] = a[i]*b[i];
            DFT(a, n, -1);
            for(int i = 0; i <= n; i++) c[i] = (a[i].a/n + 0.5);
        }
    } fft;
    
    int main() {
        cin>>la>>lb>>k;
        T['A'] = 0; T['C'] = 1; T['G'] = 2; T['T'] = 3;
        scanf("%s", str);
        for(int i = 0; i < la; i++) {
            int ty = T[str[i]];
            A[ty][i] = 1;
            for(int j = i+1; j <= min(i+k, la-1); j++) {
                if(ty == T[str[j]]) break;
                A[ty][j] = 1;
            }
            for(int j = i-1; j >= max(0, i-k); j--) {
                if(A[ty][j]) break;
                A[ty][j] = 1;
            }
        }
        scanf("%s", str);
        for(int i = 0; i < lb; i++) {
            int ty = T[str[i]];
            B[ty][lb-i-1] = 1;
        }
        for(int i = lb-1; i <= la+lb-2; i++) ANS[i] = 1;
        for(int i = 0; i < 4; i++) {
            fft.init(A[i], B[i], la, lb);
            fft.calc();
            int t = 0;
            for(int j = 0; j < lb; j++) if(B[i][j]) t++;
            for(int j = lb-1; j <= la+lb-2; j++)
                ANS[j] &= (fft.c[j] == t);
        }
        int ans = 0;
        for(int i = lb-1; i <= la+lb-2; i++) ans += ANS[i];
        cout<<ans<<endl;
    }
  • 相关阅读:
    [No0000199]设计模式总结
    [No0000197]Windows用户都应该知道的运行命令
    [No000017F]如何监控注册表的修改
    [No0000196]一文读懂Java 11的ZGC为何如此高效
    [No0000195]NoSQL还是SQL?这一篇讲清楚
    [No0000194]聊聊 Chrome DevTools 中你可能不知道的调试技巧
    [No000018A]改善C#程序的建议11-20
    [No000018C]Vim清除上次的搜索高亮结果-Vim使用技巧(1)
    [No000018D]Vim快速注释/取消注释多行的几种方法-Vim使用技巧(2)
    [No000018E]Vim快速跳转任意行、任意列以及高亮显示当前行、当前列方法-Vim使用技巧(3)
  • 原文地址:https://www.cnblogs.com/Saurus/p/6909546.html
Copyright © 2011-2022 走看看