学习笔记《FFT》

用 FFT 实现多项式乘法

一点前置知识

$e$ 是一个实数，其中 $e^{\pi i}=-1$ 。

复数

指形如 $a+bi$ 的数，其中 $i=\sqrt{-1}$ ， $a,b$ 为实数。

(a+bi)+(c+di)=(a+c)+(b+d)i\\ (a+bi)-(c+di)=(a-c)+(b-d)i\\ (a+bi)\times(c+di)=(ac-bd)+(ad+bc)i

单位根

如果有一个数 $\omega_n$ ，满足 $\omega_n^n=1$ ，称这个数为 $n$ 次单位根。
根据代数基本定理， $\omega_n$ 有且只有 $n$ 个，而且显然各不相同。
如果一个单位根 $x$ ，它的 $1,2,\dots,n$ 次方分别为 $n$ 个单位根，我们称 $x$ 为 $n$ 次本原单位根。
$e^{\frac {2\pi i}n}$ 是一个 $n$ 次本原单位根。为了方便，下文中 $\omega_n = e^{\frac {2i\pi}n}$ 。
单位根的一些性质（ $0 \le k < n$ ）

$\omega_n^0=\omega_n^n=1$
$\omega_n^k=\omega_n^{k+n}=-\omega_n^{k+\frac n2}=\omega_{2n}^{2k}$
$\omega_n^k=\cos(\frac{2k\pi}{n})+\sin(\frac{2k\pi}{n})i$
$\omega_n^{-k}=\omega_n^{n-k}=\cos(\frac{2k\pi}{n})-\sin(\frac{2k\pi}{n})i$

多项式的表示方式

系数表示法

对于一个 $n$ 次多项式可以表示为 $A(x)=\sum_{i=0}^{n}a_ix^i$ 。
注意， $n$ 次多项式有 $n+1$ 项。

点值表示法

对于一个集合 $\{(x_0,A(x_0)),(x_1,A(x_1)),\dots, (x_n,A(x_n))\}$ ，可以确定一个 $n$ 次多项式 $A(x)$ 。
点值表示法虽然不适合给人看，但是可以加快乘法。
设 $n$ 次多项式 $G(x)$ 和 $m$ 次多项式 $H(x)$ 。
我们设 $F(x)=G(x)\cdot H(x)$ ，则 $(a, F(a)) = (a, G(a)\cdot H(a))$ 。
我们取 $n+m+1$ 个 $a$ 就能得出 $F(x)$ 。
而且 $a$ 可以随便选，那么我们选有一定关系的 $a$ ，就能加快运算了。

快速傅里叶变换

快速傅里叶变换（FFT）的大致流程：

将两个相乘的多项式转化成点值表示法。
将点值相乘得到答案的点值表示法。
还原出答案。

这里，我们称操作 $1$ 为 DFT，操作 $3$ 为 IDFT。
我们用分治的思想，把大问题变为小问题，将时间复杂度降到 $O(n\log n)$ 。

DFT

我们有两种分割方式，DIT 和 DIF。

DIT（Decimation in Time，按时域抽取）

令 $f(x)=\sum_{i=0}^{n-1}a_ix^i$ ，其中 $n$ 为 $2$ 的整数次幂（原因等下会讲）。
我们根据 $i$ 的奇偶性进行分组。

\begin{aligned} f(x)&=\sum_{i=0}^{n-1}a_ix^i\\ &=(\sum_{i=0}^{\frac n2-1}a_{2i}x^{2i})+(\sum_{i=0}^{\frac n2-1}a_{2i+1}x^{2i+1})\\ &=(\sum_{i=0}^{\frac n2-1}a_{2i}x^{2i})+x(\sum_{i=0}^{\frac n2-1}a_{2i+1}x^{2i})\\ \end{aligned}

令 $g(x) =\sum_{i=0}^{\frac n2-1}a_{2i}x^{i}, h(x)=\sum_{i=0}^{\frac n2-1}a_{2i+1}x^{i}$ 则 $f(x)=g(x^2)+x\times h(x^2)$

