算法学习FFT系列（1）：初习快速傅里叶变换

算法学习FFT系列（1）：初习快速傅里叶变换

引入

这个坑已经在我脑海里占了很久了，但是一直没有水平写，今天尝试着写写看FFT的算法学习。
FFT在OI中最大的作用是加速卷积。理论上背板子是没毛病的，但是仍然遇到了一些考定义的毒瘤题，所以还是理解比较好。

多项式乘法

定义

多项式$A(x)=\sum\limits_i^n a_i*x^i,B(x)=\sum\limits_i^nb_i*x^i$
多项式乘法就是

$C(x)=\sum\limits_i^{2n-1} (\sum\limits_j^i a_j*b_{i-j})x^i$

显然，多项式乘法需要$O(n^2)$的复杂度。

表示方法

普通的多项式表示方法把n阶多项式$A(x)$表示为向量$A$
这里需要引入一种全新的表示方法——点值表示法。
也就是$A(x)$表示为$A(A(x_0),A(x_1),A(x_2) \cdots A(x_n))$
n阶多项式和n个互不相同的点值表示一一对应。
这样的表示方法的优点是什么？考虑多项式乘法的过程。

$C(x)=\sum\limits_i^{2n-1} (\sum\limits_j^i a_j*b_{i-j})x^i$
$=\sum\limits_i^{2n-1} \sum\limits_j^i (a_j*x^j)*(b_{i-j}*x^{i-j})$
$=\sum\limits_j^{2n-1} (a_j*x^j)*\sum\limits_i^{2n-1} (b_i*x^i)$

转化成点值表示的两个式子，不难发现，其乘法复杂度是$O(n)$的

FFT的总路线：系数表达式->点值表达式->乘法->点值表达式->系数表达式

插值（*）

这个东西是顺便一提，FFT中系数表达式<->点值表达式过程并不是FFT的专属，但是确是FFT的关键。这一过程被称之为插值。证明插值的唯一性（也就是n阶多项式和n个互不相同的点值表示一一对应。）需要通过范德蒙德矩阵的可逆性。
$\begin{bmatrix} 1&x_0&x_0^2&\dots & x_0^n \\ 1&x_1&x_1^2&\cdots & x_1^n \\ \dots &\dots &\dots &\dots & \dots \\ 1&x_n&x_n^2&\cdots & x_n^n \end{bmatrix}\begin{bmatrix} a_0 \\ a_1 \\ \dots \\ a_n \end{bmatrix}=\begin{bmatrix} y_0 \\ y_1 \\ \dots \\ y_n \end{bmatrix}$
左边的矩阵表示为$V(x_0,x_1 \dots x_n)$ 就是范德蒙德矩阵
证明不会，自行百度。
插值相关还有这些东西，可以看看算法学习：拉格朗日插值
既然是点值表达式，最重要的就是带入什么点值，这里我们引入一个新的概念——单位复根。

单位复根

$n次单位复根$的表达式是$\omega^n=1$，更形象地，我们可以通过复数运算的几何意义(幅角相加, 模相乘)得到下面两张图

我们可以得到n次单位复根有n个，均匀分布在复平面上半径为1的圆上。

欧拉公式与单位复根表示

$e^{ix}=cosx+isinx$
证明的方法是泰勒展开。
带入$x=2\pi k$得到
$e^{2\pi k i}=cos(2 \pi k)+isin(2 \pi k)=1=\omega^n$
于是$\omega_n^k=e^{\frac{2\pi k i}{n}},k=0,1\dots n$
其实$\omega_n^k$构成了一个乘法群。。。不加以赘述。
有了表达式，我们就可以挖掘$\omega_n^k$的性质了。

各种定理

消去引理

$\omega_{dn}^{dk}=\omega_n^k$
证明：$\omega_{dn}^{dk}=e^{\frac{2\pi dk i}{dn}}=e^{\frac{2\pi k i}{n}}=\omega_n^k$

折半引理

