使用 FFT/NTT 进行字符串匹配

2025/5/20 OI 多项式

我们一般进行字符串匹配，有两种方法：KMP 和 Hash。
我们这里介绍一种新的方法，FFT。

前置知识：

FFT/NTT

不会的可以看看我的博客？这里

我们现在有两个由小写英文字母组成的字符串 \(A\) 和 \(B\)，现在 \(A\) 串是模式串，我们在 \(B\) 中进行匹配。
（\(A\) 的长度为 \(m\)，\(B\) 的长度为 \(n\)。）常见的做法是 KMP，它可以在 \(O(n+m)\) 的时间复杂度下完成这么一个问题，非常优秀。
但是我们现在的串中含有通配符 *，思考一下，会发现 KMP 是无法解决的。

普通的字符串匹配

我们定义一种差异函数 \(f(x,y)\)，表示 \(A_x\) 与 \(B_y\) 之间的“距离”。也就是对应的 ASCII 码的差值。很显然，若 \(f\) 函数的值为 \(0\)，说明这两个位置是同一个字符。
那么我们对于 \(B\) 上的某个长度为 \(m\) 的子串，可以定义一个函数 \(P(x)\)，表示以 \(B_{x}\) 结尾的长度为 \(m\) 的子串与 \(A\) 的差异值。那么有 \(P_{x}=\sum_{i=x-m+1}^{x}f(i-x+m-1,i)\)。也很显然，\(P(x)\) 如果为 \(0\)，说明 \(B\) 的子串 \(B[x-m+1,x]\) 与 \(A\) 是完全匹配的。

但是我们会发现一个问题。
对于串 AB 和串 BA，根据我们的函数计算出来，差异值为 \(0\)。
为什么？
因为正负性。
取消正负性有两种方法，一是绝对值，二是平方。显然绝对值是不好处理的。所以我们平方一下。

\[ P_{x}=\sum_{i=x-m+1}^{x}[f(i-x+m-1,i)]^2 \]

把 \(f\) 函数展开，可以得到：

\[ P_{x}=\sum_{i=x-m+1}^{x}(A_{i-x+m-1}-B_{i})^2 \]

如果我们直接将平方展开，会变成这样：

\[ P_{x}=\sum_{i=x-m+1}^{x}\left(A_{i-x+m-1}^2+B_{i}^2-2A_{i-x+m-1}B_{i}\right) \]

看不出来什么东西。

试试将求和符号丢进每一项看看。

\[ P_{x}=\sum_{i=x-m+1}^{x}A_{i-x+m-1}^2+\sum_{i=x-m+1}^{x}B_{i}^2-2\sum_{i=x-m+1}^{x}A_{i-x+m-1}B_{i} \]

发现第一项是一个定值，也就是 \(\sum_{i=0}^{m-1}A_{i}^{2}\)。
第二项对于 \(B_{i}^{2}\) 进行一个前缀和预处理，也可以快速算出。
第三项看着有点麻烦。我们发现 \(i-x+m-1-i\) 是一个定值 \(-x+m-1\)，而如果我们能将 差为定值，变成 和为定值，那就是一个卷积的形式了！
怎么变成 和为定值 呢？很简单，我们将 \(A\) 串进行翻转操作得到一个新的串 \(S\)，易得 \(S_{i}=A_{m-i-1},A_{i}=S_{m-i-1}\)，于是最后一项就变成了 \(S_{m-(i-x+m-1)-1}B_{i}=S_{x-i}B_{i}\)，这下和就是一个定值 \(x\) 了！我们便可以使用 FFT/NTT 来加速这一过程。

至此，我们使用 FFT/NTT，在 \(O(n\log n)\) 的时间复杂度下完成了普通的字符串匹配。

带通配符的字符串匹配

这才是我们要讲的重点，也就是这道题：
Luogu P4173 残缺的字符串

这道题中，我们的字符串 \(A\) 和 \(B\) 都多了一种字符：通配符。此时我们先前的差异函数 \(f(x,y)\) 显然就存在问题了，因为 * 的 ASCII 值是 \(42\)，而 \(a\) 的 ASCII 值则是 \(97\)。相减显然会出现问题。
我们知道，当 \(f\) 函数值为 \(0\) 时表示相同，那么我们不妨令 * 对应 \(0\) 这个值，其他的 a 到 b 则表示 \(1\) 到 \(26\)，那么我们新的差异函数 \(f(x,y)\) 就应该等于 \(A_{x}B_{y}(A_{x}-B_{y})^2\)。
同理：

\[ P_{x}=\sum_{i=x-m+1}^{x}A_{i-x+m-1}B_{i}(A_{i-x+m-1}-B_{i})^2 \]

我们不妨还是将 \(A\) 串翻转一下得到串 \(S\)，并且将式子展开化简。