将 $\omega_n^k$ 和 $\omega_n^{k+\frac n2}$ 分别代入，由单位根的性质得：

\begin{aligned} f(\omega_n^k)=g(\omega_n^{2k})+\omega_n^kh(\omega_n^{2k}) &=g(\omega_{\frac n2}^{k})+\omega_n^kh(\omega_{\frac n2}^k)\\ f(\omega_n^{k+\frac n2})=g(\omega_n^{2k+n})+\omega_n^{k+\frac n2}h(\omega_n^{2k+n}) &=g(\omega_{\frac n2}^{k})-\omega_n^kh(\omega_{\frac n2}^k) \end{aligned}

这样只需要代入一半的单位根的幂就可以了。对于 $h(x)$ 和 $g(x)$ 显然可以继续递归。
即我们每次将 $a_i$ 用 $i$ 的奇偶性分组， $\omega_n^k$ 前后分成两段。
因为每次都需要严格的将多项式分成相等长度的两部分，所以 $n$ 必须为 $2$ 的整次幂。

DIF（Decimation in Frequency，按频域抽取）

这次，我们对将 $a_i$ 分成前后两段， $\omega_n^k$ 按奇偶分组。

\begin{aligned} f(x)&=\sum_{i=0}^{n-1}a_ix^i\\ &=(\sum_{i=0}^{\frac n2-1}a_ix^i)+\sum_{i=0}^{\frac n2-1}a_{i+\frac n2}x^{i+\frac n2}\\ &=(\sum_{i=0}^{\frac n2-1}a_ix^i)+x^{\frac n2}\sum_{i=0}^{\frac n2-1}a_{i+\frac n2}x^i \end{aligned}

设 $2\mid k$ ，将 $\omega_n^k$ 和 $\omega_n^{k+1}$ 带入。

\begin{aligned} f(\omega_n^k)&=(\sum_{i=0}^{\frac n2-1}a_i\omega_n^{ik})+\omega_n^{\frac {kn}2}\sum_{i=0}^{\frac n2-1}a_{i+\frac n2}\omega_n^{ik}\\ &=(\sum_{i=0}^{\frac n2-1}a_i\omega_n^{ik})+\sum_{i=0}^{\frac n2-1}a_{i+\frac n2}\omega_n^{ik}\\ &=\sum_{i=0}^{\frac n2-1}(a_i+a_{i+\frac n2})\omega_n^{ik}\\ f(\omega_n^{k+1})&=(\sum_{i=0}^{\frac n2-1}a_i\omega_n^{i(k+1)})+\omega_n^{\frac {(k+1)n}2}\sum_{i=0}^{\frac n2-1}a_{i+\frac n2}\omega_n^{i(k+1)}\\ &=(\sum_{i=0}^{\frac n2-1}a_i\omega_n^{i(k+1)})-\sum_{i=0}^{\frac n2-1}a_{i+\frac n2}\omega_n^{i(k+1)}\\ &=\sum_{i=0}^{\frac n2-1}\omega_n^i(a_i-a_{i+\frac n2})\omega_n^{ik}\\ \end{aligned}

我们同样只需要带入一半即可。

IDFT

考虑怎么变回去。
设 $G(x)=\sum_{i=0}^{n-1}F(\omega_n^i)x^i$ ，其中 $F(x)$ 为最终的答案。
结论一：对 $G(x)$ 做 DFT，但用 $\omega_n^{-k}$ 代替 $\omega_n^k$ ，结果的每一项除以 $n$ 后为 $F(x)$ 。
结论二：对 $G(x)$ 做 DFT，然后将后 $n-1$ 项翻转，结果的每一项除以 $n$ 后为 $F(x)$ 。

我代码中均使用结论一。

结论一证明

