「快速傅里叶变换(FFT)」学习笔记

FFT即快速傅里叶变换，离散傅里叶变换及其逆变换的快速算法。在OI中用来优化多项式乘法。

本文主要目的是便于自己整理、复习

---

FFT的算法思路

已知两个多项式的系数表达式，要求其卷积的系数表达式。

卷积：$$c_i = sumlimits_{j=0}^{i}a_jb_{i-j}$$

先将两个多项式分别转化为点值表达式，完成点值表达式的乘法，然后转为系数表达式得到结果。

点值表达式的乘法。整体考虑：假设已知两个多项式$A(x)$和$B(x)$。如果已知当$x=x_0$时$A(x_0)$和$B(x_0)$，则其乘积一定有点值$A(x_0)*B(x_0)$。因此点值表达式的乘法复杂度$O(n)$。

问题转化为如何快速完成系数表达式与点值表达式之间的转换。其中，我们的点值表达式都默认将单位根代入。因此，系数表达式的空间大小和点值表达式的空间大小是一样的。（牢记！）

---

离散傅里叶变换(DFT)的FFT

目的：快速求得某系数表达式对应的点值表达式。

假设$n$为2的幂，不足则补0.

对于一个$n-1$次多项式的系数表达式$A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}$，我们可以按照其指数的奇偶将其分为两部分。

$A_1(x)=a_0+a_2x+...+a_{n-2}x^{frac{n}{2}-1}$

$A_2(x)=a_1+a_3x+...+a_{n-1}x^{frac{n}{2}-1}$

（注意x的指数已经修改，不再是原来的指数）

则有$A(x)=A_1(x^2)+xA_2(x^2)$

设$k<dfrac{n}{2}$，将任意$omega_n^k$代入$A(x)$可得$$A(omega_n^k)=A_1(omega_{frac{n}{2}}^{k})+omega_n^kA_2(omega_{frac{n}{2}}^k)$$

将任意$omega_n^{k+frac{n}{2}}$代入$A(x)$可得$$A(omega_n^{k+frac{n}{2}})=A_1(omega_{frac{n}{2}}^{k})-omega_n^kA_2(omega_{frac{n}{2}}^k)$$

只要我能够知道$A_1(omega_{frac{n}{2}}^{k})$和$A_2(omega_{frac{n}{2}}^k)$，那么就可以求得$A(omega_n^k)$和$A(omega_n^{k+frac{n}{2}})$。当上式中的$k$取遍$[0,frac{n}{2}-1]$时，我们就可以得出另一半$[frac{n}{2},n-1]$。因此就得到了整个区间$[0,n-1]$，也就求出了$A(x)$的点值表达式。

所以：此时我们想要知道的是当$k$为$[0,frac{n}{2}-1]$内的整数时，$A_1(omega_{frac{n}{2}}^{k})$和$A_2(omega_{frac{n}{2}}^k)$的取值。而这就是$A_1(x)$和$A_2(x)$的点值表达式。

因此我们只需要递归求点值表达式即可。（这一步是关键中的关键，这决定了为什么FFT可以分治；容易看出这样的做法和线段树非常类似，复杂度$O(nlogn)$）

---

离散傅里叶逆变换(IDFT)的FFT

目的：将点值表达式转化回系数表达式

我们已经把两个系数表达式转为了点值表达式并完成了乘法，设此结果的点值表达式为为$(y_0,y_1,...,y_{n-1})$。先在要求它对应的系数表达式，设为$(a_0,a_1,...,a_{n-1})$

我们构造一个以$(y_0,y_1,...,y_{n-1})$为系数（注意不是点值）的多项式$B(x)=y_0+y_1x+...+y_{n-1}x^{n-1}$。此时多项式$B(x)$是系数表示的。我们依次将单位根的倒数，即$omega_n^0,omega_n^{-1},...,omega_n^{1-n}$代入得到一个点值表达式$(z_0,z_1,...,z_{n-1})$。此时我们有结论（奇妙的是它竟然将虚部消去了）$$a_k=dfrac{z_k}{n}$$证明见附录。

因此我们只需要把我们得到的点值表达式看做是系数表达式，再做一次FFT（其中单位根以倒数形式代入）就能够间接我们所要的结果了。（之所以说间接，是因为还要除以$n$）

---

代码实现

FFT的代码实现依然是一个难点。

递归实现

• 参数如何设置？考虑传入数组指针：传入系数表达式，返回点值表达式（一个系数表达式的空间应当和一个点值表达式是一样的）
• 假设一个$M$次的多项式与一个$N$次的多项式相乘，那么得到的结果是$M+N$次的。也就是需要一个$M+N+1$项的点值表达式来转化。而多项式的乘法仅仅是对应的点值相乘，所以我们需要保证先前求出的两个点值表达式也是$N+M+1$项的。
• 输出结果时如果草率地输出.0f会有"-0"这样子的局面……因此要转int。
• DFT和IDFT可以共用同一个函数