如果$n>0$且n为偶数，那么n个n次单位复根的平方的集合就是$\frac{n}{2}$个$\frac{n}{2}$次单位复根的集合
证明：其实就是证明两个东西(1) $(\omega_n^k)^2=\omega_{\frac{n}{2}}^k$ (2) ) $(\omega_n^{k+\frac{n}{2}})^2=(\omega_{n}^k)^2$
$(1)(\omega_n^k)^2=\omega_n^{2k}=\omega_{\frac{n}{2}}^k$
$(2)(\omega_n^{k+\frac{n}{2}})^2=\omega_n^{2k+n}=\omega_{n}^k\omega_n^n=(\omega_{n}^k)^2$

求和引理

$\sum\limits_j^{n-1} (\omega_n^k)^j=0$
证明：$\sum\limits_{j=0}^{n-1} (\omega_n^k)^j=\frac{(w_n^k)^n-1}{w_n^k-1}=\frac{(w_n^n)^k-1}{w_n^k-1}=0(k \mod n \neq 0)$
特别地：$\sum\limits_{j=0}^{n-1} (\omega_n^k)^j=n(k \mod n = 0)$

前置技能已经get得差不多了，开始表演

离散型傅里叶变换和逆离散型傅里叶变换（DFT和IDFT）

应该还记得快速傅里叶变换要干什么吧。
快速插值。
我们要做的事情就是

(1)$已知A(x)=(a_0,a_1\dots a_n),求A(x)=(A(\omega^0),A(\omega^1)\dots A(\omega^n))$
(2)$已知A(x)=(A(\omega^0),A(\omega^1)\dots A(\omega^n)),求A(x)=(a_0,a_1\dots a_n)$

刚才介绍了这么多优秀的单位复根的性质，于是我们容易想到，把单位复根带入多项式里面。(为了方便，我们用$\omega$表示$\omega_n$)
假设序列$A和B$乘法之后得到的序列是$C$，其长度是$A和B$长度之和减一。我们找到一个最小n使得n大于C的长度并且n是2的整数次幂（为啥？之后会提到。）
我们可以把卷积变成这样的形式。
$c_r=\sum\limits_{p,q}[p+q \mod n=r]a_pb_q$
由求和引理可得$\frac{1}{n}\sum\limits_{j=0}^{n-1} (\omega^k)^j=[k \mod n = 0]$
于是$[p+q \mod n=r]=[p+q-r \mod n=0]=\frac{1}{n}\sum\limits_{k=0}^{n-1} (\omega^{p+q-r})^k=\frac{1}{n}\sum\limits_{k=0}^{n-1} \omega^{-rk}\omega^{pk}\omega^{qk}$
$c_r=\sum\limits_{p,q}[p+q-r \mod n=0]a_pb_q$
$=\sum\limits_{p,q}\frac{1}{n}\sum\limits_{k=0}^{n-1} \omega^{-rk}\omega^{pk}\omega^{qk}a_pb_q$
$=\frac{1}{n}\sum\limits_{k=0}^{n-1} \omega^{-rk}\sum\limits_p\omega^{pk}a_p\sum\limits_q\omega^{qk}b_q$
$=\frac{1}{n}\sum\limits_{k=0}^{n-1} \omega^{-rk}A(\omega^k)B(\omega^k)=\frac{1}{n}\sum\limits_{k=0}^{n-1} \omega^{-rk}C(\omega^k)$
刚才经过一波推导我们成功地找到了$C$点值表达式和$C$系数表达式的关系。
这样子的话问题转化为

(1)$求A(\omega^k)=\sum\limits_p\omega^{kp}a_p$
(2)$求a_r=\frac{1}{n}\sum\limits_r\omega^{-kr}A(\omega^k)$
前者即为DFT，后者即是IDFT

我们不难发现，两者的过程惊人的相似，其实只是多了$\frac{1}{n}$和一个-
于是我们使用同一个算法——快速傅里叶算法（FFT）来实现这玩意儿。

快速傅里叶变换（FFT）