\begin{aligned}G(\omega_n^{-k})&=\sum_{i=0}^{n-1}F(\omega_n^i)\omega_n^{-ki}\\&=\sum_{i=0}^{n-1}\omega_n^{-ki}\sum_{j=0}^{n-1}a_j\omega_n^{ij}\\&=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\omega_n^{ij-ki}\\&=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\end{aligned}

记 $S(\omega_n^a)=\sum_{i=0}^{n-1}(\omega_n^a)^i$ 。

当 $a\equiv 0 \pmod n$ 时， $S(\omega_n^a)=n$ 。
否则，我们错位相减

\begin{aligned}S(\omega_n^a)&=\sum_{i=0}^{n-1}(\omega_n^a)^i\\\omega_n^aS(\omega_n^a)&=\sum_{i=1}^{n}(\omega_n^a)^i\\S(\omega_n^a)&=\frac{(\omega_n^a)^n-(\omega_n^a)^0}{\omega_n^a-1}=0\\\end{aligned}

也就是说

S(\omega_n^a)=\begin{cases}n,&{n\mid a}\\0,&{n\nmid a}\end{cases}

那么代回原式
$G(\omega_n^{-k})=\sum_{j=0}^{n-1}a_jS(\omega_n^{j-k})=na_k$

综上所述，将 $\omega_n^k$ 换成 $\omega_n^{-k}$ 对 $G(x)$ 跑一遍 DFT，然后除以 $n$ 即可。

结论二证明

\begin{aligned}G(\omega_n^k)&=\sum_{i=0}^{n-1}F(\omega_n^i)\omega_n^{ki}\\&=\sum_{i=0}^{n-1}\omega_n^{ki}\sum_{j=0}^{n-1}a_j\omega_n^{ij}\\&=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\omega_n^{ij+ki}\\&=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j+k})^i\end{aligned}

记 $S(\omega_n^a)=\sum_{i=0}^{n-1}(\omega_n^a)^i$ 。

同上得知

S(\omega_n^a)=\begin{cases}n,&{n\mid a}\\0,&{n\nmid a}\end{cases}

带回原式得：

G(\omega_n^k)=\sum_{j=0}^{n-1}a_jS(\omega_n^{j+k})=\begin{cases}na_0,&{k=0}\\na_{n-k},&{k\neq0}\end{cases}

所以，直接对 $G(x)$ 跑一遍 DFT，然后将后 $n-1$ 位翻转，最后除以 $n$ 即可。

实现

我们用洛谷P3803 为示例。
先给出一个用递归的简单实现。
我们先以 DIT 为例。

代码

#include<bits/stdc++.h>
using namespace std;
const double Pi = acos(-1), Pi2 = 2 * Pi;
const int N = 4e6 + 10;
struct cp{
    double x, y;
    cp(double x = 0, double y = 0):x(x), y(y){}
    cp operator+(const cp&i)const{return cp(x + i.x, y + i.y);}
    cp operator-(const cp&i)const{return cp(x - i.x, y - i.y);}
    cp operator*(const cp&i)const{return cp(x * i.x - y * i.y, y * i.x + x * i.y);}
    cp&operator*=(const cp&i){return *this = *this * i;}
}tmp[N], a[N], b[N];//实现复数类，x 代表实部，y 代表虚部。
void FFT(cp*f, int n, int on){//on = 1 表示 dft，on = -1 为 idft
    if(n == 1)return;
    cp *fl = f, *fr = f + n / 2;
    for(int i = 0; i < n; i++)tmp[i] = f[i];
    for(int i = 0; i < n / 2; i++)
        fl[i] = tmp[i * 2], fr[i] = tmp[i * 2 + 1];
    FFT(fl, n / 2, on), FFT(fr, n / 2, on);
    cp wn(cos(Pi2 / n), on * sin(Pi2 / n)), w(1, 0);
    for (int i = 0; i < n / 2; i++)
        tmp[i] = fl[i] + w * fr[i],
        tmp[i + n / 2] = fl[i] - w * fr[i], w *= wn;
    for(int i = 0; i < n; i++)f[i] = tmp[i];
}
int main(){
    cin.tie(0)->sync_with_stdio(0);
    int n, m;
    cin >> n >> m;
    for(int i = 0; i <= n; i++)cin >> a[i].x;
    for(int i = 0; i <= m; i++)cin >> b[i].x;
    int len = 1;
    while(len <= n + m)len <<= 1;
    FFT(a, len, 1), FFT(b, len, 1);
    for(int i = 0; i < len; i++)a[i] *= b[i];
    FFT(a, len, -1);
    for(int i = 0; i <= n + m; i++)cout << (long long)(a[i].x / len + 0.49) << " ";
}

