快速傅里叶变换（FFT）lyh's blog

介绍快速傅里叶变换（Fast Fourier Transformation，FFT），可以用来加速多项式乘法。快速傅里叶变换是多项式以及卷积问题的根本基础，是学习多项式以及卷积的前提。

多项式的点值表示

众所周知，一个多项式可以表示成下面的形式：

$\sum_{i=0}^na_ix^i$

这时，它的系数可以组成一个向量：

$\vec F=(a_0,a_1,\cdots,a_n)$

        这种形式称为多项式的系数表示，通过原有系数向量得到的它们所对应多项式乘积的系数向量，称为两个向量的卷积。易知，如果要求两个多项式的乘积，需要用它们的系数去分别相乘，时间复杂度为$O(n^2)$，是很低效的。。这时，我们希望可以用另一种方式来表示多项式。
        这里有一个定理：任何一个最高次幂为n的多项式，可以由n+1组满足该多项式的互异点对唯一地确定出来。用这些点表示多项式的形式称为点值表示。
        下面是定理的简单证明。
        考虑点对$(x_0,y_0),(x_1,y_1),\cdots,(x_n,y_n)$，那么有：

$\begin{bmatrix} 1&x_0&x_0^2&\cdots&x_0^{n}\\ 1&x_1&x_1^2&\cdots&x_1^{n}\\ \vdots&\vdots&\vdots &\ &\vdots\\ 1&x_n&x_n^2&\cdots&x_n^{n}\\ \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_n \end{bmatrix}=\begin{bmatrix} y_0\\ y_1\\ \vdots\\ y_n \end{bmatrix}$

最左侧的方阵对应的行列式为范德蒙行列式，其值为：

$\prod_{1\leq i<j\leq n}(x_j-x_i)$

当点互异时，x互不同，行列式不为0，该方阵可逆，乘到右边就可以求出所有的系数，即唯一确定了一个多项式。证毕。
注意到，对于两个同次数多项式，给它们赋予一组相同的x值，可以得到下面两个点值表示向量：

