Skip to content

快速傅里叶变换:应用

2023-12-22

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

大整数乘法(高精度乘法)

假设有两个整数 a=537,b=721a,b 的每一位可以视为多项式 A(x),B(x) 的系数:

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

计算 C(x)=A(x)B(x)

C(x)=35x4+31x3+60x2+17x+7

提取 C(x) 的系数:

{35,31,60,17,7}

从右往左,逢 101

{3,8,7,1,7,7}

这样就能得到 537×721 的结果是 387177

多项式乘法由 FFT 进行。如果 abn 位整数,则总时间复杂度是 O(nlogn)


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

一个整数可以分解为它的位值的和。例如,整数 123 可以分解为 1×102+2×10+3

a=a2a1a0b=b2b1b0,则

a×b= (a2102+a110+a0)(b2102+b110+b0)= a2b2104+(a2b1+a1b2)103+(a2b0+a1b1+a0b2)102+(a1b0+a0b1)10+a0b0

A(x)=a2x2+a1x+a0B(x)=b2x2+b1x+b0,则

A(x)×B(x)= (a2x2+a1x+a0)(b2x2+b1x+b0)= a2b2x4+(a2b1+a1b2)x3+(a2b0+a1b1+a0b2)x2+(a1b0+a0b1)x+a0b0

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

模板

cpp
// 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;
}

朴素字符串匹配

设有两个字符串 textpattern 的长度分别为 n,m(nm),并只包含小写字母 az。问 patterntext 中出现了几次,以及所有的出现位置。

该问题可使用 FFT 在 O(nlogn) 的复杂度下解决。


约定 str|ab 表示从 str[a]str[b] 的子串。

首先将 textpattern 的所有字母都映射成数字(a1,b2,),得到两个整型数组 AB

textA:{A0,A1,,An1}patternB:{B0,B1,,Bm1}

容易发现

(AiBj)2={0,Ai=Bj>0,AiBj

进而构造

Pk=(AkB0)2+(Ak+1B1)2++(Ak+m1Bm1)2=i=0m1(Ak+iBi)2={0,text|kk+m1=pattern>0,text|kk+m1pattern

因此可以利用 Pk 判断 pattern 是否出现在 text 的第 k 位。

展开 Pk

Pk=i=0m1(Ak+iBi)2=i=0m1(Ak+i2+Bi22Ak+iBi)=i=0m1Ak+i2+i=0m1Bi22i=0m1Ak+iBi

将数组 B 反转得到 B~,有 Bi=B~mi1。代入得

Pk=i=0m1Ak+i2+i=0m1Bi22i=0m1Ak+iB~mi1

分析该式的三项:

  1. i=0m1Ak+i2 可以通过前缀和技巧实现 O(n) 预处理,O(1) 查询;
  2. i=0m1Bi2 是常数,可以直接 O(n) 预处理;
  3. i=0m1Ak+iB~mi1 是卷积式,可以通过 FFT 实现 O(nlogn) 预处理,O(1) 查询。

对于第 3 项,更具体地说:

构造多项式

A(x)=A0+A1x+A2x2++An1xn1B~(x)=B~0+B~1x+B~2x2++B~m1xm1

利用 FFT 可以 O(nlogn) 得出多项式 C(x)=A(x)B~(x) 的系数向量。易知 C(x)xk+m1 前的系数为 i=0m1Ak+iB~mi1,这正是我们想要的第 3 项。

依次计算 Pk(k=0,1,2,,nm),当得到 0 时说明 pattern 出现在 text 的第 k 位。

TIP

用 FFT 实现的字符串匹配,其效率并不是最优的,并且常数很大。

KMP 可以在 O(n+m) 的线性复杂度下解决字符串匹配问题。

模板

cpp
// 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;
}

带通配符的字符串匹配

设有两个字符串 textpattern 的长度分别为 n,m(nm),并只包含小写字母 az 和通配符 (可与任意字符匹配)。问 patterntext 中的匹配次数和所有匹配位置。

令通配符 映射到 0,则

(AiBj)2AiBj={0,Ai=Bj>0,AiBjPk=i=0m1(Ak+iBi)2Ak+iBi=i=0m1(Ak+i2+Bi22Ak+iBi)Ak+iBi=i=0m1Ak+i3Bi+i=0m1Ak+iBi32i=0m1Ak+i2Bi2=i=0m1Ak+i3B~mi1+i=0m1Ak+iB~mi132i=0m1Ak+i2B~mi12

构造多项式

Ak(x)=A0k+A1kx+A2kx2++An1kxn1B~k(x)=B~0k+B~1kx+B~2kx2++B~m1kxm1
  1. i=0m1Ak+i3B~mi1A3(x)B~1(x)xk+m1 前的系数;
  2. i=0m1Ak+iB~mi13A1(x)B~3(x)xk+m1 前的系数;
  3. i=0m1Ak+i2B~mi12A2(x)B~2(x)xk+m1 前的系数。

进行三次 FFT 即可。

模板

cpp
// 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;
}