快速傅里叶变换详解(FFT)
一、FFT的应用场景
\(A(x) \text{=} a_0 \text{+} a_1x+a_2x^2+……+a_nx^n\)
\(B(x) \text{=} b_0 \text{+} b_1x+b_2x^2+……+b_mx^m\)
\(C(x) \text{=} A(x)\times B(x)\) 求解\(C\)的各项系数
二、前置知识
1. 复数
-
复数的表示形式：\(a+bi = r(cos\theta + isin\theta) = re^{\theta i}\)，前两种高一课程应该都学过，最后请参考欧拉公式
-
复数在平面直角坐标系上的表示：
主要参考是复数的第一，二种表示形式，实部（\(a / rcos\theta\)）为横坐标，虚部（\(b / risin\theta\)）为纵坐标，也可以看成是半径为\(r\)的一个圆上的某一个点，圆上的点就可以表示出长度为\(r\)的所有复数
2. 单位根
\(x = 1\)，\(x\)有一个解为1，在平面直角坐标系中对应\((1,0)\)
\(x^2 = 1\)，\(x\)有两个解为\(1,\text{-}1\)，在平面直角坐标系中对应\((1,0)\)和\((\text{-}1,0)\)
\(x^3 = 1\)，\(x\)有三个解为\(1,\dfrac{\text{-}1\text{+}\sqrt{3}i}{2},\dfrac{\text{-}1\text{-}\sqrt{3}i}{2}\)，在平面直角坐标系中对应\((1,0),(\text{-}\frac{1}{2},\frac{\sqrt{3}}{2}),(\text{-}\frac{1}{2},\text{-}\frac{\sqrt{3}}{2})\)
……
依次类推当\(w^k = 1\)时，\(w\)有\(k\)个解，且这些解在平面直角坐标系中表示出来后把半径为1的圆平分成\(k\)等份，这样的解叫做单位根，其中\(w_{n}^{k}\)表示\(w^n = 1\)的第\(k\)个解
- 单位根的性质
根据复数的表示形式\(w_{n}^{k} = cos\theta \text{+}isin\theta = re^{\theta i}\)，其中\(\theta = \frac{2\pi k}{n}\)
通过单位根的形式变化和三角函数计算，可以推出他的一些性质：
-
\({(w_{n}^{k})}^2 = w_{n}^{2k} = w_{\frac{n}{2}}^{k}\)
-
\(w_{2n}^{2k} = w_{n}^{k}\)
-
\(w_{n}^{\frac{n}{2}+k} = \text{-}w_{n}^{k}\)
3. 矩阵乘法
给出\(A(x) = a_0+a_1x+……a_nx^n\)
\(A = \begin{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}^n\end{bmatrix}，B = \begin{bmatrix}a_0\\a_1\\\vdots\\a_n\end{bmatrix}，C = \begin{bmatrix}A(x_1)\\A(x_2)\\\vdots\\A(x_n)\end{bmatrix}\)
三个矩阵的关系显而易见：\(A\times B = C\)，已知\(A,B\)的时候自然可求出\(C\)，如果已知\(A,C\)如何求\(B\)呢？
我们引入了两个新的名词：单位矩阵和逆矩阵
-
单位矩阵\(I\)：对角线为1，其余全为0，且满足\(A\times I = A\)
-
逆矩阵：若\(A\times A^{-1} = I\)则称\(A^{-1}\)为\(A\)的逆矩阵
在等式两边同时乘\(A^{-1}\)变成：\(A\times A^{-1}\times B = C\times A^{-1} \Rightarrow I\times B = C\times A^{-1} \Rightarrow B = C\times A^{-1}\)
这样我们就可以求解\(B\)矩阵了。
4. 函数的表示方法：
-
系数法：已知所有项的系数当然可以确定一个函数
-
点值法：\(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 \text{>} \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 = \begin{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}^n\end{bmatrix},C = \begin{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} \text{+} (i\text{%}2)\times 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;
}
