[BZOJ]4503 两个串：我的第一次FFT尝试

最新推荐文章于 2018-10-09 19:38:00 发布
最新推荐文章于 2018-10-09 19:38:00 发布
阅读量665 收藏
点赞数
分类专栏： 算法导论 文章标签： fft
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/ggn_2015/article/details/78574419
版权

序言

接触FFT半年了，水平一直停留在只会写个高精度乘法那个层次，所以就试着做了一道FFT的题。

要是我没感觉错的话，这是一道权限题，然而我并没有权限号，我只是写了个暴力对拍了一下…

题目

[BZOJ]4503 两个串

Time Limit: 20 Sec；Memory Limit: 256 MB；Submit: 398 Solved: 190；

Description

兔子们在玩两个串的游戏。给定两个字符串S和T，兔子们想知道T在S中出现了几次，
分别在哪些位置出现。注意T中可能有“?”字符，这个字符可以匹配任何字符。

Input

两行两个字符串，分别代表S和T

Output

第一行一个正整数k，表示T在S中出现了几次
接下来k行正整数，分别代表T每次在S中出现的开始位置。按照从小到大的顺序输出，S下标从0开始。

Sample Input

bbabaababaaaaabaaaaaaaabaaabbbabaaabbabaabbbbabbbbbbabbaabbbababababbbbbbaaabaaabbbbbaabbbaabbbbabab

a?aba?abba

Sample Output

0

HINT

S 长度不超过 10^5， T 长度不会超过 S。 S 中只包含小写字母， T中只包含小写字母和“?”

Source

数学 FFT 字符串 脑洞题

思路

首先，因为有通配符，所以不能KMP。而又因为数据范围1e5，所以时间复杂度被锁定在O(nlogn)左右，暴力匹配显然也是不可能的，考虑通过FFT强行优化。

（下文中，l1表示S串长度，l2表示T串长度。S表示文本串，T表示模式串。）

然后，在S中的某一个位置x，T得到成功匹配的充要条件是什么。先考虑没有通配符的情况：
(c1是S中的某个字符，c2是T中的某个字符。)

c1=c2⇔(c1−c2)2=0

Sx...x+l2−1=T⇔∑i=0l2−1(Si+x−Ti)2=0

如果有通配符，只要是通配符，那么就一定可以匹配，不妨先把通配符都改成0。

c1=c2⇔c2(c1−c2)2=0

Sx...x+l2−1=T⇔∑i=0l2−1Ti(Si+x−Ti)2=0

并没发现什么问题，平方展开试一试。

Sx...x+l2−1=T⇔∑i=0l2−1Ti(S2i+x+T2i−2Si+xTi)=0

Sx...x+l2−1=T⇔∑i=0l2−1TiS2i+x+T3i−2Si+xT2i=0

请大家强行脑补，回忆一下卷积的定义是什么？

C=A∗B⇔Ck=∑i+j=kAi×Bj

在计算卷记得某一个位置k时，被使用到的i和j的和始终是一个定值。但本题中，两个相乘的数的下标始终是，i和i+x，他们的差是一个定值，但和不是。怎么办呢？为什么不把T字符串翻转一下！

令T′i=Tl2−i−1

这样，T’就是T翻转之后的结果。我们可以试着用T’去替代原式中的T,并且令其为f(x)。

f(x)=∑i=0l2−1T′l2−i−1S2i+x+T′3l2−i−1−2Si+xT′2l2−i−1=0

f(x)=(∑i=0l2−1T′l2−i−1S2i+x)+(∑i=0l2−1T′3l2−i−1)−2(∑i=0l2−1Si+xT′2l2−i−1)=0

这次是不是非常和谐了，当x一定时，被相乘的两个数的下标分别为l2-i-1和i+x，它们的和始终等于l2+x-1，是一个定值，满足卷积的性质。

把上式看成三个部分：

中间的那部分：

∑i=0l2−1T′3l2−i−1

是一个定值，可以直接求出。

左边的那部分：

∑i=0l2−1T′l2−i−1S2i+x

构造一个数组S’令它为S的平方，则有：

S′i=S2i

令S’和T的卷积为A，则有：

Ak=∑i+j=kS′jT′i

根据定义有：

Al2+x−1=∑i=0l2+x−1S′l2+x−1−iT′i=(∑i=0l2−1S′l2+x−1−iT′i)+(∑i=l2l2+x−1S′l2+x−1−iT′i)

因为T’的长度只有l2，我们默认T’在l2及以后的位置的值都是零（也就是假定都是通配符，可以和任何字符匹配）。那么这个式子的后半部分一定等于零，即：

Al2+x−1=∑i=0l2+x−1S′l2+x−1−iT′i=(∑i=0l2−1S′l2+x−1−iT′i)+(∑i=l2l2+x−1S′l2+x−1−iT′i)=∑i=0l2−1S′l2+x−1−iT′i