//其中n为本次多项式的次数(n-1)；a为数组指针，传入系数表达式，返回点值表达式；cvs的值为1或-1，表示单位根是否需要倒数 
void FFT(int n, Complex* a, int cvs){
    if(n == 1) return;//到达叶节点，此时多项式是一个常数，也就是说点值表达式应该与系数表达式相同，不需要转换，直接返回 
    Complex a1[n>>1],a2[n>>1];//左右子多项式的辅助数组(空间只需要一半) 
    for(int i = 0; i < (n>>1); ++i){
        a1[i] = a[i<<1];
        a2[i] = a[i<<1|1];
    }
    FFT(n>>1,a1,cvs);
    FFT(n>>1,a2,cvs);
    Complex Wn=Complex(1,0), w = Omega(n,cvs);//Wn表示代入的单位根，每次都需要乘上w来增加(减少)指数 
    for(int i = 0; i < (n>>1); ++i){
        a[i] = a1[i] + Wn*a2[i];
        a[i+(n>>1)] = a1[i] - Wn*a2[i];
        Wn = Wn * w;
    }
}

迭代实现

递归实现不仅常数巨大，而且我们发现它在递归的过程中需要两个辅助数组，导致空间巨大。

依然用线段树类比FFT。众所周知，一种大众的线段树建树方法是从根节点递归左右子树建，就好像递归版的FFT一样——递归左右的点值表达式，然后合并为自己的。zkw线段树是如何建树的？先确定好叶节点，然后一步一步向上合并即可。

二进制倒序的发现

那么FFT也可以这样实现。唯一的区别在于如何事先确定叶节点的值？合并左右点值表达式这一方法是建立在左右按奇偶分开的前提下的，所以我们需要知道递归到底层时的系数顺序是如何的。利用二进制去考虑奇偶的分配（有关奇偶一般都和二进制有关……）。第一层（顶层）的分配看二进制的最后一位：是0去左边，是1去右边；到了第二层以后就看的是倒数第二位；看倒数第三位；于是我们只要看这个点是沿着什么样的路径走来的就知道这是哪一个系数了。而这个叶节点究竟在什么位置恰恰和这一条路径反过来。按大小分的是和trie树一样从前往后考虑，按奇偶分的是从后往前考虑。因此系数的位置就是其初始位置的二进制倒序。例如，系数$a_6$就应该在位置$3$，因为$110$倒序得$011$。

空间处理与蝴蝶操作

我们由递归改为了迭代，也就没有了数组指针传入返回的概念了。递归版本中，转化成的点值表达式直接作为结果返回，那么迭代中怎么办？一个很自然的想法是全都用一个辅助数组存起来，但是这样的空间就又回到递归一样巨大了；进一步考虑发现当前这一层只需要利用上一层，所以可以用一个类似滚动数组的方法两个数组解决。

利用系数表达式的空间和点值表达式一样的特点，我们可以只使用一个数组。我们发现：在扫描$k$的过程中（配合代码理解），当前这一轮我们需要利用的无非两个点值$A_1(omega_{frac{n}{2}}^{k})$和$A_2(omega_{frac{n}{2}}^k)$，我们之前将他们存在a[j+k]和a[j+k+n/2]中（此时的a还是上一轮处理完了的点值）。现在我们利用他们造出了$A(omega_n^k)$，理所应当存在a[j+k]中，覆盖掉a[j+k]；同时我们推出了$A(omega_n^{k+frac{n}{2}})$，应当存在a[j+k+n/2]中，覆盖掉原来的a[j+k+n/2]。于是我们发现新一轮的点值表达式完美覆盖了原来的点值表达式。我们把这个完美覆盖的过程叫做蝴蝶操作。