每次我们都要反复复制，常数太大了。那么提前把每个元素放在最后的位置上，就好了。
我们打表可以发现，把每个位置的二进制翻转得到的数就是最后的下标。

1
2
3

for(int i = 0; i < n; i++)
   rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? n >> 1 : 0);
//递推预处理

另外 $fl_k$ 复制在 $f_k$ ， $fr_k$ 复制在 $f_{k+\frac n2}$ ，没有重叠，那么可以直接共用。
那么可以变成

递归 DIT

#include<bits/stdc++.h>
using namespace std;
const double Pi = acos(-1), Pi2 = 2 * Pi;
const int N = 4e6 + 10;
struct cp{
    double x, y;
    cp(double x = 0, double y = 0):x(x), y(y){}
    cp operator+(const cp&i)const{return cp(x + i.x, y + i.y);}
    cp operator-(const cp&i)const{return cp(x - i.x, y - i.y);}
    cp operator*(const cp&i)const{return cp(x * i.x - y * i.y, y * i.x + x * i.y);}
    cp&operator*=(const cp&i){return *this = *this * i;}
}a[N], b[N];//实现复数类，x 代表实部，y 代表虚部。
void FFT(cp*f, int n, int on){//on = 1 表示 dft，on = -1 为 idft
    if(n == 1)return;
    int half = n >> 1;
    FFT(f, half, on), FFT(f + half, half, on);
    cp wn(cos(Pi2 / n), on * sin(Pi2 / n)), w(1, 0), z;
    for (int i = 0; i < half; i++)
        z = w * f[i + half], f[i + half] = f[i] - z, f[i] = f[i] + z, w *= wn;
}
int rev[N];
int main(){
    cin.tie(0)->sync_with_stdio(0);
    int n, m;
    cin >> n >> m;
    for(int i = 0; i <= n; i++)cin >> a[i].x;
    for(int i = 0; i <= m; i++)cin >> b[i].x;
    int len = 1;
    while(len <= n + m)len <<= 1;
    for(int i = 0; i < len; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? len >> 1 : 0);
    for(int i = 0; i < len; i++)if(i < rev[i])swap(a[i], a[rev[i]]);
    for(int i = 0; i < len; i++)if(i < rev[i])swap(b[i], b[rev[i]]);
    FFT(a, len, 1), FFT(b, len, 1);
    for(int i = 0; i < len; i++)a[i] *= b[i];
    for(int i = 0; i < len; i++)if(i < rev[i])swap(a[i], a[rev[i]]);
    FFT(a, len, -1);
    for(int i = 0; i <= n + m; i++)cout << (long long)(a[i].x / len + 0.49) << " ";
}

同样的，我们可以改写为循环迭代。

迭代 DIT