尝试用l2-i-1替换i：i=0时，l2-i-1=l2-1：i=l2-1时，l2-i-1=0。

Al2+x−1=∑i=0l2−1S′l2+x−1−iT′i=∑l2−i−1=0l2−i−1=l2−1S′l2+x−1−(l2−i−1)T′(l2−i−1)=∑l2−i−1=0l2−i−1=l2−1S′l2+x−1−(l2−i−1)T′l2−i−1

=∑i=l2−1i=0T′l2−i−1S′i+x=∑i=0l2−1T′l2−i−1S2i+x

惊人的事情出现了：T’和S’的卷积A中的第l2-x+1位的值，恰好与f(x)中的左边部分相等。

同理，你也可以令T”表示T’的平方，即：

T′′i=T′2i

令T”与S的卷积为B，用同样的方法也可以证得B中的第l2-i+1位的值，恰好与f(x)中的右边部分相等。

这样的话，你只需要先预处理出中间的部分记为W，然后求一下S’和T’的卷积A，再求一下S和T”的卷积B就可以得出所谓f(x)。

f(x)=Al2+x−1−2×Bl2+x−1+W

用FFT求两次卷积，时间复杂度为O(nlogn)。其他操作的时间复杂度都为O(n)，总时间复杂度为O(nlogn)。

代码

我的代码不能保证正确性，但是通过了我自己的对拍。

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<queue>
#include<complex>
using namespace std;

typedef complex<double> cd;
const int maxl=262145;
const double PI=acos(-1.0);

int rev[maxl];

void fft(cd* a,int n,int dft){
    for(int i=0;i<n;i++){
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    }
    for(int step=1;step<n;step<<=1){
        cd wn=exp(cd(0,dft*PI/step));
        for(int j=0;j<n;j+=step<<1){
            cd wnk(1,0);
            for(int k=j;k<j+step;k++){
                cd x=a[k];
                cd y=wnk*a[k+step];
                a[k]=x+y;
                a[k+step]=x-y;
                wnk*=wn;
            }
        }
    }
    if(dft==-1){
        for(int i=0;i<n;i++)
            a[i]/=n;
    }
}

void get_rev(int bit){
    for(int i=0;i<(1<<bit);i++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }
}

char S[maxl],T[maxl];//empty

inline double sqr(double x){return x*x;}
inline double cub(double x){return x*x*x;}

double cubeB=0;
cd a[maxl],b[maxl];
cd sqrA[maxl],sqrB[maxl];
cd bSqrA[maxl],aSqrB[maxl];

int tran_str(char* s){
    for(int i=0;;i++){
        if(s[i]==0)
            return i;//return len
        s[i]=(s[i]=='?')?0:(s[i]-'a'+1);
    }
}

int l1,l2;
bool check(int x){
    return int((cubeB+bSqrA[l2+x-1]-2.0*aSqrB[l2+x-1]).real()+0.5)==0;
}

int main(){
    scanf("%s%s",S,T);
    l1=tran_str(S);l2=tran_str(T);
    reverse(T,T+l2);
    int bit=1,s=2;
    for(bit=1;(1<<bit)<l1+l2-1;bit++)s<<=1;
    get_rev(bit);
    for(int i=0;i<l1;i++){
        a[i]=S[i];sqrA[i]=sqr(S[i]);
    }
    for(int i=0;i<l2;i++){
        b[i]=T[i];sqrB[i]=sqr(T[i]);
        cubeB+=cub(T[i]);
    }
    fft(sqrA,s,1);fft(b,s,1);
    for(int i=0;i<s;i++)bSqrA[i]=sqrA[i]*b[i];
    fft(bSqrA,s,-1);//b*(a^2)
    fft(a,s,1);fft(sqrB,s,1);
    for(int i=0;i<s;i++)aSqrB[i]=a[i]*sqrB[i];
    fft(aSqrB,s,-1);//a*(b^2)
    int cnt=0;queue<int>ans;
    for(int x=0;x<=l1-l2;x++){
        if(check(x)){
            cnt++;
            ans.push(x);
        }
    }
    printf("%d\n",cnt);
    while(!ans.empty()){
        printf("%d\n",ans.front());
        ans.pop();
    }
    return 0;
}

我的数据生成器的代码：

#include<cstdio>
#include<cstdlib>
#include<ctime>
#include<algorithm>
using namespace std;

const int maxl=10001;
char S[maxl],T[maxl];

int rand(int L,int R){
    return rand()%(R-L+1)+L;
}

