序言
接触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中的某个字符。)
如果有通配符,只要是通配符,那么就一定可以匹配,不妨先把通配符都改成0。
并没发现什么问题,平方展开试一试。
请大家强行脑补,回忆一下卷积的定义是什么?
在计算卷记得某一个位置k时,被使用到的i和j的和始终是一个定值。但本题中,两个相乘的数的下标始终是,i和i+x,他们的差是一个定值,但和不是。怎么办呢?为什么不把T字符串翻转一下!
这样,T’就是T翻转之后的结果。我们可以试着用T’去替代原式中的T,并且令其为f(x)。
这次是不是非常和谐了,当x一定时,被相乘的两个数的下标分别为l2-i-1和i+x,它们的和始终等于l2+x-1,是一个定值,满足卷积的性质。
把上式看成三个部分:
中间的那部分:
是一个定值,可以直接求出。
左边的那部分:
构造一个数组S’令它为S的平方,则有:
令S’和T的卷积为A,则有:
根据定义有:
因为T’的长度只有l2,我们默认T’在l2及以后的位置的值都是零(也就是假定都是通配符,可以和任何字符匹配)。那么这个式子的后半部分一定等于零,即:
尝试用l2-i-1替换i:i=0时,l2-i-1=l2-1:i=l2-1时,l2-i-1=0。
惊人的事情出现了:T’和S’的卷积A中的第l2-x+1位的值,恰好与f(x)中的左边部分相等。
同理,你也可以令T”表示T’的平方,即:
令T”与S的卷积为B,用同样的方法也可以证得B中的第l2-i+1位的值,恰好与f(x)中的右边部分相等。
这样的话,你只需要先预处理出中间的部分记为W,然后求一下S’和T’的卷积A,再求一下S和T”的卷积B就可以得出所谓f(x)。
用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优化卷积。
上文中,如有谬误,敬请指正。