#include<bits/stdc++.h>
using namespace std;
const double Pi = acos(-1), Pi2 = 2 * Pi;
const int N = 4e6 + 10;
struct cp{
    double x, y;
    cp(double x = 0, double y = 0):x(x), y(y){}
    cp operator+(const cp&i)const{return cp(x + i.x, y + i.y);}
    cp operator-(const cp&i)const{return cp(x - i.x, y - i.y);}
    cp operator*(const cp&i)const{return cp(x * i.x - y * i.y, y * i.x + x * i.y);}
    cp&operator*=(const cp&i){return *this = *this * i;}
}a[N], b[N];//实现复数类，x 代表实部，y 代表虚部。
int rev[N];
void FFT(cp*f, int n, int on){//on = 1 表示 dft，on = -1 为 idft
    for(int i = 0; i < n; i++)if(i < rev[i])swap(f[i], f[rev[i]]);
    for(int len, p = 2; len = p >> 1, p <= n; p <<= 1){//枚举区间长度 
        cp wn = cp(cos(Pi2 / p), on * sin(Pi2 / p));
        for(int l = 0; l < n; l += p){//枚举区间左端点 
            cp w = cp(1, 0), z;
            for(int i = l; i < l + len; i++)
                z = f[i + len] * w, f[i + len] = f[i] - z, f[i] = f[i] + z, w *= wn;
        }
    }
}
int main(){
    cin.tie(0)->sync_with_stdio(0);
    int n, m;
    cin >> n >> m;
    for(int i = 0; i <= n; i++)cin >> a[i].x;
    for(int i = 0; i <= m; i++)cin >> b[i].x;
    int len = 1;
    while(len <= n + m)len <<= 1;
    for(int i = 0; i < len; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? len >> 1 : 0);
    FFT(a, len, 1), FFT(b, len, 1);
    for(int i = 0; i < len; i++)a[i] *= b[i];
    FFT(a, len, -1);
    for(int i = 0; i <= n + m; i++)cout << (long long)(a[i].x / len + 0.49) << " ";
}

由于 DIF 是按 $k$ 的奇偶性分组，所以我们要给答案做 rev。

递归 DIF

#include<bits/stdc++.h>
using namespace std;
const double Pi = acos(-1), Pi2 = 2 * Pi;
const int N = 4e6 + 10;
struct cp{
    double x, y;
    cp(double x = 0, double y = 0):x(x), y(y){}
    cp operator+(const cp&i)const{return cp(x + i.x, y + i.y);}
    cp operator-(const cp&i)const{return cp(x - i.x, y - i.y);}
    cp operator*(const cp&i)const{return cp(x * i.x - y * i.y, y * i.x + x * i.y);}
    cp&operator*=(const cp&i){return *this = *this * i;}
}a[N], b[N];//实现复数类，x 代表实部，y 代表虚部。
void FFT(cp*f, int n, int on){//on = 1 表示 dft，on = -1 为 idft
    if(n == 1)return;
    int half = n >> 1;
    cp wn(cos(Pi2 / n), on * sin(Pi2 / n)), w(1, 0), x, y;
    for (int i = 0; i < half; w *= wn, i++)
        x = f[i], y = f[i + half], f[i] = x + y, f[i + half] = (x - y) * w;
    FFT(f, half, on), FFT(f + half, half, on);
}
int rev[N];
int main(){
    cin.tie(0)->sync_with_stdio(0);
    int n, m;
    cin >> n >> m;
    for(int i = 0; i <= n; i++)cin >> a[i].x;
    for(int i = 0; i <= m; i++)cin >> b[i].x;
    int len = 1;
    while(len <= n + m)len <<= 1;
    for(int i = 0; i < len; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? len >> 1 : 0);
    FFT(a, len, 1), FFT(b, len, 1);
    for(int i = 0; i < len; i++)if(i < rev[i])swap(a[i], a[rev[i]]);
    for(int i = 0; i < len; i++)if(i < rev[i])swap(b[i], b[rev[i]]);
    for(int i = 0; i < len; i++)a[i] *= b[i];
    FFT(a, len, -1);
    for(int i = 0; i < len; i++)if(i < rev[i])swap(a[i], a[rev[i]]);
    for(int i = 0; i <= n + m; i++)cout << (long long)(a[i].x / len + 0.49) << " ";
}

迭代 DIF