int main(){
    srand(time(NULL));
    int l1=rand(100,4000),l2=rand(2,10);
    for(int i=0;i<l1;i++)S[i]=rand('a','b');
    //因为数据是随机的，如果字母的范围太大可能最后成功匹配的次数很少
    //我认为只保留两种字母并不影响程序正确性的测试
    int st=rand(0,l1-l2-1);
    for(int i=0;i<l2;i++){
        T[i]=S[st+i];
        if(!rand(0,4))T[i]='?';
    }
    printf("%s\n%s\n",S,T);
    return 0;
}

我的暴力的代码：

#include<cstdio>
#include<cstdlib>
#include<queue>
#include<algorithm>
using namespace std;

const int maxl=10001;
char S[maxl],T[maxl];

bool check(int x){
    for(int i=0;T[i];i++){
        if(!(T[i]=='?' || S[x+i]==T[i])){
            return false;
        }
    }
    return true;
}

int main(){
    scanf("%s%s",S,T);
    int cnt=0;queue<int>ans;
    for(int i=0;S[i];i++){
        if(check(i)){
            cnt++;
            ans.push(i);
        }
    }
    printf("%d\n",cnt);
    while(!ans.empty()){
        printf("%d\n",ans.front());
        ans.pop();
    }
    return 0;
}

尾声

上文中的这种 通过反转一个向量 从而 把对应相乘操作转化成卷积的方法是一种使用FFT时比较常见的套路。其主要思想是把对应相乘位置的两个下标“差一定”的关系转化成两个下标“和一定”的关系，从而利用卷积处理，然后再利用FFT优化卷积。

上文中，如有谬误，敬请指正。