/*By DennyQi 2018*/
#include <cstdio>
#include <queue>
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
const int MAXN = 4000010;
const double Pi = acos(-1);//cos(pi)=-1,所以arccos(-1)=pi 
inline int read(){
    int x = 0; int w = 1; register char c = getchar();
    for(; c ^ '-' && (c < '0' || c > '9'); c = getchar());
    if(c == '-') w = -1, c = getchar();
    for(; c >= '0' && c <= '9'; c = getchar()) x = (x<<3) + (x<<1) + c - '0'; return x * w;
}
//定义复数类，并重载运算符 
struct Complex{
    double a,b;
    Complex(){ 
        a = 0, b = 0; 
    }
    Complex(const double x, const double y){
        a = x, b = y; 
    }
    Complex operator + (const Complex u){ 
        return Complex(a + u.a,b + u.b); 
    }
    Complex operator - (const Complex u){ 
        return Complex(a - u.a,b - u.b); 
    }
    Complex operator * (const Complex u){ 
        return Complex(a*u.a - b*u.b, a*u.b + b*u.a); 
    }
};
int N,M,Lim,lg;
int r[MAXN];
Complex a[MAXN],b[MAXN];
inline void FFT(int n, Complex* a, int cvs){
    //先求最终状态 
    for(int i = 0; i < n; ++i){
        if(i < r[i]) swap(a[i], a[r[i]]);
    }
    Complex A1,A2,Wnk,w;
    for(int i = 1; i <= n; i<<=1){//共log(n)层，自底向上枚举；i同时还表示当前这一层的节点宽度
        for(int j = 0; j < n; j+=i){//一次枚举这一层的每一个节点 
            Wnk = Complex(1,0), w = Complex(cos(2*Pi/i*cvs),sin(2*Pi/i*cvs));
            for(int k = 0; k < (i>>1); ++k){//对于每一个节点，枚举代入的单位根 
                A1 = a[j+k], A2 = a[j+k+(i>>1)];//蝴蝶操作 
                a[j+k] = A1 + Wnk * A2;
                a[j+k+(i>>1)] = A1 - Wnk * A2;
                Wnk = Wnk * w;
            }
        }
    }
}
int main(){
    freopen("FFT.in","r",stdin);
    scanf("%d%d",&N,&M);
    for(int i = 0; i <= N; ++i) scanf("%lf",&a[i].a);
    for(int i = 0; i <= M; ++i) scanf("%lf",&b[i].a);
    Lim = 1;
    while(Lim < N+M) Lim<<=1, ++lg;//补零 
    for(int i = 0; i < Lim; ++i){
        r[i] = ((r[i>>1]>>1) | ((i&1)<<(lg-1)));
        //求二进制倒序。采用递推的形式：先不考虑最后一位，将前面部分翻转，再将最后一位补到前面去。 
    }
    FFT(Lim,a,1);
    FFT(Lim,b,1);
    for(int i = 0; i < Lim; ++i) a[i] = a[i] * b[i];//点值相乘，完成点值表达式的乘法 
    FFT(Lim,a,-1);
    //注意输出时要转为int，否则会出现-0等情况 
    for(int i = 0; i <= N+M; ++i) printf("%d ",(int)(a[i].a/Lim+0.5));
    return 0;
}

 实测比递归版快约3.05倍（滑稽）

附  录

多项式的表示方法

系数表示法：

即用$n+1$个系数来表示一个$n$次多项式，即$(a_0,a_1,...,a_n)$。

我们可以将$(a_0,a_1,...,a_n)$看做一个系数向量。

点值表示法：

当自变量$x$分别取$n+1$个不同的值时，会得到$n+1$个对应的结果。可以将多项式看做一个$n$次函数，这个函数被$n+1$个点$(x_0,y_0),(x_1,y_1),...,(x_n,y_n)$确定。因此我们可以用点值表达式$(A(x_0),A(x_1),...,A(x_n))$来表示多项式$A(x)$

我们可以将$(A(x_0),A(x_1),...,A(x_n))$看做一个点值向量，也就等同于$(y_0,y_1,...,y_n)$

复数

形如$a+bi$，可以理解为复平面上的向量$(a,b)$

复数的乘法法则几何意义上可以理解为模长相乘，幅角相加。代数意义上可以理解为$(a+bi)(c+di)=ac-bd+(ad+bc)i$

单位根

在复平面内以原点为圆心，1为半径作圆。称这个圆为单位圆。

以$(1,0)$为起点，将圆$n$等分。设点$(1,0)$为第0个，那么这些等分点称为$n$次单位根。记作$omega_n^k$。由于模长均为1，所以根据复数乘法法则，这$n$个等分点可以依次表示为$omega_n^0$，$omega_n$，$omega_n^2$，...，$omega_n^{n-1}$

根据三角函数，很容易推得单位根的表示方法：$omega_n^k=cosdfrac{2kpi}{n}+isindfrac{2kpi}{n}$ （欧拉公式）

性质一：$omega_{2n}^{2k}=omega_n^k$ （表示的是同一个点）

性质二：$omega_n^k=-omega_n^{k+frac{n}{2}}$ （关于原点中心对称）

$a_k=dfrac{z_k}{n}$的证明

$$ egin{align*} z_k &= sum_{i=0}^{n-1}y_i(omega_n^{-k})^i\ &= sum_{i=0}^{n-1}(sum_{j=0}^{n-1}a_j*omega_n^{ij})(omega_n^{-k})^i\ &= sum_{i=0}^{n-1}sum_{j=0}^{n-1}a_j*(omega_n^{j-k})^i\ &= sum_{i=0}^{n-1}a_j(sum_{j=0}^{n-1}(omega_n^{j-k})^i)\ end{align*} $$ 

而因为后半部分$sum_{j=0}^{n-1}(omega_n^{j-k})^i$可由等比数列求和公式得：

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

注意，要使此式有意义，需满足分母不为0。而当$j=k$时分母为0。因此分类讨论：

当$j eq k$时，$sum_{j=0}^{n-1}(omega_n^{j-k})^i=0$。

当$j = k$时，$sum_{j=0}^{n-1}(omega_n^{j-k})^i=n$

综上，$z_k = sumlimits_{i=0}^{n-1}a_j(sumlimits_{j=0}^{n-1}(omega_n^{j-k})^i) = na_k$。所以$a_k=dfrac{z_k}{n}$