#include<bits/stdc++.h>
using namespace std;
const double Pi = acos(-1), Pi2 = 2 * Pi;
const int N = 4e6 + 10;
struct cp{
    double x, y;
    cp(double x = 0, double y = 0):x(x), y(y){}
    cp operator+(const cp&i)const{return cp(x + i.x, y + i.y);}
    cp operator-(const cp&i)const{return cp(x - i.x, y - i.y);}
    cp operator*(const cp&i)const{return cp(x * i.x - y * i.y, y * i.x + x * i.y);}
    cp&operator*=(const cp&i){return *this = *this * i;}
}a[N], b[N];//实现复数类，x 代表实部，y 代表虚部。
int rev[N];
void FFT(cp*f, int n, int on){//on = 1 表示 dft，on = -1 为 idft
    for(int len, p = n; len = p >> 1, p > 1; p >>= 1){//从大到小枚举
        cp wn = cp(cos(Pi2 / p), on * sin(Pi2 / p));
        for(int l = 0; l < n; l += p){//枚举区间左端点 
            cp w = cp(1, 0), x, y;
            for(int i = l; i < l + len; w *= wn, i++)
                x = f[i], y = f[i + len], f[i] = x + y, f[i + len] = (x - y) * w;
        }
    }
    for(int i = 0; i < n; i++)if(i < rev[i])swap(f[i], f[rev[i]]);
}
int main(){
    cin.tie(0)->sync_with_stdio(0);
    int n, m;
    cin >> n >> m;
    for(int i = 0; i <= n; i++)cin >> a[i].x;
    for(int i = 0; i <= m; i++)cin >> b[i].x;
    int len = 1;
    while(len <= n + m)len <<= 1;
    for(int i = 0; i < len; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? len >> 1 : 0);
    FFT(a, len, 1), FFT(b, len, 1);
    for(int i = 0; i < len; i++)a[i] *= b[i];
    FFT(a, len, -1);
    for(int i = 0; i <= n + m; i++)cout << (long long)(a[i].x / len + 0.49) << " ";
}

优化

简单优化

众所周知 $(a+bi)^2 = a^2-b^2+2abi$ ，故我们可以将两个多项式分别放在同一多项式的实部和虚部，FFT 后再计算结果的平方，虚部的一半就是答案。

DIF 需要最后 rev，DIT 需要开始时 rev，我们用 DIF 做 DFT，DIT 做 IDFT。就可以避免用 rev。

递归的常数太大了，我们可以用模版递归，像这样。

模版递归

#include<bits/stdc++.h>
using namespace std;
const double Pi = acos(-1), Pi2 = 2 * Pi;
const int N = 4e6 + 10;
struct cp{
    double x, y;
    cp(double x = 0, double y = 0):x(x), y(y){}
    cp operator+(const cp&i)const{return cp(x + i.x, y + i.y);}
    cp operator-(const cp&i)const{return cp(x - i.x, y - i.y);}
    cp operator*(const cp&i)const{return cp(x * i.x - y * i.y, y * i.x + x * i.y);}
    cp&operator*=(const cp&i){return *this = *this * i;}
}a[N], b[N];//实现复数类，x 代表实部，y 代表虚部。
template<const int n>void DFT(cp*f){//要分开写常数小
    const int half = n >> 1;
    DFT<half>(f), DFT<half>(f + half);
    cp wn(cos(Pi2 / n), sin(Pi2 / n)), w(1, 0), z;
    for (int i = 0; i < half; i++)
        z = w * f[i + half], f[i + half] = f[i] - z, f[i] = f[i] + z, w *= wn;
}
template<>void DFT<1>(cp*f){}template<>void DFT<0>(cp*f){}//特化
#define Case(x, fu) case x: fu<x>(f);break;
#define Runfft(x) switch(n){\
Case(1<<1,x)Case(1<<2,x)Case(1<<3,x)Case(1<<4,x)\
Case(1<<5,x)Case(1<<6,x)Case(1<<7,x)Case(1<<8,x)\
Case(1<<9,x)Case(1<<10,x)Case(1<<11,x)Case(1<<12,x)\
Case(1<<13,x)Case(1<<14,x)Case(1<<15,x)Case(1<<16,x)\
Case(1<<17,x)Case(1<<18,x)Case(1<<19,x)Case(1<<20,x)Case(1<<21,x)}
int rev[N];
void rundft(cp *f,const int& n){
    for(int i = 0; i < n; i++)if(i < rev[i])swap(f[i], f[rev[i]]);
	Runfft(DFT);
}
template<const int n>void IDFT(cp*f){
    const int half = n >> 1;
    IDFT<half>(f), IDFT<half>(f + half);
    cp wn(cos(Pi2 / n), -sin(Pi2 / n)), w(1, 0), z;
    for (int i = 0; i < half; i++)
        z = w * f[i + half], f[i + half] = f[i] - z, f[i] = f[i] + z, w *= wn;
}
template<>void IDFT<1>(cp*f){}template<>void IDFT<0>(cp*f){}//特化
void runidft(cp *f,const int& n){
    for(int i = 0; i < n; i++)if(i < rev[i])swap(f[i], f[rev[i]]);
	Runfft(IDFT);
}
int main(){
    cin.tie(0)->sync_with_stdio(0);
    int n, m;
    cin >> n >> m;
    for(int i = 0; i <= n; i++)cin >> a[i].x;
    for(int i = 0; i <= m; i++)cin >> b[i].x;
    int len = 1;
    while(len <= n + m)len <<= 1;
    for(int i = 0; i < len; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? len >> 1 : 0);
    rundft(a, len), rundft(b, len);
    for(int i = 0; i < len; i++)a[i] *= b[i];
    runidft(a, len);
    for(int i = 0; i <= n + m; i++)cout << (long long)(a[i].x / len + 0.49) << " ";
}