\[ \begin{aligned} P_{x}&=\sum_{i=x-m+1}^{x}S_{x-i}B_{i}(S_{x-i}-B_{i})^2\\ &=\sum_{i=x-m+1}^{x}S_{x-i}^3B_{i}+S_{x-i}B_{i}^3-2S_{x-i}^2B_{i}^2\\ &=\sum_{i=x-m+1}^{x}S_{x-i}^3B_{i}+\sum_{i=x-m+1}^{x}S_{x-i}B_{i}^3-2\sum_{i=x-m+1}^{x}S_{x-i}^2B_{i}^2\\ \end{aligned} \]

不难发现每一项都是卷积的形式，直接 FFT/NTT 做就行了。
那么就结束啦！

使用 FFT 的参考代码：

// Problem: P4173 残缺的字符串
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P4173
// Memory Limit: 128 MB
// Time Limit: 2000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

#include<bits/stdc++.h>
using namespace std;
typedef complex<long double> fs;
const int MAXN=2e6+10;
const int mod1=1e9+7;
const int mod2=998244353;
const int inf_int=0x7f7f7f7f;
const long long inf_long=0x7f7f7f7f7f7f7f7f;
const long double pi=acos(-1);
const long double eps=1e-5;
char Buf[1<<23],*P1=Buf,*P2=Buf;
#define getchar() (P1==P2&&(P2=(P1=Buf)+fread(Buf,1,1<<23,stdin),P1==P2)?EOF:*P1++)
template<typename type>
inline void read(type &x){
    x=0;
    bool f=false;
    char ch=getchar();
    while(ch<'0'||ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+(ch^48),ch=getchar();
    if(f) x=-x;
}
template<typename type,typename... args>
inline void read(type &x,args&... y){
    read(x),read(y...);
}

int N,M,BIT,rev[MAXN],a[MAXN],b[MAXN];
fs A[MAXN],B[MAXN],C[MAXN];
string AA,BB;

void FFT(fs *s,int opt,int bit){
    int n=(1<<bit);
    for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<bit-1);
    for(int i=0;i<n;i++) if(i<rev[i]) swap(s[i],s[rev[i]]);
    for(int i=2;i<=n;i*=2){
        fs omega=fs(cos(2*pi/i),sin(opt*2*pi/i));
        for(int j=0;j<n;j+=i){
            fs z=fs(1,0);
            for(int k=j;k<j+i/2;k++){
                fs p=s[k],q=z*s[k+i/2];
                s[k]=p+q;
                s[k+i/2]=p-q;
                z*=omega;
            }
        }
    }
    if(opt==-1) for(int i=0;i<n;i++) s[i].real(s[i].real()/n);
}

signed main(){
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cin>>N>>M>>AA>>BB;
    for(int l=0,r=N-1;l<r;l++,r--) swap(AA[l],AA[r]);
    for(int i=0;i<N;i++) a[i]=(AA[i]=='*'?0:AA[i]-'a'+1);
    for(int i=0;i<M;i++) b[i]=(BB[i]=='*'?0:BB[i]-'a'+1);
    while((1<<BIT)<N+M-1) BIT++;
    for(int i=0;i<(1<<BIT);i++) A[i]=fs(a[i]*a[i]*a[i],0);
    for(int i=0;i<(1<<BIT);i++) B[i]=fs(b[i],0);
    FFT(A,1,BIT);FFT(B,1,BIT);
    for(int i=0;i<(1<<BIT);i++) C[i]+=A[i]*B[i];
    for(int i=0;i<(1<<BIT);i++) A[i]=fs(a[i],0);
    for(int i=0;i<(1<<BIT);i++) B[i]=fs(b[i]*b[i]*b[i],0);
    FFT(A,1,BIT);FFT(B,1,BIT);
    for(int i=0;i<(1<<BIT);i++) C[i]+=A[i]*B[i];
    for(int i=0;i<(1<<BIT);i++) A[i]=fs(a[i]*a[i],0);
    for(int i=0;i<(1<<BIT);i++) B[i]=fs(b[i]*b[i],0);
    FFT(A,1,BIT);FFT(B,1,BIT);
    for(int i=0;i<(1<<BIT);i++) C[i]-=fs(2,0)*A[i]*B[i];
    FFT(C,-1,BIT);
    int ans=0;
    for(int i=N-1;i<M;i++) if(fabs(C[i].real())<=eps) ans++;
    cout<<ans<<endl;
    for(int i=N-1;i<M;i++) if(fabs(C[i].real())<=eps) cout<<i-N+2<<" ";
    return 0;
}

使用 NTT 的参考代码（from：Untitled_unrevised）

#include <algorithm>
#include <bit>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

using std::cin;
using std::cout;
using std::string;
using std::swap;
using std::transform;
using std::vector;

using u32 = std::uint32_t;
using u64 = std::uint64_t;

constexpr u64 P = 998'244'353;

template<class Grp, class Exp, class Ty>
constexpr Ty pow(Grp base, Exp exp, Ty init) {
    for(; exp; exp >>= 1, base *= base)
        if(exp & 1)
            init *= base;
    return init;
}

