FFT | 2 | FFT在字符串匹配中的应用

FFT | 2 | FFT在字符串匹配中的应用

在上一篇文章里面我们介绍了 $FFT/IFFT$ 的基本原理和应用,今天我们来了解一下 $FFT$ 在字符串匹配中的神奇应用

简单的字符串匹配问题

给出两个字符串 $A$ 和 $B$ (长度分别为 $m$ 和 $n$ )

求 $A$ 在 $B$ 中出现了几次以及出现的位置

啊哈!简单字符串匹配问题,$KMP$ 可以在 $O(m+n)$ 的时间复杂度内解决

但是我们来考虑一下如何用 $FFT$ 来解决这个问题?


(下文中字符串下标均从1开始)

首先我们想到要将这个问题转换为一个数学的形式,因此我们将字符串 $A$ 和 $B$ 看成两个函数,$A(i)$ 和 $B(i)$ 分别代表字符串 $A$ 的第 $i$ 项和字符串 $B$ 的第 $i$ 项

我们再定义函数 $f(x,y)$ :

$$
f(x,y)=\left(A(x)-B(y)\right)^2
$$

因此 $f(x,y)$ 就是判断 $A$ 的第 $x$ 项是否和 $B$ 的第 $y$ 项相等,很容易得出 $f(x,y)$ 的性质:

$$ f(x,y)\begin{cases}=0&A(x)=B(y)\\>0&A(x)\ne B(y)\end{cases} $$

然后我们再定义 $g(x)$ 函数判断 $A$ 是否在 $B$ 的第 $x$ 位出现

$$
g(x)=\sum_{i=1}^mf(i,x+i-1)
$$

也就是说,只有当 $g(x)=0$ 时,$A$ 在 $B$ 的第 $x$ 位出现

我们将 $g(x)$ 进行展开

$$
\begin{align}
g(x)&=\sum_{i=1}^m\left(A(i)-B(x+i-1)\right)^2\\
&=\sum_{i=1}^mA^2(i)-\sum_{i=1}^m2A(i)B(x+i-1)+\sum_{i=1}^mB^2(x+i-1)
\end{align}
$$

令 $A’(x)=A(m-n+1)$ ,可得

$$
\begin{align}
g(x)&=\sum_{i=1}^mA^2(i)+\sum_{i=1}^mB^2(x+i-1)-\sum_{i=1}^m2A’(m-i+1)B(x+i-1)\\
&=\sum_{i=1}^mA^2(i)+\sum_{i=1}^mB^2(x+i-1)-\sum_{i+j=m+x}2A’(i)B(j)
\end{align}
$$

前两项可以通过预处理得出,后一项可以直接使用 $FFT$ 求出

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
#include <iostream>
#include <cstring>
#include <cmath>

#define nmax 100
#define P(n) cout<<#n"="<<(n)<<endl;
const double PI = acos(-1);
using namespace std;