分裂基

我们称每次一分为二为基2，相应的我们也有基4，基8。
这里只介绍 DIT 的做法，但两种代码都会给出，相信读者能想明白 DIF 的做法。
分裂基每次把基2的第 $2$ 项用同样办法拆开，变成：

h(x) = h_1(x^2) + xh_2(x^2)\\ f(x) = g(x^2) + xh_1(x^4) + x^3h_2(x^4)

我们分别将 $\omega_n^k,\omega_n^{k+\frac n4},\omega_n^{k+\frac n2},\omega_n^{k+\frac{3n}4}$ 带入。

\begin{aligned} f(\omega_n^k)&=g(\omega_n^{2k})+\omega_n^kh_1(\omega_n^{4k})+\omega_n^3kh_2(\omega_n^{4k})\\ f(\omega_n^{k+\frac n4})&=g(\omega_n^{2k+\frac n2})-(\omega_n^kh_1(\omega_n^{4k})-\omega_n^{3k}h_2(\omega_n^{4k}))i\\ f(\omega_n^{k+\frac n2})&=g(\omega_n^{2k})-\omega_n^kh_1(\omega_n^{4k})-\omega_n^{3k}h2(\omega_n^{4k})\\ f(\omega_n^{k+\frac {3n}4})&=g(\omega_n^{2k+\frac n2})+(\omega_n^kh_1(\omega_n^{4k})-\omega_n^{3k}h_2(\omega_n^{4k}))i\\ \end{aligned}

分裂基适用于序列长度为 $2^n$ 的 FFT，并且是运算次数比基8更少，并且更灵活。
运算次数的比较见 ooura 的博客。
由于迭代版本过于复杂，所以只给出递归版本。

代码

#include<bits/stdc++.h>
using namespace std;
#define Case(x, fu) case x: fu<x>(f);break;
#define Runfft(x) switch(n){\
Case(1<<1,x)Case(1<<2,x)Case(1<<3,x)Case(1<<4,x)\
Case(1<<5,x)Case(1<<6,x)Case(1<<7,x)Case(1<<8,x)\
Case(1<<9,x)Case(1<<10,x)Case(1<<11,x)Case(1<<12,x)\
Case(1<<13,x)Case(1<<14,x)Case(1<<15,x)Case(1<<16,x)\
Case(1<<17,x)Case(1<<18,x)Case(1<<19,x)Case(1<<20,x)Case(1<<21,x)}
const double Pi = acos(-1), Pi2 = 2 * Pi;
const int N = 4e6 + 10;