$\vec {F1}=(y_0,y_1,\cdots,y_n)\\ \vec {F2}=(y_0^{'},y_1^{'},\cdots,y_n^{'})$

那么它们的乘积的点值表示向量中一定有：

$\vec F=(y_0^{'}y_0,y_1^{'}y_1,\cdots,y_n^{'}y_n)$

也就是说，在点值表示形式下，可以在$O(n)$复杂度下计算两个多项式的乘积，这是我们加速多项式乘法的思想。朴素的利用系数表示得到点值表示的操作称为离散傅里叶变换（DFT），其逆过程称为离散傅里叶逆变换（IDFT）。

复数运算

复数，这是一个常用的数学工具。这里只复习一下它的运算法则。
对于两个复数$x=a+bi$和$y=c+di$，有运算：

$x+y=(a+c)+(b+d)i\\ x-y=(a-c)+(b-d)i\\ xy=(ac-bd)+(ad+bc)i$

在复平面内，复数的加法和减法符合平行四边形定则。如果用极坐标的形式来表示复数，即令$x=(\rho_1,\theta_1),y=(\rho_2,\theta_2)$，则有：

$xy=(\rho_1\rho_2,\theta_1+\theta_2)$

简记为极径相乘，极角相加。

为了下文需要，这里需要知道一个公式：

$e^{ix}=cos(x)+isin(x),x\in \mathbb{R}$

单位根

        DFT要做的是将系数化为点值，如果随意代入几个数当然是可行的，但是效率不尽人意。是否能构造一些特殊的数来更快地求解点值？答案是肯定的，这就是FFT的优化原理。
        为方便起见，下面规定所有多项式的最高次幂为n-1次，且n为2的方幂。在实际应用中，对于不满足条件的多项式，可以从最高次幂开始向上补很多0来满足条件。
        在复平面单位圆上取n个点，使它们均匀分开，这n个点分别为：

$\omega_0=(1,0)\\ \omega_1=(cos(\frac {2\pi} {n}),sin(\frac {2\pi} {n}))\\ \vdots\\ \omega_{n-1}=(cos(\frac {2\pi(n-1)} {n}),sin(\frac {2\pi(n-1)} {n}))$

即有$\omega_i=(cos(\frac {2\pi i} {n}),sin(\frac {2\pi i} {n})=cos(\frac {2\pi i} {n}))+isin(\frac {2\pi i} {n})$。我们把这里的$\omega_1$称为单位根。
根据复数的运算法则，可以知道其余n-1个数均为单位根的幂，即有：

$\omega_1^k=\omega_k$

        $\omega_n$就是$\omega_0$。由于单位根具有这样的性质，可以将其余这n个数均记为单位根的幂次形式，而单位根本身由于是在将单位圆分为n个的基础上得到的，可以记为$\omega_n$。
        这样，这n个数就是$\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}$。
        为什么选取这n个数更好呢？这是由它们的性质决定的，下面介绍三个引理。

消去引理

$\omega_{dn}^{dk}=\omega_n^k$

证明：直接套用公式即可。

$\omega_{dn}^{dk}=(cos(\frac {2\pi} {dn})+isin(\frac {2\pi} {dn}))^{dk}=e^{\frac {2\pi} {dn}dk}=e^{\frac {2\pi k} {n}}=\omega_n^k$

折半引理

如果n为大于0的偶数，那么其对应的$\omega_n^0,\omega_n^0,\cdots,\omega_n^{n-1}$平方的集合就是$\frac {n}{2}$对应的$\omega_{\frac {n}{2}},\omega_{\frac {n}{2}}^2,\cdots,\omega_{\frac {n}{2}}^{\frac {n}{2}}$的集合。
证明：由消去引理可以得到：

$\omega_n^{2k}=\omega_{\frac {n} {2}}^k$

这样就构造了一个对应关系。另外注意到：

$\omega_n^{2(k+\frac {n} {2})}=\omega_n^{2k+n}=\omega_n^{2k}$

由此引理得证，这一条是FFT复杂度的保证。

求和引理

在$n>0$且$0\leq k<n$时，有：

$\sum_{i=0}^{n-1}(\omega_n^k)^i=0$

若$0<k<n$，可知$\omega_n^k≠1$，那么可以应用求和公式：

$\sum_{i=0}^{n-1}(\omega_n^k)^i=\frac {(\omega_n^k)^n-1} {\omega_n^k-1}=\frac {(\omega_n^n)^k-1} {\omega_n^k-1}=0$

$k=0$时，则显然有：

$\sum_{i=0}^{n-1}(\omega_n^k)^i=n$

FFT

下面真正到了FFT的算法部分。考虑一个多项式：

$A(x)=a_0+a_1x+a_2x^2+\cdots +a_{n-1}x^{n-1}$

将其按照奇偶划分得到：

$A(x)=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots+a_{n-1}x^{n-1})$

将这两部分多项式化成与原多项式相同的形式，即使x折半，得到：

$A_0(x)=a_0+a_2x+\cdots+a_{n-2}x^{\frac {n} {2} -1}\\ A_1(x)=a_1+a_3x+\cdots+a_{n-1}x^{\frac {n} {2}-1}$

那么原式可以表示为：

$A(x)=A_0(x^2)+xA_1(x^2)$

现在将n个复数$\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}$代入到$A(x)$中，就相当于将它们的平方代入到$A_0(x)，A_1(x)$中。由折半引理，它们的平方就是$\frac {n} {2}$下的那些复数。
另外注意到，代入$\omega_n^k$和$\omega_n^{k+\frac {n} {2}}$时：

$A(\omega_n^k)=A_0(\omega_n^{2k})+\omega_n^kA_1(\omega_n^k)\\ A(\omega_n^{k+\frac {n} {2}})=A_0(\omega_n^{2k+n})+\omega_n^{k+\frac {n} {2}}A_1(\omega_n^{2k})=A_0(\omega_n^{2k})-\omega_n^{k}A_1(\omega_n^{2k})$

发现只有后面的符号不一样，所以可以只计算一半，另一半直接推知。
这样就把问题转化为等价的子问题，可以递归求解，从而保证了FFT的时间复杂度$O(nlogn)$。

IFFT

下面是快速傅里叶逆变换（IFFT），它可以快速根据FFT得到的点值表达转化为系数表达。
假定$A(x)$的点值表示向量为$\vec F=(p_0,p_1,\cdots,p_{n-1})$，那么可以得到一个多项式：

$B(x)=p_0+p_1x+\cdots+p_{n-1}x^{n-1}$

设该多项式在$\omega_n^0,\omega_n^{-1},\cdots,\omega_n^{-n+1}$下的点值表示为$(c_0,c_1,\cdots,c_{n-1})$，经过证明可以得到：

$c_k=\sum_{i=0}^{n-1}a_i\sum_{j=0}^{n-1}(\omega_{n}^{i-k})^j$

根据求和引理，$i≠k$时，右侧求和式为0，否则为n。因此有：

$c_k=na_k$

即：

$a_k=\frac {c_k} {n}$

这样就可以完成逆变换。求$B(x)$点值时可以套用FFT。

代码模板

下面是朴素的代码模板（不可使用）：

#include<bits/stdc++.h>

#define N (1000005<<2)
using namespace std;
const double Pi = acos(-1.0);

struct Complex {//定义复数
    double a, b;

    Complex operator+(Complex c) { return {a + c.a, b + c.b}; }

    Complex operator*(Complex c) { return {a * c.a - b * c.b, a * c.b + b * c.a}; }

    Complex operator-(Complex c) { return {a - c.a, b - c.b}; }
} a[N], b[N];

inline int read() {
    char e = getchar();
    int s = 0, g = 0;
    while (e < '-')e = getchar();
    if (e == '-')g = 1, e = getchar();
    while (e > '-')s = (s << 1) + (s << 3) + (e & 15), e = getchar();
    return g ? -s : s;
}

void FFT(int l, Complex *c, int type) {//type控制是FFT还是IFFT
    if (l == 1)return;//长度为1，不需变换
    Complex c0[l >> 1], c1[l >> 1];
    for (int i = 0; i < l; i += 2)c0[i >> 1] = c[i], c1[i >> 1] = c[i + 1];//奇偶划分
    FFT(l >> 1, c0, type), FFT(l >> 1, c1, type);//递归
    Complex wn = {cos(2 * Pi / l), type * sin(2 * Pi / l)}, w = {1.0, 0};//构造复数
    for (int i = 0; i < (l >> 1); i++, w = w * wn)c[i] = c0[i] + w * c1[i], c[i + (l >> 1)] = c0[i] - w * c1[i];//求值
}


int main() {
    int n = read(), m = read(), l = 1;
    for (int i = 0; i <= n; i++)a[i].a = read();
    for (int i = 0; i <= m; i++)b[i].a = read();
    while (l <= n + m)l <<= 1;//保证是2的方幂
    FFT(l, a, 1), FFT(l, b, 1);
    for (int i = 0; i <= l; i++)a[i] = a[i] * b[i];//乘积
    FFT(l, a, -1);
    for (int i = 0; i <= n + m; i++)printf("%d ", (int) (a[i].a / l + 0.5));//加上0.5去精度误差
    return 0;
}

上面代码不优，并且可能是由于精度误差，模板题会WA。
一种优化是蝴蝶效应，在计算时保存复数，减少复数运算。

void FFT(int l, Complex *c, int type) {
    if (l == 1)return;
    Complex c0[l >> 1], c1[l >> 1];
    for (int i = 0; i <= l; i += 2)c0[i >> 1] = c[i], c1[i >> 1] = c[i + 1];
    FFT(l >> 1, c0, type), FFT(l >> 1, c1, type);
    Complex wn = {cos(2 * Pi / l), type * sin(2 * Pi / l)}, w = {1.0, 0}, t;
    for (int i = 0; i < (l >> 1); i++, w = w * wn) {
        t = w * c1[i];//进行一步记忆化，降低常数
        c[i] = c0[i] + t, c[i + (l >> 1)] = c0[i] - t;
    }
}

这样仍不是最优的，这是因为递归过程很消耗时间，其实可以迭代进行。注意到序列和原序列下标正好是二进制翻转后的值，可以先通过一些操作来进行迭代计算，避免递归。

void FFT(int l, Complex *c, int type) {
    for (int i = 0; i < l; i++)if (i < tr[i])swap(c[i], c[tr[i]]);//交换
    for (int mid = 1; mid < l; mid <<= 1) {//mid为区间一半的长度
        Complex wn(cos(Pi / mid), type * sin(Pi / mid));//2pi/2mid就是pi/mid
        for (int len = mid << 1, j = 0; j < l; j += len) {//len是区间长度，j为现在所对应的区间首标记
            Complex w(1.0, 0);
            for (int k = 0; k < mid; k++, w = w * wn) {
                Complex x = c[j + k], y = w * c[j + mid + k];//蝴蝶效应
                c[j + k] = x + y, c[j + mid + k] = x - y;
            }
        }
    }
}


int main() {
    int n = read(), m = read(), l = 1;
    for (int i = 0; i <= n; i++)a[i].a = read();
    for (int i = 0; i <= m; i++)b[i].a = read();
    while (l <= n + m)l <<= 1, ++le;
    for (int i = 0; i < l; i++)tr[i] = (tr[i >> 1] >> 1) | ((i & 1) << (le - 1));//构造
    FFT(l, a, 1), FFT(l, b, 1);
    for (int i = 0; i <= l; i++)a[i] = a[i] * b[i];
    FFT(l, a, -1);
    for (int i = 0; i <= n + m; i++)printf("%d ", static_cast<int>(a[i].a / l + 0.5));
    return 0;
}

比较完整的代码模板：

#include<bits/stdc++.h>

#define N (60005<<2)
using namespace std;
const double Pi = acos(-1.0);
int le, tr[N], n, m, l = 1;
long long S, p[N];

struct Complex {
    double a, b;

    Complex(double a = 0, double b = 0) : a(a), b(b) {}

    Complex operator+(Complex c) { return {a + c.a, b + c.b}; }

    Complex operator*(Complex c) { return {a * c.a - b * c.b, a * c.b + b * c.a}; }

    Complex operator-(Complex c) { return {a - c.a, b - c.b}; }
} a[N], b[N];

inline void FFT(int l, Complex *c, int type) {
    for (int i = 0; i < l; i++)if (i < tr[i])swap(c[i], c[tr[i]]);
    for (int mid = 1; mid < l; mid <<= 1) {
        Complex wn(cos(Pi / mid), type * sin(Pi / mid));
        for (int len = mid << 1, j = 0; j < l; j += len) {
            Complex w(1.0, 0);
            for (int k = 0; k < mid; k++, w = w * wn) {
                Complex x = c[j + k], y = w * c[j + mid + k];
                c[j + k] = x + y, c[j + mid + k] = x - y;
            }
        }
    }
}


int main() {
    scanf("%d%d", &n, &m);
    while (l <= n + m)l <<= 1, ++le;
    for (int i = 0; i < l; i++)tr[i] = (tr[i >> 1] >> 1) | ((i & 1) << (le - 1));
    FFT(l, a, 1), FFT(l, b, 1);
    for (int i = 0; i < l; i++)a[i] = a[i] * b[i];
    FFT(l, a, -1);
    for (int i = 0; i < n; i++)a[i].a = static_cast<int>((a[i].a / l + 0.5));
    return 0;
}

多项式乘法的一种优化

在多项式乘法时，需要进行两次DFT以及一次IDFT。myy在2016年国家集训队论文《再探快速傅里叶变换》中提出了一种优化方案，可以只进行一次DFT且一次IDFT就可以完成乘法操作。关于这种做法的数学证明这里不再赘述，详见论文。
对于两个实多项式$A(x)$以及$B(x)$，构造一个多项式如下：

$C(x)=\sum_{i=0}^l(a_i+ib_i)x^i$

这里$i$是虚数单位，$l$是经过处理后的2的方幂。
对这个构造的多项式进行FFT，可以得到一个序列。在这里，经过数学上的推导，我们可以同时得到下面这个多项式经过FFT后的序列：

$D(x)=\sum_{i=0}^l(a_i-ib_i)x^i$

这个多项式的第$k$项系数满足关系：

$d_k=conj(c_{l-k})$

注意，当$k=0$时，认为$l-k=0$，即在模$l$意义下。$conj$表示取共轭复数。
这样经过一次FFT得到了两个多项式的FFT结果，那么这时$A(x)$以及$B(x)$的FFT结果就很容易推知了：

$a_k=\frac {c_k+d_k} {2}\\ b_k=-\frac {c_k-d_k} {2} i$

        注：原论文这里的$b_k$没有取相反数，笔者认为这里有误。
        做IDFT时上式同样成立，但是要注意这个优化只适用于实多项式。
        再进行一步IDFT即可完成乘法运算。

关于精度

        FFT是关于浮点数的运算，出现精度误差在所难免，但是仍然有一些方法可以很好地解决误差问题。这里介绍一个myy在他的论文里提到的一个方法：预处理单位根。
        在上面的代码中，我们的单位根是通过累乘得出的，这样得到的单位根误差较大，容易出现问题。可以在之前预处理所有的单位根，然后在FFT中直接调用，事实证明，这样的操作能够很好地解决误差问题。
        预处理单位根，我们可以在$O(nlogn)$的时间复杂度下完成预处理（忽略求三角函数的时间）：

1
2
3

for (int i = 0, j = 1; j < l; i++, j <<= 1) {
    for (int z = 1; z < j; z++)co[z][i] = cos(z * Pi / j), si[z][i] = sin(z * Pi / j);
}

然后FFT的板子就可以改成：

void FFT(int l, Complex *c, int type) {
    for (int i = 0; i < l; i++)if (i < tr[i])swap(c[i], c[tr[i]]);
    for (int mid = 1, p = 0; mid < l; mid <<= 1, p++) {
        for (int len = mid << 1, j = 0; j < l; j += len) {
            Complex w(1.0, 0);
            for (int k = 0; k < mid; k++, w = Complex(co[k][p], type * si[k][p])) {
                Complex x = c[j + k], y = w * c[j + mid + k];
                c[j + k] = x + y, c[j + mid + k] = x - y;
            }
        }
    }
}

这样的FFT在精度上有了不错的提升，在时间条件允许的情况下，推荐使用此版本的FFT。