设$A_0(x)是A(x)偶次项的和，A_1(x)是A(x)奇次项的和$
注意到$(\omega_n^m)^2=(\omega_n^{m+\frac{n}{2}})^2=\omega_{\frac{n}{2}}^m$
$A(\omega_n^m)=A_0((\omega_n^m)^2)+\omega_n^mA_1((\omega_n^m)^2)=A_0(\omega_{\frac{n}{2}}^m)+\omega_n^mA_1(\omega_{\frac{n}{2}}^m)$
$A(\omega_n^{m+\frac{n}{2}})=A_0((\omega_n^m)^2)+\omega_n^{m+\frac{n}{2}}A_1((\omega_n^m)^2)=A_0(\omega_{\frac{n}{2}}^m)-\omega_n^mA_1(\omega_{\frac{n}{2}}^m)$
这就是单位复根最英霸的地方，折半引理和消去引理可以使得它能够把插值的过程分治。上述操作被称为蝴蝶操作。
上文有提到，n是2的整数次幂，所以这个过程可以一直递推下去。不难发现，上述算法的时间复杂度递推式是。
$T(n)=2T(\frac{n}{2})+O(logn)$
由主定理可得时间复杂度为$O(nlogn)$
IDFT有两种写法，第一种是老老实实地带一个负号进去，其实还有第二种写法。
$a_r=\frac{1}{n}\sum\limits_r\omega^{-kr}A(\omega^k)=\frac{1}{n}\sum\limits_r\omega^{(n-k)r}A(\omega^k)$
我们把A数组反转一下，DFT之后除以n就好了。

代码实现

注意到$\omega_n^k=e^{\frac{2\pi ki }{n}}=cos(\frac{2\pi k}{n})+isin(\frac{2\pi k}{n})$
所以我们全程是用欧拉公式来表示$\omega$的
还有一点，蝴蝶操作本来我们是要迭代的，但是这里有一个优化常数的小技巧。
考虑原序列是$(a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7)$
模拟蝴蝶变换的过程。
$(a_0,a_2,a_4,a_6,a_1,a_3,a_5,a_7)$
$(a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7)$
如果用二进制表示原序列和变换后的序列
原$(000,001,010,011,100,101,110,111)$
后$(000,100,010,110,001,101,011,111)$
发现其实就是二进制反转了。
然后我们就可以用递推写这个东西了，具体看代码吧。

//luoguP3803 【模板】多项式乘法（FFT） 
//一个比较易于理解的版本，其实可以写得更简洁。
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
const int N = 5e6 + 10;
const double pi = acos(-1.0);
struct cp {
    double r, i;
    cp(double _r = 0, double _i = 0) : r(_r), i(_i) {}
    cp operator + (cp a) {return cp(r + a.r, i + a.i);}
    cp operator - (cp a) {return cp(r - a.r, i - a.i);}
    cp operator * (double a) {return cp(r * a, i * a);}
    cp operator * (cp a) {return cp(r * a.r - i * a.i, r * a.i + i * a.r);}
}a[N], b[N];
int r[N], m;
void FFT(cp *F, int f) {
    for(int i = 0;i < m; ++i) if(i < r[i]) swap(F[i], F[r[i]]);
    for(int i = 1; i < m; i <<= 1) {
        cp wn(cos(pi / i), f * sin(pi / i));
        for(int j = 0;j < m; j += (i << 1)) {
            cp w(1, 0);
            for(int l = 0;l < i; ++l, w = w * wn) {
                cp x = F[j + l], y = w * F[j + i + l];
                F[j + l] = x + y; F[j + i + l] = x - y;
            }
        }
    }
    if(!~f) for(int i = 0;i < m; ++i) F[i].r /= m;
}
int main() {
    int n1, n2; scanf("%d%d", &n1, &n2);
    for(int i = 0, x;i <= n1; ++i) scanf("%d", &x), a[i] = cp(x, 0);
    for(int i = 0, x;i <= n2; ++i) scanf("%d", &x), b[i] = cp(x, 0);
    int L = 0; for(m = 1; (m <<= 1) <= (n1 + n2); ++L) ;
    for(int i = 1; i < m; ++i) r[i] = (r[i >> 1] >> 1) | (i & 1) << L;
    FFT(a, 1); FFT(b, 1); for(int i = 0;i < m; ++i) a[i] = a[i] * b[i];
    FFT(a, -1); 
    for(int i = 0;i <= n1 + n2; ++i) printf("%d ", (int)(a[i].r + 0.5)); puts("");
    return 0;
}

后记

参考博文

[学习笔记] 多项式与快速傅里叶变换(FFT)基础
Pick‘s Blog 里面有各种FFT系列的东西

后续博文

算法学习FFT系列(2):快速数论变换NTT
算法学习FFT系列(3):多项式求逆详解——NTT+分治
占坑！