struct cp{
    double x, y;
    cp(double x = 0, double y = 0):x(x), y(y){}
    cp operator+(const cp&i)const{return cp(x + i.x, y + i.y);}
    cp operator-(const cp&i)const{return cp(x - i.x, y - i.y);}
    cp operator*(const cp&i)const{return cp(x * i.x - y * i.y, y * i.x + x * i.y);}
    cp&operator+=(const cp&i){return this->x += i.x, this->y += i.y, *this;}
    cp&operator-=(const cp&i){return this->x -= i.x, this->y -= i.y, *this;}
    cp&operator*=(const cp&i){return *this = *this * i;}
}a[N];//实现复数类，x 代表实部，y 代表虚部。

template<const int n>void DFT(cp*f){//DIF
    const int half = n >> 1, quar = n >> 2;
	cp w(1, 0), wn(cos(Pi2 / n), sin(Pi2 / n)), w3, x, y;
	cp *a1 = &f[0], *a2 = a1 + quar, *a3 = a2 + quar, *a4 = a3 + quar;
	for(int i = 0; i < quar; i++){
		w3 = w * w * w, x = *a1 - *a3, y = *a2 - *a4, y = cp(y.y, -y.x);
		*a1 += *a3, *a2 += *a4, *a3 = (x - y) * w, *a4 = (x + y) * w3;
		a1++, a2++, a3++, a4++, w *= wn; 
	}
	DFT<half>(f), DFT<quar>(f + half), DFT<quar>(f + half + quar);
}
template<>void DFT<2>(cp*f){cp x = f[0], y = f[1];f[0] = x + y, f[1] = x - y;}
template<>void DFT<1>(cp*f){}template<>void DFT<0>(cp*f){}//特化
void rundft(cp *f,const int& n){Runfft(DFT);}

template<const int n>void IDFT(cp*f){//DIT
    const int half = n >> 1, quar = n >> 2;
	IDFT<half>(f), IDFT<quar>(f + half), IDFT<quar>(f + half + quar);
    cp wn(cos(Pi2 / n), -sin(Pi2 / n)), w(1, 0), w3, tmp1, tmp2, x, y;
    cp *a1 = &f[0], *a2 = a1 + quar, *a3 = a2 + quar, *a4 = a3 + quar;
    for (int i = 0; i < quar; i++){
    	w3 = w * w * w, tmp1 = w * *a3, tmp2 = w3 * *a4;
    	x = tmp1 + tmp2, y = tmp1 - tmp2, y = cp(y.y, -y.x);
    	*a4 = *a2 - y, *a3 = *a1 - x, *a2 += y, *a1 += x;
    	a1++, a2++, a3++, a4++, w *= wn;
	}
}
template<>void IDFT<2>(cp*f){cp x = f[0], y = f[1];f[0] = x + y, f[1] = x - y;}
template<>void IDFT<1>(cp*f){}template<>void IDFT<0>(cp*f){}//特化 
void runidft(cp *f,const int& n){Runfft(IDFT);}

int main(){
    cin.tie(0)->sync_with_stdio(0);
    int n, m;
    cin >> n >> m;
    for(int i = 0; i <= n; i++)cin >> a[i].x;
    for(int i = 0; i <= m; i++)cin >> a[i].y;
    int len = 1;
    while(len <= n + m)len <<= 1;
    rundft(a, len);
    for(int i = 0; i < len; i++)a[i] *= a[i];
    runidft(a, len);
    double inv = 0.5 / len;
    for(int i = 0; i <= n + m; i++)cout << (long long)(a[i].y * inv + 0.49) << " ";
}

参考资料

https://oi-wiki.org/math/poly/fft/
https://www.bilibili.com/opus/785022478912061446
https://www.luogu.com.cn/article/rj58c2eb
https://charleswu.site/archives/3065