struct Complex {
double x,y;
Complex() {}
Complex(double x, double y):x(x),y(y) {}
Complex operator+(const Complex& t) const {return Complex(x+t.x,y+t.y);}
Complex operator-(const Complex& t) const {return Complex(x-t.x,y-t.y);}
Complex& operator*=(const Complex& t) {static double q;q=x*t.x-y*t.y;y=x*t.y+y*t.x;x=q;return *this;}
Complex operator*(const Complex& t) const {Complex ret=*this;return ret*=t;}
};
char a[nmax],b[nmax];
int m,n;
int fi,pr[nmax];
int len;
Complex AA[nmax],B[nmax];
int rev[nmax];
int ans[nmax];
inline void fft(Complex* a, int inv) {
static Complex w,delta,x,y;
static int i,j,k;
for (i=0;i<len;i++)
if (i<rev[i]) swap(a[i],a[rev[i]]);
for (i=1;i<len;i<<=1) {
delta.x=cos(delta.y=PI/i);
delta.y=inv*sin(delta.y);
for (j=0;j<len;j+=(i<<1)) {
w.x=1,w.y=0;
for (k=0;k<i;k++,w*=delta) {
x=a[j+k],y=a[i+j+k]*w;
a[j+k]=x+y,a[i+j+k]=x-y;
}
}
}
if (inv==-1)
for (i=0;i<len;i++) a[i].x/=len;
}
int main() {
cin>>(a+1)>>(b+1);
m=strlen(a+1);
n=strlen(b+1);
register int i,j;
for (i=1;i<=m;i++) AA[m-i+1].x=a[i]-='a',fi+=a[i]*a[i];
for (i=1;i<=n;i++) B[i].x=b[i]-='a',pr[i]=pr[i-1]+b[i]*b[i];
i=j=m+n;
do {
len=i&(-i);
i&=~len;
} while (i);
if (len!=j) len<<=1;

for (i=1;i<len;i++) {
rev[i]=rev[i>>1]>>1;
if (i&1) rev[i]^=(len>>1);
}
fft(AA,1),fft(B,1);
for (i=0;i<len;i++) AA[i]*=B[i];
fft(AA,-1);
j=n-m+1;
for (int x=1;x<=j;x++)
if (!(fi+pr[x+m-1]-pr[x-1]-int(2*AA[m+x].x))) ans[++ans[0]]=x;
cout<<ans[0]<<endl;
for (i=1;i<=ans[0];i++)
cout<<ans[i]<<' ';
cout<<endl;
return 0;
}

不过这个时间复杂度比 $KMP$ 还高的 $FFT$ 版字符串匹配有什么实际应用呢?

我们来看下面一个问题

带通配符的字符串的匹配

P4173 残缺的字符串

给出两个带有星号通配符的 $A$ 和 $B$ ,其中星号代表这里可以匹配任何字符,例如:

aa*b可以和aacb和aaqb等等匹配

求 $A$ 在 $B$ 中出现了几次以及每次出现的位置

发现了什么吗??

我们只需要把上一个问题中的 $f(x,y)$ 函数改一下就好了!

我们令通配符代表的字符为0,并定义

$$
f(x,y)=A(x)\cdot B(y)\cdot(A(x)-B(y))^2
$$

这样我们的新函数就做好了!我们的新函数满足如下性质:

$$ f(x,y)\begin{cases}=0&A(x)=0\lor B(x)=0\lor A(x)=B(x)\\ \gt 0&otherwise.\end{cases} $$

那么我们再把这个式子代入到 $g(x)$ 中展开一下

$$
\begin{align}
g(x)&=\sum_{i=1}^mA(i)\cdot B(x+i-1)\cdot\left(A(i)-B(x+i-1)\right)^2\\
&=\sum_{i=1}^mA(i)\cdot B(x+i-1)\cdot\left(A^2(i)-2\cdot A(i)\cdot B(x+i-1)+B^2(x+i-1)\right)\\
&=\sum_{i=1}^mA^3(i)\cdot B(x+i-1)-\sum_{i=1}^m2\cdot A^2(i)\cdot B^2(x+i-1)+\sum_{i=1}^mA(i)\cdot B^3(x+i-1)\\
&=\sum_{i=1}^mA’^3(m-i+1)\cdot B(x+i-1)-\sum_{i=1}^m2\cdot A’^2(m-i+1)\cdot B^2(x+i-1)+\sum_{i=1}^mA’(m-i+1)\cdot B^3(x+i-1)
\end{align}
$$

这个式子明显不那么赏心悦目…不过仔细观察三个部分都是可以用 $FFT$ 求的,直接求出即可

提示

  • 在将字符串转换成数组时不能把’a’转化为0,只有通配符能为0

  • 三个部分分别 $IFFT$ 会超时,我们将三个部分 $FFT$ 后按原式加减,最后再一次 $FFT$ 即可

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
#include <cstring>
#include <cmath>
#include <iostream>

#define nmax 1048600
const double PI = acos(-1);
using namespace std;

