快速傅里叶变换详解(FFT)

自己也看了几篇博客，但是对我这种不擅长推导小白来说还是有一点困难，所以自己也写一篇博客也为像我一样的小白提供思路。以下内容包含各种LaTeX渲染，如果哪里有错误欢迎大家评论留言，或者添加本人qq：1403482164（无事勿扰）

一、FFT的应用场景

(A(x) ext{=} a_0 ext{+} a_1x+a_2x^2+……+a_nx^n)

(B(x) ext{=} b_0 ext{+} b_1x+b_2x^2+……+b_mx^m)

(C(x) ext{=} A(x) imes B(x)) 求解(C)的各项系数

二、前置知识

1. 复数

• 复数的表示形式：(a+bi = r(cos heta + isin heta) = re^{ heta i})，前两种高一课程应该都学过，最后请参考欧拉公式

• 复数在平面直角坐标系上的表示：

主要参考是复数的第一，二种表示形式，实部（(a / rcos heta)）为横坐标，虚部（(b / risin heta)）为纵坐标，也可以看成是半径为(r)的一个圆上的某一个点，圆上的点就可以表示出长度为(r)的所有复数

2. 单位根

(x = 1)，(x)有一个解为1，在平面直角坐标系中对应((1,0))

(x^2 = 1)，(x)有两个解为(1, ext{-}1)，在平面直角坐标系中对应((1,0))和(( ext{-}1,0))

(x^3 = 1)，(x)有三个解为(1,dfrac{ ext{-}1 ext{+}sqrt{3}i}{2},dfrac{ ext{-}1 ext{-}sqrt{3}i}{2})，在平面直角坐标系中对应((1,0),( ext{-}frac{1}{2},frac{sqrt{3}}{2}),( ext{-}frac{1}{2}, ext{-}frac{sqrt{3}}{2}))

……

依次类推当(w^k = 1)时，(w)有(k)个解，且这些解在平面直角坐标系中表示出来后把半径为1的圆平分成(k)等份，这样的解叫做单位根，其中(w_{n}^{k})表示(w^n = 1)的第(k)个解

• 单位根的性质

根据复数的表示形式(w_{n}^{k} = cos heta ext{+}isin heta = re^{ heta i})，其中( heta = frac{2pi k}{n})

通过单位根的形式变化和三角函数计算，可以推出他的一些性质：

1. ({(w_{n}^{k})}^2 = w_{n}^{2k} = w_{frac{n}{2}}^{k})

2. (w_{2n}^{2k} = w_{n}^{k})

3. (w_{n}^{frac{n}{2}+k} = ext{-}w_{n}^{k})

3. 矩阵乘法

给出(A(x) = a_0+a_1x+……a_nx^n)

(A = egin{bmatrix}{x_1}^0&{x_1}^1&cdots&{x_1}^n\{x_2}^0&{x_2}^1&cdots&{x_2}^n\vdots&vdots&vdots&vdots\{x_n}^0&{x_n}^1&cdots&{x_n}^nend{bmatrix}，B = egin{bmatrix}a_0\a_1\vdots\a_nend{bmatrix}，C = egin{bmatrix}A(x_1)\A(x_2)\vdots\A(x_n)end{bmatrix})

三个矩阵的关系显而易见：(A imes B = C)，已知(A,B)的时候自然可求出(C)，如果已知(A,C)如何求(B)呢？

我们引入了两个新的名词：单位矩阵和逆矩阵

1. 单位矩阵(I)：对角线为1，其余全为0，且满足(A imes I = A)

2. 逆矩阵：若(A imes A^{-1} = I)则称(A^{-1})为(A)的逆矩阵

在等式两边同时乘(A^{-1})变成：(A imes A^{-1} imes B = C imes A^{-1} Rightarrow I imes B = C imes A^{-1} Rightarrow B = C imes A^{-1})

这样我们就可以求解(B)矩阵了。

4. 函数的表示方法：

1. 系数法：已知所有项的系数当然可以确定一个函数

2. 点值法：(n)项的函数，在平面直角坐标系中找(n+1)个点就可以确定这个函数

具体证明就不给了，毕竟oi大多数时候还是考感性理解的

三、FFT

讲了个这么多，终于把前置知识处理完毕了，接下来就是正菜FFT了，不同于其他博客，可能我不会引用专业名词像DFT等等……但是精华都是一样的，相信大家都有超前的思考能力，发现某个问题不知道如何处理，耐心看下去，也许会发现不一样的精彩。

1. 为什么用FFT，他如何优化问题

首先还是这个式子：(A(x) = a_0+a_1x^2+……a_nx^n)