template<class Grp, class Exp, class Ty>
constexpr Ty& pow_assign(Grp base, Exp exp, Ty &init) {
    for(; exp; exp >>= 1, base *= base)
        if(exp & 1)
            init *= base;
    return init;
}

struct Fp {
    
    u32 x;

    constexpr Fp(u32 value = 0) : x(value) {}

    constexpr Fp& operator+=(const Fp &other) {
        x += other.x;
        if(x >= P) x -= P;
        return *this;
    }
    constexpr Fp& operator-=(const Fp &other) {
        if(x < other.x) x += P;
        x -= other.x;
        return *this;
    }
    constexpr Fp& operator*=(const Fp &other) {
        x = u64(x) * other.x % P;
        return *this;
    }
    constexpr Fp& operator/=(const Fp &other) {
        return pow_assign(other, P - 2, *this);
    }

    constexpr friend Fp operator+(Fp lhs, const Fp &rhs) { return lhs += rhs; }
    constexpr friend Fp operator-(Fp lhs, const Fp &rhs) { return lhs -= rhs; }
    constexpr friend Fp operator*(Fp lhs, const Fp &rhs) { return lhs *= rhs; }
    constexpr friend Fp operator/(Fp lhs, const Fp &rhs) { return lhs /= rhs; }
};

constexpr Fp prim(3);
constexpr Fp co_prim = Fp(1) / prim;

template<class RanIt>
void FFT(RanIt vec, int loglen, bool type) {

    // static_assert(std::is_same_v<RanIt::velue_type, Fp>);

    size_t len = 1ull << loglen;

    vector<size_t> bitrev(len);
    for(size_t idx = 0; idx < len; ++idx) {
        bitrev[idx] = (bitrev[idx >> 1] >> 1) | ((idx & 1) << (loglen - 1));
        if(idx < bitrev[idx]) {
            swap(vec[idx], vec[bitrev[idx]]);
        }
    }

    for(int logseg = 0; logseg < loglen; ++logseg) {

        Fp omega = pow(
            type ? co_prim : prim,
            (P - 1) >> (logseg + 1),
            Fp(1)
        );
        size_t seg = 1ull << logseg;
        
        for(RanIt lspan = vec; lspan != vec + len; lspan += (seg + seg)) {
            Fp phi = Fp(1);
            RanIt rspan = lspan + seg;
            for(size_t off = 0; off < seg; ++off) {
                Fp g  = lspan[off];
                Fp hw = rspan[off] * phi;
                lspan[off] = g + hw;
                rspan[off] = g - hw;
                phi *= omega;
            }
        }
    }

    if(type) {
        Fp invN = Fp(1) / Fp(len);
        for(size_t idx = 0; idx < len; ++idx)
            vec[idx] *= invN;
    }
}

int main() {

    cin.tie(nullptr)->sync_with_stdio(false);

    size_t m, n;
    cin >> m >> n;

    string a, b;
    cin >> a >> b;

    int logext = std::bit_width(n - 1);
    size_t extn = 1ull << logext;

    vector<Fp> A, B;

    A.reserve(extn);
    for(char ch : a) {
        A.emplace_back(ch == '*' ? 0 : ch - 96);
    }
    std::reverse(A.begin(), A.end());
    A.resize(extn);

    B.reserve(extn);
    for(char ch : b) {
        B.emplace_back(ch == '*' ? 0 : ch - 96);
    }
    B.resize(extn);

    vector<Fp> Acube(extn), Asqr(extn);
    vector<Fp> Bcube(extn), Bsqr(extn);

    transform(A.begin(), A.end(), Asqr.begin(), [](const Fp &x) {
        return x * x;
    });
    transform(B.begin(), B.end(), Bsqr.begin(), [](const Fp &x) {
        return x * x;
    });
    transform(A.begin(), A.end(), Acube.begin(), [](const Fp &x) {
        return x * x * x;
    });
    transform(B.begin(), B.end(), Bcube.begin(), [](const Fp &x) {
        return x * x * x;
    });

    FFT(A.begin(), logext, false);
    FFT(Asqr.begin(), logext, false);
    FFT(Acube.begin(), logext, false);

    FFT(B.begin(), logext, false);
    FFT(Bsqr.begin(), logext, false);
    FFT(Bcube.begin(), logext, false);

    vector<Fp> C(extn);
    for(size_t i = 0; i < extn; ++i) {
        C[i] = Acube[i] * B[i] - Fp(2) * Asqr[i] * Bsqr[i] + A[i] * Bcube[i];
    }
    FFT(C.begin(), logext, true);

    vector<size_t> res;
    for(size_t i = m - 1; i < n; ++i) {
        if(C[i].x == 0) {
            res.push_back(i - (m - 1) + 1);
        }
    }

    cout << res.size() << '\n';
    for(size_t x : res)
        cout << x << ' ';
    
    return 0;
}

最后放一道课后习题：CF528D Fuzzy Search

LOADING

使用 FFT/NTT 进行字符串匹配

普通的字符串匹配

带通配符的字符串匹配