CF-528D Fuzzy Search(FFT字符串匹配)

题意：

给定一个模式串和目标串按下图方式匹配，错开位置不多于k

解题思路：

总共只有\(A C G T\)四个字符，那么我们可以按照各个字符进行匹配，比如按照\(A\)进行匹配时,当\(k=1\)时，我们将目标串

\(ACAT\)化作

\(1~0~1~0\)

模式串

\(AGCAATTCAT\)化作

\(1~1~1~1~1~1~0~1~1~1\)

同样是反置目标串

可以得到以x为匹配终点的位置的匹配函数\(p(X)=\sum_{i+j=x}A(i)B(j)\)

如此进行4次FFT，最后如果目标位置贡献等于目标串长度，则说明匹配成功

#include <bits/stdc++.h>

using namespace std;

/*    freopen("k.in", "r", stdin);

    freopen("k.out", "w", stdout); */

//clock_t c1 = clock();

//std::cerr << "Time:" << clock() - c1 <<"ms" << std::endl;

//#pragma comment(linker, "/STACK:1024000000,1024000000")

#define de(a) cout << #a << " = " << a << endl

#define rep(i, a, n) for (int i = a; i <= n; i++)

#define per(i, a, n) for (int i = n; i >= a; i--)

typedef long long ll;

typedef unsigned long long ull;

typedef pair<int, int> PII;

typedef pair<double, double> PDD;

typedef vector<int, int> VII;

#define inf 0x3f3f3f3f

const ll INF = 0x3f3f3f3f3f3f3f3f;

const ll MAXN = 1e6 + 7;

const ll MAXM = 1e6 + 7;

const ll MOD = 998244353;

const double eps = 1e-6;

const double pi = acos(-1.0);

template <class T>

inline void in(T &x)

{

    static char ch;

    static bool neg;

    for (ch = neg = 0; ch < '0' || '9' < ch; neg |= ch == '-', ch = getchar())

        ;

    for (x = 0; '0' <= ch && ch <= '9'; (x *= 10) += ch - '0', ch = getchar())

        ;

    x = neg ? -x : x;

}

struct Complex

{

    double x, y;

    Complex(double xx = 0, double yy = 0) { x = xx, y = yy; }

} a[MAXN], b[MAXN], c[MAXN], ans[MAXN];

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); } //不懂的看复数的运算那部分

int N, M;

int l, r[MAXN];

int limit = 1;

void FFT(Complex *A, int type)

{

    for (int i = 0; i < limit; i++)

        if (i < r[i])

            swap(A[i], A[r[i]]); //求出要迭代的序列

    for (int mid = 1; mid < limit; mid <<= 1)

    {                                                    //待合并区间的长度的一半

        Complex Wn(cos(pi / mid), type * sin(pi / mid)); //单位根

        for (int R = mid << 1, j = 0; j < limit; j += R)

        {                    //R是区间的长度，j表示前已经到哪个位置了

            Complex w(1, 0); //幂

            for (int k = 0; k < mid; k++, w = w * Wn)

            {                                                 //枚举左半部分

                Complex x = A[j + k], y = w * A[j + mid + k]; //蝴蝶效应

                A[j + k] = x + y;

                A[j + mid + k] = x - y;

            }

        }

    }

    /*if (type == -1)

        for (int i = 0; i < limit; ++i)

            a[i].x /= limit;//我们推过的公式里面有一个1/n这一项*/

}

char s[MAXN], t[MAXN];

void init(int N, int M)

{

    while (limit <= N + M)

        limit <<= 1, l++;

    for (int i = 0; i < limit; i++)

        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));

    memset(a, 0, sizeof(a));

    memset(b, 0, sizeof(b));

}

int change(char str)

{

    if (str == 'A')

        return 1;

    else if (str == 'T')

        return 2;

    else if (str == 'G')

        return 3;

    else

        return 4;

}

int pre[MAXN], cnt;

int main()

{

    int n, m, k;

    scanf("%d%d%d %s %s", &n, &m, &k, s, t);

    reverse(t, t + m);

    init(n, m);

    for (int ca = 1; ca <= 4; ca++)

    {

        cnt = -1;

        memset(pre, 0, sizeof(pre));

        memset(a, 0, sizeof(a));

        memset(b, 0, sizeof(b));

        for (int i = 0; i < n; i++)

        {

            if (change(s[i]) == ca)

                pre[++cnt] = i;

            a[i].x = change(s[i]) == ca ? 1 : 0, a[i].y = 0;

        }

        for (int i = 0; i < m; i++)

            b[i].x = change(t[i]) == ca ? 1 : 0, b[i].y = 0;

        int now = -1;

        for (int i = 0; i <= cnt; i++)

        {

            int L = max(pre[i] - k, 0);

            int R = min(pre[i] + k, n - 1);

            if (now > R)

                continue;

            now = max(L, now);

            for (; now <= R; now++)

                a[now].x = 1;

            now--;

        }

        FFT(a, 1);

        FFT(b, 1);

        for (int i = 0; i < limit; i++)

            a[i] = b[i] * a[i];

        FFT(a, -1);

        for (int i = 0; i < limit; i++)

            c[i] = c[i] + a[i];

    }

    int ans = 0;

    for (int i = m - 1; i < limit; i++)

        if (int(c[i].x / limit + 0.5) == m)

            ans++;

    printf("%d\n", ans);

    return 0;

}

巴特西

CF-528D Fuzzy Search(FFT字符串匹配)

题意：

解题思路：

最新文章

热门文章