GGN_2015
关注
  • 0
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
专栏目录
BZOJ_problem_rank:BZOJ问题等级按交流比排序
05-25
#BZOJ Problem Rankrank.cpp 程序文件data.dat bzoj题库数据done.dat AC过的题，初始可以把所有A过的题粘进去，正常退出的话自动维护。black.dat 黑名单。选题时会跳过。错题、神题、没题面、不想做等等。//Thank ...
两个字符串——最长公共子序列、最长公共子串——动态规划
liu_feng_zi_的专栏
08-11 1681
一、最长公共子序列 在两个字符串中，某一些字符都存在于两个串中，并且某些字符的前后顺序与在两个字符串的前后顺序相同，这些字符组成的串，就是两个字符串的公共子序列，其中最长的公共子序列，就是两个字符串的最长公共子序列（Longest Common Subsequence, LCS） 如下面的两个字符串： abdkc bdakmn ak，bdk，bd都是两个字符串的公共子序列，其中bdk就是最长的公共子序列 二、核心思想 假设两个字符串：和 求两个字符串的子串和的最长公共子序列 分为...
最长公共子串方法及优化
大爷没钱
09-13 628
最长公共子串 问题 link 牛客网 - 问题链接 对于两个字符串，请设计一个时间复杂度为 O(m*n) 的算法(这里的m和n为两串的长度)，求出两串的最长公共子串的长度。这里的最长公共子串的定义为两个序列U1,U2,..Un和V1,V2,...Vn，其中Ui + 1 == Ui+1,Vi + 1 == Vi+1，同时Ui == Vi。 给定两个字符串A和B，同时给定两串的长度n和m。...
4503: 两个串
weixin_30753873的博客
10-09 94
4503: 两个串 https://www.lydsy.com/JudgeOnline/problem.php?id=4503 题意： 　　求第二个串在第一个中出现了几次，用通配符。求出每个串的起始位置。 分析： 　　bitset。 　　一共有26个字母，求出每个字母在第一个串中出现的位置。扫一遍第二个串，ans = p[[1]] & (p[T[2]]>>1)...
bzoj4503 两个串（dp+FFT）
Coco_T的博客
03-01 393
Description兔子们在玩两个串的游戏。给定两个字符串S和T，兔子们想知道T在S中出现了几次，分别在哪些位置出现。注意T中可能有“?”字符，这个字符可以匹配任何字符。 Input两行两个字符串，分别代表S和T Output第一行一个正整数k，表示T在S中出现了几次 接下来k行正整数，分别代表T每次在S中出现的开始位置。按照从小到大的顺序输出，S下标从0开始。 Sample Input...
BZOJ 4503 两个串
依然一笑作春温。
03-08 563
FFT+神奇的思路
BZOJ4503 两个串
luoyuef的博客
04-03 201
兔子们在玩两个串的游戏。给定两个字符串S和T，兔子们想知道T在S中出现了几次， 分别在哪些位置出现。注意T中可能有“?”字符，这个字符可以匹配任何字符。 给cjr的题解跪啦~我这种菜鸡根本想不到这题怎么做呀qwq 如果没有通配符 考虑两个串匹配的条件 对于每个ai有ai==bi 即ai-bi=0 可转化为 sigma （ai-bi）²=0 考虑通配符 如果有通配符 无论ai是否等于b...
bzoj 4503 两个串
weixin_30555515的博客
04-08 101
Description 兔子们在玩两个串的游戏。给定两个字符串S和T，兔子们想知道T在S中出现了几次， 分别在哪些位置出现。注意T中可能有“?”字符，这个字符可以匹配任何字符。 Input 两行两个字符串，分别代表S和T Output 第一行一个正整数k，表示T在S中出现了几次 接下来k行正整数，分别代表T每次在S中出现的开始位置。按照从小到大的顺序输出，S下标从...
BZOJ第一部分
11-03
BZOJ原题-BZOJP1000-P2000的题目，下载后可以离线做题。
BZOJ第四部分
11-04
BZOJ原题-BZOJP4001-P4406的题目，下载后可以离线做题。
bzoj FFT 的模版
12-20
bzoj FFT 的模版
BZOJ第三部分
11-04
BZOJ原题-BZOJP3001-P4000的题目，下载后可以离线做题。
【BZOJ4503】两个串
步六孤的博客
06-20 198
【题目链接】 BZOJ4503(权限题) 【前置技能】 利用FFT/NTT进行字符串匹配 【题解】 字符串中出现了通配符，一般的字符串算法就失去效果了。 先忽略通配符的问题。思考一下，发现一种算法：令每个位置Ak=∑i=0LenT−1(Si+k−LenT+1−Ti)Ak=∑i=0LenT−1(Si+k−LenT+1−Ti)A_k=\displaystyle\sum_{i=0}^...
【刷题】BZOJ 4503 两个串
weixin_30814329的博客
06-13 160
Description 兔子们在玩两个串的游戏。给定两个字符串S和T，兔子们想知道T在S中出现了几次， 分别在哪些位置出现。注意T中可能有“?”字符，这个字符可以匹配任何字符。 Input 两行两个字符串，分别代表S和T Output 第一行一个正整数k，表示T在S中出现了几次 接下来k行正整数，分别代表T每次在S中出现的开始位置。按照从小到大的顺序输出，S下标从0开始。 Sample Input...
FFT详解
热门推荐
总理同学的编程尝试
03-31 9万+
一直想学FFT，当时由于数学基础太差，导致啥都学不懂。请教了机房里的几位学长大神，结果还是没太明白。因此下定决心写一篇关于“FFT”的文章，一篇起码我能看得懂的“FFT”。
BZOJ4503: 两个串
weixin_30443731的博客
04-14 144
Description 兔子们在玩两个串的游戏。给定两个字符串S和T，兔子们想知道T在S中出现了几次， 分别在哪些位置出现。注意T中可能有“?”字符，这个字符可以匹配任何字符。 Input 两行两个字符串，分别代表S和T Output 第一行一个正整数k，表示T在S中出现了几次 接下来k行正整数，分别代表T每次在S中出现的开始位置。按照从小到大的顺序输...
（初学者）素数计算法
总理同学的编程尝试
03-29 8176
前言记得我刚开始学习OI的时候接触到的最早的算法莫过于素数算法。当然了，素数这个知识小学的时候我们就已经有所学习了。当时我一直有一个疯狂的设想，就是编写一个程序把我能算出来的素数都算出来。
BZOJ2179: FFT快速傅立叶 FFT实现高精度乘法
最新发布
05-14
这道题目是要求用 FFT 实现高精度乘法。我们可以先将两个高精度数转换为多项式，然后对这两个多项式进行相乘，最后将结果转换为高精度数输出即可。 具体来说，将两个高精度数 $a$ 和 $b$ 分别转换为多项式 $A(x)$ 和 $B(x)$： $$ A(x) = a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1} \\ B(x) = b_0 + b_1x + b_2x^2 + \cdots + b_{n-1}x^{n-1} $$ 其中 $n$ 是 $a$ 和 $b$ 的位数，为了方便计算，我们要将 $n$ 扩展到 $2^k$，并在低位补 0。即： $$ A(x) = a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1} + 0x^n + \cdots + 0x^{2^k-1} \\ B(x) = b_0 + b_1x + b_2x^2 + \cdots + b_{n-1}x^{n-1} + 0x^n + \cdots + 0x^{2^k-1} $$ 接下来我们需要用 FFT 对 $A(x)$ 和 $B(x)$ 进行多项式乘法，得到多项式 $C(x) = A(x) \times B(x)$。然后将 $C(x)$ 转换为高精度数输出即可。 注意，由于 FFT 实现的是循环卷积，而我们要求的是线性卷积，因此在进行 FFT 前需要对 $A(x)$ 和 $B(x)$ 进行点值扩展，即令 $A'(x) = A(x) \times x^k$，$B'(x) = B(x) \times x^k$，其中 $k = 2^k - n$。这样得到的 $C'(x) = A'(x) \times B'(x)$ 在进行 FFT 后，其前 $n$ 项即为 $C(x)$ 的系数项，可以直接转换为高精度数输出。 代码实现如下：

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值