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;
}
时间: 2024-10-03 22:32:13

51nod 1565模糊搜索(FFT)的相关文章

51nod 1172 Partial Sums V2 任意模FFT

Partial Sums V2 题目连接: https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1172 Description 给出一个数组A,经过一次处理,生成一个 数组S,数组S中的每个值相当于数组A的累加,比如:A = {1 3 5 6} => S = {1 4 9 15}.如果对生成的数组S再进行一次累加操作,{1 4 9 15} => {1 5 14 29},现在给出数组A,问进行K次操作后的结果.(输出结果

51nod 1348 乘积之和 分治 + fft

给出由N个正整数组成的数组A,有Q次查询,每个查询包含一个整数K,从数组A中任选K个(K <= N)把他们乘在一起得到一个乘积.求所有不同的方案得到的乘积之和,由于结果巨大,输出Mod 100003的结果即可.例如:1 2 3,从中任选1个共3种方法,{1} {2} {3},和为6.从中任选2个共3种方法,{1 2} {1 3} {2 3},和为2 + 3 + 6 = 11. 预处理 + O(1)回答 很容易想到用分治来做,这样在分治过程中需要一个dp,这个dp递推式是O(n^2)的, 所以复杂

51Nod 快速傅里叶变换题集选刷

打开51Nod全部问题页面,在右边题目分类中找到快速傅里叶变换,然后按分值排序,就是本文的题目顺序. 1.大数乘法问题 这个……板子就算了吧. 2.美妙的序列问题 长度为n的排列,且满足从中间任意位置划分为两个非空数列后,左边的最大值>右边的最小值.问这样的排列有多少个%998244353. 多组询问,n,T<=100000. 题解:经过分析可知,不合法的排列一定存在这样一种划分: 我们考虑答案=f[i]=i!-不合法排列个数. 形如 2 1 3 4 6 5 这种排列,会有三种划分方式不合法(

ACM学习历程—51NOD1028 大数乘法V2(FFT)

题目链接:http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1028 题目大意就是求两个大数的乘法. 但是用普通的大数乘法,这个长度的大数肯定不行. 大数可以表示点值表示法,然后卷积乘法就能用FFT加速运算了. 这道题是来存模板的. 代码: #include <iostream> #include <cstdio> #include <cstdlib> #include <cmath>

XJTUOJ wmq的A&#215;B Problem FFT

wmq的A×B Problem 发布时间: 2017年4月9日 17:06   最后更新: 2017年4月9日 17:07   时间限制: 3000ms   内存限制: 512M 描述 这是一个非常简单的问题. wmq如今开始学习乘法了!他为了训练自己的乘法计算能力,写出了n个整数,并且对每两个数a,b都求出了它们的乘积a×b.现在他想知道,在求出的n(n−1)2个乘积中,除以给定的质数m余数为k(0≤k<m)的有多少个. 输入 第一行为测试数据的组数. 对于每组测试数据,第一行为2个正整数n,

51nod 1201 整数划分(dp)

题目链接:https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1201 题解:显然是一道dp,不妨设dp[i][j]表示数字i分成j个一共有几种分法. 那么转移方程式为: dp[i][j] = dp[i - 1][j] + dp[i - 1][j - 1] 表示将i - 1划分为j个数,然后j个数都+1 还是不重复,将i - 1划分为j - 1个数,然后j - 1个数都+1,再加上1这个数. 然后就是j的范围要知道1+2+

hdu 1565 方格取数(2)(网络流之最大点权独立集)

题目链接:hdu 1565 方格取数(2) 题意: 有一个n*m的方格,每个方格有一个数,现在让你选一些数.使得和最大. 选的数不能有相邻的. 题解: 我们知道对于普通二分图来说,最大独立点集 + 最小点覆盖集 = 总点数,类似的,对于有权的二分图来说,有: 最大点权独立集 + 最小点权覆盖集 = 总点权和, 这个题很明显是要求 最大点权独立集 ,现在 总点权 已知,我们只要求出来 最小点权覆盖集 就好了,我们可以这样建图, 1,对矩阵中的点进行黑白着色(相邻的点颜色不同),从源点向黑色的点连一

对AM信号FFT的matlab仿真

普通调幅波AM的频谱,大信号包络检波频谱分析 u(t)=Ucm(1+macos ?t)cos ?ct ma称为调幅系数 它的频谱由载波,上下边频组成 , 包络检波中二极管截去负半周再用电容低通滤波,可以得到基带信号,那么,截去负半周后的AM信号必定包含基带信号的频谱.我们可以通过matlab来验证. %已知基带信号为1hz,载波为64hz,调制系数ma=0.3,采样频率1024hz,FFT变换区间N为2048 clear; fs=1024; f=1; %1hz基带信号 fc=64; %64hz载

多项式FFT相关模板

自己码了一个模板...有点辛苦...常数十分大,小心使用 #include <iostream> #include <stdio.h> #include <math.h> #include <string.h> #include <time.h> #include <stdlib.h> #include <algorithm> #include <vector> using namespace std; #de