Yharim
Snow is on the ground
So NVIDIA, Fuck You🖕
数学 信息学 算法 |多项式
发布于 2023-12-22

快速傅里叶变换：应用

FFT 的所有应用都围绕着它的本质功能：计算两个多项式的乘积。

大整数乘法（高精度乘法） #

假设有两个整数 $a=537,b=721$。$a,b$ 的每一位可以视为多项式 $A(x),B(x)$ 的系数，就像这样

$$A(x)=5x^2+3x+7, B(x)=7x^2+2x+1$$

计算 $C(x)=A(x)\cdot B(x)$：

$$C(x)=35x^4+31x^3+60x^2+17x+7$$

提取 $C(x)$ 的系数：

$$\{35,31,60,17,7\}$$

从右往左，逢 $10$ 进 $1$：

$$\{3,8,7,1,7,7\}$$

这样就能得到 $537\times721$ 的结果是 $387177$。

中间的多项式乘法由 $\text{FFT}$ 进行。如果 $a$ 和 $b$ 是 $n$ 位整数，则总时间复杂度是 $O(n\log{n})$。

整数乘法和多项式乘法具有高度的相似性。

一个整数可以分解为它的位值的和。例如，整数 $123$ 可以分解为 $1\times10^2+2\times10+3$。

设 $a=\overline{a_2a_1a_0}$，$b=\overline{b_2b_1b_0}$，则

$$\begin{aligned} &a \times b \\ = \ & (a_2 \cdot 10^2 + a_1 \cdot 10 + a_0) \cdot (b_2 \cdot 10^2 + b_1 \cdot 10 + b_0) \\ = \ & a_2b_2 \cdot 10^4 + (a_2b_1 + a_1b_2) \cdot 10^3 + (a_2b_0 + a_1b_1 + a_0b_2) \cdot 10^2 + (a_1b_0 + a_0b_1) \cdot 10 + a_0b_0 \end{aligned}$$

设 $A(x)=a_2x^2+a_1x+a_0$，$B(x)=b_2x^2+b_1x+b_0$，则

$$\begin{aligned} &A(x) \times B(x)\\ = \ & (a_2x^2 + a_1x + a_0) \cdot (b_2x^2 + b_1x + b_0) \\ = \ & a_2b_2\cdot x^4 + (a_2b_1 + a_1b_2)\cdot x^3 + (a_2b_0 + a_1b_1 + a_0b_2)\cdot x^2 + (a_1b_0 + a_0b_1)\cdot x + a_0b_0 \end{aligned}$$

可以看出整数乘法和多项式乘法基本是同构的。整数乘法在最后一步需要额外处理进位。

模板 # 

// put FFT algorithm template here

vector<Comp> string_to_vector(const string& s) {
    vector<Comp> result;
    for (auto rit = s.rbegin(); rit != s.rend(); ++rit) {
        result.push_back(*rit - '0');
    }
    return result;
}

string vector_to_string(const vector<Comp>& v) {
    string result;
    for (const Comp& c : v) {
        result += to_string(int(c.real())) + " ";
    }
    return result;
}

string upgrade(vector<Comp>& a) {
    string result;
    int carry = 0;
    for (auto i : a) {
        int num = int(real(i) + 0.5) + carry;
        carry = num / 10;
        result += (num % 10) + '0';
    }
    while (!result.empty() && result.back() == '0') {
        result.pop_back();
    }
    reverse(result.begin(), result.end());
    return result.empty() ? "0" : result;
}

string multiply_big_integers(const string& s1, const string& s2) {
    vector<Comp> a = string_to_vector(s1);
    vector<Comp> b = string_to_vector(s2);
    vector<Comp> c = multiply(a, b);
    return vector_to_string(c);
}

int main() {
    string num1, num2;
    cin >> num1 >> num2;

    string product = multiply_big_integers(num1, num2);
    cout << product << endl;
    
    return 0;
}

字符串匹配 #

朴素字符串匹配 #

设有两个字符串 $\text{text}$，$\text{pattern}$ 的长度分别为 $n,m(n\ge m)$，并只包含小写字母 $a\sim z$。问 $\text{pattern}$ 在 $\text{text}$ 中出现了几次，以及所有的出现位置。

该问题可使用 $\text{FFT}$ 在 $O(n\log{n})$ 的复杂度下解决。

tip

用 $\text{FFT}$ 做字符串匹配是一个差想法。$\text{KMP}$ 可以在 $O(n+m)$ 的线性复杂度下解决字符串匹配问题，并且支持通配符。

约定 $\text{str}\Big|_a^b$ 表示从 $\text{str}[a]$ 到 $\text{str}[b]$ 的子串。