struct Complex {
double x,y;
Complex() {}
Complex(double x, double y):x(x),y(y) {}
inline Complex& operator+=(const Complex& t) {x+=t.x;y+=t.y;return *this;}
Complex operator+(const Complex& t) const {Complex ret=*this;return ret+=t;}
inline Complex& operator-=(const Complex& t) {x-=t.x;y-=t.y;return *this;}
Complex operator-(const Complex& t) const {Complex ret=*this;return ret-=t;}
inline Complex& operator*=(const Complex& t) {static double q;q=x*t.x-y*t.y;y=x*t.y+y*t.x;x=q;return *this;}
Complex operator*(const Complex& t) const {Complex ret=*this;return ret*=t;}
};
char a[nmax],ao[nmax],b[nmax]; //ao是原来的,a是倒过来的
int m,n;
int len;
Complex A[nmax],B[nmax],C[nmax];
int rev[nmax];
int ans[nmax];
inline void fft(Complex* a, int inv) {
static Complex w,delta,x,y;
static int i,j,k;
for (i=0;i<len;i++)
if (i<rev[i]) swap(a[i],a[rev[i]]);
for (i=1;i<len;i<<=1) {
delta.x=cos(delta.y=PI/i);
delta.y=inv*sin(delta.y);
for (j=0;j<len;j+=(i<<1)) {
w.x=1,w.y=0;
for (k=0;k<i;k++,w*=delta) {
x=a[j+k],y=a[i+j+k]*w;
a[j+k]=x+y,a[i+j+k]=x-y;
}
}
}
if (inv==-1)
for (i=0;i<len;i++) a[i].x/=len;
}
int main() {
cin>>m>>n>>(ao+1)>>(b+1);
register int i,j;
for (i=1,j=m;i<=m;i++,j--)
if (ao[i]=='*') a[j]=0;
else a[j]=ao[i]-'a'+1; //要保证'a'对应的是1而不是0,只有通配符能是0
for (i=1;i<=n;i++)
if (b[i]=='*') b[i]=0;
else b[i]-='a'-1;
i=j=m+n;
do {
len=i&(-i);
i&=~len;
} while (i);
if (len!=j) len<<=1;

for (i=1;i<len;i++) {
rev[i]=rev[i>>1]>>1;
if (i&1) rev[i]^=(len>>1);
}
//第一个部分
A[0].x=A[0].y=B[0].x=B[0].y=0;
for (i=1;i<=len;i++) {
A[i].x=(double)a[i]*a[i]*a[i],A[i].y=0;
B[i].x=(double)b[i],B[i].y=0;
}
fft(A,1),fft(B,1);
for (i=0;i<=len;i++) A[i]*=B[i],C[i]+=A[i];
//第二个部分
A[0].x=A[0].y=B[0].x=B[0].y=0;
for (i=1;i<=len;i++) {
A[i].x=(double)a[i]*a[i],A[i].y=0;
B[i].x=(double)b[i]*b[i],B[i].y=0;
}
fft(A,1),fft(B,1);
Complex two(2,0);
for (i=0;i<=len;i++) A[i]*=B[i],A[i]*=two,C[i]-=A[i];
//第三个部分
A[0].x=A[0].y=B[0].x=B[0].y=0;
for (i=1;i<=len;i++) {
A[i].x=(double)a[i],A[i].y=0;
B[i].x=(double)b[i]*b[i]*b[i],B[i].y=0;
}
fft(A,1),fft(B,1);
for (i=0;i<=len;i++) A[i]*=B[i],C[i]+=A[i];

//最后再IFFT
fft(C,-1);

j=n-m+1;
for (int x=1;x<=j;x++)
if (!int(C[m+x].x+0.5)) ans[++ans[0]]=x;
cout<<ans[0]<<endl;
for (i=1;i<=ans[0];i++) cout<<ans[i]<<' ';
cout<<endl;
return 0;
}

吸氧能苟过QwQ

FFT | 2 | FFT在字符串匹配中的应用

https://mivik.gitee.io/2019/intro/fft-string/

作者

Mivik

发布于

2019-09-26

更新于

2022-08-08

许可协议

评论