当我给出一个(x)，朴素求解的时间复杂度是多少呢？ 循环叠加(x)，(O(n))求解，而我们需要(n+1)个点，那么确定一个函数的复杂度就是(O(n^2))的

显然不够优秀是吧，而他的困难就在于对每个(x)求(A(x))的过程，所以我们想借助某种方法把这个过程优化到(O(logn))，就发明了FFT

2. 具体流程

1. 先想办法把式子处理一下：

奇偶分离(方便推导就假设(n)为偶数，实际写代码的时候，往后加0系数就可以了)：

(A_0(x) = a_0+a_2x^2+……a_nx^n，A_1(x) = a_1x+a_3x^3+……+a_{n+1}x^{n+1})，其中(a_{n+1} = 0)，为了让两个函数的项数相同

此时如果(A_0,A_1,A)三者的形式相同就好了，这样就可以统一求解了

我们可以发现如果我从(A_1(x))中提取一个(x)出来，那么他们之间的形式是一样的，所以我让(xA_1(x) = a_1x+a_3x^3+a_{n+1}x^{n+1}),那么(A_1(x) = a_1 +a_3x^2+……+a_{n+1}x^n)

那么现在(A_0)和(A_1)的形式相同，我们还需要把他们两个和(A)统一形式，很简单，我们把(x)替换成(y = x^2)，此时(A_0(y) = a_0+a_2y+a_4y^2+……+a_ny^{frac{n}{2}})，(A_1(y) = a_1+a_3y+……+a_{n+1}y^{frac{n}{2}})

到现在为止，我们已经把(A_0)和(A_1)从(A)中提取出来了，又要怎么合并起来呢？

其实这里面还有一个隐藏的关系：(A(x) = A_0(x^2)+xA_1(x^2))，往上面的式子里代一下，就可以得知了，这样分成两部分算，时间缩短了一半

2. 求解点值

巧妙的地方来了，除了式子的处理，我代入的点也十分的有讲究。大家也能猜到，前置知识讲的单位根肯定不是白讲的吧。

把单位根(w_{n}^{k})代入可得：

当(k leq frac{n}{2})时，(A(w_{n}^{k}) = A_0({w_{n}^{k}}^2)+w_{n}^{k}A_1({w_{n}^{k}}^2) = A_0(w_{frac{n}{2}}^{k})+w_{n}^{k}A_1(w_{frac{n}{2}}^{k}))

当(k ext{>} frac{n}{2})时,把(k+frac{n}{2})代入，(A(w_{n}^{k+frac{n}{2}}) = A_0({w_{n}^{k+frac{n}{2}}}^2)+w_{n}^{k+frac{n}{2}}A_1({w_{n}^{k+frac{n}{2}}}^2) = A_0(w_{frac{n}{2}}^{k})-w_{n}^{k}A_1(w_{frac{n}{2}}^{k}))

我们发现两个情况只有加减号不同，所以在求解前一半的时候，后一半同时可以求解，时间又缩小了一半，(O(logn))的复杂度就出来了

3. 求解最终函数

算法进行到了最最最最后一步了，现在我们通过上述的一系列计算，已知了(A(x))和(B(x))的(n+m+1)个点，(C(x) = A(x)*B(x))，自然就得出了(C(x))的(n+m+1)个点，这就要考察到我们前置内容中的矩阵乘法部分了。

我们已知(A = egin{bmatrix}{w_n^1}^0&{w_n^1}^1&cdots&{w_n^1}^n\{w_n^2}^0&{w_n^2}^1&cdots&{w_n^2}^n\vdots&vdots&vdots&vdots\{w_n^n}^0&{w_n^n}^1&cdots&{w_n^n}^nend{bmatrix},C = egin{bmatrix}C(w_n^1)\C(w_n^2)\vdots\C(w_n^n)end{bmatrix})，来求解(B)矩阵

所以问题转化成求解(A^{-1})矩阵，给出一个结论，(A^{-1})矩阵就等于把(w)的所有上角标取负再乘(frac{1}{n})

以上就是FFT的全部流程，剩下还有动态规划优化，敬请期待下回详解……

