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$ 求出

代码实现

#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$ 即可

代码

#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