首先将 $\text{text}$ 和 $\text{pattern}$ 的所有字母都映射成数字（$a\rightarrow 1,b\rightarrow 2,\cdots$），得到两个整型数组 $A$ 和 $B$。

$$\begin{aligned} \text{text}&\rightarrow A:\{A_0,A_1,\cdots,A_{n-1}\}\\ \text{pattern}&\rightarrow B:\{B_0,B_1,\cdots,B_{m-1}\} \end{aligned}$$

容易发现

$$\left(A_i-B_j\right)^2=\begin{cases} 0,&A_i=B_j\\ \\ >0,&A_i\not=B_j \end{cases}$$

进而构造

$$\begin{aligned} P_k&=(A_k-B_0)^2+(A_{k+1}-B_1)^2+\cdots+(A_{k+m-1}-B_{m-1})^2\\ &=\sum_{i=0}^{m-1}\left(A_{k+i}-B_i\right)^2\\ &=\begin{cases} 0,&\text{text}\Big|_{k}^{k+m-1}=\text{pattern}\\ \\ >0,&\text{text}\Big|_{k}^{k+m-1}\not=\text{pattern} \end{cases} \end{aligned}$$

因此可以利用 $P_k$ 判断 $\text{pattern}$ 是否出现在 $\text{text}$ 的第 $k$ 位。

展开 $P_k$：

$$\begin{aligned} P_k&=\sum_{i=0}^{m-1}\left(A_{k+i}-B_i\right)^2\\ &=\sum_{i=0}^{m-1}\left(A_{k+i}^2+B_{i}^2-2A_{k+i}B_{i}\right)\\ &=\sum_{i=0}^{m-1}A_{k+i}^2+\sum_{i=0}^{m-1}B_{i}^2-2\sum_{i=0}^{m-1}A_{k+i}B_{i} \end{aligned}$$

将数组 $B$ 反转得到 $\widetilde{B}$，有 $B_i=\widetilde{B}_{m-i-1}$。代入得

$$P_k=\sum_{i=0}^{m-1}A_{k+i}^2+\sum_{i=0}^{m-1}B_{i}^2-2\sum_{i=0}^{m-1}A_{k+i}\widetilde{B}_{m-i-1}$$

分析该式的三项：

  1. $\sum_{i=0}^{m-1}A_{k+i}^2$ 可以通过前缀和技巧实现 $O(n)$ 预处理，$O(1)$ 查询；
  2. $\sum_{i=0}^{m-1}B_{i}^2$ 是常数，可以直接 $O(n)$ 预处理；
  3. $\sum_{i=0}^{m-1}A_{k+i}\widetilde{B}_{m-i-1}$ 是卷积式，可以通过 $\text{FFT}$ 实现 $O(n\log{n})$ 预处理，$O(1)$ 查询。

对于第 3 项，更具体地说：

构造多项式

$$A(x)=A_0+A_1x+A_2x^2+\cdots+A_{n-1}x^{n-1}\\ \widetilde{B}(x)=\widetilde{B}_0+\widetilde{B}_1x+\widetilde{B}_2x^2+\cdots+\widetilde{B}_{m-1}x^{m-1}$$

利用 $\text{FFT}$ 可以 $O(n\log{n})$ 得出多项式 $C(x)=A(x)\cdot \widetilde{B}(x)$ 的系数向量。易知 $C(x)$ 中 $x^{k+m-1}$ 前的系数为 $\sum_{i=0}^{m-1}A_{k+i}\widetilde{B}_{m-i-1}$，这正是我们想要的第 3 项。

依次计算 $P_k(k=0,1,2,\cdots,n-m)$，当得到 $0$ 时说明 $\text{pattern}$ 出现在 $\text{text}$ 的第 $k$ 位。

模板 # 

// put FFT algorithm template here

vector<int> stringMatchingFFT(const string& text, const string& pattern) {
    int n = text.size(), m = pattern.size();
    
    vector<Comp> A(n), B(m);
    long long sum_AA[n], sum_BB = 0;

    for (int i = 0; i < n; i ++) {
        int value = text[i] - 'a' + 1;
        A[i] = value;
        sum_AA[i] = value * value;
        if (i) {
            sum_AA[i] += sum_AA[i - 1];
        }
    }

    for (int i = 0; i < m; i ++) {
        int value = pattern[i] - 'a' + 1;
        B[m - i - 1] = value;
        sum_BB += value * value;
    }

    vector<Comp> C = multiply(A, B);

    vector<int> matches;
    for (int i = 0; i <= n - m; i ++) {
        long long
            part1 = sum_AA[i + m - 1] - (i ? sum_AA[i - 1] : 0),
            part2 = sum_BB,
            part3 = C[i + m - 1].real();
        if (part1 + part2 - 2 * part3 == 0)
            matches.push_back(i);
    }
    
    return matches;
}

