使用FFT进行高效多项式乘法-优快云博客

本文链接：https://blog.youkuaiyun.com/Nickwzk/article/details/56316897

快速傅里叶变换（FFT）是一种用于计算离散傅里叶变换（DFT）和逆变换的算法，将原本O(n^2)的时间复杂度降低到O(nlogn)。本文详细介绍了FFT的工作原理，包括折半引理和如何通过递归或迭代实现。此外，还讨论了FFT在多项式乘法、高精度乘法和带通配符字符串匹配等领域的应用，并给出了相关注意事项，如数据范围、精度和初始化位置等。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

快速傅里叶变换

FFT是用来计算离散傅里叶变换（DFT）及其逆变换（IDFT）的快速算法。
两个n次多项式直接相乘所需的时间为O（ ${n}^{2}$ ），而FFT可以将其复杂度降低为O（ ${n}logn$ ）。
令A(x) = $\sum_{j = 0}^{n-1}a_jx^j$
B(x) = $\sum_{j = 0}^{n-1}b_jx^j$
C(x) = $\sum_{j = 0}^{2n-2}c_jx^j$
其中 $c_j=\sum_{k=0}^{j}a_kb_{j-k}$
一个n次多项式可以由n+1个点值确定，例如两点确定一条直线，三点确定一个二次函数。
那么我们通过DFT将由系数表达式确定的A，B分别转换为点值表达式以后，可以得到C的点值表达式以后，通过IDFT可以得到C的系数表达式。

例如:
A的点值表达式{ $(x_0,y_0),(x_1,y_1)···(x_{2n-1},y_{2n-1})$ }
B的点值表达式{ $(x_0,y_0^{'}),(x_1,y_1^{'})···(x_{2n-1},y_{2n-1}^{'})$ }
那么C的点值表达式为{ $(x_0,y_0y_0^{'}),(x_1,y_1y_1^{'})···(x_{2n-1},y_{2n-1}y_{2n-1}^{'})$ }

通过精心的选点，我们可以将系数表达式转为点值表达式的求值过程，与将点值表达式转为系数表达式的插值过程变为O（ $nlogn$ ）。

实现

折半引理：若n>0且n为偶数，那么n个n次单位复数根的平方的集合就是 $\frac{n}{2}$ 个 $\frac{n}{2}$ 次单位复根的集合。
我们在定义两个次数界为 $\frac{n}{2}$ 的多项式 $A^{[0]}(x)$ 和 $A^{[1]}(x)$ ：
$A^{[0]}(x)$ = $a_0 + a_2x+a_4x^2+···+a_{n-2}x^{n/2-1}$
$A^{[1]}(x)$ = $a_1 + a_2x+a_3x^2+···+a_{n-1}x^{n/2-1}$
其中 $A^{[0]}(x)$ 包括所有下标为偶数的系数， $A^{[1]}(x)$ 包括所有下标为奇数的系数，所以有A(x) = $A^{[0]}(x^2)$ + $xA^{[1]}(x^2)$
所以我们将求A(x)的点值转化为求次数界n/2处 $A^{[0]}(x)$ 和 $A^{[1]}(x)$ 的值。
我们由折半引理可以递归的将以上两个多项式在n/2个n/2次单位复根处的值求出，此时我们需要保证n是2的幂。
然后我们将其递归的过程画出来，以0~7为例，我们可以发现，从左到右，其二进制下标分别为:
000,100,010,110,001,101,011,111 将他们反转得到
000,001,010,011,100,101,110,111
所以我们可以将其转化为迭代的版本！

$DFT_n^{−1}(y)$ 的结果为：
$a_j=\frac{1}{n}\sum_{k=0}^{n-1}y_k\omega_n^{-kj}$
所以我们只需要用 $\omega_n^{-1}$ 代替 $\omega_n$ ，最后将所有结果除以n即可。

下面给出代码：

#include <cstdio>
#include <algorithm>
#include <cmath>
#include <iostream>
using namespace std;

const double pi = acos(-1.0);
const int MAXN = 300003;
struct comp {
    double x, y;
    comp(double _x = 0, double _y = 0):x(_x), y(_y) {}
    comp operator + (const comp& a) const { return comp(x + a.x, y + a.y); }
    comp operator - (const comp& a) const { return comp(x - a.x, y - a.y); }
    comp operator * (const comp& a) const { return comp(x * a.x - y * a.y, x * a.y + y * a.x); }
};

int rev[MAXN], T;
comp Sin[MAXN], tmp;
void fft(comp *a) {
    for(int i = 1; i < T; i++) if(rev[i] > i) swap(a[rev[i]], a[i]);
    for(int i = 2, mid = 1, s = (T >> 1); i <= T; mid = i, i <<= 1, s >>= 1) {
        for(int j = 0; j < T; j += i) {
            for(int k = j, cur = 0; k < j + mid; k++, cur += s) {
                tmp = a[k + mid] * Sin[cur];
                a[k + mid] = a[k] - tmp;
                a[k] = a[k] + tmp;
            }
        }
    }
}

comp A[MAXN];
int n, m;
inline void init() {
    for(T = 1; T <= n + m; T <<= 1);
    for(int i = 1; i < T; i++) {
        if(i & 1) rev[i] = (rev[i >> 1] >> 1) ^ (T >> 1);
        else rev[i] = rev[i >> 1] >> 1;
    }
    Sin[0] = comp(1, 0); Sin[1] = comp(cos(2 * pi / T), sin(2 * pi / T));
    for(int i = 2; i < (T >> 1); i++) Sin[i] = Sin[i - 1] * Sin[1];
}

int main() {
    scanf("%d%d", &n, &m);
    for(int i = 0; i <= n; i++) scanf("%lf", &A[i].x);
    for(int i = 0; i <= m; i++) scanf("%lf", &A[i].y);
    init();
    fft(A);
    for(int i = 0; i < T; i++) A[i] = A[i] * A[i];
    for(int i = 0; i < (T >> 1); i++) Sin[i].y = -Sin[i].y;
    fft(A);
    for(int i = 0; i <= n + m; i++) printf("%d%c", (int)(A[i].y / T / 2 + 0.5), i == n + m ? '\n' : ' ');
    return 0;
}