P3763 [TJOI2017]DNA 字符串匹配 快速傅里叶变换fft

题意

给你模式串和原串,在匹配字符数量<=3时认为匹配,输出原串中有多少子串在上述匹配规则下成功匹配。

分析

P4173 残缺的字符串 fft匹配字符串
采用上题类似的匹配函数,将四个字符分别计算, f f t fft fft后计算总的匹配数量+3是否 ≥ m \ge m m
字符中=A,G,C,T的置1,其他置0

#include <iostream>
#include <string>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <queue>
#include <iomanip>
#include <map>
#include <cstdio>
#include <stack>
#include <set>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int ,int > pii;
#define endl '\n'
ll gcd(ll a, ll b){
    return b == 0 ? a : gcd(b, a % b);
}
void input(){
    freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);
}
inline int read(){
	int x=0,f=1;char c=getchar();
	while(c<'0'||c>'9') {if(c=='-') f=-1;c=getchar();}
	while (c>='0'&&c<='9') x=(x<<3)+(x<<1)+(c^48),c=getchar();
	return x*f;
}
const int N = 1e6+10, M = N * 2, inf = 1e8;
const double PI = acos(-1);

int rev[N], bit, tot;
struct Complex{
    double x, y;
    Complex operator+ (const Complex& t) const { return {x + t.x, y + t.y}; }
    Complex operator- (const Complex& t) const{ return {x - t.x, y - t.y}; }
    Complex operator* (const Complex& t) const{ return {x * t.x - y * t.y, x * t.y + y * t.x}; }
}a[N], b[N];
void fft(Complex a[], int inv){
    for(int i = 0; i < tot; i++)
        if(i < rev[i]) swap(a[i], a[rev[i]]);
    for(int mid = 1; mid < tot; mid <<= 1){
        Complex w1 = Complex({cos(PI/mid), inv * sin(PI/mid)});
        for(int i = 0; i < tot; i += mid * 2){
            Complex wk = Complex({1, 0});
            for(int j = 0; j < mid; j ++, wk = wk * w1){
                Complex x = a[i + j], y = wk * a[i + j + mid];
                a[i + j] = x + y, a[i + j + mid] = x - y;
            }
        }
    }
}
int n, m, t, ans[N];
double A[N], B[N];
char s1[N], s2[N];
char DNA[] = {'A','G','C','T'};
int main(){
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    cin>>t;
    while(t--){
        memset(ans, 0, sizeof(ans));
        cin>>s2>>s1;
        int m = strlen(s1), n = strlen(s2);
        while((1 << bit) < n + m + 1) bit++;
        tot = 1 << bit;
        for(int i = 0; i < tot; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
        for(int cnt = 0; cnt < 4; cnt++){
            char ch = DNA[cnt];
            for(int i = 0; i < tot; i++) A[i] = B[i] = 0, a[i] = b[i] = {0, 0};
            for(int i = 0; i < m; i++) if(s1[i] == ch) A[i] = 1;
            for(int i = 0; i < n; i++) if(s2[i] == ch) B[i] = 1;
            reverse(A, A + m);
            for(int i = 0; i < m; i++) a[i] = {A[i], 0};
            for(int i = 0; i < n; i++) b[i] = {B[i], 0};
            fft(a, 1), fft(b, 1);
            for(int i = 0; i < tot; i++) a[i] = a[i] * b[i];
            fft(a, -1);
            for(int i = m-1; i < n; i++) ans[i] += (int)(a[i].x/tot + 0.5);
        }
        int res = 0;
        for(int i = m-1; i < n; i++){
            if(ans[i] + 3 >= m) res++;
        }
        cout<<res<<endl;
    }
    return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值