int main() {
    string text, pattern;
    cin >> text >> pattern;
    
    vector<int> matches = stringMatchingFFT(text, pattern);

    cout << matches.size() << "\n";
    for (int i : matches) cout << i << " ";
    cout << "\n";

    return 0;
}

带通配符的字符串匹配 #

设有两个字符串 $\text{text}$，$\text{pattern}$ 的长度分别为 $n,m(n\ge m)$，并只包含小写字母 $a\sim z$ 和通配符 $*$（可与任意字符匹配）。问 $\text{pattern}$ 在 $\text{text}$ 中的匹配次数和所有匹配位置。

令通配符 $*$ 映射到 $0$，则

$$\left(A_i-B_j\right)^2\cdot A_i\cdot B_j=\begin{cases} 0,&A_i=B_j\\ \\ >0,&A_i\not=B_j \end{cases}$$

$$\begin{aligned} P_k&=\sum_{i=0}^{m-1}\left(A_{k+i}-B_i\right)^2\cdot A_{k+i}\cdot B_i\\ &=\sum_{i=0}^{m-1}\left(A_{k+i}^2+B_{i}^2-2A_{k+i}B_{i}\right)\cdot A_{k+i}\cdot B_i\\ &=\sum_{i=0}^{m-1}A_{k+i}^3B_i+\sum_{i=0}^{m-1}A_{k+i}B_{i}^3-2\sum_{i=0}^{m-1}A_{k+i}^2B_{i}^2\\ &=\sum_{i=0}^{m-1}A_{k+i}^3\widetilde{B}_{m-i-1}+\sum_{i=0}^{m-1}A_{k+i}\widetilde{B}_{m-i-1}^3-2\sum_{i=0}^{m-1}A_{k+i}^2\widetilde{B}_{m-i-1}^2 \end{aligned}$$

构造多项式

$$A_k(x)=A_0^k+A_1^kx+A_2^kx^2+\cdots+A_{n-1}^kx^{n-1}\\ \widetilde{B}_k(x)=\widetilde{B}_0^k+\widetilde{B}_1^kx+\widetilde{B}_2^kx^2+\cdots+\widetilde{B}_{m-1}^kx^{m-1}$$

  1. $\sum_{i=0}^{m-1}A_{k+i}^3\widetilde{B}_{m-i-1}$ 是 $A_3(x)\cdot \widetilde{B}_1(x)$ 中 $x^{k+m-1}$ 前的系数；
  2. $\sum_{i=0}^{m-1}A_{k+i}\widetilde{B}_{m-i-1}^3$ 是 $A_1(x)\cdot \widetilde{B}_3(x)$ 中 $x^{k+m-1}$ 前的系数；
  3. $\sum_{i=0}^{m-1}A_{k+i}^2\widetilde{B}_{m-i-1}^2$ 是 $A_2(x)\cdot \widetilde{B}_2(x)$ 中 $x^{k+m-1}$ 前的系数。

进行三次 $\text{FFT}$ 即可。

模板 # 

// put FFT algorithm template here

vector<int> stringMatchingFFT(const string& text, string pattern) {
    int n = text.size(), m = pattern.size();
    
    vector<Comp> A1, A2, A3, B1, B2, B3;

    for (auto i : text) {
        int value = i == '*' ? 0 : i - 'a' + 1;
        A1.push_back(value);
        A2.push_back(value * value);
        A3.push_back(value * value * value);
    }

    reverse(pattern.begin(), pattern.end());

    for (auto i : pattern) {
        int value = i == '*' ? 0 : i - 'a' + 1;
        B1.push_back(value);
        B2.push_back(value * value);
        B3.push_back(value * value * value);
    }

    vector<Comp> part1 = multiply(A3, B1);
    vector<Comp> part2 = multiply(A1, B3);
    vector<Comp> part3 = multiply(A2, B2);

    vector<int> matches;
    for (int i = 0; i <= n - m; i ++) {
        if (part1[i + m - 1].real() + part2[i + m - 1].real() - part3[i + m - 1].real() * 2 == 0)
            matches.push_back(i);
    }
    
    return matches;
}

int main() {
    string text, pattern;
    cin >> text >> pattern;
    
    vector<int> matches = stringMatchingFFT(text, pattern);

    cout << matches.size() << "\n";
    for (int i : matches) cout << i << " ";
    cout << "\n";

    return 0;
}