四、递归版代码

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
const double pi = acos(-1.0);const int maxn = 1e7+10;
int read(){
	int x = 1,a = 0;char ch = getchar();
	while (ch < '0'||ch > '9'){if (ch == '-') x = -1;ch = getchar();}
	while (ch >= '0'&&ch <= '9'){a = a*10+ch-'0';ch = getchar();}
	return x*a;
}
struct node{
    double x,y;
}a[maxn], b[maxn];
node operator + (node a,node b){return node{a.x+b.x,a.y+b.y};}
node operator - (node a,node b){return node{a.x-b.x,a.y-b.y};}
node operator * (node a,node b){return node{a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y};}
int n,m;
void fft(int len, node *a, int op){
    if (len == 1) return;
    node a0[(len>>1)+3],a1[(len>>1)+3];
    for (int i = 0;i <= len;i += 2) a0[i>>1] = a[i],a1[i>>1] = a[i+1]; 
    fft(len>>1,a0,op);fft(len>>1,a1,op); 
    node wn = node{cos(2*pi/len),op*sin(2*pi/len)},w0 = node{1,0};
    for(int i = 0;i < (len >> 1);i++,w0 = w0*wn){
        a[i] = a0[i]+w0*a1[i];
        a[i+ (len>>1)] = a0[i]- w0*a1[i];
    }
}
int main(){
	//freopen("in.in","r",stdin);
	//freopen("out.out","w",stdout);
	n = read(),m = read();
	for (int i = 0;i <= n;i++) a[i].x = read();
	for (int i = 0;i <= m;i++) b[i].x = read();
	int len = 1;
	while (len <= n+m) len <<= 1;
	fft(len,a,1);fft(len,b,1);
	for (int i = 0;i <= len;i++) a[i] = a[i]*b[i];
	fft(len,a,-1);
    for (int i = 0;i <= n+m;i++) printf("%.0lf ",a[i].x/len);
	return 0;
}

被WC摧残的第一天，回来更新了……

首先来模拟一下递归的过程（盗用wiki）：

规律： 其实就是原来的序列，每个数用二进制表示，然后把二进制翻转对称一下，就是最终那个位置的下标。比如(x)是001，翻转是100，也就是4，即(x)最后的位置。我们称这个变换为位逆序置换（bit-reversal permutation，国内也称蝴蝶变换）。（把二进制写出来就很好发现了，就是想不到写二进制呀）

进而我们还可以发现，每一层对应的两个数都是由他下面一层对应的两个数更新而来，这样问题就转化成了求最后一层数的顺序了，也就是求每个数的位逆序……

因为每个数都要求一遍，所以想能到dp转移

状态定义：(dp_i)表示(i)的位逆序

状态转移：(dp_i = frac{dp_{i/2}}{2} ext{+} (i ext{%}2) imes 2^{l-1})(自己模拟一下过程，还是挺显而易见的)

五、非递归版代码(话说luogu卡精度来着)

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
const double pi = acos(-1.0);const int maxn = 1e7+10;
int read(){
	int x = 1,a = 0;char ch = getchar();
	while (ch < '0'||ch > '9'){if (ch == '-') x = -1;ch = getchar();}
	while (ch >= '0'&&ch <= '9'){a = a*10+ch-'0';ch = getchar();}
	return x*a;
}
struct node{
    double x,y;
}a[maxn], b[maxn],a0[maxn],a1[maxn];
node operator + (node a,node b){node x = {a.x+b.x,a.y+b.y};return x;}
node operator - (node a,node b){node x = {a.x-b.x,a.y-b.y};return x;}
node operator * (node a,node b){node x = {a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y};return x;}
int n,m,dp[maxn];
double coss[maxn],sinn[maxn];
void fft(int len,node *a,int op){
	for (int i = 0;i <= len;i++){
		if (i < dp[i]) swap(a[i],a[dp[i]]);
	}
    for (int l = 1;l < len;l<<=1){
        node wn = {coss[l],op*sinn[l]};
        for (int i = 0;i < len;i+=(l << 1)){
            node w0 = {1,0};
            for (int j = 0;j < l;j++,w0 = w0*wn){
                node x = a[i+j],y = w0*a[i+j+l];
                a[i+j] = x+y,a[i+j+l] = x-y;
            }
        }
    }

}
int main(){
	//freopen("in.in","r",stdin);
	//freopen("out.out","w",stdout);
	n = read(),m = read();
	for (int i = 0;i <= n;i++) a[i].x = read();
	for (int i = 0;i <= m;i++) b[i].x = read();
	int len = 1,num = 0;
	while (len <= n+m) len <<= 1,num++;
	for (int i = 0;i <= len;i++) dp[i] = (dp[i>>1]>>1)|((i&1)<<(num-1));
	for (int i = 1;i <= len;i <<= 1) coss[i] = cos(pi/i),sinn[i] = sin(pi/i);
	fft(len,a,1);fft(len,b,1);
	for (int i = 0;i <= len;i++) a[i] = a[i]*b[i];
	fft(len,a,-1);
    for (int i = 0;i <= n+m;i++) printf("%d ",(int)(a[i].x/len+0.5));
	return 0;
}